(sec:relativistic.detailed)=
# Relativistic Options
The relativistic methods in ORCA are implemented in a fairly
straightforward way but do require some caution from the user. The
options are controlled through the `%rel` block which features the
following variables:
```orca
%rel
#----------------------------------------------------
# Basic scalar relativistic method
#----------------------------------------------------
method DKH # Douglas-Kroll-Hess
ZORA # ZORA (numerical integration)
IORA # IORA (numerical integration)
IORAmm # IORA with van Wuellens
# modified metric
X2C # Exact two-component
# ---------------------------------------------------
# Choice of the model potential for ALL methods
# ---------------------------------------------------
ModelPot VeN, VC, VXa, VLDA, VPC
# Flags for terms in the model potential
# =0 not included =1 included
# WARNING: default is currently 1,1,1,1 for ZORA and IORA and
# VeN = nuclear attraction term
# VC = model Coulomb potential (ZORA/IORA only)
# VXa = model Xalpha potential (ZORA/IORA only)
# VLDA= VWN-5 local correlation model pot. (ZORA/IORA only)
# VPC = external point charges (X2C only)
Xalpha 0.7 # default value for the X-Alpha potential,
# only has an effect when VXa is part of the model potential
# --------------------------------------------------
# This variable determines the type of fitted atomic
# density that enters the Coulomb potential part of the
# model potential (has no effect when using DKH):
# --------------------------------------------------
ModelDens rhoDKH # DKH4 model densities (default)
rhoZORA # ZORA model densities
rhoHF # Hartree-Fock model densities
# --------------------------------------------------
# This flag controls whether only one center terms
# are retained. If this is true an approximate treat-
# ment of relativistic effects is the result, but
# geometry optimizations CAN BE PERFORMED WITH ALL
# METHODS AND MODEL POTENTIALS
# In addition one gets NO gauge noninvariance
# errors in ZORA or IORA
# --------------------------------------------------
OneCenter false # default value
# --------------------------------------------------
# Flag for the diagonal approximation to the unitary
# decoupling matrix (DLU) in X2C. Mutually exclusive
# with OneCenter. See section "DLU approximation".
# --------------------------------------------------
DLU false # default value
# --------------------------------------------------
# Specify the speed of light used in relativistic
# calculations
# --------------------------------------------------
C 137.0359895 # speed of light used (137.0359895 is the default value)
# synonyms for C are VELIT, VELOCITY
# --------------------------------------------------
# Picture change for properties
# ---------------------------------------------------
PictureChange 0 # (or false): no picturechange (default)
1 # (or true): include picturechange
2 # for DKH: use second-order DKH transformation
# (see section "Picture-Change Effects")
2 # for X2C: include the response of the unitary
# decoupling transformation (see section
# "Exact Two-Component Theory")
# ---------------------------------------------------
# Order of DKH treatment (this has no effect on ZORA calculations)
# ---------------------------------------------------
order 1 # first-order DKH Hamiltonian
2 # second-order DKH Hamiltonian
# ---------------------------------------------------
# Kind of Foldy-Wouthuysen transformation for picturechange effects
# in g tensors (see section "Picture-Change Effects")
# ---------------------------------------------------
fpFWtrafo true # do not include vector potential into momentum (default)
false # include vector potential
# ---------------------------------------------------
# Finite Nucleus Model: (see section "Finite Nucleus Model")
# ---------------------------------------------------
FiniteNuc false # Use point-charge nuclei (default)
true # Use finite nucleus model
# ---------------------------------------------------
# X2C intermediate storage level: (see section "X2C derivatives and properties")
# ---------------------------------------------------
StorageLevel 0-5 # default: 2
end
```
```{note}
It is important to recognize that in the one-center approximation
(`OneCenter true`) ALL methods can be used for geometry
optimization. Several papers in the literature show that this
approximation is fairly accurate for the calculation of structural
parameters and vibrational frequencies. Since this approximation is
associated with negligible computational effort relative to the
nonrelativistic calculation it is a recommended procedure.
```
(sec:relativistic.Hamiltonians.detailed)=
## Approximate Relativistic Hamiltonians
In the relativistic domain, calculations are based on the one-electron,
stationary Dirac equation in atomic units (rest mass subtracted)
$$h_{D} \Psi =\left({ \left({ \beta -1} \right)c^{2}+c \boldsymbol\alpha \cdot \boldsymbol p +V} \right)\Psi =E\Psi .
$$ (eqn:161)
The spinor $\Psi$ can be decomposed in its so-called large and small
components
$$\Psi =
\begin{pmatrix}
\Psi_L \\
\Psi_S
\end{pmatrix}
$$ (eqn:162)
These are obviously coupled through the Dirac equation. More precisely,
upon solving for $\Psi_{S}$, the following relation is obtained:
$$\Psi_{S} =\frac{1}{2c}\left(1+\frac{E-V}{2c^{2} } \right)^{-1}\boldsymbol\sigma\cdot\boldsymbol p \Psi_{L} =R\Psi_{L}
$$ (eqn:exact-relativistic)
Through the unitary transformation
$U=
\begin{pmatrix}
\Omega_+ & - R^+\Omega \\
R\Omega_+ & \Omega_-
\end{pmatrix}$ with
$\Omega_{+} =\frac{1}{\sqrt{1+R^{+}R} },\Omega_{-} =\frac{1}{\sqrt{ 1+RR^{+} } }$,
the Hamiltonian can be brought into block-diagonal form
$$U^{+}h_{D} U=
\begin{pmatrix}
\tilde h_{++} & 0 \\
0 & \tilde h_{--}
\end{pmatrix}
$$ (eqn:164)
The (electronic) large component thus has to satisfy the following
relation
$$h_{++} \Psi_{L} =\Omega_{+} \left({ h_{++} +h_{\pm } R+R^{+}\left({ h_{\mp
} +h_{--} R} \right)} \right)\Omega_{+} \Psi_{L} =E_{+} \Psi_{L} .
$$ (eqn:exact-relativistic-2c)
The approximate relativistic schemes implemented in ORCA use different
methods to substitute the exact relation
{eq}`eqn:exact-relativistic-2c` with approximate ones.
Two approximation schemes are available in ORCA: the regular
approximation and the Douglas-Kroll-Hess (DKH) approach.
(sec:relativistic.regularapprox.detailed)=
## The Regular Approximation
In the regular approximation,
{eq}`eqn:exact-relativistic-2c` is approximated by
$$R=\frac{c}{2c^{2}-V} \boldsymbol\sigma \cdot \boldsymbol p.
$$ (eqn:166)
At the zeroth-order level (ZORA), $\Omega_{\pm } =1$, so that the ZORA
transformation is simply
$$U_{\mathrm{ZORA} } =
\begin{pmatrix}
1 & -R^+ \\
R & 1
\end{pmatrix}
$$ (eqn:167)
and the corresponding Hamiltonian given by
$$\tilde{{h} }_{++}^{\mathrm{ZORA} } =V+c \boldsymbol\sigma \cdot \boldsymbol p \frac{1}{2c^{2}-V}c\boldsymbol\sigma \cdot \boldsymbol p.
$$ (eqn:168)
At the infinite-order level (IORA), $\Omega_{\pm }$is taken into
account, so that
$$U_{\mathrm{IORA} } =U_{\mathrm{ZORA} }
\begin{pmatrix}
\Omega_+ & 0 \\
0 & \Omega_-
\end{pmatrix}
$$ (eqn:169)
and
$$\tilde{{h} }_{++}^{\mathrm{IORA} } =\Omega_{+} \left({ V+c \boldsymbol\sigma \cdot \boldsymbol p \frac{1}{2c^{2}-V}c \boldsymbol\sigma \cdot \boldsymbol p} \right)\Omega_{+}
$$ (eqn:170)
is the corresponding Hamiltonian. Note that despite the name --
*infinite-order* regular approximation -- this is still *not* exact.
In ORCA, the spin-free (scalar-relativistic) variant of ZORA and IORA
are implemented. These are obtained from those above through the
replacement
$$\boldsymbol\sigma \cdot \boldsymbol p \frac{1}{2c^{2}-V} \boldsymbol\sigma \cdot \boldsymbol p \to \boldsymbol p \frac{1}{2c^{2}-V} \boldsymbol p.
$$ (eqn:171)
The regular Hamiltonians contain only part of the Darwin term and no
mass-velocity term. A problem with relations
{eq}`eqn:170` and
{eq}`eqn:168` is
that due to the non-linear dependence of the resulting regular
Hamiltonians on $V$, a constant change of $V$, which in the Dirac and
Schrödinger equations will result in a corresponding change of energy
$$E\to E+\text{const}$$
does not so in the regular approximation. Several attempts have been
made to circumvent this problem. The scaled ZORA variant is one such
procedure. Another one is given through the introduction of model
potentials replacing $V$. Both approaches are available in ORCA.
### The scaled ZORA variant
This variant goes back to van Lenthe et al. {cite}`vanlenthe1994jcp`. The
central observation is that the Hamiltonian
$$h_{\mathrm{scaledZORA} } =\frac{h_{\mathrm{ZORA} } }{1+\left\langle { \Psi_{L} \left|{c \boldsymbol\sigma \cdot \boldsymbol p \frac{1}{\left({ 2c^{2}-V}
\right)^{2} }c \boldsymbol\sigma \cdot \boldsymbol p} \right|\Psi_{L} } \right\rangle}
$$ (eqn:172)
produces constant energy-shifts $E\to E+\text{const}$ when the potential
$V$ is changed by a constant -- *for hydrogenic ions*. For many-electron
systems, the scaled-ZORA Hamiltonian still does not yield simple,
constant energy shift for $V\to
V+\text{const}$. But it produces the exact Dirac energy for
hydrogen-like atoms and performs better than the first-order regular
approximation for atomic ionization energies.
(sec:relativistic.zora.detailed)=
### The regular approximation with model potential
The idea of this approach goes back to Van Wüllen {cite}`vanwuellen1998jcp`,
who suggested the procedure for DFT. However we also use it for other
methods. The scalar relativistic ZORA self-consistent field equation is
in our implementation (in atomic units):
$$\left[{ \mathrm{\mathbf{p} }\frac{c^{2} }{2c^{2}-V }\mathrm{\mathbf{p} }+V_{\mathrm{eff} } } \right]\psi_{i} =\varepsilon_{i} \psi_{i}
$$ (eqn:173)
where $c$ is the speed of light. It looks like the normal
nonrelativistic Kohn--Sham equation with the KS potential
$V_{\mathrm{eff} }$:
$$V_{\mathrm{eff} } \left({ \mathrm{\mathbf{r} }} \right)=-\sum\limits_A { \frac{Z_{A} }{\left|{\mathrm{\mathbf{r} }-\mathrm{\mathbf{R} }_{A} } \right|} } +\int \frac{\rho \left( \mathrm{\mathbf{r}'} \right)}{\left| \mathrm{\mathbf{r} } - \mathrm{\mathbf{r}'} \right|}
d\mathrm{\mathbf{r}'} + V_{\text{xc} } \left[ \rho \right](\mathbf r)
$$ (eqn:174)
($Z_{A}$ is the charge of nucleus $A$ and $R_{A}$ is its position;
$\rho (r)$ is the total electron density and
$V_{xc} \left[ \rho \right]$ the exchange-correlation potential -- the
functional derivative of the exchange-correlation energy with respect to
the density). The kinetic energy operator $T=-\frac{1}{2}\nabla^{2}$ of
the nonrelativistic treatment is simply replaced by the ZORA kinetic
energy operator:
$$T^{\mathrm{ZORA} }=\mathrm{\mathbf{p} }\frac{c^{2} }{2c^{2}-V }\mathrm{\mathbf{p} }
$$ (eqn:175)
Clearly, in the regions where the potential $V$ is small compared to
$c^{2}$, this operator reduces to the nonrelativistic kinetic energy.
$V$ could be the actual KS potential. However, this would require to
solve the ZORA equations in a special way which demands recalculation of
the kinetic energy in every SCF cycle. This becomes expensive and is
also undesirable since the ZORA method is not gauge invariant and one
obtains fairly large errors from such a procedure unless special
precaution is taken. Van Wüllen {cite}`vanwuellen1998jcp` has therefore
argued that it is a reasonable approximation to replace the potential
$V$ with a model potential $\tilde{{V} }_{\mathrm{model} }$ which is
constructed as follows:
$$
\tilde{{V} }_{\mathrm{model} } =-\sum\limits_A { \frac{Z_{A} }{\left|{ \mathrm{\mathbf{r} }
-\mathrm{\mathbf{R} }_{A} } \right|} } +\int \frac{\rho_{\mathrm{model} } \left( \mathrm{\mathbf{r}'} \right)}{\left| \mathrm{\mathbf{r} }
-\mathrm{\mathbf{r}'} \right|}
d\mathrm{\mathbf{r}'} +V_{\text{xc} }^{\text{LDA} } \left[{ \rho^{\mathrm{model} }} \right](\mathbf r)
$$ (eqn:176)
The model density is constructed as a sum over spherically symmetric
(neutral) atomic densities:
$$\rho_{\mathrm{model} } \left({ \mathrm{\mathbf{r} }} \right)=\sum\limits_A { \rho
^{A}\left({ \mathrm{\mathbf{r} }} \right)}
$$ (eqn:177)
Thus, this density neither has the correct number of electrons (for
charged species) nor any spin polarization. Yet, in the regions close to
the nucleus, where the relativistic effects matter, it is a reasonable
approximation. The atomic density is expanded in a sum of s-type
Gaussian functions like:
$$\rho^{A}\left({ \mathrm{\mathbf{r} }} \right)=\sum\limits_i { d_{i} \exp \left({-\alpha_{i} \left|{ \mathrm{\mathbf{r} }-\mathrm{\mathbf{R} }_{A} } \right|^{2} } \right)}
$$ (eqn:178)
The fit coefficients were determined in three different ways by near
basis set limit scalar relativistic atomic HF calculations and are
stored as a library in the program.
The fitting densities are available for elements up to Rn, as well as the actinoids.
Through the variable `ModelDens`
(*vide supra*) the user can choose between these fits and study the
dependence of the results in this choice (it should be fairly small
except, perhaps, with the heavier elements and the HF densities which
are not recommended). The individual components of the model potential
(eq. {eq}`eqn:176`)
can be turned on or off through the use of the variable `ModelPot`
(*vide supra*).
Van Wüllen has also shown that the calculation of analytical gradients
with this approximation becomes close to trivial and therefore scalar
relativistic all electron geometry optimizations become easily feasible
within the ZORA approach. However, since $T^{\mathrm{ZORA} }$ is
constructed by numerical integration it is very important that the user
takes appropriate precaution in the use of a suitable integration grid
and also the use of appropriate basis sets! In the case of
`OneCenter true` the numerical integration is done accurately along the
radial coordinate and analytically along the angular variables such that
too large grids are not necessary unless your basis set is highly
decontracted and contains very steep functions.
(sec:relativistic.dkh.detailed)=
## The Douglas-Kroll-Hess Method
The Douglas-Kroll-Hess (DKH) method expands the exact relation
{eq}`eqn:exact-relativistic-2c` in the external potential V. In
ORCA the first- and second-order DKH methods are implemented. The
first-order DKH Hamiltonian is given by
$$\tilde{h}_{++}^{(1) } =E_{p} +A_{p} VA_{p} +B_{p} V^{(p) }B_{p} ,
$$ (eqn:179)
with
$$E_{P} =\sqrt{ c^{4}+c^{2}p^{2} }, \, A_{p} =\sqrt{\frac{E_{p} +c^{2} }{2E_{p} }}, \, B_{p} =\frac{c}{\sqrt{ 2E_{p} (E_{p} +c^{2}) } }
$$ (eqn:180)
At second order, it reads
$$\tilde{{h} }_{++}^{(2) } =\tilde{{h} }_{++}^{(1) } +\frac{1}{2}\left[{ W_{p} ,O} \right]$$ (eqn:181)
where
$$\left\{{ W_{p} ,E_{p} } \right\}=\beta O, \, O=A_{p} \left[{ R_{p} ,V}
\right]A_{p}, \, R_{p} =\frac{c\sigma p}{E_{p} +c^{2} }
$$ (eqn:182)
define the second-order contribution. In ORCA, the spin-free part of
$\tilde{{h} }_{++}^{(2) }$ is implemented.
The occurrence of the relativistic kinetic energy, $E_{P}$, which is not
well-defined in position space, makes a transformation to the
$p^{2}$-eigenspace necessary. Thus any DKH calculation will start with a
decontraction of the basis set, to ensure a good resolution of the
identity. Then the non-relativistic kinetic energy is diagonalized and
the $E_{P}$-dependent operators calculated in that space. The potential
$V$ and $V^{(p) }$ are transformed to $p^{2}$-eigenspace. After all
contributions are multiplied to yield the (first- or second-order)
Hamiltonian, the transformation back to AO space is carried out and the
basis is recontracted.
The (spin-free) DKH-Hamiltonians contain all spin-free, relativistic
correction terms, e.g. the mass-velocity and Darwin terms. As the
potential enters linearly, no scaling or model potential is necessary to
introduce the correct behaviour of the energy under a change
$$V\to V+\text{const}.$$
In all these respects the DKH Hamiltonians are much cleaner than the
regular Hamiltonians.
(sec:relativistic.picturechange.detailed)=
## Picture-Change Effects
Irrespective of which Hamiltonian has been used in the determination of
the wave function, the calculation of properties requires some special
care. This can be understood in two ways: First of all, we changed from
the ordinary Schrödinger Hamiltonian to a more complicated Hamiltonian.
As properties are defined as derivatives of the energy, it is clear that
a new Hamiltonian will yield a new expression for the energy and thus a
new and different expression for the property in question. Another way
of seeing this is that through the transformation $U$, we changed not
only the Hamiltonian but also the wave function. To obtain the property
at hand as the expectation value of the property operator with the wave
function, we have to make sure that property operator and wave function
are actually given in the same space. This is done through a
transformation of either the property operator or the wave function.
In any case, the difference between the non-relativistic and (quasi)
relativistic property operator evaluated between the (quasi)
relativistic wave function is called the picture-change effect. From
what was said above, this is clearly not a physical effect. It describes
how consistent the quasi relativistic calculation is carried out. A
fully consistent calculation requires the determination of the wave
function on the (quasi) relativistic level as well as the use of the
(quasi) relativistic property operator. This is obtained through the
choice
```orca
%rel PictureChange 1 end # or 2 - see below
```
in the `%rel` block. It may be that the (quasi) relativistic and
non-relativistic property operator do produce similar results. In this
case, a calculation with picture-change turned off
(`PictureChange=0`) may be a good approximation. This is, however,
not the rule and cannot be predicted before carrying out the
calculation. It is therefore highly recommended to turn on picture-change
in all (quasi) relativistic property calculations!
For DKH2, the fully consistent picture-change effects are obtained
using the same transformation order for the property operator
as for the one-electron Hamiltonian, i.e. setting
```orca
%rel PictureChange 2 end
```
while with `PictureChange=1`
only first-order changes on the property operators are taken into
account, which reduces the computational cost. However, since this
is in no way a significant reduction, this choice is not recommended.
A similar argument applies to X2C (see {ref}`sec:relativistic.x2c.props.detailed`).
For magnetic properties, the DKH transformation and consequently the DKH
Hamiltonian and the corresponding property operators are not unique.
Depending on whether the magnetic field is included in the free-particle
Foldy--Wouthuysen (fpFW) transformation carried out in the first step of
the DKH protocol or not, two different Hamiltonians result. If the
magnetic field is included in the fpFW transformation, the resulting
Hamiltonian is a function of the gauge invariant momentum
$$\boldsymbol{\pi}=\textbf{p}+\textbf{A}.$$
It is therefore gauge invariant under gauge transformations of the
magnetic vector potential $\mathbf A$ and thus are the property
operators derived from it. This is referred to as f$\pi$FW DKH
Hamiltonian. If the magnetic field is not included in the FW
transformation, the resulting Hamiltonian is a function of the kinetic
momentum $\textbf{p}$ only and thus is not gauge invariant. The latter
Hamiltonian is referred to as fpFW DKH Hamiltonian. A comparison of both
Hamiltonians is given in Table {numref}`table:compdkhham`.
(table:compdkhham)=
:::{table} Comparison of the properties of the fpFW and fπFW DKH Hamiltonians. For details see Ref. {cite}`sandhoefer2012jcp`.
| Criterion | fpFW Hamiltonian | fπFW Hamiltonian |
|:-----------------------------|:----------------:|:--------------------:|
| Convergence of Eigenvalues to Dirac Eigenvalues | ? | yes |
| 1st order is bounded | no | yes |
| Reproduces Pauli Hamiltonian | no | yes |
| Gauge invariance | no | yes |
| Lorentz invariance | no | no |
:::
From this Table, it becomes clear that the f$\pi$FW DKH Hamiltonian is
clearly preferred over the fpFW Hamiltonian. To obtain the property
operators, it is however necessary to take the derivatives of these
Hamiltonians. It turns out that in the case of the hyperfine-coupling
tensor, the necessary derivatives produce divergent property operators
in the case of the f$\pi$FW DKH Hamiltonian. This may be due to the
unphysical assumption of a point-dipole as a source of the magnetic
field of the nucleus. As a physical description of the magnetization
distribution of the nucleus is not available due to a lack of
experimental data, the magnetization distribution is assumed to be the
same as the charge distribution of the nucleus, see
Section {ref}`sec:relativistic.finnuc.detailed`. This is unphysical as the
magnetization is caused by the one unpaired nucleon in the nucleus
whereas the charge distribution is generated by the protons in the
nucleus. So, physically, the magnetization should occupy a larger volume
in space than the charge. This might also be the reason why the
resulting finite-nucleus model is insufficient to remedy the
divergences in the f$\pi$FW hyperfine-coupling tensor. Consequently,
the hyperfine-coupling tensor is only implemented in the version
resulting from the fpFW DKH Hamiltonian. In the case of the g-tensor
both versions are implemented and accessible via the keyword
```orca
%rel fpFWtrafo true/false end
```
By default, this keyword is set to `true`. A detailed form of the
property operators used for the g-tensor and hyperfine-tensors can be
found in Ref. {cite}`sandhoefer2012jcp`.
(sec:relativistic.finnuc.detailed)=
## Finite Nucleus Model
Composite particles like nuclei have, as opposed to elementary
particles, a certain spatial extent. While the point-charge
approximation for nuclei is in general very good in nonrelativistic
calculations, in relativistic calculations it might lead to
non-negligible errors. A finite-nucleus model is available for all
calculations in the ORCA program package. It is accessible from the
`%rel` block via
```orca
%rel FiniteNuc true/false end
```
By default, this keyword is set to `false`. If the keyword is set to
`true`, finite-nucleus effects are considered in the following
integrals:
- nucleus potential V
- DKH-integral $V^{(p) }$
- one-electron spin-orbit integrals SOC (also in one-electron part of
SOMF)
- electric-field gradient EFG (and thus, as a consequence in the
Fermi-contact and spin-dipole terms of the HFC tensor)
- nucleus-orbit integral NUC
- angular-momentum integral l
The finite-nucleus model implemented in ORCA is the Gaussian nucleus
model of Ref. {cite}`visscher1997atom`.
(sec:relativistic.x2c.detailed)=
## Exact Two-Component Theory (X2C)
The X2C implementation in ORCA closely follows that of Franzke, Weigend,
and their coworkers, described in the references: {cite}`Peng2013` (energy),
{cite}`Franzke2018` (gradient), {cite}`Franzke2022` (EPR hyperfine coupling),
{cite}`Franzke2021a` (NMR spin--spin coupling), {cite}`Franzke2019` (NMR
shielding). These are also consistent with the work of Gauss and
coworkers in refs: {cite}`Cheng2011` (gradient, electric properties),
{cite}`Cheng2011a` (Hessian), {cite}`Cheng2013` (NMR shielding). However, despite
the name, only a scalar relativistic (spin-free) version is available at
present, resulting effectively in a one-component method (more aptly
called "SF-X2C-1e"), very similar to the DKH and ZORA approaches
described above. The main difference to the latter two is that the
decoupling of the one-electron Dirac Hamiltonian is exact, rather than
approximate. SF-X2C-1e is implemented in ORCA for energies, gradients,
and various properties (see below) with both a point- and finite nucleus
model.
We briefly describe the working equations here, using the notation of
the aforementioned references, which differs somewhat from that in the
previous sections. The one-electron Dirac equation is solved directly:
$$\mathbb{D C} = \mathbb{M C E},
\quad
\mathbb{C} =
\begin{pmatrix}
\mathbf C^{\mathrm L}_{+} & \mathbf C^{\mathrm L}_{-} \\
\mathbf C^{\mathrm S}_{+} & \mathbf C^{\mathrm S}_{-} \\
\end{pmatrix}
$$ (eq:x2c-DCSCE)
to obtain the unitary transformation matrix:
$$\begin{aligned}
\mathbb{U} &=
\begin{pmatrix}
\mathbf U_{++} & \mathbf U_{+-} \\
\mathbf U_{-+} & \mathbf U_{--} \\
\end{pmatrix}
=
\begin{pmatrix}
\mathbf 1 & -\mathbf X^\dagger \\
\mathbf X & \mathbf 1 \\
\end{pmatrix}
\begin{pmatrix}
\mathbf R & \mathbf 0 \\
\mathbf 0 & \mathbf R' \\
\end{pmatrix}
\\
\mathbf X &= \mathbf C^{\mathrm S}_{+} \left( \mathbf C^{\mathrm L}_{+} \right)^{-1}\\
\mathbf R &= \mathbf S^{-\frac{1}{2} } \left(\mathbf S^{-\frac{1}{2} }\tilde{\mathbf S} \mathbf S^{-\frac{1}{2} }\right)^{-\frac{1}{2} } \mathbf S^{\frac{1}{2} }\\
\tilde{\mathbf S} &= \mathbf S + \frac{1}{2c^2}\mathbf X^\dagger \mathbf T \mathbf X
\end{aligned}$$
which exactly block-diagonalizes the Hamiltonian:
$$\mathbb{U^\dagger D U} =
\begin{pmatrix}
\mathbf h^{+} & \mathbf 0 \\
\mathbf 0 & \mathbf h^{-} \\
\end{pmatrix}
$$ (eq:x2c-UDU)
$$\mathbf h^{+} =
\mathbf R^\dagger \mathbf L \mathbf R
$$ (eq:x2c-h)
$$ \mathbf L = \mathbf V + \mathbf T \mathbf X + \mathbf X^\dagger \mathbf T
+\mathbf X^\dagger \left(\frac{1}{4c^2} \mathbf W - \mathbf T\right) \mathbf X
$$ (eq:x2c-L)
where $\mathbf V$, $\mathbf T$, $\mathbf S$, and $\mathbf W$ are the
potential, kinetic, overlap, and relativistic potential
($\hat p\,\hat V\hat p$) integral matrices. $\mathbf h^{+}$ is thus the
matrix form of the relativistically-corrected one-electron Hamiltonian
used for the rest of the calculation.
Note that the potential operator $\hat V$ used in the above equations
only includes the electron--nuclear Coulomb interaction. External point
charges may optionally be included via `%rel ModelPot[4]=1` (default 0).
This option affects energy and NMR shielding calculations but it is
presently not available for gradients or Hessians.
(sec:relativistic.x2c.dlu.detailed)=
### DLU approximation
The diagonal local approximation to the unitary transformation matrix
(DLU), as introduced by Peng and Reiher,{cite}`Peng2012` reduces the
computational cost of the X2C transformation by approximating
$\mathbb U$ (i.e., $\mathbf R$ and $\mathbf X$) as an
atomic-block-diagonal matrix.
$$\begin{aligned}
\mathbf R &\approx \bigoplus_A \mathbf R_A
&
\mathbf X &\approx \bigoplus_A \mathbf X_A\end{aligned}$$
where $\bigoplus_A$ denotes a direct sum of atomic diagonal blocks.
Eqs {eq}`eq:x2c-DCSCE`--{eq}`eq:x2c-L` can then be solved independently for diagonal
atomic blocks $\mathbf h^+_{AA}$. Off-diagonal blocks $\mathbf h^+_{AB}$
are obtained as:
$$\begin{aligned}
\mathbf h^+_{AB} &= \mathbf R^\dagger_A \mathbf L_{AB} \mathbf R_B
\\
\mathbf L_{AB} &= \mathbf V_{AB} + \mathbf T_{AB} \mathbf X_B + \mathbf X^\dagger_A \mathbf T_{AB} + \mathbf X^\dagger_A \left(\frac{1}{4c^2}\mathbf W_{AB} - \mathbf T_{AB}\right) \mathbf X_B\end{aligned}$$
It is also possible to approximate light atoms $a,b$ (below a certain
atomic number) as non-relativistic:
$$\begin{aligned}
\mathbf X_a &= \mathbf R_a = \mathbf I
\\
\mathbf W_{ab} &= \mathbf W_{aB} = \mathbf 0
\\
\mathbf h^+_{ab} &= \mathbf V_{ab} + \mathbf T_{ab}
\\
\mathbf h^+_{aB} &= \left(\mathbf V_{aB} + \mathbf T_{aB} \mathbf X_B\right) \mathbf R_B\end{aligned}$$
Unlike the one center approximation (which is also available for X2C),
all nuclei are included in the potential operator $\hat V$. Use of the
DLU approximation is controlled via:
```orca
%rel
DLU false # default - not used
true # turn DLU on
LightAtomThresh 0 # (default) highest atomic number treated as non-relativistic
end
```
(sec:relativistic.x2c.props.detailed)=
### X2C derivatives and properties
As discussed in
section {ref}`sec:relativistic.picturechange.detailed`, when computing
properties with relativistic methods, derivatives need to be taken of
the correct Hamiltonian, namely $\mathbf h^+$
(eq {eq}`eq:x2c-h`) in the X2C case. These include contributions due
to the derivatives of $\mathbf R$ and $\mathbf X$, which are often
small. Therefore, it is possible to neglect them and save some
computational time via the following option:
```orca
%rel
PictureChange 2 # compute the full relativistic Hamiltonian derivative
1 # neglect derivatives of R and X
0 # (default) use the non-relativistic operator
end
```
The same setting is applied to all properties for which the X2C
correction is implemented. Currently, these are: geometric gradients,
electric dipoles, quadrupoles, and polarizabilities, electric field
gradients, EPR hyperfine couplings, and NMR shieldings and spin--spin couplings.
For the Hessian, the X2C correction is implemented in a semi-numeric fashion.
The DLU approximation is applied throughout, if requested,
and reduces the computational effort dramatically. Note that for
magnetic properties, the restricted magnetic balance (RMB) for the small
component basis functions is used whenever GIAOs are requested (e.g. NMR
shielding), while the restricted kinetic balance (RKB) is used otherwise
(e.g. NMR coupling). The spin--orbit coupling integrals used for various
properties only include a relativistic correction to the one-electron
term, as is the case for DKH.
For second derivative properties whose number is proportional to the
number of atoms, e.g. the DSO term of NMR couplings, first derivatives
of various intermediate quantities required for the full X2C second
derivative are stored on disk. The storage requirements can be reduced
via the `StorageLevel` keyword and the missing intermediates will be
recomputed on-the-fly, which of course increases the computation time.
```orca
%rel
StorageLevel 0 # do not store anything - not always available
1 # store integral derivatives
2 # 1 + derivatives of R and X (default, minimum for DLU)
3 # 2 + derivatives of L
4 # 3 + derivatives of S-tilde
5 # 4 + derivatives of D and M
end
```
(relativistic.basissets.detailed)=
## Basis Sets in Relativistic Calculations
For relativistic calculations, special basis sets have been designed,
both as DKH and ZORA recontractions of the non-relativistic Ahlrichs
basis sets (in their all-electron versions) for elements up to Kr, and
as purpose-built segmented all-electron relativistically contracted
(SARC) basis sets for elements beyond Kr
{cite}`Pantazis2008,Buehl2008,Pantazis2009,Pantazis2011,Pantazis2012,Rolfes2020`.
Their names are "ZORA-" or "DKH-" followed by the conventional basis set
name. For X2C calculations, the "x2c-XVPall" basis sets and their
variants are available.{cite}`Weigend2017-x2c-XVP,Weigend2019-x2c-XVP-s`
See section
{ref}`sec:basisset.detailed` for a complete list of basis sets.