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.

In this Letter, we study a paramagnetic center and a sufficiently far nucleus with magnetic moment MK. 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 (ℏ = me = 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 = χ −1/2(χ + χ) and Δχrh = χ – χ 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. Equation was described by McConnell and Robertson[21] and by Kurland and McGarvey[22] in the spin Hamiltonian (SH) approximation.

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 = −kBT ln Z (where Z = ∑e–β):The different operators that appear in this expression are derivatives of the full electronic Hamiltonian H = H(B, MK) around B = 0 and MK = 0. The summations in eq  are over eigenstates |Ψ(0)⟩ (with eigenvalues E(0)) of the field-free Hamiltonian H0 = H(B = 0, MK = 0). The indices n, m, ... label energy levels, and μ, ν, ... are additional labels for degenerate states with the same energy. Z0 is the zeroth-order canonical partition function, and β = 1/kBT.

If all relevant states can be described by a SH of the formeq  turns into[20,24−27]where the “spin dyadic” ⟨SST⟩ is a matrix that depends on the zero-field splitting tensor D and the temperature:The zeroth-order eigenfunctions and energies in eq  are derived from the zeroth-order Hamiltonian in the SH approximation, Hspin(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 ,[28,30−33] all of the contributions to the hyperfine coupling (HFC) tensor A except for the Fermi contact (FC) and spin-dipolar (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  by χ = χ′·gT/ge. The appearance of this unsymmetrical tensor was considered a necessary consequence of the more rigorous QC treatment compared with the “semiempirical” eq .[26,28,30] The two approaches are ostensibly different.[26] Several recent papers prefer the use of eq  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 . Therefore, as the use of eq  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  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 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 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]RO is a position in the region where the paramagnetic center is located, and RK is the position of the nucleus. We use RO also as the gauge origin for the homogeneous external magnetic field. Furthermore, we define r = r – RO and r = r – RK, where r 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 can be assumed to be small compared with R = RO – RK if the latter becomes large enough, i.e., r = r + R ≈ R. Following McConnell, inverse powers of r can then be expanded in a Taylor series around the point r = 0, which gives[38]Inserting these expressions into the equations for the Hamiltonian derivatives occurring in eq  that depend on the nuclear position and neglecting all terms with inverse powers of R larger than 3 (PDA), one obtainsandIn eq , ∂H/∂MK contains FC, SD, paramagnetic spin–orbit (PSO), and gauge contributions, while ∂H/∂B contains spin, orbital, and gauge contributions. The commutator terms in eq  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 H0 in eq  lead to a cancellation with the energy difference in the denominator of the sum-over-states term (third term) of eq . The resulting R–2 term exactly cancels the corresponding term originating from eq , 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 and 12 inserted into eq  exactly reproduce the “semiempirical” eq  with the susceptibility tensor defined via the equationThis 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 , which was originally derived in a semiempirical way. Our derivation therefore demonstrates that eq  has a solid theoretical foundation. This is the main result of the current Letter. The alternative expression based on the unsymmetric tensor χ′ (eq ) 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 ), our particular choice of electronic Hamiltonian, and the PDA. We expect that it will be possible to derive eq  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  that we just established, together with the SH expressions eq  and eq , suggests that in the long-distance limit one can approximate the HFC tensor asThe 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 .[26] Furthermore, a nonrigorous motivation for the validity of eq  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]where the |Ψ⟩ are the principal components (M = S) of the eigenstates of the nonrelativistic zeroth-order Hamiltonian H0nonrel = H0 – HSOC with eigenvalues E(0),nonrel. Equations 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 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 , one can express the long-range limit of the PSO operator asInserting this into eq  and noticing as before that the commutators with H0nonrel can be used to cancel the energy difference in the denominator, one obtainsThe second term (the expectation value of the commutator) prevents APSO/SOC and gorb/SOC from fulfilling eq . 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 ), 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  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  numerically. We performed calculations on the radicals H2O+, CO+, H2CO+, and NF3+ and the transition metal complexes Cu(NO3)2, Ni(CO)3H, and TiF3. 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 . 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 APSO/SOC with gorb/SOC inserted into eq , we denote the relative difference as ΔrelSOC; when we compare the two sides of eq , we denote it as ΔrelSOC+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 . We expect that with larger distance, where the PDA (eq ) becomes better, the relative deviation Δrel should get smaller. This is also what we observe for ΔrelSOC+gauge according to eq , 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 TiF3, the convergence is much faster (see the Supporting Information). This suggests that eq  is only a valid long-distance 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  holds irrespective of the basis set size. In contrast to ΔrelSOC+gauge, ΔrelSOCincreases with increasing distance for all basis set sizes. This indicates that eq  is not valid when only the orbital contribution is considered without the gauge contribution, as expected from eq . The presence of the R–2 scaling in the unphysical second term of this equation means that it will increasingly dominate APSO/SOC, which explains why the PDA for this contribution alone becomes increasingly worse.

