Local relativistic exact decoupling

The consideration of relativistic effects is essential to the proper understanding of the chemistry of heavy-element-containing molecules.^1,2 The Dirac equation provides the relativistic quantum mechanical description of a single electron in the presence of external electromagnetic potentials. As the Dirac operator comprises four-dimensional operators, its eigenfunctions possess four components. The Coulomb interaction turned out to be sufficiently accurate for the description of chemical phenomena and the corresponding many-electron Hamiltonian is called the Dirac–Coulomb Hamiltonian.

A relativistic method based on the Dirac–Coulomb Hamiltonian is called a four-component method. It yields negative-energy solutions, which are pathologic and of no use in chemistry. Also, a large number of basis functions are required to properly describe the negative-energy states. These are two drawbacks of four-component methods that motivated the development of two-component methods, which remove the negative-energy solutions and can, in principle, exactly reproduce the results of four-component calculations.² 

Several relativistic two-component methods were developed in the past decades. One of the widely used approaches is the second-order Douglas–Kroll–Hess method (DKH2).^3,4 It employs the free-particle Foldy–Wouthuysen (fpFW) (Ref. 5) transformation as well as sequential Douglas–Kroll⁶ transformations to decouple the four-component operator. Higher order^7–10 and even arbitrary-order^11–16 DKH methods have been developed. The zeroth-order regular approximation (ZORA) (Refs. 17–19) is another highly successful relativistic two-component method. Within the ZORA framework, it is particularly easy to implement the calculation of molecular properties; see Refs. 20 and 21 for very recent examples and Ref. 22 for a review.

In recent years it was shown that, if the goal is to exactly decouple the four-component Hamiltonian matrix instead of the Hamiltonian operator, the formulas of exact decoupling can be directly obtained. The Barysz–Sadlej–Snijders (BSS) (Refs. 23–26) method aims at exact decoupling of the free-particle FW-transformed four-component Hamiltonian matrix by a matrix operator of the form derived in Ref. 27. The key matrix operator, which is used to construct the decoupling transformation matrix, is obtained by solving an iterative equation. However, invoking the free-particle FW transformation turns out to be not necessary for the construction of the exact-decoupling transformation in a one-step protocol. The pioneering work of this one-step transformation was provided by Dyall^28–32 in the form of the so-called normalized elimination of the small component (NESC) approach. Later it was generalized to the so-called exact two-component (X2C) method by several groups.^33–44 Since the DKH method also employs matrix techniques to evaluate the Douglas–Kroll transformation, it is also able to exactly decouple the four-component Hamiltonian matrix. This has been shown within the arbitrary-order approach.^{11–14,16,45} For reviews of these developments see Refs. 46–52.

Even within a scalar-relativistic approximation, where the spin-dependent operators are neglected, the construction of the relativistic Hamiltonian matrix dominates the step of the one-electron integral evaluation. This is the case because the relativistic transformation is done within an uncontracted basis set and involves many matrix manipulations such as multiplication and diagonalization. For a large molecule including only one or a few heavy atoms, this procedure would waste computational resources since the relativistic effects are highly local and only significant for heavy atoms. It is therefore clear that a local relativistic method is desirable. Even for a molecule containing only heavy atoms, such a local relativistic approach would be preferable as we shall see.

The relativistic Hamiltonian operator is universal in coordinate space and a local approach is unfeasible at the operator level. However, since most quantum chemical calculations apply a linear combination of atom-centered basis functions, we may employ these basis functions to construct atomic projectors as it is done in charge and spin population analysis.⁵³ For those Hamiltonian matrix blocks for which relativistic corrections are important, we shall derive a local relativistic Hamiltonian matrix to evaluate it. The locality is then exploited in the basis-function space instead of coordinate space. A crucial aspect to determine will be which block of the Hamiltonian matrix requires a relativistic description and how to evaluate it. Undoubtedly, the heavy-atom diagonal blocks will require a relativistic description, while the relativistic description of heavy-atom off-diagonal blocks depends on their contribution to physical observables and on the feasibility of applying the relativistic description. For relativistic transformation techniques, such as DKH and X2C, it is impossible to apply the relativistic transformation to an off-diagonal block alone without information about other blocks. It requires the information of a full square matrix to construct the relativistic transformation. Such a difficulty disappears if one employs an operator-based relativistic method such as ZORA. The off-diagonal block of the Hamiltonian matrix can be directly computed from relativistic integrals. However, the ZORA method provides a poor description of atomic inner shell orbitals, which carry the largest part of the relativistic effect. In this sense, the exact-decoupling methods are the best candidate for applying a local relativistic scheme since they do not neglect important relativistic effects.

