Scalar Relativistic Effects with Multiwavelets: Implementation and Benchmark

The importance of relativistic effects in quantum chemistry is widely recognized, not only for heavier elements but throughout the periodic table. At the same time, relativistic effects are strongest in the nuclear region, where the description of electrons through a linear combination of atomic orbitals becomes more challenging. Furthermore, the choice of basis sets for heavier elements is limited compared with lighter elements where precise basis sets are available. Thanks to the framework of multiresolution analysis, multiwavelets provide an appealing alternative to overcoming this challenge: they lead to robust error control and adaptive algorithms that automatically refine the basis set description until the desired precision is reached. This allows one to achieve a proper description of the nuclear region. In this work, we extended the multiwavelet-based code MRChem to the scalar zero-order regular approximation framework. We validated our implementation by comparing the total energies for a small set of elements and molecules. To confirm the validity of our implementation, we compared both against a radial numerical code for atoms and the plane-wave-based code EXCITING.


Introduction
In his famous Nobel Price Lecture, P.A.M. Dirac stated that with the advent of quantum mechanics all fundamental problems of chemistry were in principle solved 1 .It is however interesting that he did not seem to realize that relativity would also play a role, despite his own fundamental contribution to combining relativity and quantum mechanics 2,3 .It is now widely accepted that the correct description of the electronic structure of atoms and molecules requires the inclusion of relativistic effects.This is particularly relevant for the core-and the innermost valence shell-electrons of heavier elements, which move at sizeable fractions of the speed of light.A wide number of chemical properties and phenomena, such as the yellow colour of gold 4,5 , the liquid state of mercury 5,6 , the functioning of lead-acid batteries 5,7 , the catalytic behavior of cobalt 5,[8][9][10] , the difficulties to describe fluorinated compounds with DFT [11][12][13] all depend on relativistic effects.There is a wide spectrum of methods to include relativistic effects in electronic structure calculations.They range from the use of effective core potential (ECP) [14][15][16][17] , also known as pseudopotentials in the solid state physics community 18 , to the solution of the 4-component Dirac equation.In between these two extremes are a number of methods with different degrees of accuracy and complexity 19 .
The ECP method collapses a given number of core-electrons and the nucleus to provide a core potential for the valence electron.This has the dual advantage of reducing the number of explicit electrons and at the same time including relativistic effects implicitly.On the other hand, if core electronic properties such as high pressure chemistry, core-electron spectroscopy 20 and nuclear magnetic resonance (NMR) shielding constants 9 are considered, ECP methods fall short, and all-electron calculations are necessary.The four-component Dirac equation 2 constitutes the starting point for the treatment of relativistic effects, and the full Breit Hamiltonian [21][22][23][24][25] is the the most complete treatment of relativistic effects for a many-electron system.On the other hand, it is also the most computationally expensive: spinorbitals are 4-component complex functions.Spin is no longer a good quantum number in a relativistic framework and the simplified picture of two electrons with opposite spin shar-ing the same orbital is no longer valid 19,26,27 .Most operators couple the different components of a spinorbital, leading roughly to a factor 100 in computational cost, because coupling four complex functions requires 8×8 matrices instead of a real scalar [28][29][30][31] .Although progress has been made to make 4-component calculations faster and easily available [28][29][30][31] , it is still convenient to attempt approximations that promise great reduction in the computational cost at the price of reduced accuracy.The first step in such a hierarchy of approximations is the elimination of the small components through a Foldy-Wouthuysens (FWs) transformation 32,33 .The choice of transformation leads to different kinds of methods, such as the Pauli Hamiltonian [34][35][36] , the Regular Approximations (RAs) [37][38][39] , the Douglas-Kroll-Hess (DKH) Hamiltonian [40][41][42][43][44] , or the exact two-component (X2C) Hamiltonian method 19,[45][46][47][48][49][50][51][52][53] .Further reduction to a scalar method is also possible for the Pauli Hamiltonian and the RAs [37][38][39] .
In particular, the RAs has two interesting features: (1) the decoupling part of the FWs transformation can be expanded in a convergent series to recover the exact elimination, and (2) the renormalization part can be exactly incorporated in the wave function.The ZORA 34,37,[54][55][56] is the simplest form of Hamiltonian keeping only the zeroth order in both parts.The reduction to a scalar method is carried out by applying the Dirac identity and discarding the spin-orbit term.The advantage of ZORA is to keep most of the standard algorithms of quantum chemistry in their original form just re-scaling the kinetic energy by a function that includes the potential energy (see below Theory section) 34,37,[54][55][56] .The numerical treatment of this rescaling function can be challenging, because it displays a cusp at each nucleus 34,37,[54][55][56] .Standard approaches, based on atomic orbital expansions, can struggle to get an accurate description in this region.The issue is further aggravated for heavier nuclei, where such a correction is important, and at the same time the number of basis sets available is more limited and less is known about their true precision 57,58 .Some efforts to assess the precision of ZORA for all-electron calculations of heavier elements have been undertaken using Gaussian type orbitals (GTOs), but it is anyway challenging to assess the precision without an external reference 57,58 .One such option for hydrogen-like ions (He + , Ne 9+ , Ar 17+ , an so on) is constituted by the scaling properties of the ZORA Hamiltonian in a two-component framework, which yields the (scaled) exact Dirac energies for such a system 59 .However, for many-electron systems, and/or for a scalar relativistic approach, assessing the true precision of a GTO basis is challenging and several basis sets are developed to include relativistic effects.Examples of all-electron relativistic basis sets include the universal Gaussian basis set (UGBS) 60 , the atomic natural orbitals basis sets [61][62][63][64][65] , the X2C basis sets 66 , and the segmented all-electron relativistic contracted (SARC) basis sets [67][68][69][70][71] .Importantly, such basis sets have to be fitted to the chosen Hamiltonian (ZORA, DKH).This comes on top of the required fitting of the basis set to a given electronic structure method, and if relevant, to a particular property 71 .Although some of the known all-electron relativistic basis sets may provide good results for a particular method and property, their transferability is limited.The large number of available GTO basis sets is an indication that no single basis set is good enough to describe all properties of interest to sufficient precision 58 .
In recent years, multiwavelets (MWs) 72 have emerged as a powerful alternative to traditional local basis sets 73 .Their foundation based on multiresolution analysis 74 leads to a basis set that is not empirically parameterized.Robust error control [75][76][77] means that the user can set a finite but arbitrary target precision, and adaptive algorithms [78][79][80] ensure that the representation of molecular orbitals is automatically refined until the required precision is reached.
Instead of a plethora of bases to choose from, the user only needs to set the requested precision and the polynomial order of the basis, providing in practice a robust black-box, where straightforward numerical considerations guide the user's choice.MWs have proven reliable and robust in providing Hartree-Fock (HF) and density functional theory (DFT) benchmark results for energies 73 and properties, 81 and have ventured both towards post-HF methods 82 and relativistic treatments. 83 this work we present a MW implementation of the ZORA method in the MRChem code.
To verify the correctness of the implementation, we consider two alternative methods: a radial atomic solver implemented as a separate package 84 and linearized augmented plane wave (LAPW) implemented in the electronic-structure code exciting 85 .Both of these approaches provide means for systematic improvement of the precision and are capable of yielding total energies approaching the complete basis set (CBS) limit as demonstrated in Refs.84 and 86.

