Electron correlation and vibrational eﬀects in predictions of paramagnetic NMR shifts †

Electronic structure calculations are fundamentally important for the interpretation of nuclear magnetic resonance (NMR) spectra from paramagnetic systems that include organometallic and inorganic compounds, catalysts, or metal-binding sites in proteins. Prediction of induced paramagnetic NMR shifts requires knowledge of electron paramagnetic resonance (EPR) parameters: the electronic g tensor, zero-field splitting D tensor, and hyperfine A tensor. The isotropic part of A , called the hyperfine coupling constant (HFCC), is one of the most troublesome properties for quantum chemistry calculations. Yet, even relatively small errors in calculations of HFCC tend to propagate into large errors in the predicted NMR shifts. The poor quality of A tensors that are currently calculated using density functional theory (DFT) constitutes a bottleneck in improving the reliability of interpretation of the NMR spectra from paramagnetic systems. In this work, electron correlation eﬀects in calculations of HFCCs with a hierarchy of ab initio methods were assessed, and the applicability of diﬀerent levels of DFT approximations and the coupled cluster singles and doubles (CCSD) method was tested. These assessments were performed for the set of selected test systems comprising an organic radical, and complexes with transition metal and rare-earth ions, for which experimental data are available. Severe deficiencies of DFT were revealed but the CCSD method was able to deliver good agreement with experimental data for all systems considered, however, at substantial computational costs. We proposed a more computationally tractable alternative, where the A was computed with the coupled cluster theory exploiting locality of electron correlation. This alternative is based on the domain-based local pair natural orbital coupled cluster singles and doubles (DLPNO-CCSD) method. In this way the robustness and reliability of the coupled cluster theory were incorporated into the modern formalism for the prediction of induced paramagnetic NMR shifts, and became applicable to systems of chemical interest. This approach was verified for the bis(cyclopentadienyl)vanadium( II ) complex (Cp 2 V; vanadocene), and the metal-binding site of the Zn 2+ - Co 2+ substituted superoxide dismutase (SOD) metalloprotein. Excellent agreement with experimental NMR shifts was achieved, which represented a substantial improvement over previous theoretical attempts. The eﬀects of vibrational corrections to orbital shielding and hyperfine tensor were evaluated and discussed within the second-order vibrational perturbation theory (VPT2) framework.