All current exact-decoupling methods employ the transformation technique. It is therefore not trivial to apply the relativistic description to heavy-atom off-diagonal Hamiltonian matrix elements. However, in this work we propose that if the local approximation is applied to the decoupling transformation instead of the Hamiltonian matrix, one could obtain results which include the relativistic description to heavy-atom off-diagonal blocks. We also note that the main effort for the construction of the relativistic Hamiltonian matrix with exact-decoupling methods is due to the evaluation of the decoupling transformation matrices. Therefore, a local approximation to the decoupling transformation would also lead to a significant reduction of the computational cost.

The following notation will be used throughout this paper. Upper case labels A, B, C denote heavy atoms which require a relativistic description, while lower case labels a, b, c denote light atoms for which a (standard) nonrelativistic (NR) description can be assumed to be accurate enough based on numerical evidence compiled within computational chemistry in the past decades. To be more precise, the former labels refer to nuclei with high nuclear charge numbers, while the latter refer to those with low nuclear charge numbers (typically with nuclear charges of less than two dozen protons). Of course, the distinction whether a nucleus is considered heavy or light bears some sort of arbitrariness and it will be made on the basis of the accuracy required in a quantum chemical calculation. In general, there are then five different types of Hamiltonian matrix blocks H_AA, H_AB, H_Aa, H_aa, and H_ab that need to be considered.

The evaluation of the relativistic Hamiltonian matrix with exact-decoupling approaches requires a set of square matrices. They are the nonrelativistic kinetic energy matrix T,

the external potential matrix V,

and the overlap matrix S,

as well as an additional relativistic matrix W,

(note that all expressions are given in Hartree atomic units, in which the rest mass of the electron takes a numerical value of one). In the above equations, λ_i denotes the ith 2-spinor atom-centered basis function,

$V$ the external potential, c the speed of light, and σ the vector of Pauli spin matrices. The basis functions λ_i may be grouped according to the atomic nucleus to which they are assigned, which is key to the local approaches discussed in the following. The matrix representation of the relativistic Hamiltonian is then evaluated as a function of the above-mentioned matrices

A straightforward idea for a trivial local approximation is to apply the evaluation of the relativistic Hamiltonian matrix only to atomic (i.e., diagonal) blocks

(instead of to the full matrix, but where A may also represent a group of atoms) and to ignore the relativistic correction to all off-diagonal blocks

Adding more atoms to a group A will improve the accuracy, but the computational advances will be lost when a group A becomes very large. This approach was employed in Refs. 54 and 55 and explored in detail in Ref. 56 for the low-order DKH method. We refer to it as the diagonal local approximation to the Hamiltonian (DLH) matrix. It is clear that the DLH approximation can also be applied to relativistic exact-decoupling approaches. The DLH approximation works well at large interatomic distances, but the difference to a reference energy increases with shorter distance.

However, relativistic corrections to off-diagonal blocks are also important especially when interatomic distances are short. The necessity of applying the relativistic description to off-diagonal blocks has already been discussed in Ref. 57 with a so-called two-center approximation. In Ref. 57, the DKH transformation was applied to pairs of atoms (A ⊕ B),

with

to obtain a relativistically corrected off-diagonal Hamiltonian matrix H_AB. Unfortunately, such a two-center approximation introduces an inconsistent treatment to diagonal blocks. For instance, H_AA can be obtained from either the pair (A ⊕ B) or from another pair (A ⊕ C).

So far, we have only discussed the AA and AB blocks which, by definition, require a relativistic description for both atoms. Since relativistic corrections for light atoms are very small and may be neglected, H_aa and H_ab can be approximated by the nonrelativistic Hamiltonian matrices

A relativistic scheme for the heavy–light hybrid off-diagonal blocks H_Aa should, however, be considered explicitly.

The general expression of exact-decoupling transformations can be written as

Z⁺ and Z⁻ are two-component unitary operators, and only Z⁺ is required for the evaluation of the electrons-only Hamiltonian. The X-operator generates the electronic small-component functions

$φS+$ from the large-component functions

$\varphi ^{+}_{L}$

$φL+$⁠,

Two-dimensional Hamiltonians are then obtained by applying the unitary transformation of Eq. (12) to blockdiagonalize (“bd”) the four-dimensional Dirac-based Hamiltonian D,

It yields both the electrons-only Hamiltonian H (as the upper left block in H^bd) and the one for negative-energy solutions, although the latter will be discarded.

As discussed by Dyall in Ref. 29 for the NESC method, the Breit–Pauli approximation to the X matrix, which is used to construct the relativistic transformation matrix, is

where I denotes the identity matrix. Such an approximation is not variationally stable and thus cannot be used for variational calculations. However, one could employ this approximation for light atoms only,

since the corrections to light atoms are small and may not affect variational calculations. This idea was suggested in Ref. 42 for the X2C method but no results were presented.