Theory and implementation
From the Dirac equation to the ZORA Hamiltonian We will briefly expose how the ZORA 34,37,[54][55][56] is derived starting from the 4-component Dirac equation of an electron, which describes a 4-component spinor Ψ, subject to a potential V : In the above equation p is the momentum operator, c is the speed of light, atomic units (m e = 1, ℏ = 1, e = −1) are assumed and α and β are defined as follows: .
The first step towards ZORA consists in applying the FW transformation to the Dirac Hamiltonian of Eq. ( 1): where the transformation matrix U = W 1 W 2 is a product of a decoupling matrix W 1 and a renormalization matrix W 2 : and R is the exact coupling between the large and the small components of a 4-spinorbits: The inverse potential term depends on the eigenvalue E, and this dependence can be expanded in a Taylor series as follows: .
Restricting the expansion to the zero-th order leads to RA.Once the renormalization W 2 is also considered, the following 2-component Hamiltonian is obtained: The ZORA Hamiltonian is finally obtained by Taylor-expanding the renormalization operator and retaining only the zero-order term: Using the Dirac identity we can separate the scalar-relativistic and spin-orbit contributions to the kinetic energy covered by the first and the second terms, respectively.Here, we keep only the scalarrelativistic part and obtain the following Hamiltonian: where in the last expression we have implicitly defined κ = (1 − V 2c 2 ) −1 .Given the ZORA Hamiltonian, a Kohn-Sham (KS) DFT implementation is then obtained by replacing the non-relativistic kinetic energy operator with its ZORA counterpart: It should be noted that in practical implementations the potential defining κ is usually not the full KS potential.Given the form of κ, the most important contribution is the nuclear attraction.Introducing Coulomb ( Ĵ) and exchange and correlation ( Vxc ) is possible, but the corresponding operator has to be recomputed numerically at each iteration.Additionally, operations such as function multiplications are difficult to perform in traditional linear combination atomic orbital basis representations.A common choice for GTO calculations is to use a fixed atomic potential, which is called the atomic-ZORA approximation 87 .
The ZORA kinetic operator can then be pre-computed and used throughout the calculation.
Atomic-ZORA has the advantage of being gauge-invariant 87 .

