Modelling the Effect of Zero-Field Splitting on the 1 H , 13 C and 29 Si Chemical Shifts of Lanthanide and Actinide Compounds

The prediction of paramagnetic NMR (pNMR) chemical shifts in molecules containing heavy atoms presents a significant challenge to computational quantum chemistry. The importance of meeting this challenge lies in the central role that NMR plays in the structural characterisation of chemical systems. Hence there is a need for reliable assignment and prediction of chemical shifts. In a previous study [Trends in Physical Chemistry, 17, 25–57, (2017)] we looked at the computation of pNMR chemical shifts in lanthanide and actinide complexes using a spin Hamiltonian approach. In that study we were principally concerned with molecules with S = 1/2 ground states. In the present work we extend that study by looking at the effect of zero field splitting (ZFS) for six complexes with S = 3/2 ground states. It is shown that the inclusion of ZFS can produce substantial shifts in the predicted chemical shifts. The computations presented are typically sufficient to enable assignment of experimental spectra. However for one case, in which the peaks are closely clustered, the inclusion of ZFS re-orders the chemical shifts making assignment quite difficult. We also observe, and echo, the previously reported importance of including the paramagnetic spin-orbit hyperfine interaction for 13C and 29Si atoms, when these are directly bound to a heavy element and thus subject to heavy-atom-light-atom effects. The necessary computations are very demanding, and more work is needed to find theoretical and computational approaches that simplify the evaluation of this term. We discuss the computation of each term required in the spin Hamiltonian. The systems we study in this work are restricted to a single heavy atom ion (one Nd(III) and five U(III) complexes), but typify some of the computational complexity encountered in lanthanide and actinide containing molecules.


Introduction
Interest in lanthanide and actinide chemistry continues to expand as new applications are found and existing areas of research are developed.Examples include novel synthetic techniques, such as the development of reusable reagents in the synthesis of cyclic oxycarbon compounds [1], and the investigation into the nature of actinide containing single molecule magnets [2].
Radioactivity is the most exploited property of the actinides, as they potentially offer an efficient and low carbon alternative to fossil fuels for large scale energy generation.At the opposite end of the scale, nuclear powered batteries are used for powering small devices such as heart pacemakers [3].Because of their industrial use, the reprocessing of nuclear waste to extract fissile material and reduce waste volumes is of paramount importance.This leads directly to a demand for better extractants, which in turn needs a better understanding of metal-ligand selectivity for actinide over lanthanide ions [4].This is an area that benefits from improved and reliable computational modelling techniques.
Nuclear magnetic resonance spectroscopy (NMR) is a well-established and much used analytical tool for characterising molecules and providing insight into their constituent chemical bonds.Many experimental techniques exist for obtaining NMR spectra in diverse chemical and physical situations.The analysis and interpretation of complex spectra for most diamagnetic molecules, that lack heavy nuclei, are routine.Complementary to this experimental work, the development of computational tools for predicting NMR chemical shifts has also become common [5].
NMR studies of actinide species provide significant challenges to computational techniques.Firstly, for light atoms the effects of relativity are small, and the representation of electron correlation is significantly more important.However the relativistic correction to the energy of an atomic system approximately scales as (Zc −1 ) 2 [6].For heavy elements relativistic effects have a profound influence on the ordering of energy levels within their atoms and molecules, and a concomitant effect on the molecular properties of such systems.The computational treatment of relativistic effects begins with scalar relativistic corrections and these are often sufficient for the prediction of molecular geometries.For paramagnetic molecules the presence of the unpaired electronic spin requires that the coupling between orbital and spin angular momenta are treated reliably.This can also have a significant effect on the splitting of electronic energy levels.Equally, spin-dependent properties such as the Zeeman interaction, the hyperfine interaction and the spin-spin interaction will depend on the theoretical method used for their evaluation and require calibration against reliable experimental data.These considerations provide significant challenges for the theoretical formulation of paramagnetic NMR (pNMR) chemical shifts, as well as for their computational realisation.Calculating pNMR shifts of paramagnetic molecules remains a challenging problem that is still under much investigation [7,8].
Experimental spectra are normally reported as δ shifts.The δ shift is obtained with respect to a reference chemical shift, σ ref .It is conventional to separate the experimental paramagnetic chemical shift σ exp into three terms: the orbital contribution σ orb , which is approximated as the shift that would be observed if the species was diamagnetic, and the contact σ c and pseudocontact σ pc shifts arising from the unpaired electronic spin: For simple systems, the contact term relates to the interaction of the unpaired electron spin density and the nuclear magnetic moment at the nucleus.The contact interaction is isotropic.The pseudocontact shift is the residual long-range interaction of the spin-dipoles and is anisotropic [9].Qualitatively the contact shift is traditionally interpreted as arising from the transfer of spin density through bonds and may be associated with covalent interactions.The pseudocontact shift is interpreted as acting through space and is a dipole-dipole type interaction [10].The theoretical justification for the decomposition into these terms weakens when relativistic treatments are included [11].
The theoretical framework of pNMR shifts can be formulated in a number of ways.The full pNMR shielding expression in a sum over states (SOS) formalism was derived in references [12,13].
In reference [14] multireference wave function calculations of pNMR shifts were performed, for the first time, directly employing the SOS formalism without reference to the spin Hamilitonian formalism.This study was extended to carbonate complexes of U, Np and Pu [15].In references [14,15] the wave function code used was implemented as part of the OpenMolcas code but has not yet been released.Accordingly, in this work we shall adopt the spin Hamiltonian approach, which represents each of the three paramagnetic operators with a coupling matrix and is detailed in the next section.This approach has the advantage that the necessary components are readily available in many quantum chemistry packages; however, the main disadvantage is that the computational framework that is used to calculate each component has a significant effect on the results obtained.To achieve reliable results it is necessary to assess each component, and calibrate the most effective approach against an often-limited set of experimental values for similar chemical species.Detailed overviews of the evaluation of spin Hamiltonian parameters in the context of modern computational chemistry are available from Vaara [16], Bolvin [17] and Autschbach [18].We shall mostly be concerned with studying the effect of zero field splitting (ZFS) on the predicted 1 H, 13 C, 29 Si shifts of six compounds.We will illustrate situations where the inclusion of ZFS causes a reordering in the predicted 1 H shifts, and where the inclusion of spin-orbit effects on the hyperfine term for the 13 C shifts of nuclei involved in a η-coordinated C-U bond dominate the hyperfine interaction.With the exception of the first example, our results are accurate enough to enable qualitative assignment of observed shifts.

