Relativistic DFT Calculations of Hyperfine Coupling Constants in 5d Hexafluorido Complexes: [ReF6 ]2- and [IrF6 ]2.

The performance of relativistic density functional theory (DFT) methods has been investigated for the calculation of the recently measured hyperfine coupling constants of hexafluorido complexes [ReF6 ]2- and [IrF6 ]2- . Three relativistic methods were employed at the DFT level of theory: the 2-component zeroth-order regular approximation (ZORA) method, in which the spin-orbit coupling was treated either variationally (EV ZORA) or as a perturbation (LR ZORA), and the 4-component Dirac-Kohn-Sham (DKS) method. The dependence of the results on the basis set and the choice of exchange-correlation functional was studied. Furthermore, the effect of varying the amount of Hartree-Fock exchange in the hybrid functionals was investigated. The LR ZORA and DKS methods combined with DFT led to very similar deviations (about 20 %) from the experimental values for the coupling constant of complex [ReF6 ]2- by using hybrid functionals. However, none of the methods were able to reproduce the large anisotropy of the hyperfine coupling tensor of complex [ReF6 ]2- . For [IrF6 ]2- , the EV ZORA and DKS methods reproduced the experimental tensor components with deviations of ≈10 and ≈5 % for the hybrid functionals, whereas the LR ZORA method predicted the coupling constant to be around one order of magnitude too large owing to the combination of large spin-orbit coupling and very low excitation energies.

Introduction

The magnetic properties of the seemingly simple hexafluorido complexes of ReIV and IrIV, [ReF6]2− and [IrF6]2−, have recently attracted a lot of interest because they were shown to be suitable building blocks for single‐molecule magnets.1, 2 They also exhibit an interesting electronic structure, which characterizes their single‐ion magnetic behaviour: a large zero‐field splitting (ZFS) in the case of [ReF6]2−, which is induced by strong spin–orbit coupling, and in the case of [IrF6]2−, an equally strong spin–orbit coupling without ZFS.

EPR spectroscopy is an essential tool for the structural and magnetic characterization of transition‐metal complexes; the spectra are conventionally interpreted in terms of three parameters: the electronic g‐tensor, the hyperfine coupling tensors, and the zero‐field splitting tensor. Accurate calculations of these EPR parameters are desired to predict the magnetic properties of similar compounds as well as to facilitate a detailed interpretation of the EPR experiments. In addition, quantum chemical methods can provide detailed information about the given property, such as the magnitude of the contributions from different mechanisms or the dependence on, for example, geometry or environment. In this work we have focussed on the hyperfine coupling.

Several mechanisms contribute to the hyperfine coupling constant, and they have to be treated properly to obtain results in good agreement with experimental values. In non‐relativistic theory, the only mechanism of the isotropic hyperfine coupling is the Fermi‐contact interaction, which describes the interaction between the electronic and nuclear spin when an electron is located at the position of the nucleus. This interaction requires that the spin density in the area close to the nucleus is properly described by the computational method, which affects several aspects of a quantum chemical calculation. The quality of the basis set is critical; in particular, the number of tight s‐functions is important, which has led to the development of specialized core‐property basis sets for the lighter transition metals.3, 4, 5 It has been shown that a Gaussian description of the nuclear charge (as an alternative to the standard non‐relativistic point charge model) can significantly improve the agreement with experimental results.6, 7, 8 In density functional theory calculations of hyperfine couplings, the spin polarization of the core electrons should be properly described by the exchange‐correlation functional.9 Several studies have shown that the qualitative performance of various exchange‐correlation functionals differs from system to system; therefore, the dependence on the exchange‐correlation functional of a desired property must be investigated for each class of systems.5, 9, 10