Introduction
NMR spectroscopy has evolved over the last few decades into a major instrument for probing the structure and dynamics in liquid, crystalline, and amorphous solid chemical systems.It has wide applications in chemistry, biology, and materials science.The shieldings (and therefore NMR chemical shifts) of nuclei contain accurate information about the surroundings of the corresponding atoms, including the local structure and chemical bonding.The relationships are most often complex, and as a result, theoretical NMR parameter predictions are critical for accurate signal assignment and understanding of local atomic environments in complex systems.2][3][4][5] As put by Pyykko ¨,6 ''These papers should be regarded as a whole.Together they are among the most influential ones in the quantum chemistry of the twentieth century.''Indeed, they still provide definitive answers and a comprehensive theory for closed-shell systems within the nonrelativistic limit.In addition to these, the prediction by Ramsey for the modest influence of the NMR shift on the magnetic field strength in 1970 7 has lately been proven. 8][12][13][14] Acquisition and interpretation of NMR data from paramagnetic systems, on the other hand, is still far from routine, when chemical insight into ligand atoms in close proximity to a paramagnetic transition metal (TM) or rare-earth (RE) ion is required.The total shift is measured in experimental NMR on paramagnetic systems in terms of three separate mechanisms: where the d orbital (Ramsey) chemical shift is the single shift contribution in diamagnetic systems and the other two terms relate to the unpaired electrons.The d contact term reflects the through-bond polarization and is connected with the electronnucleus hyperfine coupling constant (HFCC) between unpaired electron(s) and the nucleus.This mechanism called contact shift is operative in close vicinity to paramagnetic ions and can induce shifts as large as 10 000 ppm. 15 The pseudocontact shift (PCS, pseudocontact) is a long-range term (1/r 3 ) from the electron-nucleus dipolar coupling and magnetic anisotropy of the paramagnetic center.Paramagnetic ions can induce PCS effects for NMR active nuclei at distances of 25-40 Å. 16,17 However, at close distances to the paramagnetic center (t8 Å) both mechanisms are effective, and can shift NMR signals in either (positive or negative) direction.As a result, the theoretical prediction of induced paramagnetic NMR shifts is critical for exploring the local structure and chemical surroundings of nuclei in paramagnetic systems.Unlike the Ramsey theory for closed-shell systems, the counterpart for paramagnetic NMR shifts necessitates the evaluation of the ensemble of thermally populated electronic states, as well as the inclusion of spin-orbit (SO) and zero-field splitting (ZFS) effects.2][43] Within this formalism the induced paramagnetic NMR shift is parametrized with electron paramagnetic resonance (EPR) parameters: hyperfine A tensor, electronic g tensor, and the ZFS D tensor.
These parameters are accessible from a practical standpoint using electronic structure computations at various levels of theory.The g and D are local properties of the paramagnetic metal ion and are affected by the chemical environment only in proximity.Hence, ab initio calculations using the complete active space self-consistent field (CASSCF) 44 method followed by the second-order multi-reference perturbation theory (NEVPT2) [45][46][47][48] can be performed routinely, although within the limitation that only orbitals of the considered metal ion are included in the active space.0][51] It accounts for both static and dynamic electron correlation.On the other hand, prediction of the A tensor, and its isotropic component, HFCC, imposes substantial challenges: (i) it requires accurate quantum chemical treatment over the entire system, including all nuclei of NMR interest and their local environments, and (ii) flexible basis sets and high-quality wavefunctions to describe core-level spin polarization.Hence, accurate computations of HFCCs are very demanding.With exceptions to very small systems, multireference approaches are not viable because too many electrons and orbitals need to be included in the active space.With (ii) dynamic correlation effects should be included, and Hartree-Fock (HF) level of theory (within the unrestricted UHF formalism) provides often erratic results.3][54] DFT has been the only method chosen for routine calculations of A because it scales well with the system size and it partly accounts for electron correlation.6][57] These limitations can be related to the DFT exchange-correlation functionals and how those are affected by the Kohn-Sham self-interaction (also referred to as delocalization) error, which is critical in the prediction of HFCCs.Stepping up on the Jacob's ladder for DFT with more sophisticated descriptions of the exchange-correlation functional does not guarantee more accuracy for HFCCs, as those parameterizations are typically targeted to reproduce thermodynamic properties and not the spin polarization effects over different shells explicitly. 58,59Performance of different exchange-correlation approximations is known to be clearly system dependent. 60,61The ''parameter-free'' hybrid PBE0 functional 62 is considered the safest choice across the periodic table in this case.The Go ¨rling-Levy perturbation theory justifies using a fraction of 25% of the admixture of exact Hartree-Fock exchange (HFX) in the PBE0, 63,64 and modifying this fraction has a dramatic effect on core-level spin polarization and the predicted HFCCs.However, due to limitations of DFT approximations, the practice of making arbitrary HFX parameter adjustments to increase the agreement between calculated and experimentally observed NMR shifts in paramagnetic systems has recently evolved. 65,66Such calculations use the ''standard'' PBE0 functional and ''modified'' ones (with HFX fractions up to 50%) with the hope that the HFCCs obtained would describe the paramagnetic NMR shifts observed experimentally.Due to the nature of chemical bonding, the optimal HFX will vary even among different elements and atomic positions within the same molecule, greatly limiting the predictive power for signal assignments, etc.The DFT is essentially used as an empirical method, which is not in line with ''the right answer for the right reason'' spirit.
This journal is © the Owner Societies 2022 On the other hand, the most accurate results among the practical ab initio electronic structure theories are offered by the coupled cluster (CC) methods. 670][81][82] Especially for TM systems Kaupp 83 and Neese 84 showed that CC methods are needed for consistently accurate predictions of HFCCs.Inclusion of perturbative triples correction (CCSD(T)) in the calculations led to rather minor improvements, although descriptions of higher-order excitations as well as vibrational corrections were shown to be necessary to obtain spectroscopic accuracy. 58,85However, the applicability of the CC methods to systems of chemical interest has been limited until very recently as the computational cost scales with the size of the system for the CCSD and CCSD(T) models as N 6 and N 7 .1][92][93] Quantitative insights into challenging noncovalent interactions in endohedral fullerene complexes have for example been achieved with the DLPNO-CCSD(T) scheme in recent years. 94A complementary formulation for high-spin open-shell systems has also been developed. 95,96Moreover, within the DLPNO-CCSD scheme the use of L-equations for unrelaxed coupled cluster density provides spin densities and first-order properties. 97The extensive benchmark studies for HFCCs prediction on a large set of radicals performed by Saitow and Neese, 97 and Witwicki et al., 61 revealed that DLPNO-CCSD spin densities converge towards their canonical counterparts.[100][101][102] The main objectives of this study were to (i) demonstrate the applicability of the DLPNO-CCSD method for hyperfine tensor calculations within the current formalism for prediction of induced paramagnetic NMR shifts, and to (ii) incorporate vibrational corrections to isotropic orbital NMR shielding and hyperfine tensor.The inclusion of additional physics through electron correlation (i) and vibrational effects (ii) aimed at improving the accuracy and reliability of the existing protocol.We compared the robustness for predictions of HFCCs using the DLPNO-CCSD method against the parent CCSD model, other ab initio methods, and different DFT approximations.This analysis was performed on a variety of difficult systems including an organic radical, and complexes involving paramagnetic TM and RE ions.We evaluated the performance of the new approach for prediction of induced paramagnetic NMR shifts and the significance of vibrational corrections for the bis(cyclopentadienyl)vanadium(II) complex (vanadocene).Finally, we showed that the method worked well for the model of the metal-binding site in the Zn 2+ -Co 2+ substituted superoxide dismutase (SOD) metalloprotein.

