Department of Chemistry, UiT The Arctic University of Norway, 9037 Tromsø,
Norway Department of Chemistry, UiT The Arctic University of Norway, 9037
Tromsø, Norway First authorship is shared between AB and SRJ Department of
Chemistry, UiT The Arctic University of Norway, 9037 Tromsø, Norway Department
of Materials Science and NanoEngineering, Rice University, Houston, Texas
77005, United States of America. Department of Chemistry, UiT The Arctic
University of Norway, 9037 Tromsø, Norway Department of Chemistry, UiT The
Arctic University of Norway, 9037 Tromsø, Norway Department of Chemistry, UiT
The Arctic University of Norway, 9037 Tromsø, Norway

AO

    atomic orbital
API

    Application Programmer Interface
AUS

    Advanced User Support
BO

    Born-Oppenheimer
CBS

    complete basis set
CC

    Coupled Cluster
CTCC

    Centre for Theoretical and Computational Chemistry
CoE

    Centre of Excellence
DC

    dielectric continuum
DFT

    density functional theory
DKH

    Douglas-Kroll-Hess
ESC

    Elimination of Small Components
EFP

    effective fragment potential
ECP

    effective core potential
EU

    European Union
FW

    Foldy–Wouthuysen
GGA

    generalized gradient approximation
GPE

    Generalized Poisson Equation
GTO

    Gaussian type orbital
HF

    Hartree-Fock
HPC

    high-performance computing
HC

    Hylleraas Centre for Quantum Molecular Sciences
IEF

    Integral Equation Formalism
IGLO

    individual gauge for localized orbitals
KB

    kinetic balance
KS

    Kohn-Sham
LAO

    London atomic orbital
LAPW

    linearized augmented plane wave
LDA

    local density approximation
MAD

    mean absolute deviation
maxAD

    maximum absolute deviation
MM

    molecular mechanics
MCSCF

    multiconfiguration self consistent field
MPA

    multiphoton absorption
MRA

    multiresolution analysis
MW

    multiwavelet
NAO

    numerical atomic orbital
NeIC

    nordic e-infrastructure collaboration
NMR

    nuclear magnetic resonance
NP

    nanoparticle
OLED

    organic light emitting diode
PAW

    projector augmented wave
PBC

    Periodic Boundary Condition
PCM

    polarizable continuum model
PW

    plane wave
QC

    quantum chemistry
QM/MM

    quantum mechanics/molecular mechanics
QM

    quantum mechanics
RA

    Regular Approximation
RCN

    Research Council of Norway
RMSD

    root mean square deviation
RKB

    restricted kinetic balance
SC

    semiconductor
SCF

    Self Consistent Field
STSM

    short-term scientific mission
SAPT

    symmetry-adapted perturbation theory
SERS

    surface-enhanced raman scattering
WP1

    Work Package 1
WP2

    Work Package 2
WP3

    Work Package 3
WP

    Work Package
X2C

    exact two-component
ZORA

    zero-order regular approximation
ae

    almost everywhere
BVP

    boundary value problem
PDE

    partial differential equation
RDM

    1-body reduced density matrix

# Scalar relativistic effects with Multiwavelets: Implementation and benchmark

Anders Brakestad  Hylleraas Centre for Quantum Molecular Sciences, UiT The
Arctic University of Norway, 9037 Tromsø, Norway  Stig Rune Jensen  Hylleraas
Centre for Quantum Molecular Sciences, UiT The Arctic University of Norway,
9037 Tromsø, Norway  Christian Tantardini  Hylleraas Centre for Quantum
Molecular Sciences, UiT The Arctic University of Norway, 9037 Tromsø, Norway
Quentin Pitteloud  Hylleraas Centre for Quantum Molecular Sciences, UiT The
Arctic University of Norway, 9037 Tromsø, Norway  Peter Wind  Hylleraas Centre
for Quantum Molecular Sciences, UiT The Arctic University of Norway, 9037
Tromsø, Norway  Jānis Užulis  Department of Physics, University of Latvia,
Jelgavas iela 3, Riga, LV-1004, Latvia  Andris Gulans  Department of Physics,
University of Latvia, Jelgavas iela 3, Riga, LV-1004, Latvia  Kathrin Helen
Hopmann  Department of Chemistry, UiT The Arctic University of Norway, 9037
Tromsø, Norway  Luca Frediani  Hylleraas Centre for Quantum Molecular
Sciences, UiT The Arctic University of Norway, 9037 Tromsø, Norway
<EMAIL_ADDRESS>
###### Abstract

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 effect are strongest in the nuclear region,
where the description of electrons through linear combination of atomic
orbitals becomes more challenging. Furthermore, the choice of basis sets for
heavier elements is limited compared to lighter elements where precise basis
sets are available. Thanks to the framework of multiresolution analysis,
multiwavelets provide an appealing alternative to overcome 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 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 (ZORA) framework. We validated our
implementation 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.

##  1 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
solved1. 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 mechanics2, 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 gold4, 5, the
liquid state of mercury6, 5, the functioning of lead-acid batteries7, 5, 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 sharing the same orbital is no longer valid 26,
19, 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$\times$8 matrices instead of a real scalar 28,
29, 30, 31, 32. Although progress has been made to make 4-component
calculations faster and easily available 28, 29, 30, 31, 32, 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 transformation 33, 34. The choice of transformation leads
to different kinds of methods, such as the Pauli Hamiltonian 35, 36, 37, the
Regular Approximations 38, 39, 40, the Douglas-Kroll-Hess (DKH) Hamiltonian
41, 42, 43, 44, 45, or the exact two-component (X2C) Hamiltonian method 46,
47, 48, 49, 50, 51, 52, 19, 53, 54. Further reduction to a scalar method is
also possible for the Pauli Hamiltonian and the RAs 38, 39, 40. 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 55, 35, 56, 57, 38 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) 55, 35, 56, 57, 38. The numerical treatment of this rescaling
function can be challenging, because it displays a cusp at each nucleus55, 35,
56, 57, 38. 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 58, 59.