When considering the hyperfine coupling in complexes that involve elements with large nuclear charges, such as Re and Ir, relativistic effects have to be taken into consideration.11, 12 Relativistic contributions to hyperfine couplings are often divided into two parts: a scalar relativistic contribution and a spin–orbit coupling contribution. In the first row of the transition metals, scalar relativistic effects dominate, and therefore, more studies have considered this effect, whereas spin–orbit coupling effects were ignored or treated in an approximate manner. In the second and especially the third row of the transition metals, spin–orbit coupling effects can dominate the hyperfine coupling interaction; therefore, their accurate description is necessary.13, 14 The highest currently available level of theory with respect to relativistic effects is to start from the Dirac equation, to treat the nuclear magnetic moments as perturbations to first order and to use some approximation of the electron–electron repulsion. The most common approach is to use the Dirac–Coulomb Hamiltonian; herein, we will call this a fully relativistic 4‐component calculation. A more approximate treatment is to work only with the large component of the Dirac spinor, that is, a 2‐component treatment, which can be derived from the Dirac equation by first eliminating the small component of the four‐component wave function followed by the regular approximation approach.15 This is frequently employed in the form of the zeroth‐order regular approximation (ZORA).16 Scalar relativistic terms are always included in these treatments, whereas spin–orbit coupling can either be included variationally or by perturbation theory on top of the perturbation from the nuclear magnetic moments, which leads to a second‐order perturbation theory treatment.

In this work, we have compared the performance of three relativistic methods in the calculation of metal hyperfine coupling tensors in complexes [ReF6]2− and [IrF6]2− at the DFT level: the fully relativistic 4‐component method with electron–electron repulsion approximated by a Coulomb interaction in the form of the 4‐component Dirac–Kohn–Sham (DKS) method, the 2‐component ZORA method with variational treatment of spin–orbit coupling (EV ZORA), and the 2‐component ZORA method that treats spin–orbit coupling as another perturbation, that is, by linear response theory (LR ZORA). We have chosen the two minimally covalent 5d transition‐metal complexes [ReF6]2− and [IrF6]2− as test systems to investigate what challenges their complex electronic structure, that is, the large spin–orbit coupling in both and the large zero‐field splitting in [ReF6]2−, pose for the three theoretical methods. In addition, we review the theoretical background of the three methods in a common notation with focus on the calculation of the different contributions to the hyperfine coupling tensors. Specifically, this is the first time that the expressions for the contributions that arise in the DKS method are presented.

The article is organized as follows: initially, we have reviewed the theory of hyperfine coupling within the three relativistic methods; then in the following three sections, we have reported the details of our calculations, discussed the results, and drawn our conclusions.

A short review of the theoretical background

ZORA with variational treatment of spin–orbit coupling

The 2‐component ZORA hyperfine Hamiltonian consists of the derivative with respect to a component of the nuclear spin, I , of the 2‐component ZORA Hamiltonian16 in the presence of magnetic fields (SO‐ZORA),13, 14 which is given (here and in the following) in the Hartree system of atomic units [Eq. (1)]:

in which , σ is the vector that consists of the three Pauli spin matrices, is the momentum operator, is the position vector of the electron relative to the position of nucleus N, and is the length of r . The subscript u in Equation (1) indicates that the operator is perturbed by the u component of the nuclear spin, I . In the non‐relativistic limit, for which K→1, the first term in Equation (1) takes the form of the nuclear‐spin/electron‐orbit operator, and the last two terms reduce to the usual Fermi‐contact and spin‐dipolar operators. When evaluating the derivatives of the last two terms in Equation (1), a purely relativistic operator arises [Eq. (2)] (last two terms) as well as the ZORA analogues of the Fermi‐contact and spin‐dipolar operators (first two terms):17

In the DFT‐based EV ZORA implementation in the ADF program,18 the hyperfine coupling tensor was calculated in a restricted fashion on the basis of a Kramer's pair of the singly occupied molecular orbital (φ 1 and φ 2).13, 14 Consequently, this method applies only to doublet systems and the elements of the hyperfine coupling tensor of nucleus N, A , are calculated as described in Equations (3), (4), (5):

in which g is the nuclear g‐factor of nucleus N and μ the nuclear magneton. Notably, the current implementation does not allow for the calculation of the individual terms in the ZORA hyperfine Hamiltonian in Equation (2).

ZORA with perturbative treatment of spin–orbit coupling

