Solution of a Puzzle: High-Level Quantum-Chemical Treatment of Pseudocontact Chemical Shifts Confirms Classic Semiempirical Theory

A recently popularized approach for the calculation of pseudocontact shifts (PCSs) based on first-principles quantum chemistry (QC) leads to different results than the classic “semiempirical” equation involving the susceptibility tensor. Studies that attempted a comparison of theory and experiment led to conflicting conclusions with respect to the preferred theoretical approach. In this Letter, we show that after inclusion of previously neglected terms in the full Hamiltonian, one can deduce the semiempirical equations from a rigorous QC-based treatment. It also turns out that in the long-distance limit, one can approximate the complete A tensor in terms of the g tensor. By means of Kohn–Sham density functional theory calculations, we numerically confirm the long-distance expression for the A tensor and the theoretically predicted scaling behavior of the different terms. Our derivation suggests a computational strategy in which one calculates the susceptibility tensor and inserts it into the classic equation for the PCS.

I n this Letter, we study a paramagnetic center and a sufficiently far nucleus with magnetic moment M K . In an external magnetic field B, this nucleus experiences an NMR chemical shielding that is dipolar in origin. Thus, this shielding originates from an additional dipolar magnetic field that is created by the induced average magnetic moment of the paramagnetic center. Because of the dipolar origin, the resulting "pseudocontact shifts" (PCSs) have an inverse cubic dependence on the distance of the nucleus from the paramagnetic center. Consequently, their analysis provides a wealth of structural information that was heavily exploited during the last decades in chemistry 1−6 and structural biology. 7−15 PCSs are also gaining momentum in the characterization of the electronic structures of single ion magnets. 6,16−18 In the following, we use Hartree atomic units (ℏ = m e = e = 4πϵ 0 = 1) throughout. Traditionally, PCSs are often interpreted by the "semiempirical" (SE) expression for the chemical shielding tensor σ: where R is the distance vector between the nucleus and the paramagnetic center and χ is the paramagnetic susceptibility tensor. The PCS is then given by 19,20 where Δχ ax = χ zz − 1 / 2 (χ xx + χ yy ) and Δχ rh = χ xx − χ yy are the axial and rhombic susceptibility parameters and θ and ϕ are the polar and azimuthal angles of the direction of the nucleus in the principal axis system of the susceptibility tensor. In recent years, alternative treatments of PCSs in terms of rigorous first-principles quantum chemistry (QC) have been developed. A very general approach was given by Van den Heuvel and Soncini 23 and leads to an expression for the chemical shielding tensor as a mixed second partial derivative of the electronic Helmholtz free energy F = −k B T ln Z (where Z = ∑ i e −βE i ): The zeroth-order eigenfunctions and energies in eq 6 are derived from the zeroth-order Hamiltonian in the SH approximation, H spin (0) = S·D·S. The analogous SH expression for the susceptibility tensor is 28,29 In some previous studies on chemical shieldings based on eq 5, 28,30−33 all of the contributions to the hyperfine coupling (HFC) tensor A except for the Fermi contact (FC) and spindipolar (SD) ones were neglected. Under this assumption, the chemical shift in the point-dipole approximation (PDA) turns out to be 26,28,30,34 with the unsymmetrical tensor χ′ related to χ in eq 7 by χ = χ′· g T /g e . The appearance of this unsymmetrical tensor was considered a necessary consequence of the more rigorous QC treatment compared with the "semiempirical" eq 1. 26,28,30 The two approaches are ostensibly different. 26 Several recent papers prefer the use of eq 8 because it seemingly gives better agreement with the experimental data. 28,32,33 Even recent and authoritative reviews 29 indirectly contribute to endorsing the use of eq 8. Therefore, as the use of eq 8 is becoming increasingly attractive and the two equations can lead to significantly different results, 26,35 it is important to establish as soon as possible and beyond doubt which of the two expressions should be used for the calculation of PCSs. Comparison between the calculations and the experimental data does not resolve the issue in full: paramagnetic NMR calculations from first principles involve the estimation of SH parameters (the g and D tensors), which are not directly accessible experimentally at room temperature in solution for most metal ions. Furthermore, the uncertainties in the calculation of these parameters can collectively be on the same order of magnitude as the expected difference between the two approaches. For this reason, a recent study 34 has investigated the performance of the two equations for a copper(II) ion with S = 1 / 2 in an axial protein environment, where the results are expected to differ by a factor of roughly 2. This comparison is particularly strong, because it does not rely on calculations of the SH parameters but rather depends on the direct experimental measurement of the g anisotropy (with uncertainty smaller than 1%) under the same conditions in which the NMR data are also measured. It was found that the "semiempirical" equation performs much better for the prediction of the PCSs than the equation derived from first principles. 34 To date, a QC-based explanation for this result is missing.
A first hint at possible problems with eq 8 comes from the insight that the same orbital currents, induced by spin−orbit coupling (SOC), are responsible for the orbital contributions to both the g tensor and the HFC tensor. Equation 8 incorporates this orbital contribution in the former but not the latter, i.e., in a very unsymmetrical way. 20,26,35 Starting from these premises, in this Letter we go beyond the previous QC-based treatments in two important aspects: (1) we include in the molecular Hamiltonian the so-called gauge correction terms, 36 which are necessary to preserve gauge invariance in the presence of SOC, and (2) we explicitly include all of the contributions (not only the SD one) to the magnetic hyperfine field and to the A tensor. With regard to the electronic Hamiltonian, we assume the usual nonrelativistic Hamiltonian with magnetic fields incorporated by adding the magnetic vector potential to the canonical momentum The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter operator, supplemented with SOC in the effective nuclear charge approximation. 37 We incorporate two contributions into the magnetic vector potential, one describing a homogeneous external magnetic field and the other describing the field created by a point magnetic dipole at the nucleus under consideration. In order to obtain gauge-invariant results, it is furthermore necessary to add gauge-correction terms 36 to the Hamiltonian. The detailed expressions for the complete Hamiltonian and its derivatives are given in the Supporting Information.
We now assume that the paramagnetic part of the system is localized in one region of space and that the nucleus of interest is at a large distance from this region. Our analysis is similar to that of McConnell 38 based on Ramsey's expression for nondegenerate ground states. 39,40 R O is a position in the region where the paramagnetic center is located, and R K is the position of the nucleus. We use R O also as the gauge origin for the homogeneous external magnetic field. Furthermore, we define r iO = r i − R O and r iK = r i − R K , where r i is the position of electron i. We assume that the only relevant contributions to the chemical shift come from electrons located at the paramagnetic center. Then r iO can be assumed to be small compared with R = R O − R K if the latter becomes large enough, i.e., r iK = r iO + R ≈ R. Following McConnell, inverse powers of r iK can then be expanded in a Taylor series around the point r iO = 0, which gives 38 Table 1. Slopes of the Log−Log Plots (as in Figure 2 for CO + ) at the Beginning and End of the Data Series for Different A Tensor Contributions a first two data points (small R) last two data points (large R) The used basis sets were decontracted cc-pV5Z for transition metal complexes and decontracted cc-pV6Z for radicals. b For NF 3 + , the calculations for the first distance (5 Å) did not converge. Hence, the first two data points correspond to 7 and 9 Å. Inserting these expressions into the equations for the Hamiltonian derivatives occurring in eq 3 that depend on the nuclear position and neglecting all terms with inverse powers of R larger than 3 (PDA), one obtains and In eq 11, ∂H/∂M K l contains FC, SD, paramagnetic spin−orbit (PSO), and gauge contributions, while ∂H/∂B m contains spin, orbital, and gauge contributions. The commutator terms in eq 11 originate only from the PSO and gauge correction operators. One can observe that these expressions contain terms that scale like R −2 , apparently violating the PDA, for which the terms with the slowest decay should scale like R −3 . This problem can be resolved by noticing that the commutators with H 0 in eq 11 lead to a cancellation with the energy difference in the denominator of the sum-overstates term (third term) of eq 3. The resulting R −2 term exactly cancels the corresponding term originating from eq 12, which has the same size but opposite sign. Such a cancellation of slowly decaying terms was already observed in the pioneering work by McConnell. 38 A similar cancellation of R −2 terms that follows the same kind of reasoning was observed by Helgaker and co-workers in the long-distance limit of indirect nuclear spin−spin coupling constants. 41 As a result, eqs 11 and 12 inserted into eq 3 exactly reproduce the "semiempirical" eq 1 with the susceptibility tensor defined via the equation This is a generalization of the well-known van Vleck equation first presented by Gerloch and McMeeking 42 including a diamagnetic term that was missing in their treatment; also see the discussion by Van den Heuvel and Soncini. 23 In summary, we have demonstrated that the use of rigorous first-principles quantum mechanics confirms eq 1, which was originally derived in a semiempirical way. Our derivation therefore demonstrates that eq 1 has a solid theoretical foundation. This is the main result of the current Letter. The alternative expression based on the unsymmetric tensor χ′ (eq 8) is flawed because of the neglect of contributions to the hyperfine field beyond the SD one, which, as our derivation demonstrates, is not justified. It is important to stress that we derived our result in a very general way using only relatively weak assumptions. Notably, we did not choose a particular electronic structure method. The only assumptions that went into our derivation are the validity of the equation for the chemical shielding given by Van den Heuvel and Soncini (eq 3), our particular choice of electronic Hamiltonian, and the PDA. We expect that it will be possible to derive eq 1 also for the more exact Dirac Hamiltonian or even without resorting to an explicit choice of Hamiltonian at all. This will be discussed in more detail in a forthcoming publication. The validity of eq 1 that we just established, together with the SH expressions eq 5 and eq 7, suggests that in the longdistance limit one can approximate the HFC tensor as The validity of this equation is well-established for the spin-only contributions to A and g. 26 It was shown that within the LS coupling approximation, the orbital contributions to A and g also fulfill eq 14. 26 Furthermore, a nonrigorous motivation for the validity of eq 14 for the full A and g tensors was given by Autschbach et al. 43 The FC contribution to A vanishes in the long-distance limit, where there should be no spin density at the position of the nucleus. At the level of second-order degenerate perturbation theory (DPT2), the remaining contributions to g and A are 36,44−46 The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter and 18 are the well-known orbital contributions to g and A, which make sizable contributions for transition metals, where SOC cannot be neglected. Equations 15 and 17 are the less well-known contributions due to gauge correction terms in the molecular Hamiltonian. While the gauge contribution to g has been investigated in a number of cases, 47−53 the gauge contribution to A has only rarely been considered in the existing literature. 54,55 Similar to eq 11, one can express the long-range limit of the PSO operator as Inserting this into eq 18 and noticing as before that the commutators with H 0 nonrel can be used to cancel the energy difference in the denominator, one obtains The second term (the expectation value of the commutator) prevents A PSO/SOC and g orb/SOC from fulfilling eq 14. It can also be seen that one of the terms in the commutator scales again as R −2 , in contrast to the expected R −3 . However, when considering the long-range limit of the gauge contribution (eq 17), one notices that the same erroneous terms appear with opposite sign. Hence, these terms exactly cancel when the two contributions are added: ( ) This concludes the proof of eq 14 for the complete A and g. A more in-depth theoretical discussion of the results derived in this Letter will be presented in a forthcoming publication. In order to corroborate our theoretical results, we investigated eq 21 numerically. We performed calculations on the radicals H 2 O + , CO + , H 2 CO + , and NF 3 + and the transition metal complexes Cu(NO 3 ) 2 , Ni(CO) 3 H, and TiF 3 . A proton serving as a probe nucleus was placed at different distances from the paramagnetic center. A "naked" proton was used for radicals and an HF molecule for the metal complexes. The proton A and g tensors were then calculated with Kohn− Sham density functional theory (KS-DFT) using the BP86 exchange−correlation functional and correlation-consistent basis sets of different sizes. GGA functionals are not expected to yield extremely accurate SH parameters compared with the experiment. However, absolute accuracy is not our concern here, since our comparison is not with the experiment but entirely "internal", i.e., we compare different quantities calculated at the same level of theory. In order to provide a measure of the "size" of the calculated tensors and their individual contributions, we calculate the Frobenius norm defined as ∥ ∥ = ∑ X X ij ij 2 . This allows for the definition of a "relative deviation" Δ rel as a measure of the difference between some exact tensor and an approximation like the PDA: When we compare the exact A PSO/SOC with g orb/SOC inserted into eq 14, we denote the relative difference as Δ rel SOC ; when we compare the two sides of eq 21, we denote it as Δ rel SOC+gauge . Among all of the investigated systems, the gauge correction was found to be most important for the CO + radical, for which we present results in Figure 1. We expect that with larger distance, where the PDA (eq 14) becomes better, the relative deviation Δ rel should get smaller. This is also what we observe for Δ rel SOC+gauge according to eq 21, but only for the largest basis set (decontracted cc-pV6Z). For smaller basis sets, strong deviations from this expected behavior and slow basis set convergence can be observed, although for some systems, like TiF 3 , the convergence is much faster (see the Supporting Information). This suggests that eq 14 is only a valid longdistance approximation in the limit of a complete basis set, an issue that we will discuss further in a forthcoming publication. We note that the strong basis set dependence is solely due to the slow convergence of the A tensor, while the g tensor converges very quickly when the basis set size is increased. The observed basis set dependence is in contrast to the SD contribution to the A tensor, for which eq 14 holds irrespective of the basis set size. In contrast to Δ rel SOC+gauge , Δ rel SOC increases with increasing distance for all basis set sizes. This indicates that eq 14 is not valid when only the orbital contribution is considered without the gauge contribution, as expected from eq 20. The presence of the R −2 scaling in the unphysical second term of this equation means that it will increasingly dominate A PSO/SOC , which explains why the PDA for this contribution alone becomes increasingly worse.
We now investigate how A PSO/SOC , A gauge , and their sum scale with the distance R for the CO + radical. According to our theoretical analysis, both A PSO/SOC and A gauge have terms that scale as R −2 , and these terms cancel each other out in the sum to give the physical R −3 scaling. In Figure 2 we show the norms of these tensors as functions of the distance in a log−log plot, where the slope corresponds to the exponent of R. In order to guide the eye and emphasize deviations from the ideal R −3 scaling behavior, we have also plotted in red color a straight line with slope −3 that passes through the last data point. For A PSO/SOC , the data points proceed roughly parallel to the red line for small distances, indicating an approximately R −3 scaling. However, for larger distances, the slope is reduced, which indicates that the R −2 term that decays more slowly is gaining importance. For A gauge , the slope is already significantly different from −3 for small distances, indicating that here the R −2 term becomes important earlier. For A PSO/SOC + A gauge , the slope is very close to −3 even at large distances, which confirms eq 21 that suggests that the unphysical terms cancel each other in the sum. The slopes of the different HFC tensor contributions in the beginning (small distance) and end (large The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter distance) of the data series are summarized for all of the investigated systems in Table 1. It can be seen that CO + is the only investigated system for which the R −2 term changes the scaling of A PSO/SOC in the long-distance limit. For A gauge , the R −2 term shows up in several systems, including the metal complex Ni(CO) 3 H. In order to gain an idea of the importance of the gauge contribution, we plot the ratio of the norms of A gauge and A PSO/SOC + A gauge in Figure 3. One can observe that the gauge contribution is more important for the radicals than for the transition metal complexes. One can also see that in cases where the gauge contribution has a sizable R −2 part (CO + , H 2 CO + , and Ni(CO) 3 H; see Table 1), its importance increases with increasing distance. These results indicate that while A PSO/SOC has an unphysical part that scales as R −2 that will eventually dominate at large enough distances, the prefactor of this part might be so small that it is negligible at all distances of interest. However, the CO + example shows that the gauge contribution can be essential for a correct description even for small distances, which is why one should never neglect this contribution. Although the gauge contribution seems to be less important for metal complexes, it is definitely needed for the Ni(CO) 3 H complex from our test set, as demonstrated in Figure 4 (the analogue of Figure 1 for CO + ).
Our results have important implications for the ab initio treatment of PCSs. In principle, one can evaluate the chemical shieldings directly from eq 3 (as already demonstrated by Gendron et al. 56 ) or more commonly from the SH approximation (eq 5). Although such an approach would have the advantage of going beyond the PDA, the latter case has the disadvantage that a separate set of response equations has to be solved in order to obtain the A tensor for each individual nucleus. For systems with many nuclei, like proteins, the computational cost would quickly explode. This makes the use of the PDA attractive, as only a single electronic structure calculation is necessary to determine the PCSs for all of the nuclei in the system. In the context of the SH approximation, we have shown that eq 8 with the asymmetric χ′ tensor is fundamentally flawed because of the neglect of orbital and gauge correction contributions to A. The incorporation of these contributions shows that the correct equation to use is eq 1, which confirms the classic semiempirical treatments of PCSs. Importantly, this equation is not restricted to the SH approximation. The susceptibility tensor can also be directly calculated according to eq 13 using either ab initio or, for example, ligand field models. Such an approach has been followed already by some researchers. 4,16,17,57 We expect that this can lead to more accurate predictions than the equations involving the g tensor, especially in situations where the conditions for applying the SH approximation are not fulfilled.

