Ab initio paramagnetic NMR shifts via point-dipole approximation in a large magnetic-anisotropy Co(ii) complex

Transition metal complexes can possess a large magnetic susceptibility anisotropy, facilitating applications such as paramagnetic tags or shift agents in nuclear magnetic resonance (NMR) spectroscopy. Due to its g-shift and zero-field splitting (ZFS) we demonstrate on a Co(II) clathrochelate with an aliphatic 16-carbon chain, a modern approach for ab initio calculation of paramagnetic susceptibility. Due to its large anisotropy, large linear dimension but relatively low number of atoms, the chosen complex is especially well-suited for testing the long-range point-dipole approximation (PDA) for the pseudocontact shifts (PCSs) of paramagnetic NMR. A static structure of the complex is used to compare the limiting long-distance PDA with full first-principles quantum-mechanical calculation. A non-symmetric formula for the magnetic susceptibility tensor is necessary to be consistent with the latter. Comparison with experimental shifts is performed by conformational averaging over the chain dynamics using Monte Carlo simulation. We observe satisfactory accuracy from the rudimentary simulation and, more importantly, demonstrate the fast applicability of the ab initio PDA.


Introduction
Molecules with large magnetic anisotropy, often due to a paramagnetic metal center, find applications in chemistry, materials research, structural biology and medicine. [1][2][3][4][5] Examples of effects that can be observed due to such centers are orientation in magnetic field 6 and the unusually large range of NMR chemical shifts and rapid relaxation rate of the NMR signals. 7,8 Among the compounds of 3d transition metals, some cobalt complexes feature the combination of large zero-field splitting (ZFS) and g-tensor parameters that are commonly expected to translate to a large magnetic anisotropy. 9 Computational methods for the prediction and analysis of the magnetic properties of such paramagnetic systems have progressed during recent years. The NMR chemical shifts induced by a strongly anisotropic magnetic center have already for a long time been calculated by a longdistance PDA of the dipolar hyperfine coupling (HFC) interactions between the unpaired electrons and NMR nuclei. 1 The PDA was recently implemented in a modern computational framework that allows ab initio pNMR shift calculation of large systems such as an entire metalloprotein, 9 and which is based on the recently redressed version (ref. 10 and 11) of the classic Kurland-McGarvey theory. 12 To further investigate the performance of the method of ref. 9, we selected for the present work a large-anisotropy Co(II) complex ( Fig. 1) studied by Novikov et al., 13 denoted in their paper as complex 2. Due to the large linear dimension of this S = 3 2 (quartet spin state) clathrochelate complex, it is well-suited for testing the methods for the long-distance effect of the paramagnetic center. For a system of this size, the required  Table 1, indicated, as well as the complex including both the cage and one of the 16-carbon tails with numbering of the tail. Parallel and perpendicular directions of the cylindrically symmetric cage system are indicated as 8 and >, respectively. The cartesian X, Y, Z coordinate frame is also depicted. electron paramagnetic resonance (EPR) parameters can be calculated quite reliably using molecular quantum-chemical (QC) codes. pNMR modeling of 3d transition metal clusters has been by now well-probed 10,11,[14][15][16] and chemical shifts in such systems can be calculated by first-principles methods with a predictive accuracy. We retained only one of the two aliphatic 16-carbon chains of the complex. Since the presence of the chain does not significantly influence the EPR properties of the complex, as shown below, and the two chains are too far from each other to restrict their mutual conformational spaces, removing one of the chains from the model can be regarded as a safe approximation.
In the PDA, the HFC tensor appearing in the contemporary QC theory of pNMR shielding 10-12 is approximated by its longdistance limit originating from the electron-spin magnetic moment, which is taken to be a point dipole. The magnetic moment of the NMR nucleus interacts with the field generated by the point dipole of the electron spin, giving rise to the approximate expression for the shielding tensor 9 with the isotropic pseudocontact shielding constant (s PCS ) and the corresponding PCS (d PCS ) defined as In eqn (1), r IS is the length of the vector r IS between nucleus I and the paramagnetic center S, n IS the corresponding normalised direction vector. v 0 is a (generally non-symmetric) magnetic susceptibility tensor to be discussed below. In this paper we demonstrate the validity of the PDA by comparing with full firstprinciples QC shift calculations on the present Co(II) complex and relate to experiment by performing a Monte Carlo-based averaging over the chain conformations.