Although the lowest level approximation to X, Eq. (15), did not provide useful results, the atomic approximation to the X matrix gave results with very small errors. As discussed by Dyall in Refs. 30 and 31, the X matrix can be approximated as the direct sum of atomic blocks

This approximation resembles the linear combination of atomic four-spinors ansatz in four-component calculations. The same approach was also employed in Ref. 58 in the sense of an atomic approximation to the projection on electronic states. If the electronic states are further transformed to a two-component picture by a renormalization matrix, this gives rise to a local X2C method. As discussed in Refs. 38 and 42, the local approximation to the X matrix Eq. (17) works well for spectroscopic constants of diatomic molecules.

However, there are some drawbacks of the local X-matrix approximation to the X2C method. First, it only reduces the cost for the evaluation of X matrix, while the exact-decoupling transformation matrices are still of molecular dimension as the renormalization matrix is not atomic block diagonal even if the X matrix is. The computational demands for the evaluation of the relativistic approximation are thus still tremendous for large molecules. Second, the X = I approximation did not provide a satisfactory treatment of the Aa blocks, since it may suffer from variational collapse and Eqs. (10) and (11) cannot be fulfilled within this approximation.

In this paper, we suggest a local approximation to the exact-decoupling transformation. We may approximate the unitary transformation by taking only the “atomic” diagonal blocks (all off-diagonal blocks are then set to zero)

where the atomic unitary transformations U_AA are obtained from the diagonalization of the corresponding D_AA blocks of matrix operator D. We denote this approximation to U as the diagonal local approximation to the unitary decoupling transformation (DLU).

If we now substitute the approximate unitary transformation of Eq. (18) into Eq. (14), the diagonal blocks,

turn out to be the same as in the DLH approach, while the off-diagonal blocks then read

which is to be compared with the expression

where U has not been approximated (here, I and J run over all atomic blocks). Hence, the DLU approach also introduces a relativistic description to off-diagonal “interatomic” blocks when compared with the DLH approximation.

The cost for the assembly of the unitary transformation U within the DLU approach is then of order N, where N measures the system size, namely, the number of atoms. It is linear scaling since the atomic approximation is directly applied to the unitary transformation. However, the next step of applying the unitary decoupling transformation to obtain the relativistic Hamiltonian matrices is no longer linear scaling. If no local approximation was applied for the unitary decoupling transformation, according to Eq. (21), where the summation indices I, J run over all atomic blocks, the calculation of H_AB matrix will require 2N² matrix multiplications. Since the number of matrices to be calculated is of order N², the total cost of the relativistic transformation without local approximation is therefore of order N⁴. If the DLU approximation is applied, no summation will be needed and only two matrix multiplications will be required for each heavy–heavy block H_AB. The total cost is then of order N².

If the distance of two atoms A and B is sufficiently large, the relativistic description of H_AB can be neglected. Thus, we may define neighboring atomic pairs according to their distances. Then, only the Hamiltonian matrix

$HABbd$ of neighboring pairs requires the transformation

whereas all other pairs are simply taken in their nonrelativistic form. Since the number of neighboring pairs is usually a linear function of system size, the total cost is then reduced to order N.

To investigate the relativistic correction to hybrid H_Aa blocks of the electrons-only Hamiltonian H, we need to consider the explicit structure of the unitary decoupling transformation,

where only the electronic part of the exact-decoupling transformation represented by U^LL and U^LS [which are the upper part of Eq. (12)] is employed. From Eq. (23), it is easy to see that the electrons-only Hamiltonian reads

In the NR limit, where the speed of light approaches infinity, we have

If we insert them into the relativistic electrons-only Hamiltonian (24), we arrive at the nonrelativistic Hamiltonian matrix

This suggests a solution for the relativistic treatment of the hybrid heavy–light blocks H_Aa because we may set all light-atom diagonal blocks in the unitary transformation to identity matrices

To reproduce the nonrelativistic Hamiltonian matrices of light-atom-only blocks, we must set the W matrix of light atoms to zero

The W_Aa matrices are also set to zero

We denote this approximation for the treatment of hybrid blocks as the DLU(NR) approach. If we do not introduce any approximations to the W matrices, the transformation will yield a Hamiltonian that differs from the nonrelativistic Hamiltonian matrix on the order of 1/c², which corresponds to the order achieved in the Breit–Pauli approximation. We therefore denote it as the DLU(BP) approach.

The final local approximation to the unitary matrices is then

where the atomic unitary transformation matrices of heavy atoms are evaluated by the relativistic exact-decoupling approaches. The explicit expressions of Hamiltonian matrix blocks H_AA and H_AB are

If the DLU(NR) approach is employed, the unitary transformation matrices of light atoms are simply set to identity matrices and the explicit form of H_Aa, H_aa, and H_ab is