The Spin Hamiltonian Approach to Paramagnetic NMR
Applying a magnetic field B to an isolated atom with nuclear spin I nuc results in each energy level splitting by 2γ nuc |B|, where γ nuc is the magnetogyric ratio of the nucleus under investigation.In any electronic system, the circulation of electrons induced by the applied magnetic field generates an opposing magnetic field.This effect is represented by σ orb , the chemical shift tensor.For a diamagnetic molecule, the Hamiltonian corresponding to the energy shift, in atomic units, is given by However for a system with unpaired electrons, there is a permanent magnetic moment associated with the electronic spin.It is energetically favourable for this to align with an applied magnetic field and the resulting (thermally averaged) dipole must be accounted for.The corresponding paramagnetic NMR Hamiltonian is given by [19][20][21][22]: The Zeeman coupling matrix, g, represents the interaction of the electronic spin-dipole with the applied magnetic field.An effective parameter, referred to as the pseudospin S, is introduced.In the absence of spin-orbit coupling S = S, which is an often-used approximation.For molecules where the spin-orbit interaction is strong, the effective degeneracy of the molecule is reduced so a smaller effective spin S is chosen based on the splitting pattern of the electronic states [17].The interaction between the nuclear and electronic spin-dipoles is represented by the hyperfine coupling matrix, A. In a system with more than one unpaired electron, the interaction of the electron spin-dipoles is represented by the ZFS or D tensor.

The Orbital Shielding Tensor, σ orb
For a diamagnetic molecule the orbital, or Ramsey, shielding matrix [23] is the only contribution to the observed chemical shift and as such, its isotropic value is readily extractable from NMR spectra.For a paramagnetic molecule it is not possible to directly extract this term and as a result, the reliability of any computational method must be assessed on diamagnetic analogues of the species under investigation.Because magnetic fields are defined by the action of an infinitesimal rotation of a vector potential and this operation is independent of the origin of the co-ordinate system used, the calculation of magnetic properties has an artificial dependence on the choice of origin, known as gauge dependence.The gauge dependence vanishes in the limit of the exact wave function.It has been demonstrated that for calculations targeting a single nucleus, the dependence is minimised by setting the origin at the nucleus.For molecules in general, choosing the centre of electronic charge as the origin minimises gauge dependence [24], but the most reliable approach is to employ gauge-including atomic orbitals (GIAO) [25,26], which remove the gauge dependence.
For actinide compounds, modelling relativistic effects is essential and many frameworks can be used.The relativistic Dirac equation can be satisfied by a four component wave function, partitioned into pairs of electronic and positronic spinors.The presence of the positronic solution hinders a variational approach, since simply minimising the energy is no longer possible due to an unbounded and infinite set of negative positronic states.For four component methods, this issue must be addressed, e.g., by kinetic balance [27].However, since the positronic component is much smaller, efficient two component techniques have been developed to eliminate it.These techniques produce a spin-free (scalar) and a spin-dependent component.Observables in the Dirac framework are implicitly dependent on all four components, and the corresponding two component operator erroneously does not involve the full solution.This effect is called the picture change error [18] and can be addressed using a transformation of the operator equivalent to that used to transform the four component wave function to a two component wave function.
Scalar relativistic approaches such as the zero-order regular approximation (ZORA) [28], and the Douglas-Kroll Hamiltonian (DKH) [29] are implemented in a number of quantum chemistry packages.ZORA only holds for basis functions that are independent of the perturbation, which theoretically precludes using GIAOs, although in practice the error introduced is small for large basis sets [30].