Some efforts to assess the precision of ZORA for all-electron calculations of
heavier elements have been undertaken using Gaussian type orbitals, but it is
anyway challenging to assess the precision without an external reference 58,
59. 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
system60. 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) 61, the atomic natural orbitals basis sets 62, 63, 64, 65, 66, the X2C
basis sets 67, and the segmented all-electron relativistic contracted (SARC)
basis sets 68, 69, 70, 71, 72. 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 72. 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 59.

In recent years, multiwavelets 73 have emerged as a powerful alternative to
traditional local basis sets 74. Their foundation based on multiresolution
analysis 75 leads to a basis set that is not empirically parameterized. Robust
error control 76, 77, 78 means that the user can set a finite but arbitrary
target precision, and adaptive algorithms 79, 80, 81 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 energies74 and properties82, and have ventured both
towards post-HF methods 83 and relativistic treatments 84.

In 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
85 and linearized augmented plane wave (LAPW) implemented in the electronic-
structure code exciting 86. 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. 85 and 87.

##  2 Theory and implementation

###  2.1 From the Dirac equation to the ZORA Hamiltonian

We will briefly expose how the ZORA 55, 35, 56, 57, 38 is derived starting
from the 4-component Dirac equation of an electron, which describes a
4-component spinor $\Psi$, subject to a potential $V$:

$\left(\beta
c^{2}+V+c\boldsymbol{\alpha}\cdot\boldsymbol{p}\right)\Psi=E\Psi,$ (1)

In the above equation $\boldsymbol{p}$ is the momentum operator, $c$ is the
speed of light, atomic units ($m_{e}=1$, $\hbar=1$, $e=-1$) are assumed and
$\boldsymbol{\alpha}$ and $\beta$ are defined as follows:

$\beta=\begin{pmatrix}I_{2}&0\\\
0&-I_{2}\end{pmatrix},\quad\boldsymbol{\alpha}=\begin{pmatrix}0&\boldsymbol{\sigma}\\\
\boldsymbol{\sigma}&0\end{pmatrix},$

$\sigma_{x}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\quad\sigma_{y}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\quad\sigma_{z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix},$ (2)

.

The first step towards ZORA consists in applying the FW transformation to the
Dirac Hamiltonian of Eq. (1):

$H^{FW}=U^{\dagger}H^{D}U,$ (3)

where the transformation matrix $U=W_{1}W_{2}$ is a product of a decoupling
matrix $W_{1}$ and a renormalization matrix $W_{2}$:

$W_{1}=\begin{pmatrix}1&-R^{\dagger}\\\ R&1\end{pmatrix},\quad
W_{2}=\begin{pmatrix}1/\sqrt{1+R^{\dagger}R}&0\\\
0&1/\sqrt{1+RR^{\dagger}}\end{pmatrix},$ (4)

and $R$ is the exact coupling between the large and the small components of a
4-spinorbits:

$R=\frac{1}{2c^{2}-V+E}c\boldsymbol{\sigma}\cdot\boldsymbol{p}.$ (5)

The inverse potential term depends on the eigenvalue $E$, and this dependence
can be expanded in a Taylor series as follows:

$R=\frac{1}{2c^{2}-V}\left(1+\frac{E}{2c^{2}-V}\right)^{-1}\boldsymbol{\sigma}\cdot\boldsymbol{p}=\frac{1}{2c^{2}-V}\sum_{k=0}^{\infty}\frac{(-E)^{k}}{(2c^{2}-V)^{k}}\boldsymbol{\sigma}\cdot\boldsymbol{p}.$
(6)

.

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:

$H^{RA}=\frac{1}{\sqrt{1+R^{\dagger}R}}\left[V+\boldsymbol{\sigma}\cdot\boldsymbol{p}\frac{c^{2}}{2c^{2}-V}\boldsymbol{\sigma}\cdot\boldsymbol{p}\right]\frac{1}{\sqrt{1+R^{\dagger}R}}.$
(7)

The ZORA Hamiltonian is finally obtained by Taylor-expanding the
renormalization operator and retaining only the zero-order term:

$H^{ZORA}=V+\boldsymbol{\sigma}\cdot\boldsymbol{p}\frac{c^{2}}{2c^{2}-V}\boldsymbol{\sigma}\cdot\boldsymbol{p}.$
(8)

Using the Dirac identity

$(\boldsymbol{\sigma}\cdot\boldsymbol{A})(\boldsymbol{\sigma}\cdot\boldsymbol{B})=\boldsymbol{A}\cdot\boldsymbol{B}+i\boldsymbol{\sigma}\cdot(\boldsymbol{A}\times\boldsymbol{B}),$
(9)

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 scalar-relativistic part and obtain the following
Hamiltonian:

$H^{ZORA}_{SR}=V+\boldsymbol{p}\frac{c^{2}}{2c^{2}-V}\cdot\boldsymbol{p}=V+\frac{1}{2}\boldsymbol{p}\kappa\cdot\boldsymbol{p}$
(10)