In the LR ZORA formalism that was implemented in the ADF program,8, 18 spin–orbit coupling was treated as a further perturbation, and an element of the hyperfine coupling tensor was calculated as a second derivative of the perturbed energy with respect to both nuclear and electron spin [Eq. (6)]:

By applying the Hellmann–Feyman theorem to the perturbed DFT energy, one obtains an expectation value of the derivative of the hyperfine ZORA Hamiltonian in Equation (1) over the unperturbed ground state. As the unperturbed molecular orbitals were obtained with the scalar relativistic ZORA Hamiltonian, the spin degree of freedom can easily be integrated out, which allows for the spin derivative in Equation (6) to operate inside the integral. Consequently, the expression for an element of the first‐order contribution to the hyperfine coupling tensor reduces to the ZORA Fermi‐contact and spin‐dipolar contributions [Eq. (7)]:

in which is the unperturbed scalar relativistic spin density matrix, see ref. 8 for the exact expression and more details, n and n are the number of electrons with α and β spin and are the one‐electron basis functions. The scalar‐relativistic (spin‐free) ZORA (SR‐ZORA) Fermi‐contact and spin‐dipolar operator takes the following form [Eq. (8)]:

Notably, the two purely relativistic contributions (the two last terms) are considered as a part of the Fermi‐contact and spin‐dipolar contributions in the LR ZORA formalism.

The integrals of the ZORA PSO operator over the unperturbed orbitals will not contribute to the total hyperfine coupling tensor. To include this contribution, one has to let the molecular orbitals be perturbed by the spin–orbit interaction, that is, to solve the coupled perturbed Kohn–Sham equations. Details of this approach can be found in ref. 17 in the case of nuclear spin–spin coupling and ref. 19 for the implementation of molecular g‐tensors. Thus, an element of the spin–orbit contribution to the hyperfine coupling tensor can be written in terms of the perturbed spin density matrix, [Eq. (9)]:

in which the spin‐free ZORA paramagnetic spin–orbit operator is given in Equation (10):

‐Component Dirac–Kohn–Sham method

Finally, in the 4‐component Dirac–Kohn–Sham method, which is implemented in the ReSpect program,20 an element of the hyperfine coupling tensor is calculated as the derivative of the perturbed DFT energy with respect to a component of the nuclear spin [Eq. (11)]:6

in which is the effective spin of the system (1/2 for doublets, 1 for triplets etc.), J is the v component of the magnetization vector, and I is the spin of nucleus N. Applying the Hellmann–Feynman theorem led to an expectation value of the derivative of the Dirac Hamiltonian with respect to nuclear spin [Eq. (12)]:

in which denotes the occupied 4‐spinor molecular orbitals. The hyperfine Dirac–Kohn–Sham Hamiltonian [Eq. (13)] is the simple cross product between the Dirac matrix, α, which consists of the three Pauli spin matrices [Eq. (14)] and the relative vector between the electron and nuclear positions, r :

The superscript J indicates that the molecular orbitals depend on the chosen direction of magnetization, v, as described in refs. 21 and 22. In ref. 6, the use of a restricted kinetically balanced basis set is described, which leads to the following final expression [Eq. (15)]:

for which the density matrix is given in Equation (16):

and are the i′th molecular orbital coefficients of the large and small part of the 4‐spinor orbitals. consists of matrix elements over the large and small component basis functions with the off diagonal part of the operator in Equation (13). After applying restricted kinetic balance in the construction of the small component basis, these matrix elements turn out to consist of four contributions, all of which are explicitly calculated in the current implementation. These terms take the usual form of the Fermi‐contact, spin‐dipolar, and paramagnetic spin–orbit operators (FC, SD and PSO) as well as a purely relativistic contribution (here denoted REL) analogues to the last two terms of Equation (2) [Eqs. (17), (18), (19), (20)]:

has contributions from the PSO operator and the negative of the REL operator only. The resulting four contributions to an element of the hyperfine coupling tensor are consequently calculated as follows [Eqs. (21), (22), (23), (24)]:

Computational details

Geometry and spin state

