6.5 Relativistic effects

Turbomole provides two different possibilities for the treatment of relativistic effects: Via effective core potentials (ECPs) or via all-electron approaches (X2C, DKH, BSS). Both techniques can be employed in a one-component (scalar-relativistic) or two-component (including spin-orbit coupling) framework. The latter is only available in the modules ridft, rdgrad, escf, mpshift, ricc2, and riper.

6.5.1 One- and two-component relativistic methods

Incorporation of scalar-relativistic effects leads to additional contributions to the one-electron integrals (either from ECP or all-electron approach). The program structure is the same as in non-relativistic theory (all quantities are real). Two-component treatments allow for self-consistent calculations including spin-orbit interactions. These may be particularly important for compounds containing heavy elements (additionally to scalar-relativistic effects). Two-component treatments require the use of complex two-component orbitals (spinors) \[\begin{equation*} \psi_i({\bf x}) = \begin{pmatrix} \psi^{\alpha}_i({\bf r}) \\ \psi^{\beta}_i({\bf r}) \\ \end{pmatrix} \end{equation*}\] instead of real (non-complex) one-component orbitals in non-relativistic or scalar-relativistic treatments. The Hartree–Fock and Kohn–Sham equations are now spinor equations with a complex Fock operator \[\begin{equation*} \begin{pmatrix} \hat{F}^{\alpha \alpha} & \hat{F}^{\alpha \beta} \\ \hat{F}^{\beta \alpha} & \hat{F}^{\beta \beta} \\ \end{pmatrix} \begin{pmatrix} \psi^{\alpha}_i({\bf r}) \\ \psi^{\beta}_i({\bf r}) \\ \end{pmatrix} = \epsilon_i \begin{pmatrix} \psi^{\alpha}_i({\bf r}) \\ \psi^{\beta}_i({\bf r}) \\ \end{pmatrix}. \end{equation*}\] The wavefunction is no longer an eigenfunction of the spin operator, the spin vector is no longer an observable.

In case of DFT for open-shell systems, rotational invariance of the exchange-correlation energy is ensured by the non-collinear approach. In this approach, the exchange-correlation energy is a functional of the particle density and the absolute value of the spin-vector density $\vec{m}({\bf r})$ ($\vec{{\pmb \sigma}}$ are the Pauli matrices)

\[\begin{equation*} \vec{m}({\bf r}) = \sum_i \psi_i^{\dagger}({\bf x})\vec{{\pmb \sigma}}\psi_i({\bf x}). \end{equation*}\]

This quantity replaces the spin-density (difference between density of alpha and beta electrons) of non- or scalar-relativistic treatments.

For closed-shell species, the Kramers-restricted scheme, a generalization of the RHF-scheme of one component treatments, is applicable.

Effective core potentials

The most economic way to account for relativistic effects is via effective core potentials by choosing either the one- or the two-component ECP (and for the latter additionally setting $soghf in the control file or in define). The theoretical background and the implementation for the two-component SCF procedure is described in Ref. [162]. For recommendations concerning specific ECPs and corresponding basis sets see below.

Relativistic all-electron approaches (X2C, DKH, BSS)

Relativistic calculations are based on the Dirac rather than on the Schrödinger Hamiltonian. Since the Dirac Hamiltonian introduces pathological negative-energy states and requires extensive one-electron basis set expansions, methods have been devised which allow one to calculate a matrix representation of that part of the Dirac Hamiltonian, which describes electronic states only. For this, a unitary transformation is employed to block-diagonalize the Dirac Hamiltonian and thus to decouple the negative-energy states from the electronic states. For reasons of efficiency, this transformation is carried out only for the one-electron part of the full Hamiltonian (as a consequence, the two-electron interaction will then be slightly affected by a so-called picture-change effect). The resulting quantum chemical approach, “exact two-component” (X2C), was developed by several groups starting with formal work in the mid-1980s. X2C is related to the step-wise Douglas–Kroll–Hess (DKH) approach, which achieves decoupling in sequential manner. For the latter, the number of transformation steps is called the order of DKH. Infinite-order DKH yields identical results compared to X2C, but – in contrast to the latter – not feasible. Eighth order DKH usually yields results similar to X2C at similar cost. X2C is also related to the Barysz–Sadlej–Snijders (BSS) method, that first applies the free-particle Foldy-Wouthuysen transformation (which is the first mandatory step in DKH), and then constructs the one-step exact decoupling transformation of X2C. These three approaches have been reviewed and directly compared in terms of formalism and results, respectively, in Ref. [163] (see also this reference for a complete bibliography on exact-decoupling methods).