In the DLU(BP) approach, however, they read

The relativistic transformation of a property operator P^{13,14,16,45,59,60} requires an additional matrix Q and its matrix elements are given by

The relativistic corrected property matrix then reads

The X2C method employs

as the exact-decoupling unitary transformation. The unitary operators for the electronic states are then

For an actual implementation, we need their matrix representations in a finite non-orthogonal space of the atom-centered basis functions. The X matrix is evaluated by diagonalizing the following generalized eigenvalue equation:

where

$CL+$ and

$C_{S}^{+}$

$CS+$ denote the large- and small-component molecular-orbital coefficients, respectively. We obtain the X matrix by the following equation:

The matrix representation of U^LL and U^LS is then given by

with

$S̃$ defined as

In our atomic approximation to the unitary transformation U, we need to evaluate the above equation for each “atomic” block of the original Hamiltonian in order to obtain the diagonal blocks

$UAALL$ and

$U^{LS}_{AA}$

$UAALS$⁠. The decoupling transformations are then assembled according to Eqs. (30) and (31) and the relativistic Hamiltonian matrix follows from Eq. (23).

In the BSS approach, all matrices are transformed to an orthonormal basis-function space first. It is convenient to calculate the transformation matrix K by diagonalizing the nonrelativistic kinetic energy matrix

since the eigenvalues t can be used for later evaluation of the free-particle FW transformation. The eigenvector matrix K has the following properties:

The fpFW transformation features four diagonal block matrices

with

$E0=2tc2+c4$ and

$f=\sqrt{2c^2/t}$

$f=2c2/t$⁠. It is then applied to yield a transformed four-dimensional Hamiltonian matrix D₀,

Next, the exact-decoupling BSS transformation,

which has the same structure as the exact-decoupling transformation U_X2C, is applied. The R matrix is obtained by diagonalizing the free-particle FW-transformed four-component Hamiltonian matrix D₀ and employing a similar relation to the one in Eq. (46).

After the exact-decoupling BSS transformation has been carried out, the Hamiltonian matrix is back-transformed to the original (non-orthogonal) basis

Consequently, we obtain the matrix form of

$U BSS LL$ and

$U_{\rm BSS}^{LS}$

$U BSS LS$ for the BSS approach as

The expressions for the relativistic Hamiltonian matrix (and property matrices) are then as given above.

The DKH approach features the same initial transformation as the BSS approach to obtain the transformed Hamiltonian matrix D₀. Subsequent decoupling transformations are expressed as

with the generalized parametrization of the U_k,⁹ 

where W_k are anti-Hermitian matrix operators and m is related to the order of the DKH expansion. Since the expression of W_k is too lengthy to be presented here, we refer the reader to Ref. 16 for details. If m is large (strictly, if it approaches infinity), exact decoupling is achieved. Usually, a low value for m is often sufficient for calculations of relative energies and valence-shell properties (i.e., m = 1 for the original DKH2 approach).

With the complete decoupling DKH transformation written as

the matrix forms of

$U DKH LL$ and

$U_{\rm DKH}^{LS}$

$U DKH LS$ are

We should note that the finite-order DKH Hamiltonian matrices are not the result of directly applying the decoupling transformation. To give consistent results, the transformation in Eq. (23) is only used to obtain the off-diagonal blocks while the diagonal blocks are replaced by the traditional finite-order DKH Hamiltonian. However, for high-order calculations, the differences of with/without replacing the diagonal blocks are very small.

As we have shown, there exists a systematic hierarchy of several levels for the local relativistic approximation to Hamiltonian matrix elements. In this section, we study the accuracy of the different levels on a test molecule set and introduce the following notation to distinguish the results:

1. FULL: Full molecular relativistic transformation, no approximation.

2. DLU: Diagonal local approximation to decoupling transformations, U^LL and U^LS are evaluated for each atomic (diagonal) block. However, three variants are possible:

   • DLU(ALL): Calculate relativistic transformation for all atoms (all atoms are treated as “heavy” elements). That is, apply the “atomic” unitary transformation to the diagonal and to all off-diagonal blocks of the Hamiltonian.

   • DLU(NB): Calculate relativistic transformation for all atoms (all atoms are treated as “heavy” elements), but apply the relativistic transformation only to blocks of the Hamiltonian assigned to neighboring atoms (and to the diagonal blocks, of course) to achieve linear scaling.

   • DLU(ABC,BP): Distinguish heavy and light atoms, where ABC denote the heavy atoms and apply the BP approximation to the light atoms.

   • DLU(ABC,NR): Distinguish heavy and light atoms, where ABC denote the heavy atoms and apply the NR approximation to the light atoms.