Methods
Within the modern implementation of the Kurland-McGarvey theory by Vaara, the total NMR shielding tensor for the system with spin quantum number S 4 1/2 can be expressed as: 21,22,39 where m B , g, k, and T are the Bohr magneton, gyromagnetic ratio of the nucleus, Boltzmann constant, and absolute temperature, respectively.The orbital shielding tensor r orb corresponds to the Ramsey shielding theory for closed-shell systems, whereas the electronic g tensor and hyperfine A tensor operate with dyadic: where which represents the thermal average over the zero-field split ground state multiplet of |ni states (eigenfunctions) with E n energies (eigenvalues) of the ZFS SÁDÁS Hamiltonian.The isotropic shielding s and subsequently the isotropic chemical shift d are obtained according to where s ref and d ref are the shielding and chemical shifts of the reference compound.The hyperfine tensor A is decomposed as: where 1 is a 3 Â 3 unit matrix, A con denotes the contact isotropic HFCC, and A dip is the dipolar anisotropic coupling part of the hyperfine tensor, whereas A ZPV con and A ZPV dip denote the respective zero-point vibrational corrections (ZPV).The second contact A SO con2 , the second dipolar A SO dip2 , and the antisymmetric A SO as terms arise due to perturbational spin-orbit corrections.The electronic g tensor is decomposed into: where g e is the electron g-factor, and Dg iso and Dg ˜are the isotropic and anisotropic parts of the g-shift tensor (deviation from g e ).By evaluation of the respective components of A and g tensors the total NMR shift can be expressed as a sum of 17 physical contributions listed in Table 1.The antisymmetric s as term does not contribute to isotropic shielding, and only provides corrections to the shielding anisotropy.Note that to this end, vibrational effects were not considered in calculations of induced paramagnetic NMR shifts within the Kurland-McGarvey theory.Therefore, the seven shielding terms due to vibrational corrections (s ZPV orb , s ZPV con , s ZPV con3 , s ZPV dip , s ZPV dip3 , s ZPV ca , s ZPV pc ) have not been evaluated before in other studies.The shielding terms in Table 1 can be grouped together and in this notation, they approximate the nomenclature used in experimental studies.

Computational details
All calculations were performed with the ORCA code 103,104 in either version 4.2.1 (geometry optimizations, EPR parameters) or 5.0.1 (orbital shieldings, vibrational corrections) and employed a very tight self-consistent field (SCF) convergence tolerance of 1 Â 10 À9 E h .Evaluation of Coulomb and exchange integrals was accelerated with the RIJCOSX approximation 105 employing very tight grids (GridX8 in version 4.2.1 and DefGrid3 in 5.0.1) and either def2/J 106 or automatically generated (AutoAux) 107 Coulomb-fitting basis sets (for augmented or relativistic orbital basis sets).DFT calculations involved very tight integration grids (Grid7 NoFinalGrid and DefGrid3 in version 4.2.1 and 5.0.1, respectively).Geometry optimizations were converged to very tight thresholds (VeryTightOpt) and no symmetry constraints were imposed.DFT optimizations employed density dependent atom-pairwise dispersion correction (D4). 108Optimizations of transition metal systems involved the cc-pVTZ basis set [109][110][111] and the def2-TZVP basis set 112,113 was used for the optimization of the ytterbium complex.The global energy minima at the potential energy surfaces were confirmed by analytical Hessian calculations.The geometry of the staggered conformer of vanadocene was obtained from transition state calculation.The all electron CCSD(T)/aug-cc-pCVTZ geometry optimization of the cyanomethyl radical was performed with numerical gradients.DLPNO-CCSD(T) calculations for the vanadocene conformers employed TightPNO setting, the correlation-consistent (aug-)cc-p(wC)VXZ orbital basis sets 114,115 together with the corresponding auxiliary (aug-)cc-p(wC)VXZ/C basis sets. 116,117Cartesian coordinates of the models are provided in the ESI.† All calculations of properties involved all electrons (NoFrozenCore).The EPR/NMR properties of chemically equivalent atoms were averaged.The orbital shielding tensors r were calculated with the GIAO approach 9 at the PBE0 level using the pcSseg-2 basis set specifically developed to provide fast convergence towards the complete basis set limit for NMR shieldings. 118Methane was used as the NMR shift reference for 1 H and 13 C since chemical shifts measured in the gas phase are available for this molecule. 119he NMR parameters of the CH 4 reference (in ppm) are  65.Calculations of hyperfine tensors were performed with the core-property basis sets specifically developed for predictions of HFCCs: the aug-cc-pVTZ-J basis set developed by Sauer and coworkers 120,121 (with aug-cc-pwCVTZ/C auxiliary basis set) or the smaller EPR-II basis set developed by Barone and coworkers [122][123][124][125] (together with cc-pwCVDZ/C).For the ytterbium complex, the relativistic second-order Douglas-Kroll-Hess (DKH2) Hamiltonian was used with the SARC-DKH-TZVP 126 basis set for Yb, and the NMR-DKH 127,128 basis set for the ligand atoms.0][131][132] The sensitivity of the computed HFCCs to the completeness of the auxiliary basis sets with the DLPNO-CCSD method is shown in Table S1 in the ESI.† The CCSD and DLPNO-CCSD calculations of the hyperfine A tensors employed unrelaxed coupled cluster density, L-equations, and quasi-restricted orbitals (QROs; default for DLPNO-CCSD in ORCA).For all coupled cluster calculations the multi-reference T 1 diagnostic was r0.015 for both open-shell CCSD and L/Z_vector iterations.For the cyanomethyl radical and the TM/RE complexes, the HFC1 and HFC2 PNO truncation settings were used.The model of the CoSOD protein was divided into fragment 1 that included metal ion, nitrogen donors and Asp83 unit (atoms rendered with spheres in Fig. 2a) and fragment 2, which comprised the remaining atoms in the model (rendered with sticks in Fig. 2a).The HFC1 or HFC2 settings were invoked for the calculation, however, to reduce computational costs, the inter-and intrafragment settings for fragment 2 were reduced to LoosePNO using fragment specific settings.Spin-orbit corrections to A were estimated with at the PBE0 level (with basis sets corresponding to those used in DLPNO-CCSD calculations) using a mean-field/effective potential approach including 1-electron terms, Coulomb terms, exchange via one-center exact integrals including the spin-other-orbit interaction, and local DFT correlation, as implemented in ORCA.Calculations of A with the double hybrid DFT involved relaxed MP2 density.For calculations of vibrational corrections to r and A the second-order vibrational perturbation theory (VPT2) was employed, as implemented in ORCA.Hessian calculations for the anharmonic VPT2 force fields involved the PBE0-D4/cc-pVTZ level starting from the respective optimized geometries, and property calculations involved the PBE0/pcSseg-2 level for r and PBE0/aug-cc-pVTZ-J for A. For the numerical calculations of Hessian and property derivatives a step size of 0.05 was used for both the anharmonic and property displacements.The g and D tensors were calculated with a state-averaged, complete active space self-consistent field (CASSCF) method distributing n d-shell electrons (n = 3 for V and n = 7 for Co) over five 3d orbitals of a metal ion; CAS(n,5), followed by invoking an N-electron valence state second-order perturbation theory (NEVPT2).CASSCF calculation for vanadocene included 10 quartet and 15 doublet roots, whereas for the CoSOD protein model 10 quartet and 40 doublet roots were evaluated.In the CASSCF/NEVPT2 calculations the aug-cc-pwCVTZ-DK 133 basis set was used for V and Co, while the aug-cc-pVTZ-DK 134 basis set was used for the ligand atoms in vanadocene, and the cc-pVDZ-DK for the ligand atoms in the CoSOD model.Auxiliary basis sets generated with AutoAux were employed.Both scalar-relativistic (via the DKH2 approach, picture-changed operators, finite nucleus model) and spin-orbit effects were included. 44,50Induced paramagnetic NMR shifts were evaluated assuming temperatures of 298 and 280 K for the vanadocene and protein model, respectively, in accordance to the experimental conditions.Computations were performed on a cluster node equipped with the two Intel s Xeon s Gold 6126 CPUs (2.6 GHz; 12-core) and 256 GB of RAM.