Essentially, X2C methods change the one-electron Hamiltonian in a basis-set representation. The Schrödinger one-electron Hamiltonian (including the external potential of the atomic nuclei) is replaced by the transformed (upper-left block of the) Dirac Hamiltonian. Since the transformation is carried out in the fully decontracted primitive basis, all matrix operations needed for the generation of the relativistic one-electron Hamiltonian can be cumbersome and even prohibitive if the molecule is large. In order to solve this unfavorable scaling problem, a rigorous local approach, called DLU, has been devised [164] and is strongly recommended. Other local approximations available are the diagonal local approximation to the Hamiltonian (DLH), which uses the non-relativistic limit for the off-diagonal blocks, and the DLU(NB) scheme, which uses also the non-relativistic limit for atomic pairs with a large distance. That means, the unitary transformation matrices for these blocks are approximated with a unit matrix. Same is done for analytical derivative theory.

X2C, DKH, and BSS exist in full two-component (spin-(same-)orbit coupling including) and in a one-component scalar-relativistic form. Both have been implemented into the Turbomole package and all details on the efficient implementation have been described in Ref. [165]. Additionally implemented modifications like a finite nucleus model based on a Gaussian charge distribution [166] and a screened nuclear spin-orbit (SNSO) approach [167, 168, 169] are documented in Ref. [170], together with the implementation of analytical one- and two-component X2C gradients, also in their local variant. Note that we use the real atomic weight for the finite nucleus radius by default. The isotope number as proposed originally can be used if requested by the user.

Calculation of analytical energy gradients is available within the (local) X2C approach. The implementation is described in Ref. [170]. The finite nucleus model based on a Gaussian charge distribution can be used for energy and gradient calculations. To model the effect of spin-orbit coupling on the two-electron interaction a (modified) screened-nuclear-spin-orbit approximation can be applied to the spin-dependent one-electron integrals. The scalar-relativistic approach can be used with the modules grad, rdgrad, egrad, ricc2, mpgrad, and rirpa. The two-component Hamiltonian is only available in rdgrad.

In relativistic all-electron calculations, additional contributions from point charges or COSMO (see chapter 21.2) are not part of the relativistic decoupling. These contributions are evaluated after the X2C, DKH, or BSS step. Thus, they are calculated without picture-change transformation.

6.5.2 How to use

Scalar-relativistic calculations with ECPs

In case of scalar-relativistic calculations with ECPs, relativity is taken into account simply by the choice of a relativistic ECP together with a suited basis set. A reasonable choice are the Dirac–Hartree–Fock effective core potentials by the Stuttgart-Köln group labeled dhf-ecp in TURBOMOLE, together with optimized basis sets termed dhf-XVP (X = S, TZ, QZ)[171]. These ECPs and bases are available for Rb to Xe, except for f-elements. For H to Kr dhf and def2 bases are identical non-relativistic all-electron basis sets. For lanthanides and actinides at present no dhf ECPs are available. The usage of Wood-Boring type ECPs [172], for lanthanides together with def2-bases [173], for actinides together with the original Stuttgart–Köln bases [174, 175] (labeled def in TURBOMOLE) is recommended here. Note that for def2-bases the underlying ECPs are of different type (due to history). or p-elements, def2-bases are optimized for Dirac–Hartree–Fock ECPs, for the other elements for Wood–Boring ECPs.

Two-component calculations (general)

The keyword $soghf enforces the two-component calculations. Keywords for specification of the method of calculation are the same as for the one-component case. Additionally, for closed-shell species a Kramers invariant density functional formalism can be switched on with the keyword $kramers. The generalized DIIS scheme for complex Fock operators (GDIIS) is activated by default ($gdiis) and can be deactivated with the flag $nogdiis. For improvements on the SCF convergence behavior and gradients, please see Ref. [176]. It is recommended to use exact exchange for double- and triple-zeta bases instead of semi-numerical [149] or RI-$K$ approximations. Matters are different for response properties in