3. DLH: Diagonal local approximation to Hamiltonian matrix blocks. The relativistic correction for off-diagonal “interatomic” blocks are neglected.

   • DLH(ALL): Calculate a relativistic description for all diagonal “atomic” blocks.

   • DLH(ABC,NR/BP): Only heavy-atom blocks are considered for the relativistic description.

Note that this classification system does not account for all possible local relativistic approximations: (i) A could also represent a group of atoms and (ii) additional off-diagonal approximations may also be included.

We have implemented the local relativistic schemes discussed so far for the X2C, BSS, and DKH exact-decoupling approaches into the MOLCAS programme package,⁶¹ into which we have recently implemented X2C, BSS, and the polynomial-cost DKH schemes.⁵² In our scalar-relativistic calculations, the definition of the W matrix was

and λ_i now refers to scalar basis functions instead of 2-spinor basis functions. All other expressions derived in this work are the same for the scalar-relativistic X2C, BSS, and DKH variants. We should mention that we have not considered the transformation of the two-electron integrals as it is common in most exact-decoupling calculations. However, the local methods discussed here could be of value also for one of the approximate methods available to transform the two-electron integrals (note that a full-fledged transformation of the two-electron integrals would render the application of exact-decoupling methods inefficient and a four-component approach should be employed instead).

For our test molecule set, we should rely on closed-shell molecules with heavy elements in the vicinity of other heavy elements and of light elements. Four molecules were selected to test the validity of the different local schemes:

$I5+$⁠, WH₆(PH₃)₃, W(CH₃)₆, and

${\rm Pb}_9^{2+}$

$Pb 92+$⁠. Moreover, two reactions

were chosen to study reaction energies. Their structures (see Fig. 1) have been optimized with the TURBOMOLE program package⁶² employing the BP86 density functional^63,64 with a valence triple-zeta basis set with polarization functions on all atoms⁶⁵ in combination with Stuttgart effective core potentials for W, Pb, and I as implemented in TURBOMOLE.

All molecules are closed shell and for our analysis all-electron Hartree–Fock calculations were then performed with the MOLCAS programme package. Since our aim is to test the validity of local relativistic approximations instead of obtaining accurate results which can be compared to experiments, we do not need to consider electron-correlation effects in our calculations. The all-electron atomic natural orbital basis sets^66–68 were employed at a double-zeta level for all atoms. The exact-decoupling DKH calculations were performed at 35th order, since higher order results did not improve the accuracy for our purpose so that this high DKH order can be considered exact. Results of finite-order (standard) DKH2 calculations are also given for the sake of comparison.

The computation times on the Intel i7-870 central processing unit (CPU) for the different local schemes are compared at the example of

$Pb 92+$⁠. Table I presents the timings for the evaluation of the relativistic one-electron Hamiltonian matrix operator. (The picture change transformation of property integrals is not considered in the measurement of the CPU time.) For comparison, the computation times for the one-electron and for the electron repulsion integrals are 8 s and 29744 s, respectively. They are, of course, the same for the different relativistic approaches. Since the relativistic transformations are performed on each atomic block in the DLH scheme, its computational costs are much smaller than the FULL scheme. The ratio of DLH to FULL is almost a constant (around 1% in this case) for all relativistic approaches. The DLU scheme has an extra transformation step for off-diagonal blocks in addition to the DLH scheme, it therefore requires a constant time (about 10 s in this case) more than the DLH scheme, which is still negligible when compared to the FULL scheme. The computation times for the evaluation of both one-electron integrals and relativistic transformations (except for DKH35) are smaller than that for two-electron integrals. The relativistic transformation is usually not the bottleneck step. However, if density functional theory calculations—with a pure generalized-gradient-approximation functional in combination with some acceleration scheme like density fitting or fast multipole approximations—instead of Hartree–Fock calculations were performed, the computation time of the relativistic transformation would be significant and the local approximation would then be required. Of course, this also holds especially for much larger molecules with light atoms (see also discussion below). For the high-order DKH method, such as DKH35 employed here, the computational cost of the FULL scheme even exceed that of the two-electron integrals and therefore a local approximation is mandatory.

We performed DLU(ALL) and DLH(ALL) calculations for all four molecules in order to analyze the general error that is introduced by these local approximations. The differences of total electronic energies with respect to the FULL reference calculation without any local approximation are presented in Table II. Since the errors of total energies are rather small, all values are given in milli-Hartree (mH) atomic units. The errors in total energies can be understood as errors on atomization energies that one would observe when atomization energies are to be calculated with one of the local approaches compared to the FULL one.

For the