Methods
For comparison with the PDA, we calculated quantum-chemically the full HFC tensors of the present Co(II) complex encompassing both the cage and the tail. In the used ORCA code, 17 we employ nonrelativistic density-functional theory (DFT) with leading-order perturbational correction for the relativistic spin-orbit (SO) interaction. 18 The importance of this correction for NMR shieldings in open-shell molecules has been recently pointed out, for example, by Marek and coworkers. 19 The HFC tensors, together with g-and ZFS tensors, are incorporated in the expression for pNMR shielding tensors as detailed in ref. 11 and 20. A breakdown of the various terms appearing in the method of ref. 10-12 and featuring different contributions to the EPR parameters, is presented in Table S1 in the ESI. † To assess the spatial reach of the different mechanisms of hyperfine shielding originating from the paramagnetic center out to the ligand atom, the pNMR shielding tensors are broken down to scalar (rank-0) and dipolar (rank-2) terms. We use two different structures. First, we used a truncated Co(II) clathrochelate complex with both 16-carbon tails replaced by hydrogens. The cage structure was computationally optimised by Turbomole 21 and Orca software using several settings of the DFT method to obtain a range of geometries ( Table 1).
The second structure used in this work consisted of Co(II) clathrochelate complex with one 16-carbon tail replaced by a hydrogen atom, retaining one tail. This was optimised by PBE 27,28 /def2-SVP 29 combination of DFT functional and basis set, and served mainly in comparing methods for the HFC tensor calculation. If not stated otherwise, the geometry used in all calculations of this study was obtained at the DFT PBE0 27,28 /def2-TZVP 29 level with a quasirelativistic effective core potential 25 on the Co centre and the COSX approximation 26 for the Hartree-Fock exchange.
In the present work we use a combination of the multireference complete active space self-consistent field (CASSCF) 30 method, augmented by dynamical correlation treatment by N-electron valence-state perturbation theory (NEVPT2 [31][32][33], to calculate the EPR g-and ZFS tensors. In these computations, the scalar relativistic second-order Douglas-Kroll-Hess (DKH) Hamiltonian 34,35 is employed, with SO effects treated by quasidegenerate perturbation theory. 36,37 In state-averaged CASSCF wave function, the active space was selected to cover the d-orbitals of the metal center -meaning, for Co(II), 7 electrons in 5 orbitals. All the 10 quartet roots as well as 15 doublets possible within this active space were included. It has been shown 9,15,38,39 that this methodology gives good results for the EPR parameters including the critical ZFS tensor, for systems with two or more unpaired electrons. Table 1 Geometrical parameters and magnetic properties of the hydrogen-terminated cage unit optimised at various density-functional theory levels a Optimised bond lengths, see illustration of the structure in Fig. 1. b ZFS-and g-tensor calculations at the CASSCF/TZVP-DKH level.
where D ii are the eigenvalues of the ZFS-tensor. The ordering of the principal values is selected so that 0 r E/D r 1/3. 22,23 c Dg ax ¼ g 33 À 1 2 g 11 þ g 22 ð Þ , where the g ii are the g-tensor eigenvalues, see also eqn (S1)-(S3) (ESI) and the associated explanation. d D3 dispersion correction. 24 e With scalar relativistic effective core potential. 25 f COSX approximation of Hartree-Fock exchange. 26 We show that the, for relativistic purposes recontracted variant 40 of the valence triple-z + polarization (TZVP) basis set 41 (denoted here as TZVP-DKH), is for current calculations of the g-and ZFS tensors equally good as the much larger def2-QZVPP-DKH basis 40 (Table 2). The related def2-TZVP-DKH basis has been found to represent a good compromise between the computational cost and accuracy in previous pNMR computations. 11,15,20 In accordance with our previous experience, 9,20 there is a significant impact on the g-and ZFS tensors of dynamical electron correlation accounted for via the NEVPT2 method, as compared to the CASSCF results (Table 2). Both the D -parameter of the ZFS tensor and the anisotropy of the g-tensor decrease upon adopting the NEVPT2 correction. For comparison with experimental data, it is important to use the NEVPT2 results if possible. 15,20 For production results, the g-and ZFS tensors were calculated using ORCA by the CASSCF/NEPVT2 methods and the TZVP-DKH basis. These numbers, obtained from the cage-only model unless otherwise stated, were subsequently used for the PDA calculations including the 16-carbon chain.
In conformational averaging calculations, the energy difference of 3.9 kJ mol À1 between the trans and gauche conformers of the carbon chain (with the gauche conformer higher in energy) for each bond was adopted from the literature. 13,42 The bond lengths were kept at the values obtained from the present QC (PBE/def2-SVP) optimisation, approximately 1.53 Å for the C-C bonds. The temperature was set to 300 K both for the MC simulation as well as for all the pNMR calculations of this study. It was assumed that the conformations of the bonds are independent of the conformations of all the other bonds. However, conformations involving overlap of the atomic radii were excluded. Pseudocontact shifts were calculated using PDA, eqn (2) for every atom and for each cycle of the Monte Carlo simulation, with the g-and ZFS-tensors of the rigid cage-only model calculated only once.

Structure
It has been seen before 9,15,38 that the geometry of the paramagnetic center can dramatically influence the resulting ZFS and g-tensors. Table 1 shows the effect of the detailed structure as optimised using the PBE and PBE0 functionals, with particularly the PBE0 functional being not only first-principles-based but also performing well for structures of transition metal complexes. 9,43 Various other choices are elaborated in Table 1. It can be concluded that, in this case, applying the dispersion correction 44 or using a relativistic pseudopotential 25 at the cobalt centre does not change the geometry or the EPR parameters significantly. For the structure optimisation, the resulting effects are smaller than those arising from the choice of the DFT functional.

Comparison of full theory and the point-dipole approximation
Full QC HFC tensors were calculated at the PBE0/def2-TZVP level with the shielding results shown in Table 3 and Fig. 2.
In the vicinity of the paramagnetic center, the electron spin density distribution (Fig. 3) translates into, among others, the Fermi contact part of the hyperfine coupling, 45,46 which dominates the pNMR shielding in this region. This necessitates a full QC approach involving the evaluation of HFC at the NMR nuclei. At larger distances, where the spin density vanishes, the contact contribution to the HFC decays to zero. This leaves the dipolar contribution, manifested in the isotropic pNMR chemical shift as PCS, 47 as the dominant mechanism. Since in the long-distance limit, the dipolar interaction takes a form familiar from classical magnetostatics and, hence, is independent on the bonding structure, the pseudocontact shifts can equally well be calculated for the 1-dimensional C16 chain of this study and, for example, for a 3D protein structure. For a doublet system, comparison of the PDA-and QCcalculated shifts has been recently discussed in ref. 48. In ref. 9, a formulation of PDA consistent with the full QC treatment of arbitrary spin state of ref. 10-12 and 20 was elaborated. According to that, the magnetic susceptibility tensor is written as where m B , m 0 , g e , g, k, T and hSSi are, in the respective order, the Bohr magneton, vacuum permeability, free-electron g-value, g-tensor, Boltzmann constant, absolute temperature and a dyadic calculated from the effective electron spin operator S. This dyadic is constructed as a statistical average involving the eigenvalues and eigenfunctions of the ZFS Hamiltonian, H ZFS = SÁDÁS. It should be noted that the functional form of eqn (3) entails, depending on the spatial symmetry of the system, the possibility of a slightly non-symmetric form of v 0 , i.e., one possessing a contribution in the tensorial rank 1. A related, symmetric susceptibility tensor appears in a molecular energy expression that is quadratic in the magnetic field,  Table 1. b This structure, optimised at PBE0(def2-SVP) level, has been used only for this comparison. c Structure obtained by removing the 16-carbon chain and terminating with a hydrogen atom. d g-and ZFS tensors obtained using this structure and combination of method and basis set were used for the production results presented for example in Tables 3, 4 and Table S2 (ESI). is accessible in magnetometric experiments, 1 and can be approximately obtained by symmetrisation of v 0 . 9 Using the breakdown of the pNMR shielding into different mechanisms, presented in Table S1 of the ESI, † the PDA includes the sum of hyperfine term 2 -a combination of the dipolar hyperfine coupling and the free electron g-value g e , term 7 due to the isotropic g-shift and term 9 due to the anisotropy of the g-tensor. All of these contributions to the full pNMR shielding expression involve the dipolar hyperfine interaction.
It can be seen from Table 3 as well as Fig. 2 and 4 that the PDA is for the dipolar part of the hyperfine interaction very well valid already for the C1 carbon of the 16-carbon chain, where the PCS still reaches a high value of 74 ppm (Table 3 and Fig. 3). In fact, already for the boron atom of the cage, at 3.1 Å from the metal center, the PDA (À266 ppm) differs from the PCS obtained by the QC calculation (À237 ppm) only by 12%. In a singly-bonded, aliphatic hydrocarbon chain such as in the present model complex, the spin density does not extend over a large distance from the paramagnetic centre. The disappearance of the spin density far away from the paramagnetic center is discussed further for example in ref. 9 and 49. According to Table 3, at the C3 position of the 16-carbon chain, the contact term still amounts to 2.32 ppm, which represents 12% of the total shielding constant. The contact term becomes negligible at around C4 with 0.3 ppm contribution remaining in the total s (Fig. 2). It is worth noting that the shielding terms 2 and 9 are of a very similar magnitude (Fig. S1 in the ESI †). The nonrelativistic dipolar term 2, which would have a vanishing contribution in the Table 3 Comparison of the full quantum-chemical (QC) calculation of the isotropic pNMR shielding constants (in ppm) with the point-dipole approximation (PDA) at T = 300 K. Results for the Co(II) clathrochelate system corresponding to the static PBE/def2-SVP-optimised structure of Fig. 1 Table 2). The molecular cage structure with the C16 chain was geometry-optimised at the PBE/def2-SVP level. The QC calculations of the HFC tensors were carried out at the PBE0/def2-TZVP level using the Orca programme. For PDA, the dipolar HFC tensors (column ''PDA'') were calculated using the same structure. b Sum of the pNMR shielding terms from the full calculation, which include the rank-0 contact hyperfine coupling. Contact QC = s con + s con,2 + s con,3 + s c,aniso . See Table S1 (ESI) for the details of the breakdown of the pNMR shielding tensors to the various physical contributions and the paragraph ''Anisotropy parameters'' in the ESI for clarification of the tensorial ranks of the HFC tensor. c Relative contribution of the contact terms (in %) to the total pNMR shielding. d Sum of the pNMR shielding terms from the full quantum-chemical calculation, involving the second-rank hyperfine coupling contribution, PCS QC = s dip + s dip,2 + s dip,3 + s pc . e Relative contribution of the dipolar terms (%) to the total pNMR shielding. f Sum of the pNMR shielding terms involving the second-rank (dipolar) hyperfine coupling contribution calculated using the point-dipole approximation: PDA = s dip + s dip,3 + s pc . The same g-and ZFS-tensors were used as for other calculations (see footnote a). g Relative error calculated as (''PDA'' -''Full QC'')/ ''Full QC''. S = 1 2 case, is represented here (due to the effect of the hSSi dyadic) as equally important as the traditional pseudocontact term 9 (arising due to the g-tensor anisotropy).