Herein, we have calculated the EPR hyperfine coupling tensors of ReIV and IrIV in the [ReF6]2− (1) and [IrF6]2− (2) ions. The available experimental EPR spectra1, 2 were measured on the diamagnetic coordination polymers with alternating metal centres, Zn(viz)4[ZrF6], that were doped with ≈5 and ≈1 % of the fluorido complexes 1 and 2 respectively (Figure 1). The ZrIV chain that acts as the host for the ReIV and IrIV fluorido complexes is isostructural with the undiluted Re and Ir chains and crystallizes in the P42/n space group with the tetravalent metal centres on special positions with tetragonal symmetry. This ensures strictly axial and co‐parallel g‐ and A‐tensors. The geometries used in all calculations presented in this work are the experimental X‐ray structures of complexes 1 and 2 in the polymer Zn(viz)4[MF6] (M=Re, Ir) without further geometry optimization.

Figure 1

Illustration of the coordination polymer Zn(viz)4[MF6] (M=Zr, Re, Ir). Re, Ir: turquoise, Zn: purple, F: green, C: yellow, N: blue, H: white.

The chain structure caused a slight elongation of the M−F bonds in complexes 1 and 2 along the chain direction, which resulted in axial symmetry. As a consequence, two metal hyperfine parameters were extracted from the experimentally determined EPR spectra,1, 2 which are denoted the parallel and perpendicular components A ∥ and A ⊥. In our calculations, A ∥ and A ⊥ were determined from the principal components of the hyperfine coupling tensor. Two of the principal components will be equal and thus resemble A ⊥, and the unique component will be A ∥.

The hyperfine coupling tensor does not transform as a true tensor, and the principal values were determined as the square roots of the eigenvalues of the “squared” hyperfine tensor .23 The eigenvectors of this true tensor then make up its so‐called principal axis system. Following this procedure, the principal values were absolute values. The isotropic hyperfine coupling constant was calculated by using Equation (25) :

The individual contributions to the isotropic hyperfine coupling constant should, in principle, also be calculated by following the above described procedure. However, as the relative signs of the contributions would be lost, these values were based on the principal values of the respective symmetrized tensors, .

For both compounds, the expected spin states have previously been established by magnetometry;1, 2 specifically, the ground states of [ReF6]2− and [IrF6]2− are S=3/2 and S=1/2, respectively. These spin states have been assumed in the computations, and the most abundant isotopes with I≠0 were used: 187Re (I=5/2; 62.7%) and 193Ir (I=3/2; 62.7%).

Programs and basis sets

The 2‐component ZORA calculations that treat spin–orbit coupling either variationally (EV ZORA) or as a perturbation (LR ZORA) were carried out with the ADF program8, 14, 18 at the DFT level. The EV ZORA method was implemented in a spin‐restricted fashion, and accordingly, spin polarization effects were not taken into account. In contrast, the LR ZORA method employs the unrestricted scheme and uses the collinear description of spin density. For consistency, the collinear option was also employed in the EV ZORA calculations. The Slater‐type orbitals (STO) basis sets, DZ, TZ2P, and QZ4P, which were optimized for ZORA calculations, were employed. Convergence in all 2‐component ZORA calculations was obtained with a threshold of 4×10−3 for the linear dependence of basis functions.

The unrestricted 4‐component calculations were carried out within the Dirac–Kohn–Sham (DKS) framework, as implemented in the ReSpect program with restricted kinetic balanced basis sets and non‐collinear spin density.6, 20 Dyall's Gaussian‐type orbital (GTO) valence basis sets, vdz, vtz, and vqz, were used.24, 25

Both the ZORA and the DKS results were obtained with the GGA exchange‐correlation functional BP8626, 27 as well as the hybrid functionals B3LYP28 and PBE029 with 20 and 25 % Hartree–Fock exchange, respectively. In all self‐consistent‐field (SCF) calculations, a Gaussian‐type nucleus model was employed along with an SCF convergence criteria of 10−6.

Results and discussion

Basis set dependence