$I5+$ molecule, which possesses a pseudo-one-dimensional structure with iodine atoms connected to at most two other iodine atoms, the DLU(ALL) scheme produces an error of about 0.012 mH for both DKH and BSS, while X2C has a slightly larger error of 0.013 mH. Interestingly, it is smallest for DKH2, i.e., only 0.008 mH. The DLH(ALL) local scheme produces relatively large errors on the total electronic energies, they are 0.24 mH for DKH and BSS, 0.26 mH for X2C, and 0.29 mH for DKH2. Clearly, DLU(ALL) provides more accurate total energies than DLH(ALL), as it includes the relativistic form of the atom–other-atom off-diagonal blocks, while DLH(ALL) only takes into account the atom–same-atom diagonal blocks. The errors of DLH(ALL) are about 20 times larger than those of DLU(ALL) for the

${\rm I}_{5}^{+}$

$I5+$ molecule. This indicates that the relativistic corrections to off-diagonal blocks are quite important. Apparently, DLU(ALL) provides an accurate scheme for taking the off-diagonal relativistic description into account.

The

$Pb 92+$ molecule features a more compact structure where more than two heavy atoms bind to a central lead atom. Moreover, lead atoms are heavier than iodine atoms and relativistic effects are more pronounced. As we can see from Table II, errors of

${\rm Pb}_9^{2+}$

$Pb 92+$ are all larger than those observed for the

${\rm I}_{5}^{+}$

$I5+$ molecule. The deviations of DLU(ALL) electronic energies from the FULL reference energy are 0.038 mH for BSS and DKH, while X2C is still slightly larger with a value of 0.052 mH (again DKH2 features the smallest deviation of only 0.012 mH). However, deviations of DLH(ALL) results from the reference have become significantly larger, they are 2.27, 2.25, and 2.84 mH for DKH, BSS, and X2C, respectively. Hence, the errors of DLU(ALL) are roughly 60 times smaller than those observed for DLH(ALL).

Considering now molecules that also contain light atoms, the relative accuracy of DLU(ALL) vs. DLH(ALL) changes dramatically. First of all, we note that our DLU(ALL) scheme preserves its accuracy and the deviations from the FULL reference electronic energy are in between those obtained for

$I5+$ and

${\rm Pb}_9^{2+}$

$Pb 92+$ with the same trends as already discussed. The errors of DLH(ALL) energies obtained for the WH₆(PH₃)₃ molecule are, however, larger than 30 mH. That is, they are significantly larger than those observed for

${\rm Pb}_9^{2+}$

$Pb 92+$⁠. For W(CH₃)₆, the deviations from the reference are around 6 mH, which are also still larger than for

${\rm Pb}_9^{2+}$

$Pb 92+$ and certainly non-negligible. Note that the DLH(ALL) errors obtained for WH₆(PH₃)₃ and W(CH₃)₆ are not all below those of the

${\rm Pb}_9^{2+}$

$Pb 92+$ molecule although the nuclear charge number of tungsten is smaller than that of lead. This is in contrast to the results obtained for our DLU(ALL) scheme. Since DLH neglects any relativistic correction to off-diagonal “interatomic” blocks, the large errors must be due to the strong coupling between tungsten and its surrounding lighter atoms (phosphorus and carbon). Apparently, if a heavy-element-containing molecule has several strong bonds to this heavy element, the off-diagonal blocks of neighboring atom pairs can have significant contributions to the total electronic energy. Therefore, the DLH scheme may then introduce non-negligible errors.

Since DLU(ALL) gives a balanced description for all blocks, it will still work in those cases where DLH(ALL) fails. DLU can therefore even be recommended for calculations of molecular systems in different geometries as, for example, required for the calculation of potential energy surfaces. DLH would be unfit for this case as it will introduce large errors when the interatomic distances become short.

DLU(ALL) also works well for the finite-order (standard) DKH2 method. The errors are even smaller compared to exact-decoupling approaches. However, note that the reference energies of the different relativistic approaches are not the same. Since the DLU(ALL) scheme only discards the off-diagonal blocks of the decoupling transformation compared to the FULL reference scheme, the small errors in the total electronic energies indicate that neglecting off-diagonal terms in the decoupling transformation produces negligible errors for the total energies. This conclusion holds for both exact-decoupling and finite-order DKH approaches. With the DLH(ALL) scheme, however, the errors in the DKH2 calculation that are due to the DLH approximation are of the same order as those introduced by the finite, i.e., second order when compared to the exact-decoupling approaches. This observation leads to the very important conclusion that none of the exact-decoupling approaches is superior to a finite-order DKH approach if one employs a DLH scheme as this would render the accuracy gained by exact decoupling meaningless. In fact, if DLH is for some reason the local approximation of choice, a standard DKH2 approach would yield a more efficient relativistic scheme and nothing would be gained by achieving exact decoupling.

Besides total energies, we also study the effects of local approximations on molecular properties. We calculated the picture change corrected quadrupole moments. The property matrices are transformed according to Eq. (41). The values for the components of the quadrupole moment depend on the orientation of the molecule. Therefore, we present the spherically averaged quadrupole moments defined as