The Zeeman Coupling Matrix, g
The Zeeman coupling matrix, g, more commonly referred to as the "g tensor", (since it relates two quantities defined in space and spin coordinate frames respectively it is not a true tensor.)represents the effect of the magnetic field on the average electronic spin-dipole.With the addition of spin-orbit coupling, the electronic spin of a molecule is no longer well defined and the pseudospin parameter is introduced to describe the observed splitting of spin multiplets [19].The theory and computation of spin-orbit coupling has a detailed review by Marian [31].Because spin-orbit coupling has such a large effect on the effective spin-dipole, the reliable evaluation of the g tensor is directly dependent on how well this effect is modelled.
For cases of a weak spin-orbit interaction, density functional theory (DFT) and linear response theory (DFT-LRT) can successfully be used [32].However DFT-LRT is unable to produce reliable results where degeneracies or low-lying excited states exist, or where the ground state is multiconfigurational in nature.Instead we have used a multiconfigurational approach that can describe the ground and excited states equally well, i.e., state-averaged complete active-space self-consistent field theory (SA-CASSCF).It is possible to include spin-orbit effects at a variational level.This type of approach is available in the SPOCK-CI program [33], which is a spin-orbit configuration interaction (CI) package.While this would be the preferred approach it is computationally demanding and perturbation techniques can be used to represent the interaction more efficiently.
One multiconfigurational approach is to use double perturbation theory and a SOS formalism to generate the g tensor directly.The double perturbation theory approach can be modified for degenerate systems, but for nearly degenerate systems, individual contributions to the SOS depend inversely on the energy difference between states, if this energy difference is small the resulting perturbation theory can be divergent.The spin-orbit term can also be applied using a quasi-degenerate perturbation theory (QDPT) approach which constructs a matrix over the manifold of spin states.The diagonal elements are taken as the energies of the spin-free states ( e.g., from a SA-CASSCF calculation), and the off-diagonal elements couple the spin-free states under the action of the spin-orbit operator.Diagonalisation gives the spin-orbit energy levels and wave functions.
In the case of a molecule with spin S = 1 2 , the lowest energy pair of spin-orbit wave functions must be degenerate in the absence of an applied magnetic field due to Kramers theorem [34], which states that the degeneracy of a state with half-integer spin must be at least two in the absence of a magnetic field.The energy splitting ∆E Zeeman upon application of a magnetic field is given by the difference in expectation values of the Kramers pair, represented by |Φ ±1/2 , the subscript indicating the sign of the associated pseudospin, that is, (5) In the above equation, the magnetic moment operator µ is related to the magnetogyric ratio for the electron, g e = 2.0023193 and the spin, S, and angular momentum, L, operators by µ = −µ B (L + g e S).Squaring Equations ( 5) and (6) gives a relationship between the related tensor G = gg T , called the Abragam-Bleaney tensor [35], and the matrix elements of µ in the Kramers pair basis.The principal values of the Zeeman matrix are the square roots of the eigenvalues of the Abragam-Bleaney tensor [36].While this relationship is only valid for S = 1 2 , similar ones have been derived for S > 1 2 [37].This process has two drawbacks, the resulting g matrix is necessarily symmetric, and because of the use of a square root precursor matrix, there is a sign indeterminacy for the principal values of g.For weak spin coupling, the principal values must be close to the free electron value, and therefore positive.Otherwise, the signs can be inferred by exploiting relationships between the sign of combinations of the principal values and the sign of combinations of the eigenvalues of the µ κ matrices (where κ represents a cartesian direction) [37].
The final consideration is the treatment of the spin-orbit operator.It is possible to use the full two-electron form of the operator, but as with two-electron repulsion integrals, this results in a considerable computational cost.It is often preferable to use a one-electron operator.Three operators are in common use, the first an empirical operator based on an effective nuclear charge calibrated against calculations on small molecules as developed by Kosekii [38][39][40].The second is the spin-orbit mean field (SOMF) operator which forms an effective one-electron operator, similar to the Fock operator in Hartree-Fock theory.Finally the atomic mean-field approximation (AMFI) to the spin-orbit operator forms a similar one-electron operator, but discards the multicentre two-electron spin-orbit integrals, exploiting the short range of the two-electron spin-orbit operator while retaining the simpler one-centre, two-electron integrals [41,42].

The ZFS Tensor, D
The ZFS interaction primarily represents the interelectronic dipolar interaction, and can be split into a first order term, directly corresponding to spin-spin interaction and a second-order term representing corrections due to spin-orbit coupling [43].When represented over a manifold of spin states S, M S , the first order term is diagonal.It is evaluated over the ground state eigenfunction.The second-order term includes off-diagonal elements with selection rule −1 ∆S, ∆M S 1.This term must be evaluated over ground and excited states.Available computational approaches are similar to those used for the Zeeman matrix.Single reference wave function methods include a modified DFT coupled-perturbed approach to handle the spin-orbit perturbation [44].In multireference approaches, the spin-spin contribution can be obtained from first-order perturbation theory, while the spin-orbit contribution requires second-order perturbation theory.However, as for the Zeeman matrix, near degeneracy often produces spurious results.A QDPT effective Hamiltonian approach gives more reliable ZFS tensors in such cases [45], and it is the latter approach that we have adopted.We shall outline below for each molecule the choice of active space used and states averaged in producing the ZFS tensor.

The Hyperfine Coupling Matrix, A
The hyperfine coupling constant (HFCC) matrix represents the interaction between the nuclear and the electronic spin-dipoles, and like the Zeeman matrix, the matrix for nucleus K can be related to an operator with matrix F K,κ = S κ • A K .For a Kramers doublet it is given by [17] where g K is the g-factor for the target nucleus, and M p the mass of the proton.A wide range of approaches are available, such as DFT linear response and SOS double perturbation methods.However, the hyperfine interaction has strong inverse dependence on electron-nucleus separation and is strongly influenced by relativistic effects [6].Relativistic corrections cannot be ignored, at least scalar relativistic treatments are needed for a reliable approach.However additional corrections for spin-orbit coupling are computationally demanding, since three sets of coupled-perturbed equations must be solved for each nucleus.In principle, it is possible to avoid the solution of the CP-SCF equations by adopting a 4-component approach [46,47], although this has an associated computational overhead in performing the 4-component calculation.For a light target nucleus the corrections are small and often can be safely ignored, except when the target nucleus is bonded to a heavy nucleus; the latter can have a substantial effect on the electronic structure of its neighbour, a phenomenon referred to as the "heavy atom-light atom" (HALA) effect [48], and the spin-orbit correction to the hyperfine interaction must be included for reliable results.
A remaining issue in predicting isotropic chemical shifts is due to our calculations being based on a single geometry, ignoring vibrational and rotational motion.This is a valid approach for electronic transitions which operate on a much faster timescale than significant nuclear motion, but nuclear transitions are slow enough for the motion to be significant.The pNMR measurement is an average over nuclear motion.It is necessary to calculate the orbital shielding and hyperfine matrices for all nuclei and then average them in appropriate groupings to generate the averaged hyperfine matrices for each group of equivalent atoms.