Figure 1

Comparison of the exact A tensor contributions to an approximation calculated with eq  for the CO+ radical. The presented numbers are the relative deviations (see eq ) in %.

We now investigate how APSO/SOC, Agauge, and their sum scale with the distance R for the CO+ radical. According to our theoretical analysis, both APSO/SOC and Agauge 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 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 APSO/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 Agauge, the slope is already significantly different from −3 for small distances, indicating that here the R–2 term becomes important earlier. For APSO/SOC + Agauge, the slope is very close to −3 even at large distances, which confirms eq  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 distance) of the data series are summarized for all of the investigated systems in Table . It can be seen that CO+ is the only investigated system for which the R–2 term changes the scaling of APSO/SOC in the long-distance limit. For Agauge, the R–2 term shows up in several systems, including the metal complex Ni(CO)3H.

Figure 2

Log–log plots of the norm of the A tensor contributions versus R, the distance of the probe nucleus from the center of the CO+ radical, with the decontracted cc-pV6Z basis set. The computed data are shown with diamond symbols for only the PSO/SOC part (top), only the gauge correction part (middle), and their sum (bottom). Shown in red are straight lines with slope = −3 passing through the last data points.

Table 1

Slopes of the Log–Log Plots (as in Figure for CO+) at the Beginning and End of the Data Series for Different A Tensor Contributionsa

		first two data points (small R)	last two data points (large R)
CO+	APSO/SOC	–3.09	–2.46
	Agauge	–2.16	–2.02
	APSO/SOC + Agauge	–3.17	–3.02
H2CO+	APSO/SOC	–3.08	–3.00
	Agauge	–2.73	–2.12
	APSO/SOC + Agauge	–3.07	–3.01
H2O+	APSO/SOC	–3.02	–3.02
	Agauge	–3.01	–2.85
	APSO/SOC + Agauge	–3.02	–3.02
NF3+b	APSO/SOC	–2.98	–3.00
	Agauge	–3.00	–3.00
	APSO/SOC + Agauge	–2.98	–3.00
Cu(NO3)2	APSO/SOC	–3.02	–3.00
	Agauge	–2.82	–3.00
	APSO/SOC + Agauge	–3.02	–3.00
Ni(CO)3H	APSO/SOC	–3.05	–3.00
	Agauge	–3.01	–2.18
	APSO/SOC + Agauge	–3.05	–3.00
TiF3	APSO/SOC	–2.98	–3.00
	Agauge	–2.98	–3.00
	APSO/SOC + Agauge	–2.98	–3.00

The used basis sets were decontracted cc-pV5Z for transition metal complexes and decontracted cc-pV6Z for radicals.

For NF3+, the calculations for the first distance (5 Å) did not converge. Hence, the first two data points correspond to 7 and 9 Å.

In order to gain an idea of the importance of the gauge contribution, we plot the ratio of the norms of Agauge and APSO/SOC + Agauge in Figure . 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+, H2CO+, and Ni(CO)3H; see Table ), its importance increases with increasing distance. These results indicate that while APSO/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)3H complex from our test set, as demonstrated in Figure (the analogue of Figure for CO+).

Figure 3

Ratios of the norms of the gauge contribution to the A tensor and the sum of the PSO/SOC and gauge contributions. All of the values were calculated with the decontracted cc-pV6Z basis set for the radicals and with the decontracted cc-pV5Z basis set for the transition metal complexes.

Figure 4

Comparison of the exact A tensor contributions to an approximation calculated with eq  for the Ni(CO)3H complex. The presented numbers are the relative deviations (see eq ) in %.

Our results have important implications for the ab initio treatment of PCSs. In principle, one can evaluate the chemical shieldings directly from eq  (as already demonstrated by Gendron et al.[56]) or more commonly from the SH approximation (eq ). 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  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 , 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  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 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 ) 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 (RO) along the arbitrarily chosen unit vector . We chose distances between 5 and 45 Å in steps of 2 Å. The centers were defined as O for H2O+, the midpoint between C and O for CO+, C for H2CO+, N for NF3+, 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.