where in the last expression we have implicitly defined
$\kappa=(1-\frac{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:

$\left(\frac{1}{2}\boldsymbol{p}\cdot\kappa\boldsymbol{p}+V\right)\varphi_{i}=\epsilon_{i}\varphi_{i}$
(11)

It should be noted that in practical implementations the potential defining
$\kappa$ is usually _not_ the full KS potential. Given the form of $\kappa$,
the most important contribution is the nuclear attraction. Introducing Coulomb
($\hat{J}$) and exchange and correlation ($\hat{V}_{xc}$) 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 88. The ZORA kinetic
operator can then be pre-computed and used throughout the calculation. Atomic-
ZORA has the advantage of being gauge-invariant88.

###  2.2 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 76, 78. The standard KS equations can be
concisely written as follows:

$\hat{F}\varphi_{i}=\sum_{j}F_{ij}\varphi_{j},$ (12)

where $\hat{F}$ is the Fock operator, $\varphi_{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 $\hat{T}$, the nuclear attraction $\hat{V}_{nuc}$, the Coulomb
repulsion $\hat{J}$, the Hartree-Fock exchange $\hat{K}$ scaled by some
numerical factor $\lambda\in[0,1]$, and the exchange and correlation potential
$\hat{V}_{xc}$:

$\hat{F}=\hat{T}+\hat{V}_{nuc}+\hat{J}-\lambda\hat{K}+\hat{V}_{xc}.$ (13)

In the non-relativistic domain, the coupled KS differential equations (12) can
be rewritten in integral form,89, 90 by making use of the bound-state
Helmholtz kernel

$\displaystyle\varphi_{i}$
$\displaystyle=-2\hat{G}_{\mu_{i}}\left(\hat{V}\varphi_{i}-\sum_{j\neq
i}F_{ij}\varphi_{j}\right),$ (14)

where $\hat{G}_{\mu_{i}}=\big{[}-\nabla^{2}-\mu_{i}^{2}\big{]}^{-1}$ is the
integral convolution operator associated with the bound-state Helmholtz kernel
$G(r)=e^{-{\mu_{i}}r}/r$, using $\mu_{i}=\sqrt{-2F_{ii}}$ (i.e. the diagonal
elements of the Fock matrix).

In the ZORA Hamiltonian, the kinetic energy operator becomes:

$\hat{T}=\boldsymbol{p}\cdot\frac{c^{2}}{2c^{2}-V_{Z}}\boldsymbol{p}=-\frac{c^{2}}{2c^{2}-V_{Z}}\nabla^{2}-\nabla\frac{c^{2}}{2c^{2}-V_{Z}}\cdot\nabla=-\frac{1}{2}\kappa\nabla^{2}-\frac{1}{2}\nabla\kappa\cdot\nabla.$
(15)

where $V_{Z}=\hat{V}_{nuc}+\hat{J}+\hat{V}_{xc}$. Including $\hat{J}$ and
$\hat{V}_{xc}$ 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. 91.

Inserting the ZORA kinetic operator into Eq. (12) we obtain:

$\left(-\frac{1}{2}\kappa\nabla^{2}-\frac{1}{2}\nabla\kappa\cdot\nabla+V\right)\varphi_{i}=\sum_{j}F_{ij}\varphi_{j}.$
(16)

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 $\hat{G}=\left(\nabla^{2}+2F_{ii}\right)^{-1}$.
This is achieved first by division by $\kappa$, recalling that
$\kappa^{-1}=1-V_{Z}/2c$. The following integral equation is obtained:

$\varphi_{i}=-2\hat{G}_{\mu_{i}}\left[-\frac{1}{2}\frac{\nabla\kappa}{\kappa}\cdot\nabla\varphi_{i}+\left(\frac{V}{\kappa}+\frac{V_{Z}}{2c^{2}}F_{ii}\right)\varphi_{i}-\frac{1}{\kappa}\sum_{j\neq
i}F_{ij}\varphi_{j}\right].$ (17)

When $c\rightarrow\infty$, $\kappa\rightarrow 1$ and $\nabla\kappa\rightarrow
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 $\kappa\simeq 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\rightarrow 0$. To overcome the
convergence issues for heavier nuclei (5th row of the periodic table) we have
therefore introduced a finite nucleus model, as described in Section 2.5. For
elements beyond the 5th row, one has to keep in mind that ZORA becomes
questionable92.

An alternative approach to obtain the desired
$\left(\nabla^{2}+2F_{ii}\right)$ term from Eq. (16) is to add a Laplacian
term ($\frac{1}{2}\nabla^{2}$) directly, thus avoiding division by $\kappa$.
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.

###  2.3 Implementation with MW

An implementation of the ZORA method as outlined above, is currently in a
development version of the MRChem package93, 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 $\kappa$ function is computed as a point-wise map of the chosen ZORA
potential $V_{Z}$ through

$\kappa(r)=\frac{1}{1-V_{Z}(r)/2c^{2}},$ (18)

and similarly for its inverse

$\kappa^{-1}(r)=1-V_{Z}(r)/2c^{2}.$ (19)

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).

###  2.4 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 approximation76:

$u(r)=\frac{\text{erf}(r)}{r}+\frac{1}{3\sqrt{\pi}}\left(e^{-r^{2}}+16e^{-4r^{2}}\right),$
(20)

and then parameterized as $u(r/s)/s$, where $s$ is a scalar smoothing
parameter. The smoothing parameter is chosen from the nuclear charge $Z$ and
the desired precision $\epsilon$ as 76

$s=\left(\frac{0.00435\epsilon}{Z^{5}}\right)^{1/3}.$ (21)

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 94, 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.

###  2.5 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.
94. 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

$\rho^{G}(R)=\rho_{0}^{G}e^{-\xi R^{2}}$ (22)

with the normalisation prefactor
$\rho_{0}^{G}=eZ\Bigg{(}\frac{\xi}{\pi}\Bigg{)}^{3/2}$ and the parameter $\xi$
related to the the root-mean-square radius of the nucleus via the expression
$\xi=\frac{3}{\langle R^{2}\rangle}$. To determine $\langle R^{2}\rangle$, we
apply the empirical formula for Ref. 95:

$\sqrt{\langle R^{2}\rangle}=(0.836A^{1/3}+0.570)fm$ (23)

where $A$ is the nuclear mass number and $fm$ 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 $\sqrt{\langle R^{2}\rangle}$
provided in Ref. 94.

Although the potential for the Gaussian charge distribution cannot be
represented in analytical form, a sufficiently close representation for such
narrow Gaussians has been introduced by Visscher et al. 94:

$V^{G}(R)=-\frac{eZ}{R}\mathrm{erf}\big{(}\sqrt{\xi}R)\big{)}.$ (24)

###  2.6 Methods for verification

To verify correctness of the ZORA implementation in MRChem, we use two other
codes: a numerical atomic solver 85 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 $\varphi_{n\ell
m}(\mathbf{r})=u_{n\ell}(r)Y_{\ell m}(\hat{r})$, where $u_{n\ell}(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:

$\varphi_{n\ell
m}(\mathbf{r})=-2\hat{G}_{\epsilon_{n\ell}}\left(\hat{V}\varphi_{n\ell
m}(\mathbf{r})\right).$ (25)

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
$\epsilon_{n\ell}$ and setting $F_{ij}=0$ for $i\neq 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)96, 97 and (ii) local-orbitals (LOs)98. 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. 87, 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. 2.5 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:

$E_{\mathrm{M}}=\int\rho_{\mathrm{n}}(\mathbf{r})V_{\mathrm{C}}(\mathbf{r})d\mathbf{r}.$
(26)

If a nucleus at the site $\mathbf{R}_{0}$ is defined as a point charge, its
density distribution is
$\rho_{\mathrm{n}}(\mathbf{R})=Z\delta(\mathrm{R}-\mathbf{R}_{0})$ leading to
$E_{\mathrm{M}}=ZV_{\mathrm{C}}(\mathbf{R}_{0})$. For smeared nuclei with
$\rho_{\mathrm{n}}(\mathbf{R})=\rho^{G}(\mathbf{R})$ (see Eq. 24), Eq. 26 is
evaluated as an integral on a radial grid.

##  3 Computational details

All calculations were performed with the PBE functional 99 employing the
parameters $\beta=0.06672455060314922$ and $\mu=0.066725\frac{\pi^{2}}{3}$.

For all MW calculations, a development version of MRChem has been employed.
Interpolating polynomials of 9th order were used with a simulation box size of
$128\,\mathrm{a_{0}}$. The numerical precision thresholds used are summarized
in Table 1.

Table 1: MRChem precision parameters. 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\,\mathrm{a_{0}}$ for atoms and $30\,\mathrm{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. 82. 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. The truncation of the Coulomb interaction
makes the potential unique and thus removes this ambiguity completely. The
canonical orbitals were expressed in terms of local orbitals and LAPWs with
the cutoff $R_{\mathrm{MT}}K_{\mathrm{max}}$ sufficient to ensure a few
$\mu\mathrm{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.

Table 2: LAPW Parameters and structural data used in the calculations of the considered molecules and atoms. $R_{\mathrm{MT}}^{\mathrm{min}}G_{\mathrm{max}}$ is the product of the smallest muffin-tin radius and the largest reciprocal lattice vector.  Material | $R_{MT}^{min}G_{max}$ | $R_{MT}[a_{0}]$ | Bond length / Å
---|---|---|---
CaO | 11 | 1.63 / 1.40 | 1.8221
CuH | 10 | 1.40 / 1.00 | 1.4626
SrO | 13 | 1.63 / 1.40 | 1.9050
Cu2 | 14 | 1.8 | 2.2197
AgH | 11 | 1.68 / 1.20 | 1.6180
I2 | 14 | 1.6 | 2.6630
He | 11 | 2 |
Ne | 12 | 2 |
Ar | 13 | 2 |
Kr | 13 | 2 |
Xe | 15 | 2 |

The calculations with the atomic solver employed a 5th order polynomial on a
radial grid with the innermost and outermost points
$r_{\mathrm{min}}=10^{-8}\,\mathrm{a_{0}}$ and
$r_{\mathrm{max}}=35\,\mathrm{a_{0}}$. The number of radial points was set to
5000, which is fully sufficient to guarantee a sub-$\mu\mathrm{Ha}$ precision.

##  4 Results and discussion

###  4.1 Relative contributions of the ZORA terms

The current MRChem implementation is able to include all local contributions
of the potential into $\hat{V}_{Z}$. 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 $\hat{V}_{nuc}$ is the largest contribution, as expected. It is
followed by the Coulomb term $\hat{J}$, and the exchange and correlation
potential $\hat{V}_{xc}$ is the smallest one. Moreover the relativistic
correction increases roughly with the fourth power of the nuclear charge as
expected5, and the nuclear term becomes progressively more dominant for
heavier atoms.

(a) Helium

(b) Neon

(c) Argon

(d) Krypton

(e) Xenon

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.

###  4.2 Validation

To validate the ZORA implemetation 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 2.4 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\cdot 10^{-8}$, and between MWs and atomic solver equal to
$8.07\cdot 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\cdot 10^{-8}$, and between MWs and atomic solver equal to
$6.93\cdot 10^{-9}$, in line with what has been observed for the non-
relativistic regime, thus confirming the validity of the implementation. 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 I2) 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).

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. |  |  |  | Relative difference
---|---|---|---|---
Species | AtomSolver | exciting | MRChem |  MRChem vs |  MRChem vs
|  |  |  | AtomSolver | exciting
He | -2.892935 | -2.892935 | -2.892935 | -4.1e-10 | 4.4e-09
Ne | -128.866434 | -128.866434 | -128.866433 | -3.1e-10 | 1.9e-09
Ar | -527.346141 | -527.346138 | -527.346140 | -1.7e-09 | 3.7e-09
Kr | -2753.416138 | -2753.416137 | -2753.416138 | -1.5e-10 | 2.1e-10
Xe | -7234.233259 | -7234.233259 | -7234.233259 | -3.4e-11 | 9.7e-11
CaO | - | -752.562050 | -752.562058 | - | 1.0e-08
CuH | - | -1640.901717 | -1640.901727 | - | 6.1e-09
SrO | - | -3208.161712 | -3208.161722 | - | 3.1e-09
Cu2 | - | -3280.675840 | -3280.675846 | - | 1.7e-09
AgH | - | -5200.162245 | -5200.162259 | - | 2.8e-09
I2 | - | -13840.158307 | -13840.158328 | - | 1.5e-09
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. |  |  |  | Relative difference
---|---|---|---|---
Species | AtomSolver | exciting | MRChem |  MRChem vs |  MRChem vs
|  |  |  | AtomSolver | exciting
He | -2.893119 | -2.893119 | -2.893119 | -2.2e-09 | -1.3e-07
Ne | -129.089335 | -129.089335 | -129.089335 | -2.1e-09 | -3.4e-10
Ar | -530.224039 | -530.224036 | -530.224038 | -2.6e-09 | 3.3e-09
Kr | -2810.049764 | -2810.049762 | -2810.049755 | -3.1e-09 | -2.5e-09
Xe | -7564.596665 | -7564.596662 | -7564.596554 | -1.5e-08 | -1.4e-08
CaO | - | -757.195492 | -757.195491 | - | -1.2e-09
CuH | - | -1663.238493 | -1663.238504 | - | 6.4e-09
SrO | - | -3279.826587 | -3279.826596 | - | 2.9e-09
Cu2 | - | -3325.342851 | -3325.342843 | - | -2.3e-09
AgH | - | -5380.042000 | -5380.042003 | - | 6.1e-10
I2 | - | -14448.747787 | -14448.747761 | - | -1.8e-09

##  5 Conclusions

We have formulated the ZORA method in a form compatible with multiwavelets and
implemented in the MRChem program. This implementation has been used to
compute total energies for small set of atoms and diatomic molecules to cover
different blocks of the periodic table of elements up to the $5^{th}$ row. The
implementation was validated by comparing the results from the MRChem
implementation against both a radial numerical code for closed-shell atoms and
exciting for atoms and small diatomics. The precision obtained with the ZORA
implementation is similar to the one obtained with the non-relativistic one.

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.

## References

  * Dirac 1933 Dirac, P. A. Theory of electrons and positrons. _Nobel Lecture_ 1933, _12_ , 320–325
  * Dirac 1928 Dirac, P. A. M. The Quantum Theory of the Electron. _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ 1928, _117_ , 610–624
  * Simões 2002 Simões, A. Dirac’s Claim and the Chemists. _Physics in Perspective_ 2002, _4_ , 253–266
  * Bartlett 1998 Bartlett, N. Relativistic effects and the chemistry of gold. _Gold Bulletin_ 1998, _31_ , 22–25
  * Pyykkö 2012 Pyykkö, P. Relativistic Effects in Chemistry: More Common Than You Thought. _Annual Review of Physical Chemistry_ 2012, _63_ , 45–64
  * Keeton and Loucks 1966 Keeton, S.; Loucks, T. Electronic structure of mercury. _Physical Review_ 1966, _152_ , 548
  * Ahuja et al. 2011 Ahuja, R.; Blomqvist, A.; Larsson, P.; Pyykkö, P.; Zaleski-Ejgierd, P. Relativity and the lead-acid battery. _Physical review letters_ 2011, _106_ , 018301
  * 8 Demissie, T. B.; Garabato, B. D.; Ruud, K.; Kozlowski, P. M. Mercury Methylation by Cobalt Corrinoids: Relativistic Effects Dictate the Reaction Mechanism. _Angewandte Chemie International Edition_ _55_ , 11503–11506
  * Vícha et al. 2016 Vícha, J.; Marek, R.; Straka, M. High-Frequency 13C and 29Si NMR Chemical Shifts in Diamagnetic Low-Valence Compounds of TlI and PbII: Decisive Role of Relativistic Effects. _Inorganic Chemistry_ 2016, _55_ , 1770–1781
  * Lein et al. 2010 Lein, M.; Rudolph, M.; Hashmi, S. K.; Schwerdtfeger, P. Homogeneous Gold Catalysis: Mechanism and Relativistic Effects of the Addition of Water to Propyne. _Organometallics_ 2010, _29_ , 2206–2210
  * Boualili et al. 2021 Boualili, F. Z.; Nemouchi, M.; Godefroid, M.; Jönsson, P. Weak Correlation and Strong Relativistic Effects on the Hyperfine Interaction in Fluorine. 2021
  * Mattsson et al. 2019 Mattsson, S.; Paulus, B.; Redeker, F. A.; Beckers, H.; Riedel, S.; Müller, C. The Crystal Structure of $\alpha$-F2: Solving a 50 Year Old Puzzle Computationally. _Chemistry – A European Journal_ 2019, _25_ , 3318–3324
  * Tantardini et al. 2022 Tantardini, C.; Jalolov, F. N.; Kvashnin, A. G. Crystal Structure Evolution of Fluorine under High Pressure. _The Journal of Physical Chemistry C_ 2022, _126_ , 11358–11364
  * Hellmann 1935 Hellmann, H. A new approximation method in the problem of many electrons. _The Journal of Chemical Physics_ 1935, _3_ , 61–61
  * Gombás 1935 Gombás, P. Über die metallische Bindung. _Zeitschrift für Physik_ 1935, _94_ , 473–488
  * Xu and Truhlar 2011 Xu, X.; Truhlar, D. G. Accuracy of Effective Core Potentials and Basis Sets for Density Functional Calculations, Including Relativistic Effects, As Illustrated by Calculations on Arsenic Compounds. _Journal of Chemical Theory and Computation_ 2011, _7_ , 2766–2779
  * Stoll et al. 2002 Stoll, H.; Metz, B.; Dolg, M. Relativistic energy-consistent pseudopotentials—Recent developments. _Journal of computational chemistry_ 2002, _23_ , 767–778
  * Heine 1970 Heine, V. In _The Pseudopotential Concept_ ; Ehrenreich, H., Seitz, F., Turnbull, D., Eds.; Solid State Physics; Academic Press, 1970; Vol. 24; pp 1–36
  * Saue 2011 Saue, T. Relativistic Hamiltonians for chemistry: A primer. _ChemPhysChem_ 2011, _12_ , 3077–3094
  * Rudek et al. 2018 Rudek, B.; Toyota, K.; Foucar, L.; Erk, B.; Boll, R.; Bomme, C.; Correa, J.; Carron, S.; Boutet, S.; Williams, G. J.; Ferguson, K. R.; Alonso-Mori, R.; Koglin, J. E.; Gorkhover, T.; Bucher, M.; Lehmann, C. S.; Krässig, B.; Southworth, S. H.; Young, L.; Bostedt, C.; Ueda, K.; Marchenko, T.; Simon, M.; Jurek, Z.; Santra, R.; Rudenko, A.; Son, S.-K.; Rolles, D. Relativistic and resonant effects in the ionization of heavy atoms by ultra-intense hard X-rays. _Nature Communications_ 2018, _9_
  * Breit 1928 Breit, G. An Interpretation of Dirac’s Theory of the Electron. _Proceedings of the National Academy of Sciences_ 1928, _14_ , 553–559
  * Breit 1932 Breit, G. Dirac’s Equation and the Spin-Spin Interactions of Two Electrons. _Phys. Rev._ 1932, _39_ , 616–624
  * Moss 2012 Moss, R. _Advanced molecular quantum mechanics: an introduction to relativistic quantum mechanics and the quantum theory of radiation_ ; Springer Science & Business Media, 2012
  * Dyall and Fægri Jr 2007 Dyall, K. G.; Fægri Jr, K. _Introduction to relativistic quantum chemistry_ ; Oxford University Press, 2007
  * Helgaker et al. 2012 Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Recent Advances in Wave Function-Based Methods of Molecular-Property Calculations. _Chemical Reviews_ 2012, _112_ , 543–631
  * Jacob and Reiher 2012 Jacob, C. R.; Reiher, M. Spin in density-functional theory. _International Journal of Quantum Chemistry_ 2012, _112_ , 3661–3684
  * Marian 2012 Marian, C. M. Spin–orbit coupling and intersystem crossing in molecules. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2012, _2_ , 187–203
  * Repisky et al. 2020 Repisky, M.; Komorovsky, S.; Kadek, M.; Konecny, L.; Ekström, U.; Malkin, E.; Kaupp, M.; Ruud, K.; Malkina, O. L.; Malkin, V. G. ReSpect: Relativistic spectroscopy DFT program package. _The Journal of Chemical Physics_ 2020, _152_
  * Repisky et al. 2020 Repisky, M.; Komorovsky, S.; Kadek, M.; Konecny, L.; Ekström, U.; Malkin, E.; Kaupp, M.; Ruud, K.; Malkina, O. L.; Malkin, V. G. ReSpect: Relativistic spectroscopy DFT program package. _The Journal of Chemical Physics_ 2020, _152_ , 184101
  * Saue et al. 2020 Saue, T.; Bast, R.; Gomes, A. S. P.; Jensen, H. J. A.; Visscher, L.; Aucar, I. A.; Di Remigio, R.; Dyall, K. G.; Eliav, E.; Fasshauer, E.; Fleig, T.; Halbert, L.; Hedegård, E. D.; Helmich-Paris, B.; Iliaš, M.; Jacob, C. R.; Knecht, S.; Laerdahl, J. K.; Vidal, M. L.; Nayak, M. K.; Olejniczak, M.; Olsen, J. M. H.; Pernpointner, M.; Senjean, B.; Shee, A.; Sunaga, A.; van Stralen, J. N. P. The DIRAC code for relativistic molecular calculations. _The Journal of Chemical Physics_ 2020, _152_ , 204104
  * Zhang et al. 2020 Zhang, Y.; Suo, B.; Wang, Z.; Zhang, N.; Li, Z.; Lei, Y.; Zou, W.; Gao, J.; Peng, D.; Pu, Z.; Xiao, Y.; Sun, Q.; Wang, F.; Ma, Y.; Wang, X.; Guo, Y.; Liu, W. BDF: A relativistic electronic structure program package. _The Journal of Chemical Physics_ 2020, _152_ , 064113
  * Belpassi et al. 2020 Belpassi, L.; De Santis, M.; Quiney, H. M.; Tarantelli, F.; Storchi, L. BERTHA: Implementation of a four-component Dirac–Kohn–Sham relativistic framework. _The Journal of Chemical Physics_ 2020, _152_ , 164118
  * Foldy and Wouthuysen 1950 Foldy, L. L.; Wouthuysen, S. A. On the Dirac Theory of Spin 1/2 Particles and Its Non-Relativistic Limit. _Phys. Rev._ 1950, _78_ , 29–36
  * Foldy 1952 Foldy, L. L. The Electromagnetic Properties of Dirac Particles. _Phys. Rev._ 1952, _87_ , 688–693
  * Chang et al. 1986 Chang, C.; Pelissier, M.; Durand, P. Regular two-component Pauli-like effective Hamiltonians in Dirac theory. _Physica Scripta_ 1986, _34_ , 394
  * Kang and Uhlmann 2004 Kang, H.; Uhlmann, G. Inverse problems for the Pauli Hamiltonian in two dimensions. _Journal of Fourier Analysis and Applications_ 2004, _10_ , 201–215
  * Goldman 1977 Goldman, T. Gauge invariance, time-dependent Foldy-Wouthuysen transformations, and the Pauli Hamiltonian. _Physical Review D_ 1977, _15_ , 1063
  * van Lenthe et al. 1994 van Lenthe, E.; Baerends, E.-J.; Snijders, J. G. Relativistic total energy using regular approximations. _The Journal of chemical physics_ 1994, _101_ , 9783–9792
  * Dyall and van Lenthe 1999 Dyall, K. G.; van Lenthe, E. Relativistic regular approximations revisited: An infinite-order relativistic approximation. _The Journal of chemical physics_ 1999, _111_ , 1366–1372
  * Mohri and Nederhof 2001 Mohri, M.; Nederhof, M.-J. _Robustness in language and speech technology_ ; Springer, 2001; pp 153–163
  * Douglas and Kroll 1974 Douglas, M.; Kroll, N. M. Quantum electrodynamical corrections to the fine structure of helium. _Annals of Physics_ 1974, _82_ , 89–155
  * Hess 1985 Hess, B. A. Applicability of the no-pair equation with free-particle projection operators to atomic and molecular structure calculations. _Physical Review A_ 1985, _32_ , 756
  * Hess 1986 Hess, B. A. Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators. _Physical Review A_ 1986, _33_ , 3742
  * Jansen and Heß 1989 Jansen, G.; Heß, B. A. Revision of the Douglas-Kroll transformation. _Physical Review A_ 1989, _39_ , 6016
  * Reiher 2006 Reiher, M. Douglas–Kroll–Hess Theory: a relativistic electrons-only theory for chemistry. _Theoretical Chemistry Accounts_ 2006, _116_ , 241–252
  * Kutzelnigg and Liu 2005 Kutzelnigg, W.; Liu, W. Quasirelativistic theory equivalent to fully relativistic theory. _The Journal of chemical physics_ 2005, _123_
  * Liu and Peng 2006 Liu, W.; Peng, D. Infinite-order quasirelativistic density functional method based on the exact matrix quasirelativistic theory. _The Journal of chemical physics_ 2006, _125_
  * Liu and Peng 2009 Liu, W.; Peng, D. Exact two-component Hamiltonians revisited. _The Journal of chemical physics_ 2009, _131_
  * Iliaš and Saue 2007 Iliaš, M.; Saue, T. An infinite-order two-component relativistic Hamiltonian by a simple one-step transformation. _The Journal of chemical physics_ 2007, _126_
  * Peng et al. 2007 Peng, D.; Liu, W.; Xiao, Y.; Cheng, L. Making four-and two-component relativistic density functional methods fully equivalent based on the idea of “from atoms to molecule”. _The Journal of chemical physics_ 2007, _127_
  * Peng et al. 2013 Peng, D.; Middendorf, N.; Weigend, F.; Reiher, M. An efficient implementation of two-component relativistic exact-decoupling methods for large molecules. _The Journal of chemical physics_ 2013, _138_
  * Liu 2010 Liu, W. Ideas of relativistic quantum chemistry. _Molecular Physics_ 2010, _108_ , 1679–1706
  * Liu 2014 Liu, W. Advances in relativistic molecular quantum mechanics. _Physics Reports_ 2014, _537_ , 59–89
  * Liu 2016 Liu, W. Big picture of relativistic molecular quantum mechanics. _National Science Review_ 2016, _3_ , 204–221
  * Filatov and Cremer 2003 Filatov, M.; Cremer, D. On the physical meaning of the ZORA Hamiltonian. _Molecular Physics_ 2003, _101_ , 2295–2302
  * Heully et al. 1986 Heully, J.-L.; Lindgren, I.; Lindroth, E.; Lundqvist, S.; Martensson-Pendrill, A.-M. Diagonalisation of the Dirac Hamiltonian as a basis for a relativistic many-body procedure. _Journal of Physics B: Atomic and Molecular Physics_ 1986, _19_ , 2799
  * Lenthe et al. 1993 Lenthe, E. v.; Baerends, E.-J.; Snijders, J. G. Relativistic regular two-component Hamiltonians. _The Journal of chemical physics_ 1993, _99_ , 4597–4610
  * Pantazis and Neese 2014 Pantazis, D. A.; Neese, F. All-electron basis sets for heavy elements. _WIREs Computational Molecular Science_ 2014, _4_ , 363–374
  * Güell et al. 2008 Güell, M.; Luis, J. M.; Sola, M.; Swart, M. Importance of the basis set for the spin-state energetics of iron complexes. _J. Phys. Chem. A_ 2008, _112_ , 6384–6391
  * Leeuwen et al. 1994 Leeuwen, R. v.; Lenthe, E. v.; Baerends, E. J.; Snijders, J. G. Exact solutions of regular approximate relativistic wave equations for hydrogen-like atoms. _The Journal of Chemical Physics_ 1994, _101_ , 1272–1281
  * de Castro and Jorge 1998 de Castro, E. V. R.; Jorge, F. E. Accurate universal Gaussian basis set for all atoms of the Periodic Table. _The Journal of Chemical Physics_ 1998, _108_ , 5225–5229
  * Zobel et al. 2020 Zobel, J. P.; Widmark, P.-O.; Veryazov, V. The ANO-R Basis Set. _Journal of Chemical Theory and Computation_ 2020, _16_ , 278–294
  * Zobel et al. 2021 Zobel, J. P.; Widmark, P.-O.; Veryazov, V. Correction to “The ANO-R Basis Set”. _Journal of Chemical Theory and Computation_ 2021, _17_ , 3233–3234
  * Roos et al. 2004 Roos, B.; Veryazov, V.; Widmark, P.-O. Relativistic atomic natural orbital type basis sets for the alkaline and alkaline-earth atoms applied to the ground-state potentials for the corresponding dimers. _Theoretical Chemistry Accounts_ 2004, _111_ , 345–351
  * Roos et al. 2004 Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. Main Group Atoms and Dimers Studied with a New Relativistic ANO Basis Set. _The Journal of Physical Chemistry A_ 2004, _108_ , 2851–2858
  * Roos et al. 2005 Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. New Relativistic ANO Basis Sets for Transition Metal Atoms. _The Journal of Physical Chemistry A_ 2005, _109_ , 6575–6579
  * Pollak and Weigend 2017 Pollak, P.; Weigend, F. Segmented Contracted Error-Consistent Basis Sets of Double- and Triple-$\zeta$ Valence Quality for One- and Two-Component Relativistic All-Electron Calculations. _Journal of Chemical Theory and Computation_ 2017, _13_ , 3696–3705
  * Pantazis and Neese 2009 Pantazis, D. A.; Neese, F. All-electron scalar relativistic basis sets for the lanthanides. _Journal of chemical theory and computation_ 2009, _5_ , 2229–2238
  * Pantazis and Neese 2011 Pantazis, D. A.; Neese, F. All-electron scalar relativistic basis sets for the actinides. _Journal of chemical theory and computation_ 2011, _7_ , 677–684
  * Pantazis and Neese 2012 Pantazis, D. A.; Neese, F. All-electron scalar relativistic basis sets for the 6 p elements. _Theoretical Chemistry Accounts_ 2012, _131_ , 1–7
  * Aravena et al. 2016 Aravena, D.; Neese, F.; Pantazis, D. A. Improved segmented all-electron relativistically contracted basis sets for the lanthanides. _Journal of chemical theory and computation_ 2016, _12_ , 1148–1156
  * Rolfes et al. 2020 Rolfes, J. D.; Neese, F.; Pantazis, D. A. All-electron scalar relativistic basis sets for the elements Rb–Xe. _Journal of Computational Chemistry_ 2020, _41_ , 1842–1849
  * Alpert 1993 Alpert, B. K. A Class of Bases in $L^{2}$ for the Sparse Representation of Integral Operators. _SIAM Journal on Mathematical Analysis_ 1993, _24_ , 246–262
  * Jensen et al. 2017 Jensen, S. R.; Saha, S.; Flores-Livas, J. A.; Huhn, W.; Blum, V.; Goedecker, S.; Frediani, L. The elephant in the room of density functional theory calculations. _The journal of physical chemistry letters_ 2017, _8_ , 1449–1457
  * Keinert 2009 Keinert, F. In _Encyclopedia of Complexity and Systems Science_ ; Meyers, R. A., Ed.; Springer New York, 2009; pp 5841–5858
  * Harrison et al. 2004 Harrison, R. J.; Fann, G. I.; Yanai, T.; Gan, Z.; Beylkin, G. Multiresolution quantum chemistry: Basic theory and initial applications. _The Journal of Chemical Physics_ 2004, _121_ , 11587–11598
  * Kato et al. 2012 Kato, T.; Yokoi, Y.; Sekino, H. Basis set limit computation of dynamic polarizability at near-resonance region. _International Journal of Quantum Chemistry_ 2012, _113_ , 286–289
  * Jensen et al. 2017 Jensen, S. R.; Saha, S.; Flores-Livas, J. A.; Huhn, W.; Blum, V.; Goedecker, S.; Frediani, L. The Elephant in the Room of Density Functional Theory Calculations. _The Journal of Physical Chemistry Letters_ 2017, _8_ , 1449–1457
  * Alpert et al. 2002 Alpert, B.; Beylkin, G.; Gines, D.; Vozovoi, L. Adaptive Solution of Partial Differential Equations in Multiwavelet Bases. _Journal of Computational Physics_ 2002, _182_ , 149–190
  * Beylkin et al. 2008 Beylkin, G.; Cheruvu, V.; Pérez, F. Fast adaptive algorithms in the non-standard form for multidimensional problems. _Applied and Computational Harmonic Analysis_ 2008, _24_ , 354–377
  * Frediani et al. 2013 Frediani, L.; Fossgaard, E.; Flå, T.; Ruud, K. Fully adaptive algorithms for multivariate integral equations using the non-standard form and multiwavelets with applications to the Poisson and bound-state Helmholtz kernels in three dimensions. _Molecular Physics_ 2013, _111_ , 1143–1160
  * Jensen et al. 2016 Jensen, S. R.; Flå, T.; Jonsson, D.; Monstad, R. S.; Ruud, K.; Frediani, L. Magnetic properties with multiwavelets and DFT: the complete basis set limit achieved. _Physical Chemistry Chemical Physics_ 2016, _18_ , 21145–21161
  * Bischoff et al. 2012 Bischoff, F. A.; Harrison, R. J.; Valeev, E. F. Computing many-body wave functions with guaranteed precision: The first-order Møller-Plesset wave function for the ground state of helium atom. _The Journal of Chemical Physics_ 2012, _137_ , 104103
  * Anderson et al. 2019 Anderson, J.; Sundahl, B.; Harrison, R.; Beylkin, G. Dirac-Fock calculations on molecules in an adaptive multiwavelet basis. _The Journal of Chemical Physics_ 2019, _151_ , 234112
  * Užulis and Gulans 2022 Užulis, J.; Gulans, A. Radial Kohn–Sham problem via integral-equation approach. _Journal of Physics Communications_ 2022, _6_ , 085002
  * Gulans et al. 2014 Gulans, A.; Kontur, S.; Meisenbichler, C.; Nabok, D.; Pavone, P.; Rigamonti, S.; Sagmeister, S.; Werner, U.; Draxl, C. exciting: a full-potential all-electron package implementing density-functional theory and many-body perturbation theory. _Journal of Physics: Condensed Matter_ 2014, _26_ , 363202
  * Gulans et al. 2018 Gulans, A.; Kozhevnikov, A.; Draxl, C. Microhartree precision in density functional theory calculations. _Phys. Rev. B_ 2018, _97_ , 161105
  * Knuth et al. 2015 Knuth, F.; Carbogno, C.; Atalla, V.; Blum, V.; Scheffler, M. All-electron formalism for total energy strain derivatives and stress tensor components for numeric atom-centered orbitals. _Computer Physics Communications_ 2015, _190_ , 33–50
  * Kalos 1962 Kalos, M. H. Monte Carlo Calculations of the Ground State of Three- and Four-Body Nuclei. _Phys. Rev._ 1962, _128_ , 1791–1795
  * Beylkin and Mohlenkamp 2002 Beylkin, G.; Mohlenkamp, M. J. Numerical operator calculus in higher dimensions. _Proceedings of the National Academy of Sciences_ 2002, _99_ , 10246–10251
  * Faas et al. 2000 Faas, S.; van Lenthe, J. H.; Hennum, A. C.; Snijders, J. G. An ab initio two-component relativistic method including spin–orbit coupling using the regular approximation. _The Journal of Chemical Physics_ 2000, _113_ , 4052–4059
  * Visscher and Saue 2000 Visscher, L.; Saue, T. Approximate relativistic electronic structure methods based on the quaternion modified Dirac equation. _The Journal of Chemical Physics_ 2000, _113_ , 3996–4002
  * 93 MRChem program package 2020. DOI: 10.5281/zenodo.3786996
  * Visscher and Dyall 1997 Visscher, L.; Dyall, K. DIRAC–FOCK ATOMIC ELECTRONIC STRUCTURE CALCULATIONS USING DIFFERENT NUCLEAR CHARGE DISTRIBUTIONS. _Atomic Data and Nuclear Data Tables_ 1997, _67_ , 207–224
  * Johnson and Soff 1985 Johnson, W.; Soff, G. The lamb shift in hydrogen-like atoms, 1 $\leqslant$ Z $\leqslant$ 110\. _Atomic Data and Nuclear Data Tables_ 1985, _33_ , 405–446
  * Slater 1937 Slater, J. C. Wave functions in a periodic potential. _Physical Review_ 1937, _51_ , 846
  * Andersen 1975 Andersen, O. K. Linear methods in band theory. _Physical Review B_ 1975, _12_ , 3060
  * Sjöstedt et al. 2000 Sjöstedt, E.; Nordström, L.; Singh, D. J. An alternative way of linearizing the augmented plane-wave method. _Solid State Communications_ 2000, _114_ , 15–20
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Physical Review Letters_ 1996, _77_ , 3865–3868

  *[ZORA]: zero-order regular approximation
  *[ECP]: effective core potential
  *[NMR]: nuclear magnetic resonance
  *[DKH]: Douglas-Kroll-Hess
  *[X2C]: exact two-component
  *[RAs]: Regular Approximation
  *[FWs]: Foldy–Wouthuysen
  *[GTO]: Gaussian type orbital
  *[MWs]: multiwavelet
  *[HF]: Hartree-Fock
  *[DFT]: density functional theory
  *[MW]: multiwavelet
  *[LAPW]: linearized augmented plane wave
  *[CBS]: complete basis set
  *[FW]: Foldy–Wouthuysen
  *[RA]: Regular Approximation
  *[KS]: Kohn-Sham
  *[SCF]: Self Consistent Field
  *[LAPWs]: linearized augmented plane wave
  *[RMSD]: root mean square deviation