Zero-Field Spitting-Theory
Ramsey [12,23,49] posited that because nuclear motion is much slower than electronic motion, the changes in the nuclear magnetic moment µ nuc changes much more slowly than the electronic magnetic moment µ, and can be treated as time-independent, effectively as a parameter.Assuming equilibrium is reached within experimental time scales, the observed NMR shielding tensor can be shown to be given by the second derivative of the electronic energy with respect to the components of the magnetic field and the nuclear magnetic moment at the limit of zero field and magnetic moment.
The resulting energy levels are often close enough to have significant occupation above the ground state, so Boltzmann averaging of observables is required.The Boltzmann average A 0 of a set of observables {A n } with n energy levels, {W n } is given by: Following Moon and Patchkovskii's summary [50] the approach we use is based on Kurland and McGarvey's work [51] which allows the pNMR Hamiltonian to be split into Zeeman and hyperfine terms.The field-independent form of the shielding tensor [50] can be decomposed into diamagnetic and paramagnetic shielding tensors, σ orb and σ p respectively: Taking g and A as the following partial derivatives it can be shown that, in the low spin-orbit coupling limit [10,11,13] where Z is a 3 × 3 coupling matrix.In the absence of zero-field splitting, Z κυ = 1 3 S(S + 1)δ κυ , where δ κυ is the Kronecker delta and the resulting matrix is diagonal.
Pennanen and Vaara [22] formulated a practical approach to modelling ZFS in the context of modern quantum chemistry, the "no coupling approximation."This requires the formation of the (2S + 1) × (2S + 1) matrix, S • D • S, which is then diagonalised to obtain eigenvalues E n and eigenstates |n , the latter being linear combinations of pure spin states |SM S .Soncini and van den Heuvel [13] presented a fully coupled approach, which can be expressed as: with q nm being a Boltzmann weighting term Pierre Curie observed that the magnetisation of a paramagnetic material is approximately inversely proportional to the temperature.If the occupation of the ground state ( i.e., its Boltzmann weight) is nearly one, this is observed for a doublet and in the no coupling approximation.In the coupling described in Equations ( 12) and ( 13), the collection of terms corresponding to E m = E n have approximate 1/T dependence and the combined contribution is referred to as the Curie term.The remaining terms have an approximate 1/T 2 dependence.
The derivation used for both approaches assumes weak spin-orbit coupling.Van den Heuvel and Soncini presented an extended approach [52] which is valid for strong spin-orbit coupling.Because of its potential as an improved treatment, an outline of the approach is presented in Appendix A. The approach depends on using Equation A1 to form irreducible tensor operator (ITO) representations of the µ and F K matrices introduced in Sections 2.2 and 2.4 respectively.

Contact and Pseudocontact Shifts
The separation of the paramagnetic chemical shift into contact and pseudocontact shifts provides insight into to the chemical environment of the target nucleus; however, as additional effects are included in the model, the separation becomes increasingly difficult and debatable.The simplest approach is to decompose the Zeeman and hyperfine coupling matrices into isotropic and anisotropic components, and evaluate their nature and relative sizes in terms of their scaling with respect to the fine structure constant α [22,53].
In the decomposition of the Zeeman matrix, the only zeroth-order term corresponds to the magnetogyric ratio of the free electron, the remaining terms relate the g-shift, ∆g = g − g e I 3 which is decomposed into an isotropic scalar g-shift, ∆g iso = 1  3 trace (∆g), and the remaining anisotropic matrix, ∆ g = ∆g − ∆g iso I 3 .The hyperfine matrix has no first order terms, but is entirely second-order in the absence of relativistic corrections and we will denote the scalar relativistic (strictly the scalar relativistic corrections implicitly include O(α 4 ) terms, but to simplify the analysis we will consider A SR to be a O(α2 ) term.) matrix as A SR .This is decomposed into its isotropic scalar A iso =1 3 trace A SR and an anisotropic spin-dipolar matrix A SD = A SR − A iso I 3 , which is symmetric.The relativistic corrections are entirely fourth order, and are separated into an isotropic scalar A ∆SO iso and an anisotropic A ∆SO aniso matrix.The anisotropic component is further split into a symmetric matrix A ∆SO,SYM = 1 2 A ∆SO aniso + A ∆SO aniso T and an asymmetric matrix A ∆SO,AS = A ∆SO aniso − A ∆SO,SYM .These decompositions are outlined in Table 1.
The product of the matrices, g • A, has fifteen terms, six of which are sixth-order in α.In this work we use two methods to qualitatively divide the paramagnetic chemical shift into contact and pseudocontact contributions [50].The first method is to relate the terms derived from isotropic scalars of the hyperfine matrix to the contact shift and those from the anisotropic matrices to the pseudocontact shift [10,54].The second is to truncate the expansion to fourth order, resulting in nine remaining terms.The remaining difference is that all contributions with an anisotropic component, most notably the term A iso ∆ g (a contact term in method 1), are considered to relate to the pseudospin [11,55,56].We will refer to these decompositions as methods 1 and 2 respectively, as outlined in Table 2. Method 2 is elaborated on by Kaupp [7], who presents a separation of terms for a doublet.In EPR terminology, A ∆SO iso is sometimes referred to as the pseudocontact term.The anisotropic matrices are traceless, so without ZFS, Z is diagonal, so only the isotropic hyperfine/isotropic Zeeman and the anisotropic hyperfine/anisotropic Zeeman terms produce an interaction with a non-zero trace.In the absence of spin-orbit coupling these interactions can unambiguously be related to the contact and pseudocontact terms, respectively.However the nature of the fourth order spin-orbit corrections to the hyperfine is unclear and the application of zero-field splitting results in non-diagonal elements in Z which couple isotropic and anisotropic terms.As a result, the separation of magnetic chemical shifts into contact, bond-mediated effects and pseudocontact, dipolar-interaction effects is not strictly consistent.With this in mind, we will use both methods to analyse our results.
The truncation to fourth order terms used in the second method is based on the assumption that the order in α determines the size of the contributions.This implies that the spin-orbit corrections to the hyperfine are small compared to the scalar relativistic hyperfine.If these, computationally demanding, corrections are not included, no sixth-order terms are present, so the final shift is the same for both decompositions.However in cases of hyperfine calculations on heavy nuclei or on nuclei that are bonded to heavy nuclei, this assumption is not valid: spin-orbit corrections typically scale in terms of the atomic number to the fourth power (Z 4 ) making the spin-orbit correction considerably larger.In such cases, the isotropic spin-orbit correction may be larger than the isotropic value of the uncorrected hyperfine matrix [17].