Calculations of magnetic properties are known to put extra requirements on the quality and size of basis sets compared with simple calculations of energies.30 For example, previous studies showed that, in ZORA calculations of nuclear magnetic shielding constants, it was important to use the large QZ4P basis sets.31, 32 Therefore, we also studied the basis set dependence of both the ZORA and DKS results for the hyperfine coupling constants in both ions. Table 1 shows the LR ZORA and EV ZORA results of the isotropic hyperfine coupling constant of ReIV and IrIV complexes, respectively, which were obtained with the special ZORA STO basis sets DZ, TZ2P, and QZ4P, as well as DKS 4‐component results, which were obtained with Dyalls vdz, vtz, and vqz GTO basis sets.

Table 1

Basis set dependence of the isotropic metal hyperfine coupling constant A (MHz) in complexes [ReF6]2− and [IrF6]2− calculated at the ZORA (LR ZORA for [ReF6]2− and EV ZORA for [IrF6]2−) and DKS levels by using the PBE0 functional.

Method	Basis set	[ReF6]2−	[IrF6]2−
ZORA	DZ	1231	81.4
	TZ2P	1224	86.0
	QZ4P	1324	82.8
DKS	vdz	1283	97.1
	vtz	1303	97.9
	vqz	1303	98.1

Evidently, the isotropic hyperfine coupling constant changed only slightly but in a consistent manner when the basis set was increased by using the DKS method and Dyall's basis sets, whereas we observed the same irregular, no‐monotonic behaviour of the ZORA basis sets, as seen previously, for nuclear magnetic shielding.31, 32 It was shown there, that the increasing cardinal number and the increasing amount of polarization functions have opposite effects, which lead to an oscillating and not yet converging behaviour. Therefore, the change from the TZ2P to the QZ4P basis set in the ZORA calculations still led to a change in the isotropic hyperfine coupling constants of 8 and 3 % for [ReF6]2− and [IrF6]2−, respectively. Therefore, in the following sections, the ZORA results are all with the QZ4P basis set, whereas the DKS results are with the vdz or vtz basis set, which is stated in each section, as we see only a small change between the results with the vdz and vqz basis set.

Effect of exchange‐correlation functional and amount of Hartree–Fock exchange

In Table 2, the metal isotropic hyperfine coupling constants calculated at the 2‐ and 4‐component levels with the GGA functional, BP86, and the hybrid functionals, B3LYP and PBE0, are shown along with the experimental values. It was not possible to obtain acceptable convergence with the BP86 functional in combination with the DKS method.

Table 2

Exchange‐correlation functional dependence of the isotropic metal hyperfine coupling constant A (MHz) in complexes [ReF6]2− and [IrF6]2− calculated at the ZORA (LR ZORA for [ReF6]2− and EV ZORA for [IrF6]2−) and DKS levels by using the QZ4P and vtz basis sets, respectively. Relative deviations from the experimental values calculated as are given in parentheses.

	ZORA				DKS				Exp.
	BP86	B3LYP	PBE0		BP86	B3LYP	PBE0
[ReF6]2−	971 (40 %)	1240 (23 %)	1324 (18 %)		940 (42 %)	1193 (26 %)	1303 (19 %)		1607
[IrF6]2−	79.5 (14 %)	83.3 (10 %)	82.8 (10 %)		–	97.1 (5 %)	97.9 (6 %)		92.4

Three conclusions can be drawn when the relative deviations from the experimental values are considered. Firstly, for both relativistic methods and all three functionals, the deviations from the experimental values were considerably smaller for complex [IrF6]2− than for [ReF6]2−. Secondly, for complex [IrF6]2−, the 4‐component DKS results were in better agreement with the experimental values than the EV ZORA results were. For the [ReF6]2− system, the approximate 2‐component LR ZORA results seemed to be in marginally better agreement with the experimental values than the 4‐component DKS results. However, neither of the two methods in combination with DFT performed very well for this complex. Thirdly, for both complexes and both relativistic methods, the hybrid functionals reduced the deviation from the experimental values compared with the pure GGA functional.