$(Q XX 2+Q YY 2+Q ZZ 2)/3$⁠, in Table III, where Q_XX denotes the XX (diagonal) component of the quadrupole moment. The spherically averaged quadrupole moments are isotropic and thus independent of molecular orientation.

The relative errors of DLU(ALL) are all below 0.001%, this upper limit can be further reduced to 0.0001% if we do not include

$Pb 92+$⁠. Hence, DLU(ALL) is not only a good approximation for total electronic energies, but also for molecular properties such as quadrupole moments. The errors introduced by DLU(ALL) are all negligible, whereas DLH(ALL) shows large errors for quadrupole moments. While for total electronic energies, the largest error of the DLH(ALL) scheme shows up in WH₆(PH₃)₃ and the smallest in

${\rm I}_{5}^{+}$

$I5+$⁠, it is, for the quadrupole moments, largest for

${\rm Pb}_9^{2+}$

$Pb 92+$ and smallest for W(CH₃)₆. This different behavior of DLH(ALL) indicates again the important contribution of off-diagonal “interatomic” blocks, whose importance must be different for electronic energies and for properties. In our DLU(ALL) scheme,

${\rm I}_{5}^{+}$

$I5+$ has the smallest error and

${\rm Pb}_9^{2+}$

$Pb 92+$ the largest one for both, energies and properties. Therefore, DLU(ALL) turns out to give a consistent description of relativistic effects, while DLH(ALL) yields an unbalanced relativistic treatment. We should stress that we have also investigated other properties (radial moments, electric field gradients, contact densities) drawing basically the same conclusions (which is the reason why we do not report those data here).

For both, energies and quadrupole moments, in the DLU(ALL) scheme, DKH and BSS yield almost the same deviations from the reference and X2C always has larger ones than DKH and BSS. This similarity of DKH and BSS can be well understood since both of them employ the free-particle FW transformation as well as the subsequent exact-decoupling transformation. The close results of DKH and BSS were also observed in Ref. 52. The slightly larger errors of X2C indicate that the off-diagonal blocks of the X2C exact-decoupling transformation have slightly larger contributions compared to DKH and BSS. However, such differences of errors are too small to be considered a serious drawback for the X2C approach; they are smaller than the differences of the total values. We may use WH₆(PH₃)₃ as an example. The difference of total isotropic quadrupole moments (57.61285−57.61240) is 0.45 milli atomic units, while the difference of errors by DLU(ALL) (0.00029−0.00018) is 0.11 milli atomic units. The errors introduced by the DLU(ALL) approximation are smaller than the differences between exact-decoupling approaches. This is also the case for the differences between DKH2 and exact-decoupling approaches. As we can see, the difference between DKH2 and BSS of total isotropic quadrupole moments (57.61393–57.61285) is 1.08 milli atomic units, which are much larger than the DLU(ALL) errors.

For chemical purposes, comparisons of total electronic energies are meaningless because chemical reaction kinetics and thermodynamics are governed by energy differences. We therefore calculated the electronic reaction energies of the two reactions mentioned above. The results are given in Tables IV and V. The errors of relative energies are less than those of absolute total energies. For example, with the DKH approach, the error of

$I5+$ is reduced from 0.012 mH to 0.002 mH for DLU(ALL), 0.24 mH to 0.12 mH for DLH(ALL). For the WH₆(PH₃)₃ molecule, the error is decreased from 0.0126 mH to 0.0006 mH for DLU(ALL), 30.9 mH to 19.4 mH for DLH(ALL). This also holds for other relativistic approaches.

Next, we study the neighboring-atomic-block (NB) approximation at the example of the

$I5+$ molecule. As shown in Figure 1, it is composed of a chain of iodine atoms and thus suitable for applying the neighboring approximation. Only the pairs of atoms which are connected with a bond displayed in Figure 1 are considered as neighbors. The neighboring groups could also be determined by introducing a cut-off distance parameter. Then, the pair of atoms which has shorter distance than the given parameter is counted as a neighboring pair. The relativistic transformations are now only applied to such neighboring pairs according to Eq. (22) within the DLU scheme. As we can see from Table IV, the neighboring approximation DLU(NB) gives very small errors for relative energies. Therefore, it is a very good (additional) approximation within our DLU scheme. Hence, with DLU(NB) the computational cost will become linear scaling, while the DLU(ALL) scheme has an order-N² scaling, which may be a bottleneck for calculations on very large molecules. But be aware that the success of the NB approximation depends on how one defines the neighboring atoms. If the cut-off distance is too small, the DLU scheme will be reduced to the DLH scheme since no neighbors would be found.