Prediction of hyperfine coupling: the challenge
In this section the impact of electron correlation on the prediction of HFCCs was evaluated by comparing results at the HF, MP2 (different variants), DLPNO-CCSD, and canonical CCSD levels of theory with those obtained at four different levels of DFT approximations.There were three types of systems considered (with available experimental data).First, the HFCC for the 14 N in cyanomethyl radical NRC-CH 2 was calculated.The nitrogen atom is in general challenging for calculations of EPR/NMR parameters.For example, the CCSD(T) level of theory was needed to provide quantitative agreement with experimental estimates for the HFCC of the nitrogen atom in its ground state 80 or NMR shielding in the NRN molecule. 13evious DFT studies on nitrogen containing organic radicals have revealed that the results depend on the particular choice of the exchange-correlation functional and the basis set, which often have led to fortuitous cancellation of errors. .Although the spin-orbit contributions to the hyperfine tensor cannot be currently assessed at the DLPNO-CCSD level, these are expected to be small for high-spin d 3 and d 5 configurations of the V 2+ , Cr 3+ , and Mn 2+ ions. 84,97As the last, yet the most challenging system, a RE metal ion complex was considered: [Yb(H 2 O) 8 ] 3+ .
The results for the prediction of hyperfine and quadrupolar coupling constants for the 14 N nucleus of the cyanomethyl radical are presented in Table 2.The 14 N hyperfine coupling constant (A con ) in the cyanomethyl radical calculated with the simplest DFT method for open-shell systems, local spin density approximation (LSDA), 141 was severely underestimated compared to experimental values, but the predicted 14 N quadrupolar coupling constant (w) was quite close to the experiment.As one could assume, predicting the spin density at the nucleus is substantially more challenging than predicting the electric field gradients.Incorporation of electron density gradient corrections to the LSDA in the generalized gradient approximation (GGA) PBE functional 142 significantly improved on the predictions of both A con and w.Nonetheless, the calculated A con was still significantly underestimated.DFT underestimate the core-level spin polarization, whereas the Hartree-Fock method overestimates it.Incorporating 25% of the exact Hartree-Fock exchange energy component in the hybrid-GGA PBE0 exchange-correlation functional 62 improved prediction of A con significantly.However, further increase of the HF exchange admixture to 50% (PBE0 (50% HFX) ) severely deteriorated the predictions of both A con and w when compared to the original PBE0 formulation.The last step on the Jacob's ladder for DFT approximations consists of double hybrid, ''perturbatively corrected'' exchange-correlation functionals that incorporate an exact HF exchange component and a perturbative second-order correlation contribution. 145The double hybrid DSD-PBEP86 functional 146 that incorporates HF exchange, PBE and P86 exchange and correlation potentials, as well as spin-component-scaled (SCS) MP2 mixing, constitutes one of the most accurate DFT methods for general chemistry applications. 1479][150][151][152] However, as shown in Table 2, for the case of cyanomethyl radicals the DSD-PBEP86 functional suffered from instabilities that resulted in erratic spin density on the nitrogen atom.Special care has to be taken when applying double hybrid functionals to open-shell systems. 60This also indicated that in general, no systematic improvements for prediction of the hyperfine coupling constants were feasible within the current DFT framework. 61The Hartree-Fock and MP2/SCS-MP2 methods provided essentially erratic results, and values obtained with MP2 and its SCS-MP2 variant were similar.Orbital optimization (OO-MP2/OO-SCS-MP2) remarkably improved the results by removing spin contamination of the reference wavefunction, however, the predicted HFCCs were still far from experimental ones.The robustness of the CCSD level of theory was clearly necessary to provide good agreement with experimental data for the NRC-CH 2 .We noted however that the CCSD seemed to slightly overestimate the 14 N hyperfine coupling in cyanomethyl radical.This was attributed to the lack of perturbative triple correction (CCSD(T)) that was shown to slightly lower the HFCCs predictions compared to CCSD. 82Before we discuss the predictions with the DLPNO-CCSD method, we add a remark.To facilitate the accuracy and performance control of the DLPNO-CCSD calculations, five levels of predefined PNO truncation settings have been developed: LoosePNO, NormalPNO, TightPNO, HFC1, and HFC2.
These settings provide convergence towards the canonical limit at increasing computational cost, 90,97 by invoking refined thresholds for the electron pair correlation energies evaluation.Essential pairs are selected for accurate coupled cluster treatment, and those pairs with a less critical contribution to the total correlation energy are subjected to a much more computationally efficient second-order perturbation theory.The LoosePNO, NormalPNO, and TightPNO thresholds were designed for general chemistry applications, i.e. evaluations of total energies or intermolecular interactions.The tighter HFC1 and HFC2 settings were specifically developed for calculations of hyperfine coupling constants.As shown in Table 2, for predictions of the 14 N quadrupolar coupling constants both the HFC1 and HFC2 settings provided converged results of CCSD quality.They were clearly better than those from hybrid DFT.However, the convergence for the prediction of the 14 N hyperfine coupling constants was more difficult to achieve.With the HFC1 setup full convergence was not achieved and the more conservative (and expensive) HFC2 setting was needed to provide results close to the CCSD reference.This in turn indicated that electron correlation effects in NRC-CH 2 are substantial.The sum of first-order components of 51 V, 53 Cr, and 55 Mn hyperfine tensors (A xx/yy/zz con+dip ) and the isotropic components of the 17 2+ , and [Yb(H 2 O) 8 ] 3+ complexes are shown in Table 3.For predictions related to 51 V, 53 Cr, and 55 Mn similar patterns were observed as for the cyanomethyl radical.Calculations with the two lowest levels of DFT approximations, LSDA and GGA, severely underestimated the magnitudes of hyperfine couplings.Although the hybrid DFT method provided considerable improvements, the predicted values were still significantly underestimated.In contrast to the cyanomethyl radical case, where the increase of Hartree-Fock exchange admixture to 50% deteriorated the results severely, the predictions of the hyperfine couplings were improved for 51 V, 53 Cr, and 55 Mn.However, the hyperfine couplings obtained with the PBE0 and PBE0 (50% HFX) approximations were lower than the experimental values.The double hybrid DSD-PBEP86 approximation did offer improvements, and the predictions with this functional were generally closer to the experimental values than those with hybrid DFT.Considering the ab initio methods, HF was expected to overestimate hyperfine couplings, although this effect was much smaller for the TM complexes than for the cyanomethyl radical.The predictions with MP2 and orbital-optimized MP2 were quite close to the experimental values.The values were notably closer than those calculated with DFT, with the exception of 55 Mn.Finally, both the CCSD and DLPNO-CCSD predictions were in good agreement with the experimental results for 51 V and 53 Cr, with somewhat larger deviations for 55 Mn, which might be attributed to the lack of triples correction.Although the results for the radical and metal hyperfine couplings provide an assessment of the applicability of the respective computational methods, the accuracy for the prediction of EPR and NMR properties of the ligand atoms in the TM and RE complexes is the most important, from the point of view of practical calculations for the prediction of induced paramagnetic shifts.Let us focus first on the 17 O isotropic hyperfine coupling constant (A con ) in [Mn(H 2 O) 6 ] 2+ , for which experimental result is available.Both the CCSD and DLPNO-CCSD methods provided excellent predictions that were almost exact matches with experimental data.The MP2 and OO-MP2 results were also close, whereas the HF values deviated.When it comes to DFT, predictions with the LSDA and PBE approximations resulted in severe errors.Also, the calculations with the PBE0 and PBE0 (50% HFX) functionals were not able to approach the experimental result.In fact, considering the results for ligand atoms for all systems in Table 3, the predictions with the PBE0 and PBE0 (50% HFX) functionals never encompassed experimental results or the CCSD predictions.Furthermore, the effect of increasing the HF exchange admixture was not always beneficial.It deteriorated the predictions for [V(H 2 O) 6 ] 2+ and [Cr(H 2 O) 6 ] 3+ as compared to the genuine PBE0 formulation.The predictions as presented in Table 3 suggested that it is not possible to assess the intrinsic DFT errors in the prediction of HFCCs by performing computations with PBE0 and PBE0 (50% HFX) approximations.Predictions of properties with hybrid DFT for the ligand atoms were inconclusive as shown in Table 3, and for the RE complex, it was not even possible to establish the sign of the coupling constant.In contrast, the predictions performed with the double hybrid DSD-PBEP86 functional consistently improved on the hybrid DFT, although the deviations from experimental values or reference CCSD calculations were still significant.The results for [Yb(H 2 O) 8 ] 3+ displayed an even more complicated pattern.The predictions performed by LSDA, PBE, and PBE0 provided erratic values.For this complex, the DFTbased calculations benefited from high HF exchange admixtures.PBE0 (50% HFX) and DSD-PBEP86 (which incorporates 72% HFX) improved on the predictions.With HF the calculated HFCCs were overestimated, but with MP2 improvements were achieved.Predictions with DLPNO-CCSD using HFC1 and HFC2 thresholds revealed convergence towards the experimental value, which was encouraging in the perspective of future studies involving RE systems.The OO-MP2 calculation did not converge even with significantly increased number of iterations, whereas the computational cost of each canonical CCSD iteration was so high that it was not expected to converge within a reasonable time frame.
To summarize, results, as presented in Tables 2 and 3, confirmed outstanding reliability of the CCSD method for the prediction of HFCCs for all systems considered.Most importantly, DLPNO-CCSD spin densities were shown to converge towards the canonical CCSD counterparts.Given the substantially reduced computational cost of the DLPNO-CCSD method compared to CCSD, predictions of coupled cluster quality can be achieved for large systems of chemical interest.Notably, although the openshell DLPNO-CCSD method in its current implementation exhibits very favorable scaling (linear to quadratic) when considering the iterative part, integral transformations on larger systems add substantial prefactors.In fact, the memory and disk space requirements of these integral transformations limit the feasibility of the method on large systems rather than long computation times.Yet, the open-shell formulation of the DLPNO-CCSD method is twice more demanding compared to its closed-shell counterpart.Therefore, although the HFC2 setup produced excellent results, very close to canonical CCSD, it is not a practical choice to employ such tight thresholds for large systems, which was indicated by Saitow and Neese.The HFC1 setup produces HFCCs ''of desirable accuracy while maximizing the balance between acceptable computational cost and high accuracy.'' 97Thus, in the rest of this current study, the HFC1 setup is assumed when discussing the DLPNO-CCSD results unless otherwise specified.

