## Abstract

This article outlines some basic concepts of relativistic quantum chemistry and recent developments of relativistic methods for the calculation of the molecular properties that define the basic parameters of magnetic resonance spectroscopic techniques, i.e. nuclear magnetic resonance shielding, indirect nuclear spin–spin coupling and electric field gradients (nuclear quadrupole coupling), as well as with electron paramagnetic resonance *g*-factors and electron–nucleus hyperfine coupling. Density functional theory (DFT) has been very successful in molecular property calculations, despite a number of problems related to approximations in the functionals. In particular, for heavy-element systems, the large electron count and the need for a relativistic treatment often render the application of correlated wave function *ab initio* methods impracticable. Selected applications of DFT in relativistic calculation of magnetic resonance parameters are reviewed.

## 1. Introduction: where the theory of relativity and the quantum theory for NMR and EPR parameters intersect

The purpose of this review is to highlight some recent theoretical developments for calculations of magnetic resonance spectroscopy parameters in heavy-element systems from first principles. Selected applications of such methods will be discussed, with emphasis on applications of density functional theory (DFT). The aim is to illustrate rather than to elaborate. It is hoped that the reader will find this article suitable as an introduction to the background (which takes up most of the space) and a pointer to recent and on-going developments. A rather comprehensive overview of relativistic effects on atomic and molecular properties can be found in a recent review by Ilias *et al.* [1]. Molecular properties are more generally dealt with in [2], where the reader can also find information about relativistic corrections. There are many journal articles, reviews and textbooks available on the topic of relativistic quantum chemistry. For recent reviews, see [3–9]. Issues more specific to DFT are discussed in [10]. Several textbooks on the topic of relativistic quantum chemistry are also available; a selection can be found in [11–16]. Various older review articles are also recommended [17–25]. Much of the material in this article has been compiled, adapted, updated and extended (in some cases shortened) from reviews and book chapters written by me in the past few years [9,26,27].

The observed world is relativistic, meaning that the speed of light has a fixed (presumably) and finite value that represents a universal speed limit. However, non-relativistic (NR) physical theories can represent such good approximations that relativistic effects can often be neglected. Dirac originally thought that relativity would be of no importance for chemistry [28]. As it turned out, relativistic effects rapidly gain importance in the quantum theory of molecules as the nuclear charges in a system become larger—even for valence shells and ‘chemical’ properties [23]. The reader is referred to one of the cited reviews and textbooks on relativistic quantum chemistry for details as to why this is the case. For the purpose of this introduction, it suffices to state that relativistic effects on the energy and on calculated properties rise, relative to NR quantities, roughly as (*Zα*)^{2} in leading order, and in higher powers of (*Zα*)^{2} beyond leading order. Here, *Z* is the charge number of the heaviest nucleus (typically) in the system, and *α*≈1/137.036=7.29735×10^{−3} is the fine structure constant. In dimensionless Hartree atomic units (au), *α*=1/*c*, where *c*=137.036 au is the speed of light. One may therefore use (*Z*/*c*)^{2} in place of (*Zα*)^{2}. The NR limit is defined conveniently as , such that in this limit velocities can be added up to arbitrarily large values just as it is assumed in a Galilei-relativistic (=NR) physical theory. In molecules where there are typically many different nuclear charges present, it is convenient simply to use *c*^{−2} as a relativistic perturbation parameter in first and higher orders, with *c*^{−2}→0 corresponding to the NR limit. Sometimes, ultra-relativistic cases with very small values of *c* are also of interest.