An Example of a Spin-Orbit Description of a Molecular Ground State
Our calculations allow for mixing of electronic states under the action of the spin-orbit operator.The number of electronic states included in each calculation was motivated by the quality of the Zeeman matrix produced, as discussed in our previous paper [57].As an example, for our first actinide molecule (Section 4.2) we performed state-averaged CASSCF calculations over 14 quartet and 25 doublet states (i.e., being a subset of the 35 quartet and 112 doublet states implicitly included in the CASSCF(7,3) configuration state function (CSF) expansion.)to generate the Zeeman and ZFS matrices.The resulting ground state wave function, corrected for spin-orbit coupling is detailed in Table 3.
Table 3.One of the Kramers pair wave functions of U(MeBTP) 3+  3 , spin-orbit corrected, SA-CASSCF(3,7) wave function, with 14 quartet and 25 doublet states.The second state can be inferred, via Kramers degeneracy, exchanging the weights of each +M S state with that of the corresponding −M S state.States with weight less than 0.01 are omitted.Chibotaru states that the spin Hamiltonian approach is applicable when a QDPT-SOC calculation is convergent for a manifold of spin states [37], and later he discusses the differences between the fictitious pseudospin and the spin.In the limiting cases of weak SOC the pseudospin can be set equal to the spin.In the limiting case of strong SOC, the degeneracy of the states is broken in low symmetry molecules, and the pseudospin mirrors a doublet, S = 1/2 [17].Chibotoaru considers actinide ionic complexes to fall in the intermediate range [37], where the model Hamiltonian cannot be defined by the total angular momentum because of the mixture of states that we describe above; however, he does still associate the pseudospin with a manifold of electronic states, and so this is the approach we employ.Furthermore, in systems for which S 3/2 the calculation ideally should include third order tensorial spin, Zeeman and Hyperfine matrices [52], as we summarise in the Appendix A. However, it must be noted that to implement this scheme would require the generalised hyperfine tensors in Equation (7), that are not currently generally available.We are not aware of the use of the generalised hyperfine operators by any other researchers.
In Table 4 we list the relative energies of the first five Kramers pairs for all molecules studied in this work.

Previous Work and the ZFS Tensor
In a recent study we looked at the pNMR shifts of 14 f-element compounds [57].The present work studies the effects of zero field splitting on the six f 3 compounds from that study.One compound contains neodymium as a sample lanthanide, detailed in Section 4.1 and shown in Figure 1.The remaining five are uranium compounds detailed in Sections 4.2-4.6.Equation (4) shows that the pNMR shift depends explicitly on the four magnetic interactions (σ orb , g, D and A) and implicitly on the geometry at which the various terms are evaluated.We summarise the findings of our previous study in justifying the approach we have adopted in evaluating each of the necessary terms.More details can be found in the aforementioned paper.

Geometry Optimisation, R
Geometry optimisation was performed using GAUSSIAN 09 [58].Since geometries largely depend on the ground state electronic structure, it was considered reasonable to use DFT with the PBE0 [59] functional.Due to the presence of uranium, relativistic effects are important, but it is sufficient to only consider scalar relativistic effects.We prefer an all-electron approach and previously created a program in which implements the ZORA [27,60], and applies the correction to the one-electron Hamiltonian in atomic blocks, removing the gauge dependence issue.This program was interfaced to the GAUSSIAN programs and used in all geometry optimisations.We used the segmented all-electron relativistically corrected (SARC) [61,62] basis set for the heavy atoms with Def2-SVP [63] for the remaining atoms.This approach showed reliable agreement against a test set of crystal structure data for a set of 11 molecules, producing correlation coefficients of R 2 = 1.00 for bond length and R 2 = 0.99 for bond angles [57].

Orbital Shielding, σ orb
The orbital components of the chemical shifts were obtained using GAUSSIAN 09, DFT linear response theory (PBE0), the ZORA relativistic treatment outlined above, and GIAO [25,64].A Def2-TZVP basis was used for 1 H chemical shifts, and an uncontracted Def2-SVP basis for the 13 C, otherwise the basis used was as in the previous section.This approach was calibrated against 20 1 H and 32 13 C chemical shifts of 4 diamagnetic molecules, producing correlation coefficients for 1 H and C 13 of R 2 = 0.99 and R 2 = 0.98 respectively [57].Deviations were noted when a light atom under investigation was directly bonded to a heavy atom.These cases are subject to the HALA effect working through the chemical bond and as such, need a more elaborate relativistic treatment [48].

Hyperfine Coupling Matrix, A
The reported A tensors were calculated using the ORCA 3.0.3program [65] , DFT (PBE0), the Douglas-Kroll-Hess transformation to second-order (DKH2) [66,67].Analysis of the effect that picture change has on these hyperfine coupling coefficients (without spin-orbit correction) was performed in our previous paper [57], and except for 1 H calculations, picture change effects were found to be significant.Therefore they were included as implemented in ORCA and documented by Sandhoefer et al. [68].The second-order DKH transformation is applied to the Fermi contact, spin-dipole and Zeeman operators.The transformation of the paramagnetic spin-orbit term was treated to first order only.Relevant experimental results are extremely scarce, but the approach was compared against [U(C 7 H 7 ) 2 ] − , showing reasonable accuracy.
For the nuclei of interest, a Def2-TZVP basis was used with three additional uncontracted s orbitals with large exponents added to allow improved description of the electronic environment close to the nucleus, otherwise the basis used was as in the previous section.The exponents were derived using a geometric progression starting from the largest s function exponents of the existing basis, following the requirements of basis set design for hyperfine coupling calculations outlined by Chipman [69].Because of the high computational cost, the calculation generally did not include spin-orbit corrections to the hyperfine calculation.We indicate in the text the two molecules for which this term was evaluated.In these cases the calculation used the same methods as detailed above for the other hyperfine terms but included solution of the relevant coupled-perturbed Kohn-Sham equations to obtain the derivative of the spin density, followed by contraction with the nuclear (spin)-electron (orbit) integrals.