ZORA equations in a Multiwavelet framework
To obtain a MW implementation of the ZORA eigenvalue problem, it is necessary to transform the differential equation into an integral equation, in analogy with the non-relativistic case 75,77 .The standard KS equations can be concisely written as follows: where F is the Fock operator, φ i refer to an occupied molecular orbital, and F ij are the matrix elements of the Fock operator between two occupied orbitals, assuming a general non-canonical (non-diagonal) form.
Within the framework of KS-DFT, the Fock operator consists of the kinetic energy T , the nuclear attraction Vnuc , the Coulomb repulsion Ĵ, the Hartree-Fock exchange K scaled by some numerical factor λ ∈ [0, 1], and the exchange and correlation potential Vxc : In the non-relativistic domain, the coupled KS differential equations ( 12) can be rewritten in integral form, 88,89 by making use of the bound-state Helmholtz kernel where Ĝµ i = − ∇ 2 − µ 2 i −1 is the integral convolution operator associated with the boundstate Helmholtz kernel G(r) = e −µ i r /r, using µ i = √ −2F ii (i.e. the diagonal elements of the Fock matrix).
In the ZORA Hamiltonian, the kinetic energy operator becomes: where V Z = Vnuc + Ĵ + Vxc .Including Ĵ and Vxc does not pose any issue in the MW framework, other than the computational overhead of having to update the potential at every iteration, because all potentials are anyway treated on an equal footing using a numerical grid.The non-local Hartree-Fock exchange, in turn, does not seem to contribute significantly to V Z as shown in Ref. 90.
Inserting the ZORA kinetic operator into Eq.( 12) we obtain: In order to make use of the same framework as the non-relativistic implementation of Eq. ( 14), it is necessary to isolate the Laplacian and the diagonal element of the sum on the right-hand side, which together make up the bound-state Helmholtz operator Ĝ = . This is achieved first by division by κ, recalling that κ −1 = 1 − V Z /2c.The following integral equation is obtained: When c → ∞, κ → 1 and ∇κ → 0, and the non-relativistic form as in Eq. ( 14) is recovered.Eq.17 can therefore be seen as a level-shifted version of its non-relativistic counterpart.
Although it cannot be expected that the iterative solution of Eq. ( 17) will work for arbitrary shifts, the approach is justified by recalling that κ ≃ 1 almost everywhere, except close to the nuclei.Our tests indicate that it becomes more difficult to converge the above equation when the ZORA contribution becomes larger, either physically by going down the periodic table, or artificially by letting c → 0. To overcome the convergence issues for heavier nuclei (5 th row of the periodic table) we have therefore introduced a finite nucleus model, as described in Section .For elements beyond the 5 th row, one has to keep in mind that ZORA becomes questionable 91 .
An alternative approach to obtain the desired (∇ 2 + 2F ii ) term from Eq. ( 16) is to add a Laplacian term ( 1 2 ∇ 2 ) directly, thus avoiding division by κ.Our tests indicate that the strategy presented above works better, despite the additional singularity introduced.The main reason seems to be that the former strategy removes the Laplacian altogether, which is ill-conditioned in the discontinuous MW basis, whereas the latter keeps part of it on the right-hand side.