Susceptibility
Using eqn (3) and NEVPT2 calculation with the selected computational parameters listed in Table 2   In an alternative approach 50 one uses the symmetric relation where the full g-tensor (instead of the isotropic g e factor) occurs twice causing a quadratic dependence of Dw ax on Dg ax and g iso . As discussed in ref. 9 with the axiality parameter Dw ax = 27.3 Â 10 À32 m 3 . For comparison, Dw ax = 10 Â 10 À32 m 3 is the experimental value obtained by Novikov et al. 13 using superconducting quantum interference device (SQUID) magnetometry. The parameters of ZFS and g-tensors obtained by fitting the temperature-dependent SQUID data by these authors are  Table 3. The different panels (top, middle, down) represent the same data in increasing magnification of the vertical axis.  of the ZFS-tensor, as well as on g iso and Dg ax , according to eqn (3). The plot assumes cylindrical symmetry for both the g-and ZFS tensors for this quartet-state molecule. Furthermore, it assumes that the unique axes of these tensors coincide. These assumptions are valid for an idealised geometry of the studied cage complex. The same dependence of Dw ax 0 on the different magnetic parameters is shown in separate graphs in Fig. 6, where only one variable is changed at a time while the others are kept at the experimental or computational results. Dw ax 0 depends linearly on Dg ax and g iso , as well as (in the presented range of parameters) slightly non-linearly on the D parameter of the ZFS-tensor. This shows that there is no simple dependence of Dw ax 0 on any single magnetic parameter. Instead, the resulting Dw ax 0 depends principally on the combination of the tensors hSSi 8 g 8 À hSSi > g > , where the component along the cylindrical symmetry axis is denoted with 8 and the perpendicular component with >.
Even though our Dg ax is larger than that used by Novikov et al. by a factor of five, the corresponding, presently calculated Dw ax 0 is only roughly twice as large as the one obtained using their parameters. It is also worth noting that the vanishing of magnitude of Dw ax 0 (corresponding to nodes of the PCS surface) occurs when the condition hSSi 8 g 8 = hSSi > g > holds. The interplay of the value of D and temperature in forming the hSSi matrix has been elaborated in ref. 11.