The root mean square (RMS) momentum expectation value of the electron in the hydrogen atom ground state is 1 au. Given that the electron mass is 1 au, even in this light atom the speed of the electron is on average just under 1% of the speed of light, which is extremely fast compared with speeds of macroscopic objects that we observe in everyday life. Yet, (*Z*/*c*)^{2} for the H atom is only about 5×10^{−5}. Relativistic effects for this and other atoms in the upper part of the periodic table can be neglected in all but the most accurate calculations. For Hg, however, the value of (*Z*/*c*)^{2} is about 0.34, and even the next order (*Z*/*c*)^{4} reaches 0.12. Dirac's argument that relativity should not be important for chemistry was based on the assumption that valence shell orbitals in heavy atoms are kept away from the nucleus by inner shells. The order of magnitude for relativistic effects would then correspond to (*Z*_{eff}/*c*)^{2} with *Z*_{eff} being a small effective screened nuclear charge. However, the empirical finding is that it is the full nuclear charge that matters [29]. The main reason is that orbitals without orbital angular momentum (ℓ=0) have nuclear ‘tails’ that reach inside the K shell all the way to the nucleus. This causes a cascade of relativistic effects on these and other orbitals. Consequently, relativistic effects are so important for heavy elements that they may reverse periodic trends found for lighter elements in the same group, and they may even dominate the chemical and physical behaviour of these elements. Well-known examples of this type are the favoured low oxidation states of heavy p-block elements, for instance, Pb(II) versus Pb(IV), or the smaller radii of sixth-row elements, such as Cs, Ba, Au and Hg, compared with their fifth-row counterparts in the same group. The inert electron pair effect that defines much of the chemistry of heavy p-block elements, as well as the small radii of heavy elements of groups 1, 2, 11 and 12, in particular, stem from a strong relativistic stabilization and contraction of the 6s valence orbital. These and many other examples for the importance of relativistic effects in the chemical and physical behaviour of heavy elements have been compiled by Pyykkö [23,30,31]. For a brief summary of relativistic effects in atoms and molecules, see [9].

Magnetic resonance parameters fall into the class of *response* properties [2,26,32,33]. These are derivative properties that represent the response of the electronic structure of a molecule (or atom) to the presence of external or internal electromagnetic fields and other types of perturbations (e.g. nuclear displacements). The starting point of a calculation is usually the field-free (unperturbed) electronic ground state. This zeroth-order step is a single-point quantum chemical calculation with a chosen basis set and electronic structure model such as DFT or a wave function method. This step is then followed by a calculation of the desired derivative property. In the simplest case, a first-order property can be calculated simply as an expectation value using variational methods. The electric field gradient (EFG) is an example of such a property. Electron paramagnetic resonance (EPR) parameters can also be obtained from ground-state orbitals or wave functions. In the case of nuclear magnetic resonance (NMR) chemical shifts and indirect nuclear spin–spin coupling (*J*-coupling), a set of linear response (LR) equations is usually solved after the ground-state calculation in order to obtain these properties. EPR parameters have also been calculated in this fashion. A note on the nomenclature used here is in order: a first-derivative property such as the EFG can be obtained without solving response equations; the unperturbed ground-state electronic structure parameters are sufficient. However, in non-variational methods it is necessary to solve additional equations for the electronic structure parameters in order to obtain the energy derivative correctly, and additional terms need to be calculated in either case if the basis set depends on the derivative parameter(s). A second-derivative (double perturbation) property such as the NMR shielding requires solving at least one set of LR equations. By virtue of the 2*n*+1 theorem of perturbation theory, which states that an *n*th-order response of the electronic structure gives access to energy derivatives up to order 2*n*+1, one can also calculate a third-derivative property from LR quantities. Unless the perturbations are identical, one needs to solve LR equations for several perturbations. Alternatively, a quadratic response equation can be solved to gain access to a third derivative. Similar techniques exist for frequency-dependent derivative properties.

At this point, it is illustrative to introduce expressions for static first- and second-order derivative properties. One can expand the electronic energy in terms of perturbations by external and internal fields characterized by perturbation parameters *a*,*b*, etc., around the field-free case:
1.1The superscripts in parentheses indicate derivatives of the indicated order with respect to the perturbation parameters, taken at *a*=0,*b*=0,…, for the energy, and below for the Hamiltonian and the wave functions. Assuming knowledge of the exact unperturbed wave functions for the ground state (0) and excited states (*j*) with energies *E*_{0} and *E*_{j}, the energy derivatives in first and second order read
1.2aand
1.2bEquation (1.2a) expresses the Hellmann–Feynman theorem, which states that the wave function derivatives are not needed to obtain the energy perturbation in first order. For EFG calculations, the formal derivative parameter would be a point quadrupole at some position; usually the EFG is evaluated at the position of a quadrupolar nucleus [34,35]. For EPR calculations, the derivative parameter would be the component of an external magnetic field or a nuclear spin or the nuclear spin magnetic moment [36,37]. Equation (1.2b) is an example for a ‘sum over states’ (SOS) expression. The SOS enters the expression from the derivative of the ground-state wave function with respect to one of the perturbations, or , which can formally be expressed in the basis of the excited-state wave functions. The SOS is largely symbolic; usually, in approximate electronic structure calculations, response equations are solved to obtain or without explicit calculations of excited-state wave functions. Such LR techniques are also used in DFT, where excited-state densities are typically not available. With *a*,*b* representing the component of an external magnetic field and of a nuclear spin magnetic moment, respectively, *E*^{(a,b)} represents a component of the NMR shielding tensor. If *a*,*b* are components of nuclear magnetic moment vectors for two different nuclei, *E*^{(a,b)} represents a component of the *J*-coupling tensor. If one of the perturbations in *E*^{(a,b)} is spin–orbit (SO) coupling, second-derivative techniques give access to EPR parameter calculations based on spin-free (scalar) relativistic references.

For heavy-element systems, relativistic effects must be included in a response calculation. For example, some one-bond *J*-coupling constants involving Hg have been shown to increase relativistically by a factor of 2 or more [38–40]. NR theory is simply wrong in this case. By using *c*^{−2} as a perturbation parameter in first and higher orders, ‘relativity’ can be treated as an additional perturbation in response calculations based on an NR ground-state electronic structure. Relativistic effects on magnetic resonance parameters are then obtained as response properties of at least one order higher than what would be required in the NR case. In the lowest relativistic order, *c*^{−2}, in particular within a scalar framework where SO coupling is neglected, this approach offers a fast and convenient way to obtain relativistic corrections to a spectroscopic property of interest. For example, for NMR shielding, one can obtain the leading-order relativistic correction as a third-derivative quantity, where the derivative parameters are a component of the nuclear spin magnetic moment vector *m*_{N} of a selected nucleus *N*, a component of the external field ** B**, and

*c*

^{−2}. In higher orders of

*c*

^{−2}, the approach becomes more involved particularly when SO coupling is included. Ground-state calculations with a relativistic Hamiltonian have become rather fast, and computational resources are abundant, such that magnetic resonance parameter calculations based on a relativistic zeroth-order electronic structure are now commonplace. The focus of the present article on such methods where relativistic effects on magnetic resonance parameters are then automatically included—ideally at a fully relativistic level (see below). Perturbation terms in order

*c*

^{−2}may be extracted from such calculations if they are carried out with varying values for the speed of light, which is sometimes useful for purposes of analysis.

Relativistic quantum chemistry is approaching maturity, but a lot of research still needs to be done. For a long time, ‘fully relativistic’ molecular calculations based on the four-component (4c) Dirac Hamiltonian for the kinetic and potential energy (with or without suitable relativistic corrections to the electron–electron (e–e) interaction) were considered too demanding on computational resources for all but the smallest systems. Tremendous effort has consequently been put into the development of increasingly accurate hierarchies of approximate variationally stable two-component (2c) Hamiltonians (quasi-relativistic methods), and into the development of response calculations based thereupon. References are provided in §§2–7. These efforts continue, in particular to give access to the vast number of highly important spectroscopic parameters that researchers use in chemistry, physics and related disciplines. Quasi-relativistic Hamiltonians such as low-order Douglas–Kroll–Hess (DKH) or the zeroth-order regular approximation (ZORA) are now available in many program codes, which makes these methods an attractive framework for the development of response methods. In recent years, a convergence has been achieved in the sense that, on the one hand, fast four-component machinery has been developed, and, on the other, robust ways have been devised by which to construct matrix representations of exact two-component Hamiltonians (some of which are collectively referred to as X2C) as an alternative way to perform fully relativistic calculations. ‘Fully relativistic’ means here that no approximations are made in the one-electron part of the Hamiltonian, with the Dirac equation (DE) for a single electron in some chosen basis and its X2C counterpart meant to deliver the full range of relativistic effects. Relativistic corrections to the e–e interaction and corrections from quantum electrodynamics (QED) [3] may still have to be added. The development of response methods within the X2C framework and other formally exact two-component methods is fully under way. Many new research articles on implementations and tests of such methods can be expected in the coming years in addition to new developments taking place within the available four-component and approximate two-component frameworks. There are also further developments needed for a robust treatment of the e–e interaction in a relativistic framework, both in wave function and in DFT. The single-reference problem in Kohn–Sham DFT is also a source of problems when dealing with open-shell d- and f-element compounds.

The following sections contain a brief account of relativistic quantum chemistry methods, a description of how electric and magnetic perturbation operators are obtained, a discussion of picture-change (PC) effects, recent developments for magnetic resonance parameters, and selected examples for calculations of magnetic resonance parameters for heavy-element systems using recent developments, with emphasis on applications of DFT to molecules. Unless stated otherwise, dimensionless atomic units [41] are used in many of the equations, such that *m*_{e}=1, , *e*=1, 4*πε*_{0}=1 and *c*=137.03599976(50)=*α*^{−1}. The discussion is restricted to electronic motion, and the clamped nucleus approximation is adopted.

## 2. Relativistic quantum chemistry: from Schrödinger to Dirac

The time-independent Schrödinger equation (SE) for an electron in a potential *V* , in the absence of other electromagnetic fields, reads
2.1In the previous equation, is the momentum operator for the electron. The SE corresponds to the NR limit of quantum theory. The electron spin does not appear in this equation. Spin is introduced in the SE for many-electron systems in an ad-hoc manner, for instance, in order to account for the antisymmetry of wave functions with respect to electron permutations. There are also electron-spin-dependent properties, including NMR *J*-coupling and EPR hyperfine coupling, that have non-vanishing NR limits. A modification of the one-electron SE that includes the electron spin can be made, which then allows one to derive equations for such properties in the limit .

The relativistic description of an electron is furnished by the Dirac equation (DE)
2.2where ** α** is some kind of three-vector such that a dot product with the momentum operator

**can be defined. Dirac obtained the Hamiltonian in (2.2) by proposing a linearization of the Einstein energy–momentum expression 2.3followed by quantization. The energy**

*p**W*

^{rel}includes the rest-mass energy

*m*

_{e}

*c*

^{2}of the electron, whereas

*E*

^{nrel}in the SE does not, hence the different symbol. In order to establish the equality in (2.3),

*β*and the components of

**must not be just numbers. The standard representation of the Dirac matrices, as they are usually referred to, is 2.4The**

*α**α*matrices can be written in terms of the Pauli spin matrices. Let

**be a three-vector of the 2×2 Pauli spin matrices, with components 2.5A ‘split’ 2×2 block notation for the Dirac matrices and related operators, indicated by square brackets, is then 2.6where the zeros and ones represent zero and unit matrices of dimension 2. In this notation, the DE reads 2.7In order to facilitate the connection with the SE, the origin of the energy scale can be shifted to +**

*σ**m*

_{e}

*c*

^{2}. Replacing

*W*

^{rel}by

*E*

^{rel}=

*W*

^{rel}−

*m*

_{e}

*c*

^{2}and subtracting

*m*

_{e}

*c*

^{2}on the diagonal of the Dirac Hamiltonian accordingly, the DE can be written in the following form: 2.8where 2.9A useful equation for manipulation of spin operators is the following: 2.10where

**and**

*a***are vector operators (such as the momentum) acting on spatial electron degrees of freedom.**

*b*The DE has some interesting features that led to the idea of the existence of an antiparticle of the electron, the positron. The eigenvalue spectrum, in terms of *W*^{rel}, for unbound particles ranges from negative infinity to −*m*_{e}*c*^{2} and from +*m*_{e}*c*^{2} to positive infinity. Dirac had wondered whether the infinitely many negative-energy states (NESs) are all occupied (the ‘Dirac sea’ of electrons) such that the observable positive-energy electrons would be prevented from collapsing into a lower NES. What may seem like a strange idea at first glance led Dirac to realize that an amount of energy of 2*m*_{e}*c*^{2} or more would excite one of the negative-energy electrons to a level above +*m*_{e}*c*^{2}, leaving a positively charged hole behind. The process would appear as the creation of a pair of particles with opposite charge, the electron with positive energy and the hole (initially, Dirac identified the positively charged particle with the proton). Lo and behold, the prediction qualitatively explains the creation of particles and antiparticles. Particle creation and annihilation is properly dealt with by quantum field theories. This topic shall not be considered any further. Electronic states are considered by using a charge of −*e* when defining the potential energy *V* in the DE. Positive-energy as well as negative-energy solutions of the DE then describe electronic states. The DE can also be used to describe positronic states, which requires charge conjugation [5]. In chemistry, one is interested in positive-energy states of electrons. In the *E*^{rel} energy scale, energies below zero and above −2*m*_{e}*c*^{2} indicate bound electronic states of atoms and molecules.

Another important aspect of the DE is that it incorporates the electron spin from the onset. This is already indicated by the fact that the Pauli spin matrices can be used in order to construct the *α* matrices in the Dirac Hamiltonian, which means that electron spin degrees of freedom are incorporated in its eigenfunctions. The wave function in the SE, equation (2.1), is a scalar function. On the contrary, the 4×4 matrix structure of the Dirac Hamiltonian means that the solutions of the DE are four-component objects called four-spinors:
2.11Each of the spinor components *ψ*_{i} is a function of the electron coordinates. The square brackets indicate a split notation corresponding to the block notation introduced in equation (2.6). The index U is used here to group the upper two components of the four-spinor, *ψ*_{1} and *ψ*_{2}, into a two-spinor, and the index L groups the lower two components, *ψ*_{3} and *ψ*_{4}, into a second two-spinor, with the index ordering defined by how the Dirac Hamiltonian is written in equations (2.7) and (2.8). For bound electronic states in atoms and molecules, the lower components are much smaller than the upper components outside of the atomic cores. Therefore *ψ*^{U} and *ψ*^{L} are usually referred to as the ‘large’ and the ‘small’ component, a terminology that is used here interchangeably with ‘upper’ and ‘lower’ in this article. As shown below, *ψ*^{L} and *ψ*^{U} are coupled and therefore they do not vary independently. Furthermore, the components in each two-spinor have certain transformation properties under rotations in spatial and spin coordinates that differ from those of vectors of the same dimension, hence the term spinor. Early attempts to solve self-consistent field-type equations with a four-component one-electron operator in a basis set suffered from variational collapse, because the components of the four-spinors were allowed to vary independently [42]. Further comments on this issue are provided below. The eigenvalues of the Dirac Hamiltonian are stationary but obey a minimax condition [43] rather than a boundedness from below. The latter feature renders a variational treatment of the Schrödinger Hamiltonian comparatively much easier.

## 3. From four to two components

As pointed out in §2, the components of the Dirac four-spinor are coupled. From equation (2.8), one obtains explicitly
3.1aand
3.1bThe superscript ‘rel’ for the relativistic energy is dropped from here on, and are set in atomic units. The second equation gives
3.2where
3.3The fact that *k* in the previous equation depends on the energy *E*, which is not known prior to solving the DE, means that the operator *X* or its inverse cannot be constructed up-front and that it is state-dependent. In the NR limit, *X*=(1/2*c*)** σ**⋅

**. In calculations with basis functions, it has been demonstrated that, when using a basis set {**

*p**χ*} for the upper components and a basis for the lower components that includes the function set {

**⋅**

*σ*

*p**χ*}, the aforementioned variational collapse in four-component relativistic calculations can be kept under control. This is referred to as ‘kinetic balance’ [44]. Restricting the lower-component basis to only functions {

**⋅**

*σ*

*p**χ*} is referred to as ‘restricted kinetic balance’ (RKB); other schemes exist as well. RKB still allows the wave function coefficients for the upper and lower components to be different. They must become equal in the NR limit. The additional factor of

*k*in the operator

*X*requires flexibility in the lower-component basis, in particular in the atomic cores, where

*k*may differ significantly from unity. However, the construction of exact two-component Hamiltonians (see below) is facilitated by the use of RKB. Some researchers have therefore strongly advocated the use of RKB for this and various other reasons [6]. For a discussion of various kinetic balance schemes, see [45].

The fact that one of the components can be calculated once the respective other component is known suggests that it may be possible to get rid of one of the components in the formalism altogether. This would lead to a two-component relativistic framework that does not necessarily introduce any new approximations. In quantum chemistry, one seeks an electrons-only positive-energy formalism where one gets rid of the small component *ψ*^{L}. The general framework, referred to as Foldy–Wouthuysen (FW) transformation, assumes the existence of a transformation of the four-component Hamiltonian such that [46]
3.4With such a Hamiltonian, the equations for *ψ*^{U} and *ψ*^{L} are decoupled. The eigenfunction of the (still four-component) FW Hamiltonian is
3.5where the goal is to produce a wave function of the form
i.e. the lower components of the transformed wave function should vanish. This would lead to a two-component relativistic quantum equation for electrons:
3.6There are two points to consider: first, similar to the relationship between the components of the Dirac wave function, equation (3.2), a decoupling operator *X* can be used to eliminate *ψ*^{U}. Second, if one wants to get rid of one of the components of the Dirac wave function, the resulting function should be normalized properly. The Dirac wave function *ψ*^{D} is supposed to be normalized, as in
3.7with the exact coupling *X* defined in equation (3.2). However, the idea is to find a transformation that is not explicitly energy-dependent [47]. One may define normalization operators
3.8(Note that *X* and *X*^{†} do not generally commute.) The FW transformation *U* can then be written as a product of a decoupling step *U*^{d} and a renormalization step , as [5,47,48]
3.9This transformation produces a wave function
3.10As *ψ*^{L}=*Xψ*^{U} per equation (3.2), the lower component indeed vanishes. By comparison with equation (3.7), it is seen that the upper component *ψ*^{+} is correctly normalized.

For many years, researchers have tried to construct a two-component relativistic Hamiltonian from an FW transformation. For a free electron, the FW transformation is known exactly, but for bound electrons the knowledge of *X* needed to construct the transformation for general situations would require one to solve the DE. Solving the DE, however, is what one tries to avoid in two-component methods. The approach traditionally has therefore been to obtain the Hamiltonian in the form of approximations that can be made exact in some infinite limit.

For the one-electron operators in relativistic calculations, it is now possible to employ formally exact two-component operators [5,6,49–53]. One-step decoupling schemes are referred to as X2C [54]. There have been other elimination schemes devised to construct formally exact two-component relativistic operators, for instance, infinite-order methods by Barysz and co-workers (BSS=Barysz–Sadlej–Snijders, IOTC=infinite-order two-component) [55–58] and normalized elimination of the small component (NESC) methods [59–66]. In its recent incarnations, the X2C technique rests on a representation of the four-component equation in an RKB basis. The two-component operators are not derived explicitly because in the end one seeks a matrix representation in the chosen basis set anyway for quantum chemical calculations. Thus, the operators are directly constructed in matrix form. The following simplified description has been adapted from [67]; an algorithmic recipe has been given by Saue [5]. The notation for matrices used in this section is as follows. Symbols formatted as *M* indicate AO matrices in the basis of a set {*χ*} of SF basis functions (‘atomic orbitals’, AOs). Formatting as ** M** indicates 2×2 super-matrices in two-component spinor space. For matrix representations of four-component operators, a split notation for the components is used where it is convenient. To keep the notation in some places more compact, a notation like is used to indicate a 4×4 super-matrix representation of a four-component operator or transformation matrix. In this notation, the generalized four-component (Dirac) matrix eigenvalue problem reads
3.11which corresponds to the numerical solution of (2.8) in the basis {

*χ*}. The eigenvectors are collected in a matrix , and contains the eigenvalues. The AO overlap matrix is and its indices run over the basis set for the upper and lower components. Instead, one wishes to solve a two-component equation, 3.12where

**is a matrix representation of the operator**

*H**h*

^{+}of equation (3.6), and ensure that it produces only positive-energy solutions. The exact decoupling transformation in matrix form reads 3.13For an electrons-only Hamiltonian, only

*U*_{UU}and

*U*_{LU}are required. It has been realized that the so-called modified Dirac equation (mDE) 3.14is equivalent to the four-component DE when an RKB basis set is used [68–71]. The mDE is the target equation to be exactly decoupled by a matrix transformation. Here,

*C*_{U}and

*C*_{L}are separate coefficient matrices for the positive-energy solutions for the upper and lower components, respectively, which have two spin components each. The matrices 3.15are constructed from the overlap (

*S*), NR kinetic energy (

*T*) and one-electron potential energy (

*V*) AO matrices that are typically available in quantum chemistry programs. Furthermore, 3.16is a matrix with elements 3.17Using equation (2.10) with

**=**

*a*

*p**V*,

**=**

*b***, one can write 3.18(a 2×2 unit matrix is not written explicitly). In equation (3.16), W**

*p*^{0}is the AO matrix of the scalar operator

*p**V*⋅

**, and the W**

*p*^{u}(

*u*∈{

*x*,

*y*,

*z*}) are the AO matrices of the Cartesian components of the vector operator

*p**V*×

**. The calculation proceeds with the solution of the generalized one-electron matrix eigenvalue equation (3.14). From the positive-energy solutions, one can then construct**

*p***, which is the matrix that produces the lower components of the positive-energy eigenfunctions of the Dirac one-electron Hamiltonian from the upper components, 3.19The matrix**

*X***can therefore be calculated as 3.20The components of can be obtained as follows [52,67]: 3.21where 3.22A X2C two-component Hamiltonian represented in the RKB basis can then be constructed as 3.23This technique, or variants thereof, offers a convenient way to implement two-component relativistic calculations that are equivalent to calculations using the four-component Dirac operator for the one-electron part of the Hamiltonian. Regarding e–e terms, see §6. A concern may be the memory requirements to solve the four-component equation (3.14), as one of the motivations for performing two-component calculations is the lower demand for calculating and storing wave function or density matrix parameters. Recent works [6,72,73] have demonstrated that an atomic and nearest-neighbour (NN) partitioning of decoupling transformations is feasible, for the reason that the regions in space where the Dirac lower component differs significantly from (1/2**

*X**c*)

**⋅**

*σ*

*p**ψ*

^{U}are confined to the deep atomic cores, close to the nuclei. See also a review by Saue for an illustration [5]. Therefore, the memory-intensive parts of a decoupling transformation can be carried out in a reduced basis, accumulating contributions from different parts of a molecule. This approach also eliminates a potential bottleneck from the diagonalization step of equation (3.14) in large-scale computations.

## 4. Picture change

Suppose there is a four-component operator *A* for some molecular property other than the energy. If the calculation is switched to a two-component framework, either formally exact or with an approximate relativistic Hamiltonian, all other operators must be transformed accordingly, i.e.
4.1The upper left block of the transformed operator is then used together with the two-component wave function to calculate the property. This change of the operator is usually referred to as picture change (PC) [74]. The PC also applies to the position and density operators. For instance, for one electron,
is the four-component electron number density associated with the electron charge distribution. In a Hartree–Fock (HF) or DFT calculation, the density is given by similar contributions from all orbitals. The last expression is formally generated via an expectation value
where *δ*(** r**−

**′) is a density operator. In the two-component formalism, one needs to calculate where [**

*r**U*

^{†}

*δ*(

**−**

*r***′)**

*r**U*]

_{UU}is the picture-changed density operator in the two-component framework. This result is different from the two-component density The PC is a relativistic effect of the same leading order,

*c*

^{−2}, as other relativistic effects. Therefore, it can lead to large differences in calculations of properties in a two-component framework.

It is noted that the difference between *ρ*^{+} and *ρ*^{2c} is significant mainly deep in the atomic cores and to some extent in the outer atomic core shells of heavy atoms [5,75,76]. In the valence regions of light and heavy atoms, the four-component and the two-component densities are very similar. This finding has the consequence that, for example, a calculation of the dipole moment of a molecule is weakly affected by PC because the corresponding one-electron dipole operator (here in four-component form)
4.2gives strong weight to the valence regions. On the contrary, the four-component operator for the electronic EFG at some position ** R** is given by
4.3The EFG is typically evaluated at the position of an atomic nucleus, which means that the overall 1/

*R*

^{3}distance behaviour of the operator gives a strong weight to inner and outer atomic core regions. Correspondingly, PC effects can be large for the EFG [77–81]. Selected examples are discussed in §8.

In the X2C framework, property transformations are relatively convenient to apply, because the transformation is constructed explicitly. Because of the implicit RKB basis for the lower components, the operator matrices must be calculated in the full basis prior to the transformation. The dipole or EFG operators do not have spin-dependent terms. Using *A* as a general one-electron operator of this type, it is then sufficient to calculate
4.4aand
4.4bwhere the matrix elements represent the operator in the small-(lower-)component RKB basis. One then constructs
4.5The matrix elements *A*^{L,{0,x,y,z}}_{rs} have the same structure as the matrix elements in (3.17) with the potential replaced by the operator *A*. The PC-transformed operator matrix is then
4.6with the transformation matrices defined in equation (3.21). Contraction of ** A** with the two-component density matrix in the AO basis gives the PC-transformed property, whereas contraction of

*A*_{UU}with the density matrix gives the property without PC correction.

Cheng & Gauss [82] pointed out that in calculations on many-electron systems the PC-corrected property operator should also contain contributions from the decoupling transformation to first order in the field. For EFGs, the contributions were calculated to be small (scalar X2C) when compared with the magnitude of relativistic and PC effects themselves.

The PC transformation of a general four-component operator with non-vanishing UL and LU blocks requires additionally matrix elements of the type 〈*χ*_{r}|*A*(** σ**⋅

**)|**

*p**χ*

_{s}〉 and 〈

*χ*

_{r}|(

**⋅**

*σ***)**

*p**A*|

*χ*

_{s}〉 representing the operator in the mixed upper/lower-component basis sets. Of prime interest would be the operators representing the perturbations due to an external magnetic field or a nuclear spin magnetic moment. This approach requires consideration of the transformation correct to first order in the field as well. The approach is outlined in §7.

## 5. Quasi-relativistic operators

The term ‘quasi-relativistic’ is usually used for approximate two-component operators. For instance, one may envision a construction of the decoupling transformation (3.4) in powers of *c*^{−2}, which is also hinted at by the series expansion of *k* in (3.3). In the lowest order, one obtains the Pauli operator (or Breit–Pauli operator, if two-electron terms are also transformed). Using the approximations in the FW transformation,
5.1which are correct to order *c*^{−2}, the following terms are obtained from the transformation of the one-electron Dirac operator:
5.2It is important to note that in the NR limit, *c*^{−2}→0, the NR Schrödinger Hamiltonian is recovered. The derivation involves the application of some operator identities in order to bring the relativistic terms into their usual form. The first term on the right-hand side, in brackets, corresponds to *h*^{nrel} in (2.1). Using equation (2.10) with ** a**=

**=**

*b***and**

*p***×**

*p***=0, one finds that (**

*p***⋅**

*σ***)(**

*p***⋅**

*σ***)=**

*p***⋅**

*p***times a unit matrix. This identity has been used to obtain**

*p**h*

^{nrel}in (5.2) and to eliminate spin-dependent terms in some of the relativistic operators. The 2×2 unit matrix is not written explicitly, but implies that even in the NR limit the wave function carries a spin-dependent term (i.e. a spin-up or spin-down spin function or a linear combination thereof). The first term after

*h*

^{nrel}is called the mass–velocity operator, indicating its relation to the relativistic mass increase of the electron at high velocities. The operator can also be obtained from an expansion of the square root in equation (2.3) in orders of

*c*

^{−2}followed by quantization. The second term in (5.2) is called the Darwin term, representing a correction to the potential. The third term, containing the Pauli spin matrices, is the one-electron SO coupling operator. SO coupling is responsible, for example, for splitting the 2(2ℓ+1)-fold degeneracy of AOs with non-zero angular momentum ℓ (the factor of 2 counting both spin-up and spin-down orbitals) into 2ℓ degenerate levels with total angular momentum and 2ℓ+2 levels with . Unfortunately, the construction of the FW transformation in orders of

*c*

^{−2}creates operators with pathological features. For instance, for a point nucleus

*p*

^{2}

*V*gives a Dirac delta (

*δ*) distribution, and the

*p*

^{4}term as well as the SO operator are also numerically very problematic. The situation gets worse in higher orders [83,84]. The Pauli operator cannot be used in variational calculations because it leads to variational collapse. The Pauli operator is suitable for perturbation calculations in order

*c*

^{−2}. This has been done extensively, but because of the undesirable properties of this operator, and of the perturbation terms obtained in the presence of external fields, a more in-depth discussion is skipped. See [11,37,85,86] for details. The operator remains in use, even for calculations of heavy-atom chemical shifts [87].

‘Direct’ or ‘Dirac’ four-component perturbation theory (DPT) avoids the problematic features of the Pauli Hamiltonian. After a change of the metric between upper and lower components in the DE, the resulting four-component equation is straightforward to expand in powers of *c*^{−2}. The approach leads to non-singular first- and higher-order expressions for the energy and the wave function [88–90]. A treatment of magnetic properties within the DPT framework has been formulated [91].

Another approach constructs the transformation in orders of the potential to arrive at a variationally stable Hamiltonian. The DKH sequence [92,93] starts with the free-particle FW transformation and decouples the wave function components in a sequence (DKH1, DKH2, etc.) of increasing accuracy in the presence of a potential. Regular approximations (RAs) [94] are motivated by the series expansion of *k* in equation (3.3) but attempt expansions such that the series is convergent and the resulting operator is variationally stable. For a point nucleus, for instance, *V* in equation (3.3) can adopt arbitrarily large magnitudes close to the nucleus, causing the series to diverge. Ultimately, this is related to the undesirable properties of the Pauli operator.

The procedure leading to the DKH2 Hamiltonian can be outlined as follows. The transformation is written as *U*=*U*_{0}*U*_{1}*U*_{2}…. Application of the free-particle FW transformation *U*_{0} to the Dirac Hamiltonian with an external potential gives [5,15]
5.3with
5.4For a vanishing potential, *V* =0, the second term on the right-hand side of (5.3) vanishes. In this case, the Hamiltonian is indeed block diagonal. For non-vanishing *V* , the goal is then to block-diagonalize the second term as well. The operators and in (5.3) are
5.5and
5.6The notation indicates an ‘even’ and an ‘odd’ operator. Traditionally, in the DKH framework *U*_{1} has been parametrized as [15]
5.7where the anti-Hermitian operator *W*_{1} is to be determined. Other types of expansions, such an exponential parametrization for a unitary matrix, can also be used [16]. The next-order DKH-transformed Hamiltonian is [15]
5.8*E*_{p} is of order *c*^{2}, which is the highest power of *c* in the expression. The term [*W*_{1},*βE*_{p}] is selected to cancel , which gives an equation for *W*_{1} as
5.9After elimination of , the remaining terms in *h*_{2} are of second or higher order in the potential. Truncation to terms of second order in the potential gives the second-order Hamiltonian,
5.10The transformed operator is then used as an approximate relativistic operator in two-component calculations. Note that because of (5.9), *W*_{1} is odd and therefore the last term of (5.10) is even, i.e. the Hamiltonian has the desired decoupled form (apart from the neglected terms of higher order in *V*). The DKH transformation involves momentum-dependent terms with square roots that are not straightforward to evaluate in position representation. Even in momentum space the integrations would be very tedious. Hess made the important realization that because the terms of concern (see equation (5.4)) contain *p*^{2}, they can be evaluated easily in matrix form in a basis that diagonalizes *p*^{2} [95–97]. Afterwards, the matrices can be transformed back to the original AO basis. Similar techniques have also been used for elimination of the small component, such as the approach developed by Nakajima & Hirao [98]. The construction of the DKH operator has been carried out essentially to arbitrary order [99,100].

Another variationally stable quasi-relativistic operator that is in widespread use is obtained by rewriting *k* of equation (3.3) as
5.11The term in parentheses can then be expanded in a power series in *E*/(2*c*^{2}−*V*), which has the advantage that the expansion parameter is small as long as *E* is small compared with 2*c*^{2}, even close to a point nucleus. The ZORA [94,101] represents, as the name says, the zeroth-order term of this expansion. Using *X*≈(1/2*c*)[2*c*^{2}/(2*c*^{2}−*V*)]** σ**⋅

**in the FW transformation, and ignoring the renormalization (i.e. using ), the resulting two-component operator reads 5.12a 5.12b with 5.13The NR limit is obtained from . The ZORA operator is dependent on a change in the origin of the energy scale (gauge dependence), via the potential**

*p**V*in of equation (5.13), such that

*V*+Δ does not translate into an eigenvalue

*E*+Δ. The available implementations usually make use of model potentials in (for instance, sums of atomic potentials, or the potential from a sum of atomic densities plus the sum of nuclear potentials) [102–104], which bypasses many of the gauge-dependence problems in practical applications and simplifies the calculation of the matrix elements of the operator.

The ZORA Hamiltonian misses some contributions in order *c*^{−2}. Going back to the expression for *X* in equation (3.2), ZORA amounts to neglecting the constant *E* in *k*. The approximation is justified for valence orbitals in heavy atoms, for which *V* can be extremely large locally but *E* is small. For core orbitals with large *E*, the error inherent in ZORA is substantial. In the very many applications of ZORA to date, it has been shown that for valence-shell properties in heavy-element systems ZORA is quite accurate. This includes ‘nuclear’ properties such as NMR chemical shifts or *J*-coupling, but not the absolute shielding [105]. A convenient feature of ZORA is that it is straightforward to derive magnetic perturbation terms at the operator level and calculate the corresponding matrix elements by numerical integration or a mix of numerical and analytical techniques. Consequently, many important magnetic properties have been implemented and applied with ZORA for a long time [26].

A ‘scaled’ ZORA method [106] eliminates most of the gauge-dependence errors and gives much improved energies for the core orbitals in heavy atoms. With the approximations *X*≈*X*_{0}=(1/2*c*)[2*c*^{2}/(2*c*^{2}−*V*)]** σ**⋅

**and , one obtains a first-order regular approximation (FORA). The FORA expressions have been used to derive the scaled-ZORA equations more rigorously. Dyall & van Lenthe [107] devised an infinite-order regular approximation (IORA), which takes the square normalization operators in of equation (3.9) fully into account, but leaves the ZORA approximation intact in the operator**

*p**X*. Therefore, IORA is still an approximate relativistic Hamiltonian. However, the IORA Hamiltonian includes all relativistic terms in order

*c*

^{−2}as well as some higher-order terms. The technique has been applied successfully in magnetic property calculations [40].

## 6. Many-electron systems

A many-electron Hamiltonian can be written as
6.1where *k*,*l* are electron labels. In the NR case the Hamiltonian includes the instantaneous electron–electron (e–e) repulsion,
6.2This instantaneous repulsion is not consistent with a relativistic framework where particle interactions are retarded according to the finite speed at which the interactions can take place. Consequently, there are relativistic corrections to the e–e interaction. The relativistic interaction between two classical charged particles (see the textbooks cited in the Introduction for details) can be quantized by using the relativistic velocity operator *c*** α**. In Coulomb gauge, the four-component operator reads
6.3The second term on the right-hand side is a magnetic interaction known as the Gaunt term. The sum of the Gaunt term and the third term is referred to as the frequency-independent Breit operator. The operator may be considered as an approximation of the interaction between two electrons calculated from QED [108]. There is a more general frequency-dependent version as well.

For two-component relativistic calculations, the e–e interaction needs to be transformed to two-component form. In order *c*^{−2}, the corresponding terms are well known [2,109]. The leading-order relativistic corrections to can be written as
6.4a
6.4b
6.4c
6.4d
6.4e
When added together with to the one-electron terms in the Pauli Hamiltonian (5.2), the operator is referred to as the (external field-free) Breit–Pauli Hamiltonian. Operators (6.4a) and (6.4b) are two-electron analogues of the Darwin and SO one-electron operators in (5.2), with the electron–nucleus potential replaced by the NR potential for the e–e interaction. These operators arise from the transformation of to two-component form. The operator (6.4c) is the spin—other orbit (SOO) two-electron SO operator (as opposed to the spin—same orbit (SSO) term in (6.4b)) and has its origin in the relativistic corrections to the e–e interaction in equation (6.3). The operator (6.4d) represents a magnetic electron-spin–electron-spin dipolar interaction that contributes to zero-field splitting in open-shell systems but is often neglected. The last term, (6.4e), is a two-electron orbit–orbit interaction.

For transformations of two-component form (such as X2C, DKH or ZORA) that depend on the potential *V* , the e–e interaction represents a difficulty because the potential interaction from the electrons should be included from the beginning. A relatively straightforward approximation is the use of a model potential [51,110,111] for the e–e interaction, for instance, from relativistic calculations on atomic fragments of a molecule. As pointed out by van Wüllen [110], the potential must be calculated from the correct electron charge density and therefore needs to correspond to the four-component density of an atomic fragment. Otherwise, additional PC corrections would be required. If atomic nuclear and Coulomb potentials are used to set up the model potential, one obtains mean-field approximations to the 2eDar and SSO terms (in their respective X2C, DKH, ZORA, etc., forms). The SOO contribution in Breit–Pauli form is often subsequently added to the two-component Hamiltonian. Atomic mean-field integrals (AMFI) are in widespread use for calculations of two-electron SO terms [112].

The four-component Hamiltonian obtained with *h*^{D} for the one-electron part is the Dirac–Coulomb–Breit (DCB) Hamiltonian, or Dirac–Coulomb (DC) Hamiltonian if the NR e–e interaction is used. DCB serves as a starting point for correlated electronic structure calculations. Additional QED corrections may be added [3]. There has been some debate whether these and the Breit term [17,113,114] may be used in variational calculations. In order to account for effects due to the Breit term on atomic and molecular properties other than total energies and orbital energies, a variational treatment would be preferred over a perturbative scheme [113]. The DCB Hamiltonian is in general not Lorentz-invariant, but it has been said to provide an excellent approximation to the full theory [22].

When correlated wave function methods known from NR theory are transferred to the relativistic realm, some new problems arise. When using a basis of Slater determinants for a correlated wave function, in the DC(B) framework there is the possibility of having positive and negative continuum one-particle states in an infinite number of determinants that are degenerate with a determinant formed from bound states. As a result, the DCB Hamiltonian has only continuum solutions (Brown–Ravenhall disease) [115]. A ‘no-pair’ Hamiltonian has been introduced as a cure for the disease, where one-particle NESs are excluded. As discussed by Saue [5], there are important considerations as far as a full treatment of the correlation is concerned, as the no-pair technique may eliminate some of the flexibility contained in the one-particle basis. There is also the issue that the NESs depend on the potential. Liu [116] has pointed out that neither a configuration nor a Fock space formulation is fully satisfactory for a treatment of NESs in the context of correlation. He advocated for a no-pair plus QED approach as the proper tool for relativistic molecular quantum mechanics. Intrinsic errors of the no-pair approximation may be removed by accounting for some QED one-body terms involving NESs, which would give a potential-independent no-pair approximation. For related discussions, see some earlier studies [4,117].

DFT is often called the workhorse of numerical quantum chemistry. DFT offers a comparatively efficient way of introducing electron correlation into electronic structure calculations, along with various unwanted side effects (an example of the old adage that there is no such thing as a free lunch) that tend to dampen the enthusiasm about DFT among practitioners. In the NR domain, the energy functionals are functionals of the electron charge density (*ϱ*=−*eρ*, with *ρ* being the electron number density associated with the charge distribution). For open-shell systems, approximate functionals are dependent on the spin densities *ρ*^{α} and *ρ*^{β}, for practical reasons (i.e. approximate density-only functionals do not work sufficiently well in this case). There are relativistic versions of the Hohenberg–Kohn theorem and the Kohn–Sham approach [118–121]. In relativistic theory, electromagnetic interactions are related to a four-current
where ** j** is the current density. The electromagnetic four-potential is
with

*ϕ*being the electric potential, and

**the vector potential. The corresponding scalar electromagnetic interaction is**

*a**ϕϱ*−

**⋅**

*a***. The four-current has a spin-dependent part that can be identified with a spin density in the NR limit. Therefore, one may seek to develop relativistic exchange–correlation (XC) functionals that depend on the four-current. However, such functionals are not available. Instead, NR spin-density functionals are routinely employed in relativistic DFT calculations. The electron and spin densities used in the functionals are constructed from relativistic four-component or two-component orbitals. The reader is referred to a recent book chapter by van Wüllen where the principles of constructing relativistic XC functionals are explained and major issues are highlighted [10]. The reader is also pointed to an article describing the implementation of a four-component DFT code (‘BERTHA’) that touches upon the issue of four-current functionals [122], and literature on suitable definitions of spin density in relativistic DFT methods [123,124]. One point worth emphasizing is that relativistic corrections to the e–e interaction, i.e. the Breit term of equation (6.3), are supposed to come from relativity included in the functional. Another point is that, within the framework where electron motion is treated in a fixed nuclear reference frame, a relativistic XC functional should only depend on the charge density. However, for open-shell systems, one would then run into similar problems as in NR DFT, namely that information on the distribution of spins should be used in an approximate functional in order to obtain reasonable energetics. In the absence of external vector potentials, the neglect of functionals that depend on the current has been argued to follow from neglecting the Breit interaction [125].**

*j*## 7. Magnetic field-perturbed operators

Electromagnetic fields are usually treated semi-classically in molecular quantum mechanics (i.e. the fields are not quantized). ‘Minimal substitution’
7.1is a way to incorporate electric potentials *ϕ* and magnetic vector potentials ** A** in a quantized equation for a particle with charge

*q*. For an electron,

**→**

*p***+**

*p***when atomic units are used. For an external static homogeneous magnetic field, 7.2aand for the field from a point-nucleus spin magnetic moment**

*A*

*m*_{N}, 7.2bIn the previous equations,

*R*_{N}is the position of the nucleus, and

*r*_{N}=

**−**

*r*

*R*_{N}. Furthermore,

**is a chosen gauge origin for the vector potential of the external field, which does not need to coincide with the origin of the laboratory coordinate system. For**

*G*

*A*^{N}the position of the nucleus is a natural gauge origin. The corresponding magnetic fields are obtained from

**=∇×**

*B***. The definition leaves a choice of gauge for the vector potential since**

*A***+∇**

*A**f*(

**) with an arbitrary scalar function**

*r**f*gives the same magnetic field. Typically, one uses Coulomb gauge: ∇⋅

**=0, which is satisfied by (7.2a) and (7.2b). However, this does not fix the gauge origin**

*A***. The gauge freedom has the unpleasant consequence that quantum chemical calculations of magnetic properties with incomplete AO basis sets lead to an unphysical origin dependence of the calculated results. The problem is usually circumvented by the use of distributed gauge-origin methods such as the adoption of gauge-including AO (GIAO) basis sets.**

*G*The magnetic operators arising from the substitution (7.1) in an NR one-electron Hamiltonian with spin (using instead of ) read
7.3The notation {⋯ } indicates that the derivative (from ** p**) is taken only in the operator, not with functions that the operator acts upon. The magnetic-field-dependent part of the NR Hamiltonian affords a ‘diamagnetic’ term proportional to

*A*

^{2}, and ‘paramagnetic’ terms that are linear in

**. In Coulomb gauge, the electron-spin-free paramagnetic operator works out to be . For the vector potential of a homogeneous external field, equation (7.2a), this operator gives the usual orbital Zeeman interaction.**

*A*In the one-electron Dirac Hamiltonian, minimal substitution gives the following four-component magnetic field term:
7.4A diamagnetic contribution to a second-order derivative property such as NMR shielding or *J*-coupling is implicit in standard formulations of four-component theory despite the apparent lack of a diamagnetic operator. Originally, diamagnetic contributions were traced back to contributions from negative-energy eigenfunctions of *h*^{D} if one were to express the shielding or *J*-coupling in an SOS expansion [126]. The role of the NESs for diamagnetic shielding has also been pointed out by Pyykkö [127]. The four-component nature of the magnetic perturbation mixes upper and lower components of occupied and unoccupied states in an SOS expression of NMR parameters. The NESs can produce large contributions despite the large energy gaps in the denominators in the SOS because their ‘small’ components are of order *c*^{0}. For additional references and further discussion, see also references [15,128,129]. When a matrix formulation or suitable transformations at the operator level are employed, diamagnetic terms can be made to appear explicitly in a four-component formalism, which is the preferred way for implementations at this level of theory [130–132]. A recent article has summarized various approaches and the inclusion of magnetic balance schemes in relativistic calculations of NMR parameters with HF and DFT, along with a discussion of how to obtain these properties from correlated MP2 calculations [133].

In many-electron systems, relativistic corrections to the e–e interaction afford momentum-dependent terms, as can be seen from equations (6.4). Inclusion of magnetic-field-dependent terms in the formalism via (7.1) in the DCB Hamiltonian correspondingly gives rise to a number of field-dependent two-electron terms. The operators in order *c*^{−2} can be found, for instance, in the book by Harriman [37].

In four-component calculations of magnetic-field-dependent properties with basis sets, it is important to consider the concept of ‘magnetic balance’ [126]. As mentioned above, it is important to reflect in the basis set the coupling between upper and lower components of the Dirac wave function resulting from the *c* ** α**⋅

**term in the Hamiltonian, which is taken care of by kinetic balance. In the presence of magnetic fields, the lower-component basis set should consequently contain {**

*p***⋅**

*σ*

*A**χ*} functions in addition to the {

**⋅**

*σ*

*p**χ*} set [91]. As in the case of kinetic balance, the term ‘restricted’ indicates a lower-component basis that consists of just this set of functions, whereas ‘unrestricted’ indicates that additional functions may be present in the basis. Implementations of four-component NMR shielding tensor calculations with restricted magnetic balance (RMB) were reported not long ago [134–137], including extensions to using gauge-including (magnetic-field-dependent) basis functions [135,136]. Olejniczak

*et al.*[138] have come up with a simplified magnetic balance scheme for four-component magnetic property calculations, in which the flexibility of unrestricted magnetic balance can be combined with the advantages of a formalism that assumes restricted magnetic balance.

The ZORA framework has one-electron magnetic perturbation operators that can be written to resemble those of the NR framework:
7.5The NR limit (7.3) is obtained by setting . It is important to note that, despite the apparent similarity to (7.3), the operators differ in crucial ways. For example, for *A*^{N} the derivative in the electron-spin-dependent term in (7.3) generates a Dirac delta distribution; the corresponding operator is the famous Fermi contact term (FC). The presence of in the ZORA relativistic version suppresses the contact term and causes the operator to sample regions with very small *r*_{N} instead.

The treatment of magnetic terms in the DKH operator has taken comparatively long to develop, compared with ZORA. The expressions in second order (DKH2) are not particularly compact and therefore are not repeated here. The reader is referred to references [129,139–146] for details. The formalism has recently been extended to infinite order [147], whereby a recursive equation is solved for the decoupling operator following the free-particle FW transformation (§5). As the DKH transformations involve operators that contain functions of *p*^{2}, the transformations should in principle include the vector potential from the onset [148] (because of (7.1)). At the DKH1 level [141], an expression has been given for the one-electron magnetic perturbation operator that structurally resembles that of the ZORA framework:
7.6As in other two-component relativistic Hamiltonians, the DKH perturbation terms contain ‘kinematic factors’ that modify the operator, when compared with its NR counterpart, mainly in the core regions and near the nucleus of a heavy atom. This is also true in higher orders of the DKH expansion. A diamagnetic operator appears at the DKH2 level [141].

The developments for constructions of formally exactly decoupled two-component Hamiltonians (X2C) have led to formulations of how to construct X2C magnetic-field perturbed operator matrices, along with a pilot application to NMR shielding [149,150]. See also previous work on NMR parameters with formally exactly decoupled two-component methods [143,151], including the infinite-order DKH mentioned above [147]. The X2C procedure is only sketched here. The reader can find the corresponding matrix expressions in the relevant references. The calculation starts from an X2C ground-state calculation as outlined in §3. The field-free X2C Hamiltonian matrix can be determined from the solution of the generalized matrix eigenvalue equation (3.14), which provides the matrix that couples the upper and lower components of the wave function or of a set of molecular orbitals (in HF theory and in Kohn–Sham DFT). Details depend on the chosen kinetic (and magnetic) balance scheme. The procedure suggests that one needs to solve a first-order perturbation equation for (3.14) in order to obtain the matrix elements for the perturbation Hamiltonian. In a nutshell, that is how the calculation can be set up.

Sun *et al*. [149,150] specifically considered two perturbations, as in equation (1.2b), related to magnetic fields. The perturbations are described by vector potentials
as in equations (7.2a) and (7.2b). The perturbation parameters *a* and *b* are then vector components of the external magnetic field and of the nuclear spin magnetic moment, respectively. The energy *E*^{(a,b)} in second order represents a component of the NMR shielding, *J*-coupling or the magnetizability tensor. The ingredients needed to calculate *E*^{(a,b)} in a DFT framework are as follows: *E*^{(a,b)} is written as a sum of contributions from each occupied orbital. One needs to solve coupled-perturbed (CP) equations for one of the perturbations, say perturbation *a*, for each orbital. A suitable recipe might be the following:

Solve the unperturbed equation (3.12) within the X2C framework.

Solve CP equations for

*C*^{(a,0)}. The CP equations contain a perturbed Fock matrix*F*^{(1,0)}, which, in turn, contains the perturbed X2C operator*H*^{(a,0)}. One also needs, among other quantities, , which is the perturbation of of equation (3.22).*H*^{(a,0)}is constructed from the magnetic perturbations of the Dirac operator, equation (7.4), subject to transformation to two-component form, from various unperturbed quantities, and from*X*^{(a,0)}. The latter quantity is also needed to construct and can be obtained from a perturbation expansion of the equation that defines.*X*

It has been pointed out [149] that there is an interdependence of *X*^{(a,0)} and *C*^{(a,0)}, which may in principle cause the X2C calculation to be more expensive than a four-component one. However, it was found that there are a number of very good approximations under which *X*^{(a,0)} can be calculated at the outset without further iteration [149]. The working equations depend, for instance, on how a magnetic balance criterion is worked into the formalism. Liu *et al.* recently set up a flexible scheme for working kinetic (k) and magnetic (m) balance as well as GIAO gauge factors (g) attached to the basis functions into the formalism, by writing the Dirac orbitals in the AO basis {*χ*} as
7.7where each of the *Z* operators has a 2×2 block structure,
7.8The orbital coefficients and *C*^{L}_{ri} have two spin components each, as in equation (3.14). The scheme allowed for NMR shielding calculations using three set-ups to achieve the decoupling transformation, two kinetic balance schemes, and five flavours of magnetic balance [149].

## 8. Selected examples

The abbreviations 4c and 2c for four-component and two-component relativistic calculations, respectively, are used in this section. Before turning to NMR and EPR parameters, it is interesting to mention a few general applications of relativistic calculations. For instance, it is now textbook knowledge that, because of relativity, gold is yellow coloured and mercury is liquid at room temperature. The low melting point of mercury has for a long time been attributed to relativistic effects [17,152]. The argument makes intuitive sense because the Hg 6s valence shell is chemically rather inert owing to its strong relativistic stabilization. It is therefore likely that the low melting point of Hg metal, −39^{°}C, compared with 321^{°}C for Cd, is a consequence of relativity. Recent developments of Hg interatomic potentials [153] may soon demonstrate this via first-principles simulations of the metal's melting point. The closed-shell interactions between Hg atoms are difficult to treat accurately by theoretical methods, and DFT falls short in terms of reliable performance. It is possible that recent developments on functionals incorporating long-range correlation [154] or range-separation techniques with long-range correlation from wave function approaches [155] will help DFT to catch up with correlated wave function methods for this challenging case.

DFT calculations were recently used to show that relativistic effects contribute 1.7–1.8 V out of the 2.0 V of a standard lead–acid battery cell [156]; the article has received much attention. The effect has been traced back to the oxidative power of solid PbO_{2}, which, in turn, is caused by the exceptional relativistic stabilization of the Pb 6s shell.

For scalar (SF) X2C and NESC techniques, calculations of analytic nuclear gradients and second derivatives, along with other derivative properties, have been implemented [82,157–160]. These methods therefore give access to the traditional exploration of potential energy surfaces in studies of reaction mechanisms that make up a vast portion of applications of computational chemistry methods.

The examples in the remainder of this section focus on applications to EFGs and magnetic resonance parameters. A recent review by Helgaker *et al.* discusses wave function-based methods and applications thereof for a large variety of molecular properties, including NMR and EPR parameters [2]. Furthermore, as stated above, the examples are meant to highlight recent developments. A few applications of established methods are also included. For references to older literature, the review articles cited in the Introduction are recommended (among the author's publications see references [9,26,27]). Reviews covering NMR and EPR calculations can also be found in earlier studies [161–167]. The articles cited in the following subsections of course also contain relevant references to prior art.

### (a) Electric field gradients

Electric properties in 4c calculations can be calculated from the 4c density without PC errors. PC errors in quasi-relativistic calculations of properties are more difficult to avoid for methods where the transformation to two-component form is not explicitly available. A comparison has been made between different exactly decoupled two-component methods, *viz*. high-order DKH, BSS and X2C [54]. It is noted in passing that the different 2c methods worked equally well and required comparable computational effort. It was shown that PC corrections are mandatory for inverse radii expectation values and contact densities in heavy atoms and small molecules with heavy elements. This section emphasizes EFGs because they are subject to sizable relativistic effects in general [168], and to substantial PC effects in particular. In recent years, a considerable number of papers have been published where PC-corrected EFGs were obtained with two-component methods.

The canonical set of molecules used for testing new EFG implementation is the series of hydrogen halides, HX, with X=F, Cl, Br, I, At. Relativistic and PC effects of the EFG at the X nucleus grow rapidly as X becomes heavier. Mastalerz *et al.* have studied relativistic and PC effects on EFGs in molecules using DKH transformations of increasing order in a scalar relativistic (spin-free, SF) framework [81]. Results for the HX series are provided in table 1. The EFGs are in reasonable agreement with four-component HF data, which include SO coupling. However, for the heaviest element, At, the deviations approach −0.8 au. The basis sets used for the calculations are close to saturation and have therefore been adopted in several more recent studies.

Cheng & Gauss [82] implemented a scalar X2C approach and techniques for first- and second-order derivative properties [82,157]. Representative EFG data for the HX series are listed in table 1 along with SF 4c data, EFGs that are correct to order *c*^{−2}, and NR results that were reported in the same paper. EFGs have also been calculated within a scalar NESC framework [169] for X in the HX series and for Hg in 11 complexes with mercury, along with Hg contact densities for the latter. Excellent agreement between SF four-component, DKH7, X2C and NESC calculations is demonstrated for the popular HX benchmark set, as seen from the results provided in table 1. The NESC approach was subsequently extended to study finite-nuclear-volume effects on nuclear quadrupole coupling constants (NQCCs) [170]. The ratio of NQCCs obtained from different isotopomers of a molecule may differ from the ratio of the corresponding bare-nucleus quadrupole moments, which has been dubbed ‘quadrupole anomaly’. The calculations based on a Gaussian nuclear model predicted that the effect due to different isotope radii may reach magnitudes of 0.1% in some Au compounds, which may render it experimentally detectable.

We have recently implemented X2C with SO coupling and applied it to calculations of EFGs [67]. Representative data for the HX series are included in table 1. The SO data are in reasonable agreement with 4c reference values and improve notably upon the scalar relativistic results for HAt. In addition to X2C transformations using only the nuclear potential, calculations with an NR model potential (mp) for the two-electron terms were also implemented. For HF to HI, the inclusion of the model potential gives an improvement towards the 4c reference data. For HAt a slight overestimation with respect to DC is in part due to the use of finite nuclei for the latter. Overall, in comparison with the NR data, it is seen that relativistic effects on heavy elements range from significant to essential. For comparison, the dataset in table 1 also includes X2C results that were obtained without a PC transformation of the EFG operator (see §4). The EFGs for heavy elements are clearly not reliable. For uranyl, UO^{2+}_{2}, and a uranyl–carbonate complex, it has been shown that the PC corrections can be as large in magnitude as the EFG itself [67].

Most of the PC errors in order *c*^{−2} in ZORA EFG calculations can be eliminated by explicitly reconstructing an approximate 4c density [171,172]. The technique has been dubbed ZORA-4 [173]. When the EFG is instead obtained as a derivative of the scaled-ZORA DFT energy expression, one ends up with the same ZORA-4 working expressions [171]. While the operator transformation and the reconstruction of the 4c density should formally give equivalent results, the latter tends to be more convenient in conjunction with numerical integration (which is typically used in the ZORA framework), while the former tends to be more convenient for analytical techniques such as high-order DKH, BSS, NESC or X2C.

SO effects on the EFG tend to be small, with some exceptions. The iodine EFG in the diatomic molecule Tl–I exhibits a relatively large SO effect when compared with the overall magnitude of the EFG and its relativistic effects. ZORA-4 CAM-B3LYP DFT calculations yielded 2.06 (nrel), 2.37 (scalar) and 2.61 (SO) au [171]. The experimentally derived EFG is 2.68(5) au. The iodine EFGs in group-13 iodides were also calculated with scalar and SO X2C and the B3LYP functional [67]. For TlI, the scalar and SO EFGs obtained with a model potential in the decoupling transformation were 2.34 and 2.76 au, respectively. The remaining deviations from experiment are probably due to approximations in the density functional and basis set incompleteness. The PC-uncorrected SO iodine EFG for TlI was 2.96 au.

### (b) NMR nuclear magnetic shielding and *J*-coupling

Sun *et al.* [149] recently provided NMR shielding benchmark data comparing four-component DFT with various flavours of X2C, demonstrating excellent agreement. Representative data from the study for HI, I_{2}, HAt and At_{2} are provided in table 2. The small deviations are probably within the numerical precision of the calculations. A large collection of benchmark data for the HX series published up to 2008 can be found in Autschbach & Zheng [27] for comparison.

A 4c technique for NMR shielding tensor calculations using restricted magnetic balance for the small component has been developed by Komorovsky *et al.* [134]. The method was subsequently extended for the use of GIAO basis functions in order to deal with the spurious gauge-origin dependence in finite-basis calculations discussed in §7 [135]. Results for ethene and the HX series with non-hybrid density functionals were reported in [135]. The article specifically addressed the extent of the gauge-origin problem, which, as in other works on NMR shielding constants, was found to be severe. Adopting a GIAO basis was shown to solve the problem fully, as it does in NR theory.

Aucar *et al.* [163] have put forward a relativistic theory for polarization propagators as a framework for correlated relativistic calculations of NMR shielding and *J*-coupling. The three lowest levels of the approach have been designated as PZOA (pure zeroth-order approximation), FOPPA (first-order polarization propagator approximation) or RPA (random phase approximation), and SOPPA (second-order polarization propagator approximation). The first two, PZOA and RPA, have their correspondence in uncoupled-perturbed HF theory and in coupled-perturbed HF theory, respectively. At the SOPPA stage and beyond, the scheme starts to include dynamic electron correlation. The SOPPA approach has been available for some time for NR NMR calculations [174]. In recent years, a number of articles have reported relativistic NMR data at the RPA level. An example is a study of the heavy-atom shielding in the systems SnXH_{3} (X=H, F, Cl, Br, I), SnXYH_{2} (X, Y=F, Cl, Br, I) and PbXH_{3} (X=H, F, Br, I) [175]. Comparisons were made with a linear-response ESC method (LR-ESC) [176]. The two methods agreed well for Sn shielding constants. The LR-ESC method was used to analyse various mechanisms that contribute to the shielding. A similar set of computational tools was applied by Melo *et al.* to investigate various mechanisms that contribute to relativistic effects on the carbon shielding in halide-substituted methanes [177]. An SO contribution to the carbon shielding of as a function of the nuclear charge of a heavy halogen substituent X was determined from fits. NMR shielding and *J*-coupling constants for HX and halide-substituted silanes and stannanes were investigated in [178]. Calculations on the FX series (X=F, Cl, Br, I, At) found a dependence of the SO effects on the X shielding on the nuclear charge *Z*_{X} of the X atom that can be fitted to *Z*^{y}_{X}, with the exponent *y* around 3.5–3.65 [179]. In this work, the authors also addressed the correspondence between the paramagnetic shielding and the nuclear-spin–molecular-rotation coupling constant [180]. The latter point was more thoroughly investigated in a follow-up study [181]. Within the applied approximations, it was found that the analogy between the two properties does not extend into the relativistic regime. The spin-rotation tensor is less strongly affected by relativistic corrections than the NMR shielding tensor. See also [182] regarding a proposed revision of the absolute shielding scale of ^{119}Sn.

In NMR calculations, one typically solves equations to obtain the perturbed orbitals or wave functions to first order in the external field. The advantage is then that the chemical shifts for all nuclei can be computed subsequently in one sweep (although it is just as valid to calculate the result from orbitals that are perturbed by a nuclear spin, or from both perturbed sets). The magnetic-field-perturbed orbitals can be used to calculate the induced current density, which has found chemical applications in deciding if a molecule is aromatic. Sulzer *et al.* [183] have presented a 4c method for GIAO calculations of induced ring currents and applied the method to the series C_{5}H_{5}E, where E=CH, N, P, As, Sb, Bi. It was found that all molecules should be considered as aromatic, but decreasingly so as E becomes heavier.

Many examples of interesting relativistic effects on NMR shielding constants and chemical shifts have been documented in the literature. See [27,162] for references to older literature. SO effects on chemical shifts of light ligand atoms have long been investigated in transition metal complexes. Hrobarik *et al.* have recently applied 4c DFT calculations in order to study proton chemical shifts in a number of hydride complexes [184]. The samples included the series [HMCl_{2}(PMe_{3})_{2}] with M=Co, Rh, Ir. As an example, for the Ir system, a large SO-induced proton shielding of about 30 ppm was observed. The correlation between calculated and experimental shifts for the whole range of complexes was substantially improved upon inclusion of SO effects. This ‘real life’ application demonstrates that 4c DFT NMR methods have become computationally efficient and robust to the point that routine applications are possible.

Hrobarik *et al.* [185] used ZORA DFT calculations with a PBE0 hybrid functional with 40% exact exchange to predict hydride NMR shifts in U(VI) hydride complexes. A number of complexes with transition metals and other actinides were also included in the study in order to calibrate the computational model. The authors found extremely large hydride shifts induced by SO coupling at the U(VI) centres ranging from about +40 to +150 ppm. Such shifts are unprecedented for diamagnetic compounds and help to explain why the hydride shifts have not been assigned experimentally—presumably the experimentalists did not consider the appropriate frequency range. Large SO-induced carbon shifts were also reported in this study.

A strongly increased shielding of ^{29}Si along a Ni, Pd, Pt triad of novel hypervalent silicon complexes with a direct Si–metal bond has been observed experimentally. ZORA DFT calculations have demonstrated that the effect can be attributed to SO coupling [186]. SO coupling also substantially reduces the ^{29}Si shielding tensor span for the Pt compound, as was shown by a comparison of experimental and simulated magic angle spinning (MAS) spectra for low spinning rates. In a follow-up paper [187], an analysis method for SO effects on NMR shielding has been applied to the triad and a number of related complexes. A ‘partial ZORA’ operator allowed one to switch off SO coupling or all relativistic effects for selected atoms. (For a related idea applied to magnetic anisotropy, see [188].) It was demonstrated that a complex electronic interplay determines the experimentally observed NMR trends.

It has been pointed out in §5 that the ZORA quasi-relativistic Hamiltonian is quite accurate for valence shells in heavy atoms, but not for core shells. For instance, core orbital energies [189] and hyperfine integrals of core orbitals [190] afford substantial errors. Orbital energies can be improved substantially by applying the scaled ZORA corrections, as has been shown again recently by Saue [5]. In ZORA NMR calculations, the absolute shielding is adversely affected by errors from core orbitals, as has been demonstrated in a comparative benchmark study of ^{199}Hg shielding constants [105]. However, the errors cancel when the chemical shift is evaluated, as was shown by a comparison of Hg NMR shifts of HgX_{2} (X=Cl, Br, I) obtained with the ZORA implementation in the ADF code and 4c results calculated with the Dirac package, with Hg(CH_{3})_{2} as the reference. In fact, even with the Breit–Pauli Hamiltonian of order *c*^{−2}, agreement with 4c NMR data within a few per cent has been obtained for mercury chemical shifts [87]. Notable differences for absolute metal shielding constants in transition metal cyanides calculated with 4c methods versus ZORA were also noted in [191]. The implications regarding the chemical shift are the same as noted above. For *J*-coupling constants, the agreement between ZORA and 4c theory was very close (different codes with different basis set types were used), which reflects the accuracy of ZORA hyperfine integrals calculated for valence and outer core orbitals of heavy atoms [190].

A ZORA NMR shielding code developed by Wolff *et al.* [192] has recently been extended to accommodate calculations with hybrid functionals [193]. Chemical shifts were calculated for C, Os and Pt in 5d metal carbonyls, several square planar Pt(II) complexes, and Os phosphines. For ^{195}Pt shifts, some improvements towards experiment were obtained with the B3LYP hybrid, compared with the BP and BLYP non-hybrid functionals. For the other systems, the chemical shifts obtained with hybrid and non-hybrid functionals were similar. It is noted in passing that a corresponding ZORA NMR shielding module for periodic solids, using Slater-type AO basis functions, has been implemented at the SF level [194].

A new implementation for SF ZORA NMR calculations was used to investigate the functional dependence of heavy-atom shifts in DFT calculations [195]. The selection of functionals included variants of hybrid functionals with range-separated exchange. The test set included HX (X=F, Cl, Br, I, At), H_{2}X (X=O, S, Se, Te, Po), a number of tellurium compounds, and the series UF_{6−n}Cl_{n}. In particular, the experimental ^{19}F chemical shifts in the U(VII) halide family have in the past been notoriously difficult to reproduce by calculations [196,197]. Based on calculations with a variety of different classes of functionals, it was shown again in [195] that the fluorine shifts in uranium halides are sensitive to approximations in the functionals, indicating a delicate balance of DFT self-interaction versus correlation. The results with the range-separated functionals showed no clear improvement over previous calculations with global hybrids.

Relativistic and finite-nucleus [198] effects comparable to those for electron–nucleus hyperfine coupling discussed in the next section have been reported for NMR *J*-coupling involving heavy nuclei. For instance, the reduced nuclear spin–spin coupling *K*(Hg–C) in Hg(CN)_{2} was calculated non-relativistically as 2266×10^{19} T^{2} J^{−1}, and with scalar ZORA as 4408 (BP functional, point nuclei) [38]. Similarly, large relativistic effects were subsequently obtained by Filatov & Cremer with an IORA implementation [40]. The experimental coupling is 5778 in methanol. The difference between the relativistic calculations for an isolated complex and experimental data was later shown to be due to missing solvent effects [199].

The 4c restricted magnetic balance NMR code mentioned above has been extended for calculations of *J*-coupling [200]. Benchmark data were generated for the XH_{4} series of molecules, where X=C, Ge, Si, Sn, Pb. While scalar relativistic effects dominate the one-bond X–H couplings, SO coupling gives a large correction for X=Pb. The ^{1}*J*(Pb–H) value obtained for this system with the BP non-hybrid density functional was 2345 Hz with the four-component method, compared to an experimental estimate of around 2600–2800 Hz from extrapolations of measured couplings for methyl-substituted plumbanes. For comparison, a point-nucleus result from [201] obtained with ZORA and the PBE non-hybrid functional and converted to Hz (^{1}H, ^{207}Pb) would be 2836 Hz. In [201], it was shown that the use of a hybrid functional strongly increases the Pb–H *J*-coupling in plumbanes, while finite-nucleus corrections lead to a reduction again. For the plumbane series, reasonable agreement with experiment was obtained from hybrid DFT calculations with finite-nucleus corrections [201].

Typically, with *J*-coupling, the dominant effect is from the scalar relativistic increase of the electron density near the nuclei, leading to a strong response to the perturbation by the nuclear spin magnetic moment. An example for which SO effects are particularly pronounced is the *J*-coupling in the Tl–I diatomic. (The same molecule has been mentioned in §8*a* because of the large SO effects on the iodine EFG.) Both the isotropic *J*-coupling and the coupling tensor anisotropy are dominated by SO effects [202].

The functional dependence of *J*-coupling involving heavy metals and light ligands has been investigated in [203], using a ZORA hybrid DFT implementation [202]. The test set included 47 complexes with heavy metals (W, Pt, Hg, Tl, Pb). A total of 88 one-, two-, three- and four-bond coupling constants involving one or two heavy-metal atoms were evaluated at different levels of theory: non-hybrid versus hybrid DFT, scalar versus SO ZORA, gas phase versus a continuum solvent model. The relative median deviations between various computational models and experiment ranged between 13 and 21%. The highest-level computational model performed best (hybrid density functional computations including scalar plus SO relativistic effects, a continuum solvent model and a Gaussian finite-nucleus model). Further improvements would require explicit modelling of solvent effects [199,204]. The ZORA DFT code was also used to study *J*-coupling tensors in alkali halide diatomics, reproducing the experimental trends for the isotropic coupling and the tensor anisotropy with high correlation [205]. The use of a hybrid functional was found to be highly beneficial.

Finite-nucleus effects on *J*-coupling constants between Hg and a light atom were shown to be very significant [201], reducing the magnitude of the coupling by about 10% relative to point-nucleus calculations. The effects were found to be of comparable magnitude for Pb, but below 10% for Pt. For the Hg–Hg coupling in a crown ether complex of [Hg–Hg]^{2+}, the finite-nucleus effects reached 19%, and went up to 28% for the coupling in the free [Hg–Hg]^{2+} ion [201]. The finite-nucleus effect on ‘nuclear properties’ has two aspects. One is the relativistic increase in the electron density at a heavy nucleus, which is strongly altered when switching to an extended nucleus. The other is the modification of the vector potential describing the distribution of the magnetization density (equation (7.2b) is for a point-nuclear dipole). Maldonado *et al.* [206] considered finite-nucleus corrections to *J*-couplings and NMR shielding constants in 4c HF (RPA) calculations, by using finite-nucleus nuclear potentials and point-nucleus magnetic perturbation operators. The effects on *J*-coupling were in line with those previously observed in DFT calculations. The finite-nucleus effects on heavy-atom shielding constants were similar in magnitude to those determined in a benchmark study of ^{199}Hg NMR shielding [105]. It appears that the effect on the nuclear spin perturbation operator is secondary to the effect from the modification of the electron density at the heavy nucleus due to the finite-nucleus potential. It is, however, probably best to apply a consistent set of potentials and operators. For 2c methods, this also means that one needs to consider the finite-nucleus potential in the decoupling transformation.

### (c) EPR *g*-factors, hyperfine coupling, zero-field splitting

The electron *g*-factor and the electron–nucleus hyperfine coupling (hfc) are the basic parameters of the EPR spectrum of a radical. It is common practice to refer to the EPR Zeeman coupling matrix as the *g*-tensor [36], and to the hfc matrix as the hfc tensor. For systems with multiple unpaired electrons, there is also a zero-field splitting (ZFS) that must be included in the phenomenological spin Hamiltonian. The discussion in this section is limited to *g*-factors and hfc. The reader is referred to an interesting DFT approach for ZFS calculations recently put forward by van Wüllen and co-workers [188,207]. In this method, the ZFS tensor is extracted from the magnetic anisotropy (MA) of the ground-state DFT energy in generalized-collinear 2c calculations with different orientations of the spin quantization axis. See also reference [208].

Sandhoefer & Neese [148] have claimed the first complete derivation of the *g*-tensor at the DKH2 level. A point was made that the magnetic field has in previous work not been included in the free-particle FW transformation. The authors argued that this is necessary. The set-up for the *g*-tensor was made in an LR framework where SO coupling is included as a first-order perturbation. A similar LR approach within the ZORA relativistic framework was put forward previously in [195,209]. The *g*-tensor calculations in [148] were carried out with the BP non-hybrid functional on a test set of transition metal complexes and the lanthanide series LnH_{2}. Comparisons with ‘NR’ calculations (i.e. an NR starting point for the LR calculations) demonstrated consistently better agreement of the relativistic results with experimental data.

Scalar relativistic NESC calculations have been used to evaluate ^{199}Hg hfc constants of a number of radicals containing mercury [210]. NR results were also provided, showing relativistic increases by factors of about 2–3 in HF calculations, consistent with the effects reported in the past for ^{199}Hg NMR *J*-couplings [38,40,211]. The reason is that hfc and *J*-coupling involve, in part, the same electron–nucleus contact operators, which are very sensitive to the relativistic increase of the electron density and spin density at a heavy nucleus. Results from correlated wave function methods and HF theory were reported; a comparison of CCSD with DFT data from other works is provided in table 3 for HgF and HgAg. Correlation effects were modest for the two molecules listed, but large for alkyl mercury radicals.

A number of recent articles [145,212–214] have highlighted the importance of effects from finite nuclear volumes in calculations of hfc. These effects are extremely small in NR calculations but can get very strongly amplified by relativistic effects. Therefore, finite-nucleus effects are sometimes considered as being of relativistic nature. Consideration of finite-nuclear current-density distributions affecting the hyperfine structure of heavy elements dates back to the early 1950s publications by Bohr & Weisskopf [216,217] and further back to the early days of modern quantum theory [218].

Given the strong inverse dependence of the hyperfine operators on the electron–nucleus distance, one may expect very large relativistic effects on properties that depend on these operators, along with sizable finite-nucleus effects. This is indeed the case. Consider the comparison of point-nucleus versus finite-nucleus ^{199}Hg hfc constants in the HgF and HgAg radicals listed in table 3. The data were selected from larger collections reported in [145,212–214]. There are sizable finite-nucleus corrections of the order of 10–15%. Finite-nucleus effects for the Cs and Fr atoms obtained from 4c calculations were reported in [214] to reach −2 and −12%, respectively. For various Hg radicals, the mercury HFCC finite-nucleus corrections were found to be lower than previously obtained from DKH2, and thus more in line with the ZORA data listed in table 3. Higher-order SO effects (beyond a first-order perturbation treatment) for ^{199}Hg isotropic hfc were found to be small [213] (which is reflected in the small differences between the ZORA data in table 3 for calculations where SO coupling is treated as a perturbation versus variationally).

Deviations of atomic and molecular electronic *g*-factors from the free-electron value *g*_{e}≈2.0023 may be considered a relativistic effect for orbitally non-degenerate states in the sense that such deviations vanish for vanishing SO coupling (which is a relativistic effect). The *g*-shifts (Δ*g*=*g*−*g*_{e}) are often small and reported in parts per thousand (ppt). A 4c DFT implementation using non-collinear spin densities and restricted magnetic balance has been reported recently [219]; for previous related work by the same group, see [220]. In a series of test calculations on d^{1} metal complexes, it was found that higher-order SO effects in *g*-shift calculations are larger than previously thought [221]. Therefore, theoretical approaches that treat SO coupling as a linear perturbation (LR) are not capable of producing the full magnitude of the *g*-shift for the right reasons.

We recently implemented ZORA *g*-shift calculations with GIAO basis sets and a variety of functionals, including hybrids with range-separated exchange, within an LR framework [195,209], and in the form of two approaches that include SO coupling variationally in the ground state [222]. One of them has been a new implementation of a method developed originally by van Lenthe *et al.* [223] (LWA), which does not allow for spin polarization in the calculations. The second variant affords spin polarization. It has been inspired by the aforementioned DFT-based approach for ZFS calculations proposed by van Wüllen and co-workers and is therefore referred to as MA (for magnetic anisotropy). The LWA and MA calculations reaffirmed the importance of higher-order SO effects previously highlighted in [221]. For related hfc methods, see [213].

Interesting cases with large *g*-shifts are NpF_{6} and related actinide hexahalides such as UF^{−}_{6} and UCl^{−}_{6}, which have formal 5f^{1} configurations. Representative data are collected in table 4. For instance, the accepted experimental *g*-factor of NpF_{6} is −0.6, implying a *g*-shift of the order of −2600 ppt. The signs and overall magnitudes of the experimentally derived *g*-factors have been reproduced by SO CASPT2 calculations [224,225] as well as by ZORA DFT [222], although not without significant variability. In LWA and MA calculations, a complex interplay of approximations in different classes of functionals versus the presence of lack of spin polarization emerged. Overall, when considering also results from a test set of radicals with transition metals, the MA approach for *g*-shift calculations in conjunction with hybrid functionals performed best (often in conjunction with a long-range corrected (LC) variant of PBE0).

## 9. Concluding remarks

There is currently a tremendous research activity focused on implementing molecular derivative properties such as EFGs, NMR and EPR parameters within relativistic electronic structure methods. DFT remains the ‘workhorse’ for applications. Four-component methods have become computationally very efficient, and two-component methods have evolved to a fully relativistic level, at least for the one-electron part of the operators. This is good news because the level at which relativistic effects are dealt with in one-electron operators and electric- and magnetic-field perturbations thereof is becoming irrelevant. There remains work to be done at the many-electron level, both in wave function theory and in DFT. Furthermore, the code base for relativistic molecular properties needs to be expanded much beyond the currently available implementations. Existing quasi-relativistic methods such as DKH and ZORA are likely to remain in use for the foreseeable future. In many practically relevant scenarios, these quasi-relativistic operators are accurate enough in terms of treating relativity, considering that the result of a computation is affected by many other factors such as basis set limitations, approximations in the electronic structure model, and other approximations in computational models (neglecting solvent, vibrations, temperature and so on).

## Funding

The author is grateful for financial support by the US Department of Energy (Basic Energy Sciences, Heavy Element Chemistry Program, grant no. DE-FG02-09ER16066) for method developments related to relativistic calculations of NMR and EPR parameters, for the implementation of exact two-component Hamiltonians in the NWChem quantum chemistry package, and for applications to NMR and EPR parameter calculations in heavy-element systems.

## Acknowledgements

The author acknowledges stimulating discussions with Dr Daoling Peng on the topic of two-component Hamiltonians.

## Footnotes

One contribution of 8 to a Theme Issue ‘Density functional theory across chemistry, physics and biology’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.