Implementation with Multiwavelets
An implementation of the ZORA method as outlined above, is currently in a development version of the MRChem package 92 , and is expected to appear in the next official release v1.2.
MRChem is a numerical quantum chemistry code based on a MW framework, in which all functions and operators are represented on their own fully adaptive multi-resolution numerical real-space grid.This allows for efficient all-electron treatment of medium to large molecules (hundreds of atoms) at Self Consistent Field (SCF) level of theory (both HF and DFT).
The κ function is computed as a point-wise map of the chosen ZORA potential V Z through and similarly for its inverse Both functions are represented on their own adaptive numerical grid, and subsequently treated as standard multiplicative potential operators in the SCF procedure of solving Eq. ( 17).

Point nucleus models
In a MW framework, the singularity of the nuclear potential can lead to numerical problems when V nuc is computed and used.In the non-relativistic case, the issue is circumvented by replacing the analytic 1/r potential with a smoothed approximation 75 : and then parameterized as u(r/s)/s, where s is a scalar smoothing parameter.The smoothing parameter depends on the nuclear charge Z and the desired precision ϵ as explained in Ref 75 The above prescription constitutes a numerical smoothing of the nucleus to avoid accidental infinities in the representations.It should not be confused with physical finite nucleus models 93 , which are common in relativistic methods.This numerical smoothing is much sharper than the common finite nucleus models, and it is meant to yield results which are -within the requested precision -equivalent to using a point charge.

Finite nucleus models
In order to overcome the numerical issues faced by the pointwise nuclei, we have introduced the Gaussian nuclear model, as described by Visscher et al. 93 .Not only does the model overcome the numerical problems of pointwise charges, but it is also a sounder physical description of larger nuclei.The nuclear charge is modelled as a Gaussian distribution with the normalisation prefactor ρ G 0 = eZ(ξ/π) 3/2 and the parameter ξ related to the the root-mean-square radius of the nucleus via the expression ξ = 3/⟨R 2 ⟩.To determine ⟨R 2 ⟩, we apply the empirical formula for Ref. 94: where A is the nuclear mass number and f m is a femtometer length unit (10 −15 m).We note that a choice of the nuclear mass A is depends on whether the abundance of isotopes is taken in to account.To avoid the ambiguity, we use the tabulated values of ⟨R 2 ⟩ provided in Ref. 93.
The potential for a Gaussian charge distribution can be represented in an analytic form by Visscher et al. 93 : which is the solution to the Poisson equation for the charge density defined by Eq. ( 22).
Among the different possibilities presented in Ref. 93