This difference between GGA and hybrid exchange‐correlation functionals was investigated further by following the procedure of other recent studies,10, 33 in which it was shown that raising the amount of Hartree–Fock exchange in the PBE0 functional increased the agreement with the experimental values. In Figure 2, it is shown how the perpendicular and parallel components of the hyperfine coupling tensor, as well as the isotropic hyperfine coupling constant, depend on the amount of Hartree–Fock exchange in the PBE0 functional. The results were obtained at the 4‐component DKS level with the vdz basis set (circles) and the 2‐component ZORA level with the QZ4P basis set (triangles). The left axis shows the absolute values and the right axis shows the relative deviation of the calculated isotropic results from the experimental isotropic values. The experimental values are shown as solid lines.

Figure 2

Dependence of A ⊥, A ∥, and A in [ReF6]2− (a) and [IrF6]2− (b) on the amount of Hartree–Fock exchange in the PBE0 functional as well as the percentage deviation of A from the experimental values. Calculations were carried out on the 4‐component level with the vdz basis set (circles) and the 2‐component ZORA level with the QZ4P basis set (triangles).

The [ReF6]2− results show a linear dependence of the parallel, perpendicular, and isotropic hyperfine coupling constants on the amount of Hartree–Fock exchange. The isotropic value agrees with the experimental values at around 50 % Hartree–Fock exchange. However, the experimental splitting of the parallel and perpendicular components of the hyperfine coupling tensor, as illustrated by the separation of the top and bottom grey lines, was not at all reproduced by the DFT calculations, regardless of the amount of Hartree–Fock exchange.

The results for complex [IrF6]2− did not exhibit the same linear trend as those for [ReF6]2−. From the DKS values, we can conclude that increasing the amount of Hartree–Fock exchange increases the deviation from the experimental values instead of reducing it. The opposite was true for the ZORA values. Furthermore, the calculated DKS values for the parallel and perpendicular components follow slightly different trends, which resulted in a reordering of the components around 40 % Hartree–Fock exchange. The same behaviour was seen when the ZORA method was used; however, it was not possible to obtain good convergence with a large amount of exact exchange.

For complex [IrF6]2−, the range of the deviation from the experimental values was around ±5 to ±12 %, which is considerably smaller than the deviation that was observed for complex [ReF6]2− (−20 to +15 %). The admixture of Hartree–Fock exchange in hybrid functionals is known to enhance the spin polarization of the core‐shell electrons,9 an effect that is often underestimated by GGA functionals, and one would expect this enhancement to affect the Fermi contact interaction more strongly than both spin‐dipolar and spin–orbit coupling interactions. As we will discuss in the following section, the hyperfine coupling in complex [ReF6]2− was in fact dominated by the Fermi contact term, whereas the spin–orbit coupling interaction dominated the hyperfine coupling of complex [IrF6]2−. This presumably leads to the different behaviour of the two complexes with respect to the amount of Hartree–Fock exchange, although in absolute terms, the spin–orbit contribution in complex [ReF6]2− was actually larger than in [IrF6]2−.

Importance of the treatment of the spin–orbit coupling contribution

In Table 3, the contributions to the parallel, perpendicular, and isotropic components of the hyperfine coupling tensors, which were calculated by using the three methods LR ZORA, EV ZORA, and DKS and the standard PBE0 functional, are listed. As the EV ZORA method only applies to systems with only one unpaired electron, no results are shown for complex [ReF6]2−, which has a total spin of S=3/2. Furthermore, no contributions are shown for the EV ZORA results, because they were not calculated explicitly. Above, we showed that a purely relativistic correction to the hyperfine coupling tensor arises in both ZORA and DKS theory. In the LR ZORA implementation in ADF, this contribution was considered as a part of the Fermi contact (FC) and spin‐dipolar (SD) contributions, whereas in the DKS implementation in ReSpect, this contribution was calculated explicitly (REL).

Table 3

The individual contributions to the parallel A ∥, perpendicular A ⊥, and isotropic components A of the hyperfine coupling tensor calculated with the ZORA and DKS methods by using the QZ4P and vtz basis set, respectively, and the PBE0 functional.

	LR ZORA			EV ZORA	DKS						Exp.
	FC+SD	PSO/SO	Total	Total	FC	SD	REL	F+S+R	PSO	Total