From news to everyday use: prediction of NMR shifts in vanadocene
The bis(cyclopentadienyl)vanadium(II) complex (Cp 2 V; vanadocene) was chosen as a test case for an improved method for predictions of NMR shifts.High-quality experimental and theoretical EPR/NMR data are available for comparison, and the spin-orbit contributions to hyperfine tensor (which cannot be examined at the coupled cluster level) are expected to be modest for the high-spin d 3 configuration of the V 2+ ion in the complex.Vanadocene can have two distinct conformations, eclipsed and staggered, see Fig. 1(a) and (b).The gas phase electron diffraction analysis by Gard et al. 153 found an average distance between the vanadium and cyclopentadienyl rings r(V-C 5 H 5 ) = 1.928AE 0.006 b Å and an eclipsed conformation, however, the staggered configuration could not be ruled out completely.The single crystal X-ray diffraction study by Rogers et al. 154 indicated a distance r(V-C 5 H 5 ) = 1.92 AE 0.02 Å and although the vanadium atom resided on a crystallographic center of inversion, the cyclopentadienyl rings were dynamically disordered at room temperature.This was further corroborated by a temperature-dependent X-ray study by Antipin et al. 155 Optimization of the molecular geometry of vanadocene was performed at the PBE0-D4/cc-pVTZ level of theory without symmetry constrains and resulted in the eclipsed conformer with a distance r(V-C 5 H 5 ) = 1.929Å, which agreed excellently with the experimental estimates.The energy difference between the two conformers was found to be small, which we have also corroborated by high-level DLPNO-CCSD(T) calculations, see calculation of the energy difference between eclipsed and staggered conformers (E eclipsed À E staggered ) was À0.17 kcal mol À1 .These calculations were performed using the all electron DLPNO-CCSD(T) method with the aug-cc-pwCVTZ basis set, incorporating diffuse functions and core-valence correlation.The small difference in turn indicated that the energy barrier for the rotational diffusion of the rings would be very small too.
Experimental and calculated parameters are shown in Table 5 for both the eclipsed and staggered conformers of vanadocene.7][158][159][160][161][162] The calculated EPR properties of both conformers were nearly identical.The isotropic components of the hyperfine tensor (A con ) were positive for 1 H and negative for 13 C, which was also reflected by the corresponding positive and negative regions of the spin densities plotted in Fig. 1(c)-(e).Vibrational corrections to the hyperfine coupling constant (A ZPV con ) were minor, but they were more than twice as large as the corresponding spin-orbit corrections (A SO con2 ).Both the vibration and spin-orbit corrections were more important for 13 C than for 1 H.The calculated ZFS parameters D and E/D were in very good agreement with the experiment.Note that within the CASSCF/NEVPT2 approach, only the spin-orbit coupling (SOC) contribution to the ZFS tensor was calculated, whereas with the SORCI method both the SOC and spin-spin coupling (SSC) contributions can be obtained.However, for TM complexes the SSC contributions are typically small, so the CASSCF/NEVPT2 approach is a good approximation, as was the case here.
When considering computational expenses, the DLPNO-CCSD/aug-cc-pVTZ-J calculation for A, and the DKH2-CASSCF/ NEVPT2/aug-cc-pwCVTZ-DK/aug-cc-pVTZ-DK calculation of g and D, each took less than 2 days on our hardware.Even though the coupled-perturbed equations have to be solved for all ligand nuclei in the system, the PBE0/aug-cc-pVTZ-J calculation of the spin-orbit corrections to A required only a few hours.The computations of vibrational corrections to A, with the VPT2 PBE0-D4/cc-pVTZ//PBE0/aug-cc-pVTZ-J method took 10 days, since it required 114 Hessian and property calculations.On the other hand, the computed Hessian and displaced geometries data could be used again for the VPT2 calculation of vibrational corrections to the orbital shielding (this time with the PBE0/pcSseg-2 property evaluation), avoiding repetition of the most time consuming part of the calculation.We also note Table 4 Energy difference between the eclipsed and staggered conformers (E eclipsed À E staggered ; kcal mol À1 ) of vanadocene