Reduction in terms of computational costs can also be achieved by employing the BP or NR approximation to the light atoms. It should be a good approximation if a molecule has many light atoms such as hydrogen. We study these approximations for the phosphine ligand dissociation energy of WH₆(PH₃)₃. The results are presented in Table V, where DLU(WP,NR) denotes the NR approximation applied to the H atoms, while W and P atoms are relativistically described. By contrast, DLU(W,NR) denotes application of the relativistic transformation only to the W atom, while for P and H the NR approximation applies. From Table V we can see that the results of the BP approximation are very bad, especially if the BP approximation was employed for the phosphorus atoms. The BP approximation to P and H atoms gives errors which are roughly seven times the value for the reactions energy. Hence, this approximation turns out to be completely useless. Such huge errors stem from the variational collapse of the BP approximation.

However, the NR approximation is much better. The errors on the reaction energy within the DLU(NR) scheme, which are in the range of 0.1–0.3 mH and thus less than 1 kJ/mol, are perfectly acceptable. This success is due to the correct nonrelativistic limit of the NR approximation within the DLU scheme. Not unexpectedly, the DLH(NR) results have quite larger errors because DLH(ALL) already does not provide good results for the phosphine dissociation energy from the WH₆(PH₃)₃ complex. If we study the difference between DLH(NR) and DLH(ALL), we find that the NR approximation actually gives small errors on top of DLH(ALL). However, the DLH scheme is not recommended since its accuracy is not guaranteed at all. Only in those cases for which the DLH scheme may work, we can further add the NR approximation on top of it for light atoms. As can be seen from both Tables IV and V, different relativistic approaches behave basically identical with respect to the discussion of the NB and NR approximations. Therefore, the approximations for the reduction of computational cost discussed in this paper can be applied to all relativistic transformation approaches.

In this work, we aimed at a rigorous local approximation to relativistic transformation schemes, which makes them applicable to and efficient for calculations on large molecules. We developed a systematic hierarchy of approximations that is based on the assembly of the unitary transformation from “atomic” contributions (DLU) rather than on a local approximation directly applied to the matrix representation of the Hamiltonian (DLH).

The straightforward DLH scheme, which has evolved in the development of the DKH method (see references given above), turns out to be not a very accurate local approximation since it only covers the relativistic treatment of the atom–same-atom diagonal blocks in the Hamiltonian. It may fail if interatomic distances become short so that a relativistic treatment of the off-diagonal “interatomic” blocks becomes important. As a consequence, the DLH approximation is not suitable, for example, in studies that crucially depend on accurate electronic energy surfaces. Even if the molecular geometry is fixed, the DLH approach does not provide a balanced description for different molecular properties, since the contributions of the off-diagonal blocks to total observables are quite different for different properties.

By contrast, the DLU scheme proposed by us in this paper overcomes the drawbacks of the DLH scheme. It is an excellent approach to take into account the atom–other-atom (i.e., “interatomic”) off-diagonal relativistic transformations. It does not show an instability for varying molecular structures (i.e., for electronic energy surfaces) and properties. The size of the error introduced by DLU is much smaller than that of DLH; it is typically only 1% of the latter. The errors of the DLU approach when compared to the full relativistic transformation without any local approximation turn out to be tiny. They are even smaller than the difference among the different exact-decoupling approaches, which are already very small.

The DLU scheme is valid for all exact-decoupling approaches, while the local X-matrix scheme only works for the X2C approach. Furthermore, the DLU scheme can also be applied to finite-order DKH approaches, such as DKH2, which has been very successful in computational chemistry and whose computational costs are less than those of any exact-decoupling approach.

If one has to use the DLH scheme for some reason, such as the lack of a DLU implementation or because of its linear-scaling behavior, the selection of a relativistic approaches is not decisive as the errors introduced by the DLH scheme are already higher than the difference of different relativistic approaches. Hence, the approximate, but fast DKH2 method may be safely combined with the DLH approximation, while this cannot be recommended for any exact-decoupling approach.

However, the linear-scaling behavior of DLH is no good reason since the DLU scheme can also be made linear scaling within the neighboring-atomic-blocks approximation, which produces negligible errors. One can further reduce the computational costs by employing the BP or NR approximation to all light atoms. However, the BP approximation turns out to be not suitable for variational calculations since it introduces large errors. Only the NR approximation yields acceptable results. The errors of the NR approximation will be larger if the nuclear charge of the atom for which the NR approximation is invoked increases. Therefore, the application of the NR approximation will depend on the balance of accuracy and computational cost in an actual calculation.

While this work focused on the accurate approximation of exact-decoupling methods at the self-consistent-field level with respect to total electronic energies and first-order properties of a small, but well-selected test molecule set, further investigations into the DLU scheme concerning, for example, higher order properties, spin–orbit, and electron-correlation effects, are under way in our laboratory.

This work has been financially supported by the Swiss National Science Foundation.