Conformational averaging
The PDA is crucial in modelling biological or materials systems of real applications interest, such as the pseudocontact shifts of distant ligands in a protein, as demonstrated elsewhere. 9 In the case of the present Co(II) clathrochelate complex, we demonstrate the utility of the PDA in a dynamical averaging over the conformationally flexible chain part of the molecule. An advanced analytical and numerical method to calculate pseudocontact shifts caused by a mobile paramagnetic center has been proposed, together with the solution of the inverse problem -finding the position of the paramagnetic center from the known pseudocontact shifts. 51 For our problem it was natural to adopt a simple model and sample the conformational space of the hexadecyl chain by the Metropolis Monte Carlo method 52 over the trans/gauche conformers at each of the C-C bonds. Averages over the conformers are calculated from the PCSs based on the PDA and listed in Table 4. The same averaging is also performed using the experimental g-and ZFS tensors. Fig. 4 pNMR shielding of the clathrochelate cage + tail nuclei shown as a correlation between hyperfine shielding calculated using the full QC hyperfine couplings and using only the dipolar part of hyperfine coupling calculated using the point-dipole approximation (PDA). All common NMR nuclei ( 1 H, 11 B, 13 C), are included in the plot. The shielding constant of nitrogen nuclei directly bonded to Co remain outside of the scale of the plot, however. For the H1 nucleus, the Fermi contact and the relativistic spin-orbit correction terms (see columns 1 and 3 of Table S2, ESI †) largely cancel each other. Fig. 5 Calculated susceptibility anisotropy Dw ax 0 in the point-dipole approximation for pNMR shielding in a Co(II) clathrochelate system. Dependence on the D-parameter of zero-field splitting, isotropic g-factor, g iso , and g-tensor anisotropy, Dg ax . The reddish surface is calculated with the value Dg ax = 0.98 as obtained presently for the cage complex using the NEPVT2 calculation. The blue surface is obtained with Dg ax = 0.2 used by Novikov et al. 13 The highlighted blue point corresponds to values of D = À65 cm À1 and g iso = 2.27 in ref. 13, whereas the red point corresponds to the presently calculated values of D = À85.5 cm À1 and g iso = 2.32. The highlighted red and blue lines correspond to vanishing Dw ax 0 .
Since the simulation has been started from the conformation representing the global minimum of the potential energy, there was no need for any equilibration run in the MC simulation.
The achieved error, estimated by the method described in ref. 53, remains similar in size along the length of the chain. This means that the relative error with respect to the diminishing shielding constant increases with the distance from the paramagnetic center ( Table 4). The total length of the MC simulation was 4.5 Â 10 6 cycles as compared to altogether 3 15 E 14.3 Â 10 6 possible conformers of the 16-carbon chain.
The results are depicted in Fig. 7 and listed in Table 4. The simulated PCS results with experimental Dw ax agree very well with the experimental PCSs starting at around the position C3/H3 of the 16-carbon chain. This finding is in accordance with the static structure calculation reported in Table 3, where the error of PDA compared to the full QC pNMR shielding constant decreases to around 10% at this position in the chain. The simulation results obtained with the QC-computed v 0 tensor systematically overshoot the experimental PCS values, in complete accordance with the overestimation of the computed Dw ax as compared to the experimental value. The fact that the calculations of the ZFS tensor using state-average CASSCF with the minimal active space, even when augmented with the NEVPT2 method, tend to overestimate the D parameter, has been discussed previously. 9,54,55 However, the fact that the experimental trend is reproduced by the PDA confirms that  a Fully theoretical results using the susceptibility tensor v 0 of eqn (4). Note that the results are given as contributions to chemical shifts (opposite sign as compared to shieldings) for consistency with the experimental results. All values are given in ppm. b Results obtained from the MC simulation using the experimental Dw ax = 10 Â 10 À32 m 3 . 13 c Table S3 of ref. 13. For the calculated results, the error margins were calculated using the ''data halving'' method described by Flyvbjerg and Petersen. 53 for this application, conformational averaging by the Monte Carlo simulation over the discrete trans/gauche conformers represents the conformational space of the real system well enough -regardless of its simplicity and a lack of incorporation of any solvent effects.