Method
Basis set aug-cc-pwCVTZ À0.17 b All electron, no frozen core.that the VPT2 approach is much more convenient than setting up molecular dynamics DFT (MD-DFT) simulations and performing calculations of r orb and A for tens or hundreds of system evolution trajectory frames.After collection and evaluation of the obtained A, g, D, and r orb data, the resulting 16 shielding terms contributing to the isotropic 1 H and 13 C NMR shifts in the eclipsed conformer of vanadocene are shown in Table 6.The largest contributions originated from orbital (s orb ) and contact (s con ) terms together with the respective vibrational corrections (s ZPV orb and s ZPV con ).In the case of 13 C, terms originating from the spin-orbit corrections (s con2 ) and g-shift tensor (s con3 ) made important contributions.The calculated total 1 H and 13 C NMR shifts were very close to the experimental values (within 2% of error). 1 H and 13 C NMR shifts for vanadocene have been computed by Kaupp, 36 Vaara, 20,24 and Autschbach. 31In those studies, the hyperfine tensors were obtained at the DFT level using different implementations (including four-component relativistic approach) 24 and several exchange-correlation approximations.The range of previously reported NMR shifts is shown in the last row of Table 6.By including additional physics through the inclusion of electron correlation (DLPNO-CCSD) and vibrational effects (VPT2) we were able to significantly improve the predictions of 1 H and 13 C NMR shifts in vanadocene.To illustrate better the relative importance of NMR shielding mechanisms in vanadocene, the shielding terms were grouped according to the classification in eqn (7) in Table 7.The contact mechanism had by far the largest contribution to 1 H and 13 C NMR shifts.Due to the cancellation of s dip and s pc terms the contribution from the pseudocontact mechanism was minimal (see Table 6).To provide an assessment on the necessity to include vibrational and spin-orbit corrections to the predictions of the hyperfine tensor, the summed contributions to induced paramagnetic NMR shielding in vanadocene are shown in Table 8.These shielding contributions originated from the firstorder, ZPV, and SO components of A. Components beyond the first-order were significantly more important for 13 C than for 1 H.Although the first-order components of A already accounted for 94% of the calculated induced paramagnetic 13 C NMR shielding for vanadocene, the remaining 6% could easily be imagined to be crucial in situations when many resonances with small shift differences were to be assigned in the spectrum.These contributions relate to the ZPV and SO components of A.
We expect that the remaining shortcomings and errors in the proposed protocol are dominated by the lack of corrections for the triples excitations (T) and spin-orbit effects within the coupled cluster calculations of A (these are not implemented; the SO part has to be done with DFT), as well as by deficiencies in the DFT level prediction of the orbital part of the shielding.Although the excellent agreement with experimental data for vanadocene indicated that these effects are minor, however, it is not guaranteed to be the case for other systems.The current formalism is limited to the ground-state multiplet, which could cause issues with molecular systems having low-lying excited states.Considering the calculation of r orb , unfortunately no improvements beyond the hydrid DFT level are currently feasible.This is because MP2 and double hybrid DFT (e.g.DSD-PBEP86) are unreliable when applied to calculations of r orb in transition metal complexes, as illustrated in Table S2 (ESI †) for the closed-shell ferrocene molecule.In this molecule, the orbital mechanism is the sole contributor to the 1 H and 13 C NMR shielding.