escf and mpshift, where we generally recommend to use the seminumerical exchange approximation for the response density with the $esenex flag. These keywords can be inserted into the control file manually or added in the scf section of define. By default, the RI-$J$ approximation is employed. The analytical Coulomb integrals can be used with the option $coulex in ridft, rdgrad, and mpshift. Note that two non-collinear formalisms are available, namely the canonical one (default) [162, 176, 43, 177] and the Scalmani–Frisch approach [178, 179, 180]. The latter is used with the flag $xcsf and the numerical stabilization approximation is used with $sfnumstab. Note that the stabilization approximation is not fully rotational invariant. Alternatively, a regularized version [181] is available with

$sfregstab. Here, the work on GGAs in [181] was also generalized to the meta-GGA and current-density ingredients as well as the corresponding local hybrid functionals. Further using

$sfcolstab will check additionally check for local collinearity at each grid point for the case, where the magnetization variables are non-zero, i.e. case 1 in [181].

For the two-component initial guess, Hückel theory, a superposition of atomic densities, or scalar SCF wavefunctions may be used. We recommend to use converged one-component molecular orbitals or a two-component superposition of atomic densities as starting point. The first option means that first a scalar-relativistic calculation is carried out without the keyword $soghf. Then, the keywords for a two-component calculation can be added and the two-component calculations can be started. This allows for a much smoother convergence. Here, the one-component orbitals are transformed to the two-component picture in the first iteration. By default, this is done with the Pauli matrix of the z component, meaning that the wavefunction is an eigenfunction of the z spin operator after the first iteration. Alternatively, you may specify the desired spin alignment with $sxeig, $syeig, or $szeig. The second option means that a non-collinear superposition of atomic densities with the desired spin alignment is formed. The resulting density matrix is employed to evaluate the two-electron Coulomb and HF exchange integrals followed by a diagonalization of the Fock matrix including the one-electron terms with the spin–orbit integrals. That is, this formalism uses Hartree–Fock (without the RI approximation) and adds a spin–orbit energy directly in the first iteration.

For open-shell molecules it is often helpful to increase the value for $scforbitalshift closedshell; a value of ca. 1.0 may serve as a rough recommendation. Likewise, the default value for the automatic orbitalshift can be adapted to improve the convergence. Additionally, canonical orthogonalization may be helpful to improve the convergence behavior, see Sec. 6.7 for details.

Generally, spin–orbit coupling induces a (paramagnetic) current density. Functionals such as meta-GGAs and local hybrids, which depend on the kinetic energy density, require the inclusion of this current density from a formal point of view. That is, functionals need to be constructed in the framework of current density functional theory (CDFT) [43, 177, 180]. This is applied by the keyword $curswitchengage in ridft, rdgrad, riper, mpshift, and escf. Currently, this formalism is restricted to a common gauge origin in mpshift. We stress that the canonical 2c formalism makes use of a projection onto the spin magnetization. This is avoided with the Scalmani–Frisch approach. Therefore, these two formalism lead to different results with CDFT, especially different closed-shell limits and the Kramers-restricted ansatz ($kramers) as shown in [180].

It is possible to turn off spin–orbit coupling with the keyword $noso. Then, no one-electron spin–orbit integrals are employed. Complex algebra is still employed. Likewise, the one-electron spin–orbit pVp integrals and their derivatives can be scaled with $soscal followed by the desired scaling factor. Default is one. Note that not only the respective integrals for the SCF procedure but also those for properties are rescaled. This option is available for 2c ECP and DKH/BSS/X2C calculations. Additionally, it is possible to only consider one spin-orbit component of the relativistically modified potential or the ECPs by using $socx, $socy, or $socz.

Two-component calculations with ECPs

For spin-orbit treatments, the two-component variants of ECPs (suffix -2c) are required, the use of extended basis sets accounting for the spatial splitting of inner p-shells (also suffix -2c) is recommended. ECPs dhf-ecp-2c and bases dhf-XVP-2c (X = S, TZ, QZ) [171] are available for Rb to Kr (f elements excepted); for references concerning ECPs see http://www.tc.uni-koeln.de/PP/clickpse.en.html. The two-component formalism may be most easily prepared and applied in the following way:

Run define: choose C$_1$-symmetry; select ECPs and basis sets with suffixes -2c for the respective elements. RI-$J$ and RI-$JK$ auxiliary basis sets are the same for dhf and def2 bases. Moreover, they are of sufficient flexibility for two-component treatments and provided automatically upon request. Switch on soghf and further desired options in the scf menu.
Start the two-component calculation with ridft.
At the end of the SCF procedure real and imaginary parts of spinors are written to files spinor.r and spinor.i, eigenvalues and spinor occupations are collected in the file EIGS, the total energy is added to data group $energy. The data groups $closed shells ($alpha shells and $beta shells for open shell cases) are no longer significant, but nevertheless kept in the control file; additionally the spinor occupations are deposited in data group $spinor.

One- and two-component all-electron calculations: relham data group

All-electron calculations require relativistic all-electron basis sets. It is recommended to use the x2c-type basis sets and RI-$J$ auxiliary basis sets x2c-XVPall for one-component and x2c-XVPall-2c for two-component treatments [182, 183, 184]. These are available for X = S, TZ, QZ for H to Rn. Alternatively, the Dyall basis sets can be used, which are available for H to Rn in the uncontracted form of the double, triple, and quadruple-$\zeta$ quality, e.g. dyall-vdz.

The input is controlled by the $relham data group which includes all possible keywords. The Hamiltonian is simply defined by x2c, bss, or dkh n (where n stands for the order of DKH, $n \geq 4$) are used to activate the X2C, BSS or DKH Hamiltonian. n defaults to 4. X2C is the default and is strongly recommended.

Local relativistic approximations can be used to increase the efficiency. Here, the DLU and DLH families are available. These are applied with dlu or dlh. By default, the DLU/DLH all scheme is selected, i.e. dlu all or dlh all. Especially, the DLU All scheme is highly accurate for energy differences [164, 165], gradients [170], NMR properties [42, 44, 185, 186], and EPR parameters [187, 188, 85]. Thus, it is recommended and the default. Further schemes, such as a neighboring block DLU scheme are available. In this DLU(NB) scheme (dlu nb), the atomic-off-diagonal blocks of two well separated atoms are approximated with the non-relativistic limit, i.e. the decoupling and renormalization matrices are assumed to become unit matrices. Its derivatives are assumed to vanish for property calculations. The distance threshold for DLU(NB) is set with rlocthr real (in bohr, default is 12) [189]. Similarly, the decoupling and renormalization matrices of light elements can be assumed to become unit matrices, this leads to the DLU(NR) scheme (dlu nr). The light elements are defined with the element number and the corresponding threshold is set with rlocelem integer. The default for rlocelem is 18 (argon). For completeness, a DLH(NR) scheme is also available (dlh nr), where the non-relativistic limit is formed for atomic diagonal blocks based on rlocelem integer. Local approximations can be turned off with rlocal off.

A picture-change correction for several expectation values and for the virial theorem in the proper section is available for relativistic all-electron Hamiltonians also in their local variant and applied with pccexpvals on/off. Default is on.

For NMR/EPR properties and geometry gradients, it is also possible to set the derivatives of the unitary transformation to zero. This is done with noresponse on/off. Default is off and this option should only be used for comparison with other programs.

In 2c calculations, the screened nuclear spin-orbit (SNSO) approximation may be used to estimate the effect of spin-orbit coupling on the two-electron interaction by using the keyword snso on/off. This applies the SNSO approximation to the relativistically modified potential and its derivatives. In case of snso 2e, we only apply the SNSO approximation for derivatives if there is also a two-electron integral derivative contribution to the given property. With snso all, it is applied to all derivatives. Alternatively, the keyword snso ham can be used to apply it to the spin–orbit parts of the Hamiltonian. The same way, snso pvp applies it to the spin-orbit part of the relativistically modified potential. The original set of parameters by Boettger was derived for lower-order DKH whereas the modified parameters were optimized for X2C and greatly improve the spinor energies for virtual states. The original parameters [167] can be applied by setting snsopara to 0. Additional options are the mSNSO DC and DCB (universal) [190] which can be applied by setting snsopara to 4 and 5, respectively. 6 applies the row-dependent DCB approach [190], where we assign row based on the nuclear charge. 7 applies the mSNSO approach of the Cremer group [191], which is the default option. Note that the original publication missed a square for the p-type parameters. [168, 169, 191] For completeness, this version is available with snsopara 1.

