7.23. 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:

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

7.23.1. Approximate Relativistic Hamiltonians

In the relativistic domain, calculations are based on the one-electron, stationary Dirac equation in atomic units (rest mass subtracted)

(7.162)\[h_{D} \Psi =\left({ \left({ \beta -1} \right)c^{2}+c \boldsymbol\alpha \cdot \boldsymbol p +V} \right)\Psi =E\Psi . \]

The spinor \(\Psi\) can be decomposed in its so-called large and small components

(7.163)\[\begin{split}\Psi = \begin{pmatrix} \Psi_L \\ \Psi_S \end{pmatrix} \end{split}\]

These are obviously coupled through the Dirac equation. More precisely, upon solving for \(\Psi_{S}\), the following relation is obtained:

(7.164)\[\Psi_{S} =\frac{1}{2c}\left(1+\frac{E-V}{2c^{2} } \right)^{-1}\boldsymbol\sigma\cdot\boldsymbol p \Psi_{L} =R\Psi_{L} \]

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

(7.165)\[\begin{split}U^{+}h_{D} U= \begin{pmatrix} \tilde h_{++} & 0 \\ 0 & \tilde h_{--} \end{pmatrix} \end{split}\]

The (electronic) large component thus has to satisfy the following relation

(7.166)\[h_{++} \Psi_{L} =\Omega_{+} \left({ h_{++} +h_{\pm } R+R^{+}\left({ h_{\mp } +h_{--} R} \right)} \right)\Omega_{+} \Psi_{L} =E_{+} \Psi_{L} . \]

The approximate relativistic schemes implemented in ORCA use different methods to substitute the exact relation (7.166) with approximate ones.

Two approximation schemes are available in ORCA: the regular approximation and the Douglas-Kroll-Hess (DKH) approach.

7.23.2. The Regular Approximation

In the regular approximation, (7.166) is approximated by

(7.167)\[R=\frac{c}{2c^{2}-V} \boldsymbol\sigma \cdot \boldsymbol p. \]

At the zeroth-order level (ZORA), \(\Omega_{\pm } =1\), so that the ZORA transformation is simply

(7.168)\[\begin{split}U_{\mathrm{ZORA} } = \begin{pmatrix} 1 & -R^+ \\ R & 1 \end{pmatrix} \end{split}\]

and the corresponding Hamiltonian given by

(7.169)\[\tilde{{h} }_{++}^{\mathrm{ZORA} } =V+c \boldsymbol\sigma \cdot \boldsymbol p \frac{1}{2c^{2}-V}c\boldsymbol\sigma \cdot \boldsymbol p. \]

At the infinite-order level (IORA), \(\Omega_{\pm }\)is taken into account, so that

(7.170)\[\begin{split}U_{\mathrm{IORA} } =U_{\mathrm{ZORA} } \begin{pmatrix} \Omega_+ & 0 \\ 0 & \Omega_- \end{pmatrix} \end{split}\]

and

(7.171)\[\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_{+} \]

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

(7.172)\[\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. \]

The regular Hamiltonians contain only part of the Darwin term and no mass-velocity term. A problem with relations (7.171) and (7.169) 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.

7.23.2.1. The scaled ZORA variant

This variant goes back to van Lenthe et al. [865]. The central observation is that the Hamiltonian

(7.173)\[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} \]

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.

7.23.2.2. The regular approximation with model potential

The idea of this approach goes back to Van Wüllen [867], 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):

(7.174)\[\left[{ \mathrm{\mathbf{p} }\frac{c^{2} }{2c^{2}-V }\mathrm{\mathbf{p} }+V_{\mathrm{eff} } } \right]\psi_{i} =\varepsilon_{i} \psi_{i} \]

where \(c\) is the speed of light. It looks like the normal nonrelativistic Kohn–Sham equation with the KS potential \(V_{\mathrm{eff} }\):

(7.175)\[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) \]

(\(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:

(7.176)\[T^{\mathrm{ZORA} }=\mathrm{\mathbf{p} }\frac{c^{2} }{2c^{2}-V }\mathrm{\mathbf{p} } \]

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 [867] 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:

(7.177)\[ \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) \]

The model density is constructed as a sum over spherically symmetric (neutral) atomic densities:

(7.178)\[\rho_{\mathrm{model} } \left({ \mathrm{\mathbf{r} }} \right)=\sum\limits_A { \rho ^{A}\left({ \mathrm{\mathbf{r} }} \right)} \]

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:

(7.179)\[\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)} \]

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

7.23.3. The Douglas-Kroll-Hess Method

The Douglas-Kroll-Hess (DKH) method expands the exact relation (7.166) in the external potential V. In ORCA the first- and second-order DKH methods are implemented. The first-order DKH Hamiltonian is given by

(7.180)\[\tilde{h}_{++}^{(1) } =E_{p} +A_{p} VA_{p} +B_{p} V^{(p) }B_{p} , \]

with

(7.181)\[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}) } } \]

At second order, it reads

(7.182)\[\tilde{{h} }_{++}^{(2) } =\tilde{{h} }_{++}^{(1) } +\frac{1}{2}\left[{ W_{p} ,O} \right]\]

where

(7.183)\[\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} } \]

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.

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

%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

%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 X2C derivatives and properties).

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 Table 7.19.

Table 7.19 Comparison of the properties of the fpFW and fπFW DKH Hamiltonians. For details see Ref. [746].

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 Finite Nucleus Model. 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

%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. [746].

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

%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. [872].

7.23.6. Exact Two-Component Theory (X2C)

The X2C implementation in ORCA closely follows that of Franzke, Weigend, and their coworkers, described in the references: [659] (energy), [274] (gradient), [276] (EPR hyperfine coupling), [273] (NMR spin–spin coupling), [272] (NMR shielding). These are also consistent with the work of Gauss and coworkers in refs: [165] (gradient, electric properties), [166] (Hessian), [167] (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:

(7.184)\[\begin{split}\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} \end{split}\]

to obtain the unitary transformation matrix:

\[\begin{split}\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}\end{split}\]

which exactly block-diagonalizes the Hamiltonian:

(7.185)\[\begin{split}\mathbb{U^\dagger D U} = \begin{pmatrix} \mathbf h^{+} & \mathbf 0 \\ \mathbf 0 & \mathbf h^{-} \\ \end{pmatrix} \end{split}\]
(7.186)\[\mathbf h^{+} = \mathbf R^\dagger \mathbf L \mathbf R \]
(7.187)\[ \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 \]

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.

7.23.6.1. DLU approximation

The diagonal local approximation to the unitary transformation matrix (DLU), as introduced by Peng and Reiher,[660] 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 (7.184)–(7.187) can then be solved independently for diagonal atomic blocks \(\mathbf h^+_{AA}\). Off-diagonal blocks \(\mathbf h^+_{AB}\) are obtained as:

\[\begin{split}\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}\end{split}\]

It is also possible to approximate light atoms \(a,b\) (below a certain atomic number) as non-relativistic:

\[\begin{split}\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}\end{split}\]

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:

%rel
  DLU         false # default - not used
               true # turn DLU on
  LightAtomThresh 0 # (default) highest atomic number treated as non-relativistic
end

7.23.6.2. X2C derivatives and properties

As discussed in section Picture-Change Effects, when computing properties with relativistic methods, derivatives need to be taken of the correct Hamiltonian, namely \(\mathbf h^+\) (eq (7.186)) 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:

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

%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

7.23.7. 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 [62, 641, 642, 643, 644, 729]. 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.[275, 690] See section Choice of Basis Set for a complete list of basis sets.