Application to Co(II) biding site in SOD metalloprotein
Usign Co 2+ as a spectroscopic surrogate for Zn 2+ is established in metallobiochemistry. 164 The high-spin d 7 Co 2+ ion usually adopts a coordination geometry similar to the native Zn 2+ .The Co 2+ allows EPR, paramagnetic NMR, and ENDOR spectroscopies to be used for probing local environments in enzymes and proteins.This approach was employed in the study by Bertarello et al. 65 where the geometry of metal-binding site in the thermostable mutant of human superoxide dismutase 1 (SOD) was investigated with magic-angle spinning (MAS) NMR.These models resulted from ten potential binding sites of Zn 2+ in ten protein chains derived from the X-ray crystal structure of SOD. 166After Zn 2+ -Co 2+ substitution, the geometries of the models were optimized with a few constraints to keep the overall fold of the  9 were prepared.The ranges of shifts obtained by Bertarello et al. with the PBE0/ PBE0 (50% HFX) approach for A were very wide, with a span of around 200 ppm for each model.In contrast, when the DLPNO-CCSD method was employed for the prediction of A, the accuracy and consistency of NMR shift prediction improved remarkably.Considering models E and F, the PBE0/ PBE0 (50% HFX) calculations provided rather inconclusive chemical shifts ranging from 299 to 531 ppm, whereas with the DLPNO-CCSD approach the values of 337 and 331 ppm were obtained for model E and F, respectively; both very close to the experimental result of 350 ppm.These results suggested that model E might be closer to the real coordination sphere of Co 2+ ion in CoSOD.We note however that vibrational corrections to r orb and A could not be evaluated because models derived from the X-ray structure of protein and optimized with constraints did not allow for Hessian evaluations.We checked the DLPNO-CCSD convergence for Model E, and d for 13 C b changed only marginally from 337.29 to 337.23 ppm when going from HFC1 to HFC2 thresholds and keeping all other parameters same.This indicated that results for 13 C b in the Asp83 residue in our calculations were close to the CCSD limit.Larger changes of HFCCs were observed for nuclei directly bonded to the Co 2+ ion, however, signals from those nuclei were not observed in the experiments.This presence of a ''NMR blind sphere'' around the paramagnetic Co 2+ ion in CoSOD is indicted in Fig. 2c-e where spin densities rendered with the isodensity of 0.002 e bohr À3 are observed for the donor nitrogen and oxygen atoms.Since the magnetic and EPR properties of CoSOD have been measured, 167,168 comparisons of the calculated g and D with experimental results could provide further insights.The calculated isotropic g-values of 2.228 and 2.225 for the models E and F were similar, but that from model E was slightly closer to the experimental result of 2.24.However, the calculated D-values of 5.62 and 4.25 cm À1 for model E and F, respectively, were significantly lower than the experimental estimate of 10.8 cm À1 .In general, the ZFS of tetrahedral, high-spin Co 2+ complexes is strongly dominated by the SOC interaction, 161,169 and errors in the DKH2-CASSCF/NEVPT2 approach are not expected to be significant.Moreover, reported experimental D-values in four-coordinated Co 2+ complexes involving nitrogen donors were in the range from 2.4 to 7.5 (cm À1 ), 170,171 which corresponded very well with our theoretical predictions for CoSOD.In the work by Bertarello et al. it was indicated that the experimental estimate of D obtained by fitting of temperature dependent magnetic susceptibility data with susceptibility equations derived by using isotropic g-factors and axial parameters D might not be accurate.Nonetheless, all three: the NMR shift of 13 C b (Asp83), isotropic g-value, and D parameter indicated that the coordination geometry of model E was the most accurate.