Methods for verification
To verify correctness of the ZORA implementation in MRChem, we use two other codes: a numerical atomic solver 84 and an all-electron full-potential LAPW code exciting.In this section, we briefly introduce them and provide details on the implementation of new features needed for making a direct comparison with MRChem.
The atomic solver assumes atoms with spherically symmetric densities, and the monoelectronic wave functions are represented as φ nℓm (r) = u nℓ (r)Y ℓm (r), where u nℓ (r) is a radial function defined on a one-dimensional radial grid.The SCF solver uses a similar approach as the MW framework described above.It reduces the non-relativistic Kohn-Sham equation to the integral form as follows: This equation matches Eq. ( 14), with vanishing off-diagonal terms of the Fock matrix, because the canonical representation is used.The radial solver originally supported only non-relativistic calculations: to implement ZORA, we used Eq. ( 17) in the canonical form, i.e., replacing F ii with ϵ nℓ and setting F ij = 0 for i ̸ = j.This approach avoids the evaluation of second derivatives and the corresponding numerical noise which accumulates during the self-consistency iterations.
To consider systems beyond atoms, we use exciting.In a nutshell, the code relies on partitioning the unit cell into non-overlapping atomic spheres and the interstitial region.In the atomic spheres, the wavefunctions are expressed in terms of atomic-like orbitals that are updated during the self-consistence cycle.In the interstitial region, one represents the wavefunctions with plane waves using the smoothness of the Kohn-Sham potential.Based on such an approach, we use two types of basis functions: (i) augmented plane-waves (APWs) 95,96 and (ii) local-orbitals (LOs) 97 .Each APW combines a plane wave in the interstitial region with atomic-like orbitals in the spheres, whereas each LO is a linear combination of two atomic-like orbitals in one particular sphere and strictly zero everywhere else.As shown in Ref. 86, it is possible to obtain a systematic convergence of the total energies simply by increasing the number of APWs and LOs.
The ZORA Hamiltonian is already available in the released version of the code, whereas the smeared nucleus feature described in Sec. is implemented in the present work.Aside from adopting the electrostatic potential due to a Gaussian charge density, we also revise the evaluation of the following integral: If a nucleus at the site R 0 is defined as a point charge, its density distribution is Eq. 26 is evaluated as an integral on a radial grid.

Computational details
All calculations were performed with the PBE functional 98 using the XCFun 99 library, the LibXC 100 library and the native implementation in MRChem, the atomic solver and exciting, respectively.Being aware of slight inconsistencies in the PBE parameters employed in these libraries, we set β = 0.06672455060314922 and µ = 0.066725 π 2 3 as defined in XCFun and use these values in calculations with all three codes.For all MW calculations, a development version of MRChem has been employed.Interpolating polynomials of 9 th order were used with a simulation box size of 128 a 0 .The numerical precision thresholds used are summarized in Table 1.

Parameter
Value Explanation world prec 1.0e-6 Overall numerical precision energy thrs 1.0e-6 Convergence threshold total energy orbital thrs 1.0e-4 Convergence threshold maximum orbital residual All exciting calculations were performed using large cubic unit cell with a side length of 25 a 0 for atoms and 30 a 0 for molecules.The electrostatic interaction of the periodic images was eliminated by introducing the truncation of the Coulomb interaction following the approach explained in Ref. 81.Aside from acquiring the isolated limit, this adjustment allows us to obtain ZORA energies consistent with both MRChem and the atomic solver.
The unmodified electrostatic potential for periodic densities is defined uniquely except for an additive constant which introduces an ambiguity in the ZORA energies 34,37,[54][55][56] .The truncation of the Coulomb interaction makes the potential unique and thus removes this ambiguity completely 34,37,[54][55][56] .The canonical orbitals were expressed in terms of local orbitals and LAPWs with the cutoff R MT K max sufficient to ensure a few µHa precision.The specific settings in the case of each atom and molecule are stored in Table 2 and a set of input and output data files are available in the repository.The calculations with the atomic solver employed a 5 th order polynomial on a radial grid with the innermost and outermost points r min = 10 −8 a 0 and r max = 35 a 0 .The number of radial points was set to 5000, which is fully sufficient to guarantee a sub-µHa precision.
The total energies of diatomic molecules were calculated in MRChem and exciting without geometry relaxation using the internuclear distances given in Table 2.The bond lengths in

Results and discussion
Relative contributions of the ZORA terms

Validation
To validate the ZORA implementation in MRChem, we perform total energy calculations of noble gas atoms and a small set of molecules broadly covering the first five rows of the periodic table (H-Xe).In the case of the atoms, we apply three different types of calculations: the atomic solver, LAPWs and MWs.For the diatomic molecules we compare the results obtained using MW and LAPW by means of MRChem and exciting, respectively.
To assess the potential agreement which can be achieved between the various codes, we have performed non-relativistic calculations with all three codes, using point-charge nuclei (numerically smoothed as described in Section in the case of MWs).The results are summarized in Table 3.We have then obtained a root mean square deviation (RMSD) of the relative error between MWs and LAPWs equal to 5.21 In terms of absolute errors, we find that most discrepancies in the total energies are within 10 microHartrees, and only in two cases (AgH and I 2 ) the differences are larger, yet they do not exceed 21 microHartrees.This level of agreement is well below the so called chemical accuracy threshold (1 kcal/mol).
Finally, we find that introducing the relativistic corrections in the multiwavelet formalism as expressed in Eq. 17 leads to a minor increase in MRChem runtimes.Therefore we anticipate that the analysis of the performance (runtimes and parallel scaling) given in Ref. 103 remains valid in the ZORA case.