Zeeman Coupling Matrix, g
It is possible to calculate g tensors with DFT and use linear response theory to include spin-orbit coupling as a first-order perturbation [17,70]; however, this approach is only effective for weak spin-orbit coupling and in the absence of low-lying excited states.For the molecules under investigation here neither condition is true, and the ground states are all clearly degenerate.As a result a multiconfigurational approach is necessary.ORCA 3.0.3was used to perform a series of State-Averaged Complete Active-Space Self-Consistent Field (SA-CASSCF) calculations, using an active space of 3 electrons and 7 orbitals, i.e., the f electrons of the heavy element distributed in the relevant manifold of f orbitals.Relativistic effects were added using DKH2.Spin-orbit coupling was included by utilising the SOMF operator [71], using the centre of electronic charge as the origin of the molecule to minimise gauge dependence [72].A SARC basis set was used for uranium, with a Def2-SVP for hydrogen and the remainder using a Def2-TZVP.The g tensor was evaluated within the effective Hamiltonian approach.The number of states included in the SA-CASSCF calculations were explored empirically, by increasing the number of states, until some convergence in the principal values of the g tensor was achieved [17,60].

The ZFS Tensor, D
For efficiency, the D tensor was calculated at the same level of theory as detailed in the previous section to obtain g.

1 H pNMR Shifts of Neodymium Tetrapyridyl Appended Cyclen
The first molecule under investigation is a neodymium 1,4,7,10-tetrakis(2-pyridylmethyl)-1,4,7,10-tetraazacyclododecane complex [ Nd L py ] 3+ with a triflate CF 3 SO − 3 counterion.This was synthesised and studied by Natrajan et al. [73], who produced 500Mz 1 H NMR spectra of the species in a D 2 O solvent at room temperature.The molecule is shown in Figure 1.The observed shifts were −0.2, 2.4, 4.6, 5.2, 6, 8, 8.3, 8.4, 10.8 and 12.4 ppm and were not assigned in the paper.In our previous work [57], we found that inclusion of the triflate counterion produced a significant change in the geometry of the complex and so needs to be included in modelling this structure.S1 of the ESI.The resulting isotropic g value is negative, a phenomenon observed for states with significant spin-orbit coupling as discussed in reference [37,74].
We find in the absence of ZFS corrections from our calculations a range of shifts spanning −2.56 -16.51 ppm.Inclusion of ZFS produces a wider spread of values which are compared in Figure 2. Table 5 outlines the details of the components.It is noticeable that adding ZFS has very little effect on the contact part of the shielding, but the size (and composition) of the pseudocontact term significantly changes.Examining the pseudocontact components, the term derived from the isotropic hyperfine changes only slightly and the interaction between the anisotropic Zeeman and the spin-dipole hyperfine term is negligible.The change derives from the interaction of the spin-dipole part of the hyperfine and the isotropic g term (decomposed into the free electron part and the shift).While the predicted chemical shift changes by only a few ppm, it is sufficient to shift the predicted order of the peaks and as a result assignment of all ten 1 H shifts are not possible.We have omitted the spin-orbit correction to the hyperfine coupling matrix as the calculations are extremely demanding.However it is unlikely that these 1 H shifts will be strongly influenced by this omission since the protons are well removed from the heavy atom.

4.2.
1 H pNMR Shifts of Tris(2,6-bis(5,6-dialkyl-1,2,4-triazin-3-yl)pyridine) Uranium(III), U(MeBTP) 3+ 3 2,6-Bis(5,6-dialkyl-1,2,4-triazin-3-yl)pyridine (BTP) is an effective chelating agent and has been investigated as a reagent for processing waste containing uranium(III) species [4,75], with the family of ligands showing preference for coordinating selectively with trivalent actinides over trivalent lanthanides due to the "much softer character of the N atoms in this weakly basic ligand" [4].It is thought that actinide/lanthanide selectivity results from the difference in the nature of the f orbitals, the lanthanide 4f orbitals displaying "core-like" behaviour and participating weakly in bonding, whereas the more extended 5f orbitals of actinide species participate more in covalent bonding [76].As such there is much interest in better understanding the nature of this U-N bond.The experimental results that we are using for comparison are 200 MHz 1 H NMR spectra taken at 30 • C in deuterated pyridine [4].The molecule has four chemically equivalent types of hydrogen atoms.The numbering used to refer to these sets of protons is shown in Figure 3 and the stereochemistry in Figure 4.All hydrogen atoms are relatively distant from the U-N bond, and as such we expect that spin-orbit component of the hyperfine interaction is not significant, hence they have not been included here.6 and 7.For each distinct type of hydrogen we show the effect of the ZFS correction.We list the contribution to the paramagnetic shift from each component of g • Z • A and a qualitative partition of contact and pseudocontact terms.Overall the paramagnetic shift is calculated to be slightly larger than the orbital contribution, except for the H (2) atoms.In the absence of the ZFS correction H (1) and H (3) proton chemical shifts are predicted to be larger in magnitude that the observed shifts, although the signs of the shifts are predicted correctly.The ZFS correction reduces the magnitudes in the direction of the experimental values.In this system including the ZFS produces shifts that are closer to the observed values but are still too large in magnitude.The experimental shifts span a range of about 52 ppm, while the neglect of ZFS altogether gives a span of calculated values of 117 ppm.The ZFS correction reduces the range to 90 ppm.This is shown in Figure 5.The contact and pseudocontact terms vary according to the method of partitioning them, but the contact shift is consistently predicted to be larger than the pseudocontact shift in all cases.However for H (2) and H (3) the values are much closer in magnitude.3 with and without ZFS at 303 K.For numbering of protons see Figure 3.The order is consistent in all cases.Table 6.Calculated 1 H pNMR shifts for the methyl H atoms of U(MeBTP) 3+  3 in ppm at 303K).For numbering of protons see Figure 3.The final, ZFS corrected, computational values are highlighted.