[ReF6]2−
A ∥	1087	234.7	1322	–	1672	−15.05	−591.1	1066	219.0	1285	1850
A ⊥	1087	237.6	1325	–	1674	−3.049	−591.7	1079	233.1	1312	1486
Aiso	1087	236.7	1324	–	1673	−7.948	−591.5	1074	228.4	1303	1607
[IrF6]2−
A ∥	−145.4	−87.81	−233.2	81.69	−43.62	−11.88	16.07	−39.43	137.1	97.62	93.1
A ⊥	−66.49	1616	1549	83.37	−52.43	−10.91	19.28	−44.06	142.1	98.02	92.0
Aiso	−92.79	1048	955.3	82.81	−49.50	−11.23	18.21	−42.52	140.4	97.89	92.4

The results for complex [ReF6]2− showed how the components of the hyperfine coupling tensor are indeed dominated by the Fermi contact contribution. If the DKS SD and REL contributions are added to the DKS FC contribution (Table 3, F+S+R), one can see that the LR ZORA results are very close to the DKS results, particularly the perpendicular component A ⊥. The slightly worse agreement for the parallel component A ∥ was caused by the spin‐dipolar and spin–orbit contributions, which in DKS theory, deviate slightly from the corresponding A ⊥ contributions, but in LR ZORA theory, only the PSO/SO components differ by a very small amount.

Compared with the experimental values, one observes that A ⊥ deviates for both methods by only around 10 %, whereas a much larger deviation of around 30 % was calculated for the parallel component A ∥. Furthermore, the ordering of the two components is reversed compared with the experimental results, in which A ∥≫A ⊥. The fact that neither of the two DFT‐based methods were able to reproduce the anisotropy of the hyperfine coupling tensor of complex [ReF6]2− might be because the methods did not take into account the splitting of the quartet (S=3/2) into two effective doublet states, also known as zero‐field splitting.

The components of the hyperfine coupling tensor of complex [IrF6]2− at the 4‐component DKS level were dominated by the paramagnetic spin–orbit coupling contribution, which contrary to the results for complex [ReF6]2−, had the opposite sign of the Fermi contact contribution. The EV ZORA results were similar to the DKS and experimental results, but for this method, it was not possible to obtain individual contributions. However, the LR ZORA results for complex [IrF6]2− were very different from the EV ZORA, DKS, and experimental results. The magnitude of the Fermi contact plus the spin‐dipolar contribution to A ⊥ was similar to the DKS results, whereas it was almost a factor of 3 too large for A ∥. However, the main difference originated in the spin–orbit contribution; compared to the DKS results, LR ZORA predicted a contribution to A ∥ of the wrong sign and a contribution to A ⊥ that was more than a factor of 10 too large.

The spin–orbit coupling contribution with the LR ZORA method was, as previously mentioned, calculated as a perturbation through first order, that is, by linear response theory, on top of a scalar relativistic result. This approach, which neglects higher order corrections, is only applicable if the perturbation, that is, the spin–orbit coupling, is small compared with the energy difference between the ground and excited states.8, 13, 19 Therefore, we have calculated the excitation energies in complex [IrF6]2− (and [ReF6]2− for comparison) by using the implementation of unrestricted TDDFT and the scalar relativistic ZORA Hamiltonian in ADF.34 The two lowest excitation energies for both complexes calculated at the PBE0/QZ4P level are presented in Table 4. The transitions were assigned by inspection of the corresponding eigenvectors. The symmetry labels of the molecular orbitals were assigned on the basis of the irreducible representations in the D 4 point group.

Table 4

Excitation energies calculated at the PBE0/QZ4P level with the scalar relativistic ZORA method and unrestricted TDDFT.

	Transition	Dominant orbital excitation		Rotations	ΔE [eV]
[ReF6]2−	1	α b 2g→a 1g	99.5 %		3.06
	2	α eg→a 1g	96.6 %	Rx, Ry	3.10