It is possible to further study the impact of spin–orbit coupling by scaling the spin–orbit part of the one-electron relativistically modified potential and their derivatives with the option soscal real. Default is soscal 1.0.

The X2C/BSS/DKH procedure is performed in the decontracted basis. This can easily lead to linear dependencies, which can be removed by diagonalizing the overlap matrix with sdiag on/off. By default, this is activated and eigenvectors associated with an eigenvalue of less than $1 \cdot 10^{-13}$ are removed for the transformation to the linear-independent basis. This value can be adjusted with epssmt real. The decoupling may be performed in the Cartesian atomic orbital basis (basis cao) or the spherical atomic orbital basis (basis ao). The latter is the default. Note that the CAO basis only uses snso ham and the standard snso keyword will be translated to the application on the Hamiltonian.

Point-group symmetry can be exploited in the X2C/BSS/DKH steps of the SCF and gradient runs. This is done with the option symalg on/off, default is on. Note that DLU/DLH ignore this option.

Finally, the finite nucleus model based on a Gaussian charge distribution is selected by finnuc on/off. By default, isotope number is used (finnucpara 0). The real atomic mass is used with finnucpara 1. Finnucpara is also ignored if finnuc is turned off. The finite nucleus model for the vector potential can be set with vecfinnuc on/off independently of the scalar potential. If vecfinnuc is not set, the finnuc option automatically sets it for consistency.

Taking together the default $relham data group implies

$relham
  x2c
  dlu all
  pccexpvals on
  snso
  snsopara 7
  finnuc off
  vecfinnuc off
  sdiag on
  epssmt 1.0E-13
  basis ao

You may set the group in the dedicated relham option at the end of define. 2c calculations can also be turned on there and the keyword $soghf will be added automatically and the Kramers-restricted formalism can be set ($kramers). Note that not specifing an option, will simply take the default values. Therefore, a simple line $relham is equivalent to the full specification given above.

In DFT calculations, it is further strongly recommended to use grids with an increased number of radial points [183]. These grids are selected by appending “a” to the gridsize, i.e. gridsize 4a. The modified grids such as m3 should not be used in X2C, especially for magnetic properties such as NMR or EPR.

For elements beyond Rn, there are no Hückel vectors for the initial guess. Therefore, we suggest to use the guess from a superposition of atomic densities, see Sec. 25.2.2. The basis set itself should be used for the atomic orbitals. For uncontracted basis sets or general-contracted basis sets, this option should be a very good guess. For segmented-contracted basis sets, the energy in the first iterations might not be that good.

One- and two-component all-electron calculations: Legacy Keywords

The keywords $rx2c, $rbss and $rdkh n (where n stands for the order of DKH) are used to activate the X2C, BSS or DKH Hamiltonian. n defaults to 4.

Local approximations are applied with $rlocal and the specific approximation is taken from $rlocpara n with n between 0 and 4. 0 is the default (DLU all), additionally 1 (DLH all), 2 (DLU NB), 3 (DLU NR), and 4 (DLH NR) are possible. The distance threshold for DLU(NB) is set with $rlocthr real (in bohr, default is 12) [189]. The element number for DLU/DLH NR is set by $rlocelem integer.

A picture-change correction for several expectation values and for the virial theorem in the proper section is available for relativistic all-electron Hamiltonians ($rdkh, $rbss, $rx2c, also in their local variant $rlocal) and enforced by $pcc.

Linear dependencies are removed by default based on the eigenvalues of the overlap matrix. The default value is $5 \cdot 10^{-14}$. This value can be adjusted with $x2c_lsdiag real. The transformation can be omitted by setting $x2c_lsnodiag.

For symmetric molecules, point-group symmetry is not exploited by default, but can be used in the one-component case by setting $rsym. Without setting $rsym, symmetry is only exploited in the rest of the program. Note that $rsym cannot be combined with $rlocal. The exploitation of symmetry is also supported in X2C gradient calculations [42]. Therein, nuclear symmetry is exploited without setting $rsym.

The finite nucleus model based on a Gaussian charge distribution is selected by $finnuc. By default ($finnucpara 1), the real atomic mass is used for the finite radius. The isotope number is used when $finnucpara 0 is set.

The SNSO approach is activated the keyword $snso. This applies the SNSO approximation to the relativistically modified potential and its derivatives. The parameters are selected with $snsopara n.