Component
Methyl H (1)  Methyl H (2)   A g Order  It was originally thought impossible to synthesise an organometallic complex consisting of a single atom and three η 5 -bonded pentamethylcyclopentadienyl rings due to the steric hindrance presented by the methyl groups.However a samarium complex was synthesised in 1991 followed by a uranium complex shortly after.The optimised geometry of the complex is shown in Figure 6, with the methyl carbons residing above the plane of the ring due to steric effects.The paper reporting the synthesis [77] also reported the chemical shifts observed in 1 H and 13 C spectra taken at 500 MHz, 25 • C in a deuterated benzene solvent and these act as our baseline for comparison.The SA-CASSCF(3,7) calculation of the g and D matrices used 28 quartet and 55 doublet states.The matrices are listed in Table S3 of the ESI.Our results for the methyl hydrogen atoms are shown in Table 8, again reflecting the treatment of ZFS, the individual contributions to the paramagnetic shift and the qualitative partition of contact and pseudocontact terms.Overall the paramagnetic shift is calculated to be significantly smaller (14-17%) than the orbital contribution.−0.90 In the absence of the ZFS correction, the predicted chemical shift was −3.02 ppm, with the correction the predicted chemical shift was −3.36 ppm, hence the full ZFS correction only contributes −0.34 ppm to the final computed shift.This is unlikely to be significant compared to the errors associated with the approximations used to generate the pNMR matrices, and this result remains a good match for the experimental value.Both the contact and pseudocontact shielding terms are positive, although the two decompositions used disagree with respect to the relative magnitudes of the contact and pseudocontact shifts.It must be noted that the hyperfine couplings have been evaluated using a DFT method.The choice of exchange-correlation functional can introduce a variation in the predicted hyperfine operators.It is difficult currently to deal with this property using purely wave functional techniques, which would in principle would allow a degree of systematic refinement.Given the size and complexity of the systems studied here, the approach is a suitable starting point.However, as wave function methods become more available for these operators further refinements of the computations can and should be sought [78].
Table 9 presents the shifts for 13 C pNMR shifts for the ring and methyl carbon atoms.Due to the computational expense of calculating the spin-orbit corrections to the hyperfine coupling matrix and the fact that 13 C is a light atom, we initially did not apply this correction to our calculations.However this complex contains a set of η-coordinated carbon-uranium bonds, and the HALA effects alluded to earlier are likely to be in operation.Additionally, in the work of Autschbach et al. on actinide carbonate molecules [14,15], it was noted that the spin-orbit correction dominated the hyperfine interaction.Accordingly, we computed the spin-orbit correction to the hyperfine tensor using the methodology detailed in Section 3.3, and report our findings in Table 10.The format of this table is similar to the previous ones, but the individual contributions have been expanded.The table includes a comparison with the values before the spin-orbit correction to the hyperfine matrix was added and the bottom section highlights the two decomposition methods.Table 8 shows that the predicted isotropic chemical shielding (and chemical shift) of the two decomposition methods is different.This is due to the neglect of the O(α 6 ) terms that method 2 uses, which is not valid when spin-orbit effects are significant.With reference to the method 1 decomposition, we note that the effect of the inclusion of the spin-orbit hyperfine term is less marked on the methyl carbon atoms than for the ring carbons.There is now qualitative agreement with experiment for the latter.For the methyl carbons, the effect is smaller, perhaps due to their being further removed from the heavy nucleus.Figure 7 shows the deviation from the averaged value for the isotropic paramagnetic shielding of the ring and methyl 13 C nuclei.Overall the ring 13 C range from −63 -92 ppm and the methyl 13 C range from −74 -94 ppm.Note that the experimental values of the shifts lie within both ranges.Our next molecule is a single molecule magnet which, unusually, has a trigonal planar coordination around the uranium atom when in the solid state, instead of the more common trigonal pyramidal geometry that is seen in other three-coordinate actinide compounds.This is due to the large steric hindrance introduced by the t-butyl groups.Steric hindrance is a commonly used tool to prevent the molecules oligomerising, and there is much interest in understanding how the magnetic properties of the molecule are affected by the ligand.For comparison with experiment, we use the 73 Hz 29 Si chemical shift of −296.04 ppm taken at 25 • C in a solution with deuterated benzene [2].Our calculation refers to an optimised gas phase structure [57], and is shown in Figure 8.
The SA-CASSCF(3,7) calculation of the g and D matrices used 21 quartet and 39 doublet states.The matrices are listed in Table S4 of the ESI.Table 11 shows the resulting ZFS calculations.The orbital and paramagnetic contributions are predicted to have a similar size, and the contact interaction is positive and much larger than the pseudocontact, so the paramagnetic shift is largely due to increased spin density at the nucleus.Without applying ZFS, the chemical shift predicted was 18.54 ppm less than the experimental value.Using the full coupling method, the predicted chemical shift reduced by 46.19 ppm, which is further from the experimental value.Given our findings for the 13 C shifts in U(C 5 Me 5 ) 3 , the spin-orbit correction to the hyperfine tensor cannot be neglected.Due to computational cost we have only obtained the corrected tensors for two silicon nuclei.The hyperfine matrices of the two silicon atoms considered result in very different isotropic shielding, specifically 559.96 and 371.12 ppm.We note that in our optimised structure, the uranium-silicon distance of the first atom is 3.31 Angstroms whereas the second is 3.66 Angstroms.The magnitude of the change in the isotropic shielding for the silicon atoms that we have computed clearly indicates that it is essential to carry out the average over the full set of silicon atoms to obtain a reliable estimate of the shifts.This remains to be done and we are exploring more computationally efficient ways to obtain the spin-orbit hyperfine interaction.−296.00Tris(trimethylsilyltetramethylcyclopentadienyl)uranium(lll) demonstrates even more steric crowding than the pentamethylcyclopentadienyl ring studied in Section 4.3, and was part of a systematic study of 29 Si shifts of similar compounds, focused on discovering trends [79] in the chemical shifts of uranium-silicon compounds.The SI for reference [79] includes a graph that may indicate very weak correlation between U-Si distance and 29 Si, but cautions against treating their results as definitive, citing the complexity involved when studying paramagnetic species and effects that are in opposition.That study reported an experimental value of −155 ppm (99.2 MHz spectra, taken at 298 K in deuterated benzene with an iodide counter anion).The structure is shown in Figure 9.  S5 of the ESI.Table 12 shows the resulting ZFS calculations.Our calculated value before splitting was +18.92 ppm from the observed value, applying ZFS using the full coupling approach resulted in a small shift of +1.73 ppm, so in this case ZFS has very little effect on the total shift.For silicon nuclei, the implication of our results presented in Section 4.3 is that spin-orbit effects must also be important here.Further investigation by including spin-orbit corrections to the hyperfine matrices is needed to clarify if this is the case.Although the situation is slightly different in that the silicon atoms are not directly bound to the heavy atom and this may ameliorate the relative importance of including the spin-orbit part of the hyperfine interaction.Mixed-sandwich organo-uranium complexes have been used for more environmentally friendly synthesis of oxocarbo rings from carbon dioxide and carbon monoxide with the ability to reform the complex after use.Similar uses for olefin processing exist, and there is interest in understanding this process better [1].The structure of the compound studied is shown in Figure 10, and the cited study reported an experimental value of −136.70 ppm (99.2 MHz spectra, taken at 298 K in deuterated benzene with an iodide counter anion).
The SA-CASSCF(3,7) calculation of the g and D matrices used 13 quartet and 32 doublet states.The matrices are listed in Table S6 of the ESI.Inspection of our results in Table 13 shows that the contact term is substantially larger than the pseudocontact, and as for the example in Section 4.5, we interpret this as implying a strong interaction through the C-Si bond.−136.70