[IrF6]2−	1	β eg→b 2g	98.7 %	Rx, Ry	0.10
	2	β eg→b 1g	53.4 %	Rx, Ry	2.90
		α eg→b 1g	38.3 %

In agreement with a simplistic ligand field picture, the lowest calculated excitation energy in the orbitally pseudodegenerate complex [IrF6]2− was very small compared with the lowest values for [ReF6]2−, which might cause the substantial overestimation of the spin–orbit coupling contribution for complex [IrF6]2− when calculated as a linear response function. This behaviour was previously observed for actinide systems,13 but not to this extent. Notably, this transition can only cause such an overestimation if the corresponding matrix elements over the perturbation operators are non‐vanishing. In the present implementation in ADF, it is, to our knowledge, not possible to calculate transition moments of the PSO or spin–orbit operators, but we can at least investigate whether they are symmetry allowed or forbidden.

The spin–orbit operator in Equation (10) contains the angular momentum operators ( ), which must transform like the rotation functions. Accordingly, in the D 4 point group, they transform in accordance with Equation (26):

Evaluating the direct products of the irreducible representations of the rotation functions and the orbitals involved in the transitions (Table 4) showed that the second lowest transition in complex [ReF6]2− and the two lowest transitions in complex [IrF6]2− transform like the rotation functions R and R. This means that the transition moments over the corresponding components of the angular momentum operators are symmetry allowed. Inspection of the principal values of the hyperfine coupling tensor of complex [IrF6]2− that was calculated with the LR ZORA method (Table 3) showed that it is indeed the perpendicular component, which corresponds to the x‐ and y‐directions, that is drastically overestimated.

The tetragonal symmetry imposed on the two ions by the Zn(viz)4[ZrF6] host resulted in a splitting of the t 2 level into an e and b 2 level. Assuming that this splitting is small, the 4 A ground state in complex [ReF6]2− was unaffected and the excitation from the e to b 2 level is spin forbidden in a scalar relativistic treatment. In complex [IrF6]2−, this excitation is allowed and explains the very small, first excitation energy, which is doubly degenerate.

Conclusions

We have investigated the performance of the 2‐component EV ZORA and LR ZORA methods and the 4‐component DKS method in the framework of DFT for the calculation of the metal hyperfine coupling constants in complexes [ReF6]2− and [IrF6]2−. Both the basis set dependence and the dependence on the choice of exchange‐correlation functional of the results as well as the effect of varying the amount of Hartree–Fock exchange in the hybrid functionals were studied.

Both the 2‐component EV ZORA and 4‐component DKS methods reproduced the hyperfine coupling tensor of complex [IrF6]2− in good agreement with experimental results. However, when the 2‐component LR ZORA method was applied to the [IrF6]2− ion, in which the spin–orbit coupling contribution was calculated by perturbation theory, the perpendicular component of the hyperfine coupling tensor was overestimated by around one order of magnitude and the opposite sign was predicted for the parallel component relative to the DKS method. This breakdown of a perturbation theory treatment through first order of the spin–orbit coupling is due to a very small energy difference between the ground and first excited states, which is caused by the tetragonal symmetry that is imposed on the ion in combination with the presence of large spin–orbit coupling.

For complex [ReF6]2−, the 2‐component LR ZORA and 4‐component DKS methods reproduced the perpendicular component of the hyperfine coupling tensor in good agreement with the experimental values, but they significantly underestimated the parallel component. The agreement with the experimental isotropic value for complex [ReF6]2− was improved further by increasing the amount of Hartree–Fock exchange in the PBE0 hybrid functional. However, the anisotropy of the hyperfine tensor, which was found in experimental results, was not reproduced by either the LR ZORA or the DKS method, regardless of the amount of Hartree–Fock exchange.

Overall, the LR ZORA method only applies to systems without low‐lying excited states ([ReF6]2−), and the current implementation of the EV ZORA approach can only be used for systems with one unpaired electron ([IrF6]2−). However, the DKS method produced reasonable results for both systems and allowed for an investigation of the individual contributions to the hyperfine coupling tensor.

Conflict of interest

The authors declare no conflict of interest.