Conclusions
We have formulated the ZORA method in a form compatible with multiwavelets and implemented it in the MRChem program.The validity and precision of the implementation has been tested against a radial, numerical atomic code and a plane wave code, showing excellent agreement.The current study was done with the specific idea to validate method and theory with a small benchmark: atoms and diatomics, covering broadly the periodic table up to and including 5 th -row elements.The model is also capable of dealing with all parts of the electronic potential self-consistently, with the exception of the HF exchange, which is the only one which cannot be expressed in closed form.

Acknowledgement
The current MRChem implementation is able to include all local contributions of the potential into VZ .It is therefore possible to measure their relative weights for a given atom and how their contribution changes with the nuclear charge.We have performed a series of calculations for the noble gases from helium to xenon, with all seven possibilities: one contribution only, two contributions and all three.The results are summarized in Fig.1.The nuclear potential Vnuc is the largest contribution, as expected.It is followed by the Coulomb term Ĵ, and the exchange and correlation potential Vxc is the smallest one.Moreover the relativistic correction increases roughly with the fourth power of the nuclear charge as expected,5 and the nuclear term becomes progressively more dominant for heavier atoms.

Figure 1 :
Figure 1: Relative contributions to the total energy (in Hartree) for all seven possible ZORA operators, compared to the non-relativistic total energy, for the noble gases He -Xe.Note that the scales on the y-axes are different for each sub-plot.

• 10 − 8 ,
and between MWs and atomic solver equal to 8.07 • 10 −10 .We concluded that the three methods are in very good agreement, setting the mark for what can be expected in the ZORA domain, using a Gaussian nuclear charge distribution as in Equation(24).The ZORA results are summarized in Table 4 showing a RMSD for the relative errors between MWs and LAPWs equal to 3.95 • 10 −8 , and between MWs and atomic solver equal to 6.93 • 10 −9 , in line with what has been observed for the non-relativistic regime, thus confirming the validity of the implementation.
This work was supported by the Norwegian Research Council through a Centre of Excellence grant (Hylleraas Centre 262695), a FRIPRO grant (ReMRChem 324590), by the Tromsø Research Foundation (TFS2016KHH), and by UNINETT Sigma2 through grants of computer time (nn9330k and nn4654k).Jānis Užulis acknowledges funding provided by the project "Strengthening of the capacity of doctoral studies at the University of Latvia within the framework of the new doctoral model", identification No. 8.2.2.0/20/I/006.Andris Gulans acknowledges funding provided by European Regional Development Fund via the Central Finance and Contracting Agency of Republic of Latvia under the grant agreement 1.1.1.5/21/A/004.We would like to thank Dr. Susi Lehtola from Department of Chemistry at the University of Helsinki in Finland for the useful feedback to the original draft of the manuscript.
, we have chosen the Gaussian model for the present work because it is simpler than the more realistic Fermi model and because it is implemented in all GTO codes.The goal of the present work is the validation of the implementation.It is therefore less relevant which model is used, provided that it is consistent throughout the software packages employed.

Table 2 :
LAPW Parameters and structural data used in the calculations of the considered molecules and atoms.R min MT G max is the product of the smallest muffin-tin radius and the largest reciprocal lattice vector.

Table 2 ,
are taken from the NIST 101 database with the exception of SrO and I 2 , which are private communication from not published work.

Table 3 :
Non-relativistic total energies (given in Hartrees) obtained using three different codes and their relative differences.In all cases, the point-like nucleus model is used.

Table 4 :
Scalar-relativistic ZORA total energies (given in Hartrees) obtained using three different codes and their relative differences.In all cases, the smeared nucleus model is used.