Summary
Considering our six examples as a whole, Figure 11 shows the calculated values of pNMR shifts of the five uranium molecules plotted against the experimental values.Since this is a small sample, firm conclusions are not possible but, perhaps not surprisingly, the correlation coefficients improve as the treatment of ZFS improves; from R 2 = 0.86 for the treatment without ZFS to R 2 = 0.90 with of ZFS.

Figure 1 .
Figure 1.Structure of the cyclen complex [ Nd L py ] 3+ .Labels refer to the 10 chemically distinct H atoms, with 1,3,5 representing equatorial CH 2 protons and 2,4,6 representing axial protons.Colours indicate atom type (grey C, blue N, red O, yellow S, light blue F).The calculation of the g and D matrices used 21 quartet and 40 doublet states from a SA-CASSCF(3,7) calculation.The matrices are listed in TableS1of the ESI.The resulting isotropic g value is negative, a phenomenon observed for states with significant spin-orbit coupling as discussed in reference[37,74].We find in the absence of ZFS corrections from our calculations a range of shifts spanning −2.56 -16.51 ppm.Inclusion of ZFS produces a wider spread of values which are compared in Figure2.Table5outlines the details of the components.

Figure 2 .
Figure 2. Experimental and computed pNMR spectra for [ Nd L py ] 3+ with and without ZFS at 298.15 K.For numbering of protons see Figure 1.

Figure 4 .
Figure 4. Stereochemistry of U(MeBTP) 3+ 3 , showing the three planar ligands, staggered with approximate 3-fold rotational symmetry.The SA-CASSCF(3,7) calculation of the g and D matrices used 14 quartet and 25 doublet states.The matrices are listed in Table S2 of the ESI.Our results are shown in Tables6 and 7.For each distinct type of hydrogen we show the effect of the ZFS correction.We list the contribution to the paramagnetic

Figure 5 .
Figure 5. Experimental and predicted pNMR shifts of U(MeBTP) 3+3 with and without ZFS at 303 K.For numbering of protons see Figure3.The order is consistent in all cases.

Figure 7 .
Figure 7. Difference of the predicted isotropic shielding of individual 13 C nuclei from the averaged value at 298 K.

Figure 8 .
Figure 8.(a) Chemical structure of U(N(SiMe t Bu 2 ) 2 ) 3 , (b) Stereochemistry of U(N(SiMe t Bu 2 ) 2 ) 3 with uranium and silicon nuclei shown as spheres and the remaining non-hydrogen nuclei in wireframe.Stated distances are in Angstroms.
The SA-CASSCF(3,7) calculation of the g and D matrices used 21 quartet and 39 doublet states.The matrices are listed in Table

Figure 11 .
Figure 11.Calculated results against experimental results.

Table 1 .
Decomposition of the components of the Zeeman and hyperfine coupling matrices by isotropic and anisotropic components, I 3 indicates a 3 × 3 identity matrix.

Table 2 .
Classification of components of the g • A matrix.

Table 4 .
Relative energies (cm −1 ) for the 10 lowest spin-orbit states of the six compounds studied here, showing the increased gap between first and second excited pairs compared to that of the first pair and the ground state pair.

Table 5 .
Calculated 1 H pNMR shifts of [ Nd L py ] 3+ organised in increasing δ.For numbering of protons see Figure1.The reference chemical shift is 31.58ppm (TMS) and the temperature is 298.15K.

Table 7 .
Calculated 1 H pNMR shifts for the ring H atoms of U(MeBTP) 3 in ppm at 303 K.For numbering see Figure5.The final, ZFS corrected, computational values are highlighted.H and 13 C pNMR Shifts of Tris(pentamethylcyclopentadienyl) Uranium (III), U(C 5 Me 5 ) 3 1

Table 8 .
Calculated 1 H pNMR shifts for U(C 5 Me 5 ) 3 in ppm at 298 K.The final, ZFS corrected, computational values are highlighted.

Table 9 .
Calculated13C pNMR shifts for U(C 5 Me 5 ) 3 in ppm with no spin-orbit corrections to the hyperfine term at 298 K).The final, ZFS corrected, computational values are highlighted.

Table 10 .
Calculated13C pNMR shifts for U(C 5 Me 5 ) 3 in ppm with spin-orbit corrections to the hyperfine term at 298 K).The final, ZFS corrected, computational values are highlighted.
* Rows marked A-no SOC refer to the results before the hyperfine matrix has been corrected with spin-orbit effects, A-with SOC are results including this correction.