Conclusions
We improved on the implementation of the Kurland-McGarvey theory for prediction of NMR shifts in paramagnetic systems by using the DLPNO-CCSD method to predict the critical hyperfine coupling tensor A, and by including vibrational corrections to the orbital shielding and the A tensor.The description of the additional physics involved electron correlation and vibrational effects and improved the accuracy and reliability of the calculations over existing protocols.The calculated NMR shifts for vanadocene were within 2% of the experimental result, indicating the importance of vibrational corrections, especially for 13 C.For the much larger and challenging model of the CoSOD metalloprotein, the predicted values were within 4% of the experiment.For these calculations, the vibrational corrections could not be calculated, and a smaller basis set had to be used.The advancement over earlier studies is a significant step towards reliable interpretation of NMR spectra from paramagnetic molecular systems.In terms of application of an alternative approach using hybrid DFT methods, the calculations for complexes with transition metal and rare-earth ions indicated that the hybrid DFT approach is not able to reliably predict the HFCCs for the ligand atoms.Manipulation of the Hartree-Fock exchange ratios did not lead to conclusive results, and we advise against this.On the other hand, double hybrid DFT demonstrated systematic improvements, but should be used with caution due to the risk of instabilities associated with deficiencies of the HF and MP2 components in open-shell systems.Finally, predictions of HFCCs with the DLPNO-CCSD method were shown to converge towards the canonical CCSD limit at much lower computational cost.This method exhibited outstanding reliability and consistency for all systems considered in this study.This confirmed that DLPNO-CCSD is the method of choice for calculations of A tensors within the current

Fig. 1
Fig. 1 Eclipsed (a) and staggered (b) conformers of vanadocene.Spin densities calculated at the DLPNO-CCSD/aug-cc-pVTZ-J level of theory shown with an isodensity of 0.002 (panel c), 0.0002 (d), and 0.00002 (e) e bohr À3 ; orange and cyan denote the positive and negative regions, respectively.
By probing induced paramagnetic NMR shifts of ligand nuclei, information on the coordination sphere of Co 2+ ion in Cosubstituted superoxide dismutase (CoSOD) could be obtained, see Fig.2a.Progress in MAS NMR hardware and radiofrequency irradiation (pulse) schemes together with a very fast MAS rate (100 kHz) enabled improved excitation and detection conditions compared to earlier NMR studies of the same system.165It was in this way possible to observe resonances from nuclei relatively close to the paramagnetic Co 2+ ion.To assign the signals affected by the induced paramagnetic NMR shifts, Bertarello et al. used quantum chemistry methods in a set of 10 models (86 atoms each, labeled from A to J).

Fig. 2
Fig. 2 Model of the binding site of Co 2+ ion with His63, His71, His80, and Asp83 residues in CoSOD metalloprotein (panel a).Calculated NMR shifts for 13 C b (Asp83) (panel b); DFT results for models A-J from ref. 65, predictions with DLPNO-CCSD for models E and F from this work.Spin densities

Table 1
NMR shielding terms and the respective tensorial ranks for a system with spin quantum number S 4 1/2 Dg iso-A con hS e S t i a Dg ea hS a S t i ab Dg ea A dip;ZPV bt hS a S b i s ZPVThis journal is © the Owner Societies 2022 58,135-138 Second, we considered challenging calculations of both metal and ligand-related HFCCs in transition metal complexes: [V(H 2 O) 6 ] 2+ , [Cr(H 2 O) 6 ] 3+ , and [Mn(H 2 O) 6 ] 2+

Table 2
14N hyperfine (A con ) and quadrupolar (w) coupling constants (MHz) in cyanomethyl radical NRC-CH 2 calculated with different methods using the aug-cc-pVTZ-J basis set

Table 5
Molecular geometries and EPR properties of vanadocene

Table 7
13ntributions (in ppm) to isotropic 1 H and13C NMR shielding for the eclipsed conformer of vanadocene grouped according to the classification in eqn(7)

Table 9
Calculated 13 C NMR shifts (ppm) of b-carbon in the Asp83 residue of CoSOD