Conclusion
We have investigated the point-dipole approximation (PDA) for the pseudocontact shifts of paramagnetic NMR in a Co(II) clathrochelate complex with a flexible 16-carbon chain. The PDA method was recently rendered consistent 9 with the modern formulation 10,11 of the full quantum-chemical pNMR shielding theory. 12 We have used ab initio CASSCF/NEVPT2 methods to compute the ZFS and g-tensors necessary for the calculation of, first, the susceptibility tensor using the non-symmetric formula of ref. 9 and, secondly, the resulting pseudocontact shifts. Using the example complex, we have illustrated the interdependence on the ZFS and g-tensor parameters of the resulting Dw ax and, consequently, the pseudocontact shift. The results indicate that Dw ax 0 depends principally on the difference between the products hSSi 8 g 8 and hSSi > g > . We have further confirmed the usefulness of the combination of the PDA and a simple trans/gauche conformational space averaging to simulate the dynamics of the aliphatic 16-carbon chain. The PDA makes it possible to compute a very large number of pseudocontact shifts for simulation snapshots, based on a precalculated susceptibility tensor of the metal coordination center, and the mere coordinate data of the flexible part of the model. There is great promise of the method in the prediction and analysis of pNMR shifts in complex materials and biological systems, including those with conformationally flexible parts.

Conflicts of interest
There are no conflicts to declare. Fig. 7 Pseudocontact shifts of the 16-carbon chain calculated in a Monte Carlo simulation using theoretical and experimental axiality parameter of the susceptibility tensor, Dw ax 0 and Dw ax , respectively. Experimental PCSs 13 are plotted for reference. See Table 4 for numerical values.