■ COMPUTATIONAL DETAILS
All of the calculations were performed with a development version of the ORCA electronic structure program. 58 A and g The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter tensors were calculated using KS-DFT 50,59 with the BP86 60,61 exchange−correlation functional, grid5 integration grids, and without the resolution of the identity (RI) approximation. A number of correlation-consistent basis sets of increasing accuracy 62−64 were used for the calculations. For the radicals, we used cc-pVDZ, cc-pVQZ, cc-pV6Z, and for even higher accuracy, the decontracted cc-pV6Z basis set. For transition metals, cc-pV6Z is not defined, and we used cc-pV5Z and decontracted cc-pV5Z instead. The integrals needed for the gauge correction to the A tensor (see eq 17) diverge if nuclei K and A are identical. Therefore, this contribution was omitted, as is common practice. 54,65 For computational tests of the PDA behavior, we placed probe nuclei (protons) at different distances from the center of the paramagnetic molecule (R O ) along the arbitrarily chosen unit vector − (1, 7, 4)/ 66 . We chose distances between 5 and 45 Å in steps of 2 Å. The centers were defined as O for H 2 O + , the midpoint between C and O for CO + , C for H 2 CO + , N for NF 3 + , and the metal atom for the transition metal complexes. For the radicals, bare protons were used as probe nuclei. For the metal complexes, this was not possible because of their smaller ionization energy. Therefore, the probe nucleus was chosen to be the proton of a HF molecule. The geometries were taken from an earlier g value study. 50 More details can be found in the Supporting Information.
Detailed expression for the QC Hamiltonian and its derivatives, geometries and calculation details, and plots for the remaining molecules (PDF) Raw data for all calculated contributions to the SH parameters (ZIP)