Spatial Contributions to Nuclear Magnetic Shieldings

We develop a methodology for calculating, analyzing, and visualizing nuclear magnetic shielding densities which are calculated from the current density via the Biot–Savart relation. Atomic contributions to nuclear magnetic shielding constants can be estimated within our framework with a Becke partitioning scheme. The new features have been implemented in the GIMIC program and are applied in this work to the study of the 1H and 13C nuclear magnetic shieldings in benzene (C6H6) and cyclobutadiene (C4H4). The new methodology allows a visual inspection of the spatial origins of the positive (shielding) and negative (deshielding) contributions to the nuclear magnetic shielding constant of a single nucleus, something which has not been hitherto easily accomplished. Analysis of the shielding densities shows that diatropic and paratropic current-density fluxes yield both shielding and deshielding contributions, as the shielding or deshielding is determined by the direction of the current-density flux with respect to the studied nucleus instead of the tropicity. Becke partitioning of the magnetic shieldings shows that the magnetic shielding contributions mainly originate from the studied atom and its nearest neighbors, confirming the localized character of nuclear magnetic shieldings.


INTRODUCTION
Second-order magnetic properties such as nuclear magnetic shieldings, indirect spin−spin coupling constants, and magnetizabilities are usually calculated by using the gradient theory of electronic structure calculations as the second derivative of the electronic energy with respect to the external magnetic perturbation(s) in the limit of vanishing perturbation(s). 1−4 However, the elements of the nuclear magnetic shielding and magnetizability tensors can also be obtained as second derivatives of the magnetic interaction energy, which can be written as an integral over the scalar product of a current density caused by a magnetic perturbation and the vector potential of the second magnetic perturbation. 5−7 The current density J B (r) induced by an external magnetic field Bor the current density J r ( ) m I induced by the nuclear magnetic moment m I of nucleus Iis formally defined as the real part ( ) of the mechanical momentum density where p = −i∇ is the momentum operator and Ψ(r) is a complex wave function because of the vector potential of the magnetic perturbation A B (r) or A r ( ) m I . As will be discussed in later in this work, the magnetic properties evaluated within this scheme have no reference to the magnetic gauge origin if the current density is gauge origin independent, as is the case in our GIMIC approach 8−11 as well as in the ipsocentric approach. 12−16 While the end results of the gradient-theory and the integration approaches are the same, the method based on integration can be used for providing additional information about orbital and spatial contributions to a given magnetic property. For instance, magnetizabilities, which are usually calculated using gradient theory as the second derivative of the electronic energy with respect to the external magnetic field, can also be obtained as the second derivative of the magnetic interaction energy expressed in terms of the current density induced by the magnetic field 17−19

B B
A r J r r As we have recently discussed in ref 19, eq 2 can be used to extract information about the spatial contributions to components of the magnetizability tensor.
The nuclear magnetic shielding tensor for nucleus I, in turn, is determined by the second derivative of the magnetic interaction energy with respect to the external magnetic field B and the nuclear dipole moment m I . The shielding tensor can then be calculated from the current density induced by the external magnetic field J B (r) and the vector potential of the nuclear magnetic moment A r ( ) Alternatively, the shielding tensor can be calculated from the current density induced by the nuclear magnetic moment and the vector potential of the external magnetic field A B (r) and the current density induced by the nuclear magnetic moment J r ( ) Equation 3 is typically used in computations since picking the expression with J B (r) means that the current density has to be computed only for the three components of the external magnetic field instead of the 3N components of the magnetic dipole moments of N nuclei. However, efficient algorithms have also been developed by using eq 4, as the localized nature of the current densities induced by nuclear magnetic moments allows for powerful use of screening and parallelization. 21, 22 Because the current density J B (r) induced by an external magnetic field is a function of the strength of the external magnetic field, differentiation of the magnetic interaction energy yields the first derivative of the current density with respect to the external magnetic field (∂J B (r)/∂B), which is the current-density susceptibility tensor (CDT) induced by the external magnetic field. 18,23 Analogously, the differentiation with respect to the nuclear magnetic moment acts only on the vector potential of the nuclear magnetic moment, yielding A r m ( )/ I m I ∂ ∂ , since only that term in eq 3 depends on the nuclear magnetic moment. The dot product of these two quantities, ∂J B (r)/∂B and A r m ( )/ I m I ∂ ∂ , is a scalar function known as the nuclear magnetic shielding density. 5−7 The spatial distribution of the shielding density provides detailed information about the origin of the individual elements of the nuclear magnetic shielding tensor as well as the shielding constants. 15,24−28 Further information about the magnetic shielding density can be obtained from the individual orbital contributions to the magnetic shieldings 15 and shielding functions. Dividing the magnetic shielding density into positive and negative parts as well as into orbital contributions shows the spatial origins of the shielding and deshielding contributions to the shielding tensor and the isotropic shielding constants. 27,29 Thus, calculations of magnetic shielding densities provide a rigorous physical basis for interpreting nuclear magnetic resonance (NMR) chemical shifts.
In this work, we develop a methodology for analyzing spatial contributions to nuclear magnetic shielding constants. We apply the methods to the hydrogen and carbon nuclei in benzene (C 6 H 6 ) and cyclobutadiene (C 4 H 4 ), which are test cases representing aromatic and antiaromatic hydrocarbons, respectively. Next, we present the underlying theory in section 2.1 and continue in section 2.2 with the employed numerical methods. Then, in section 2.3, we describe the computational methods. We discuss the magnetic shielding densities of the studied molecules in section 3 and summarize our study and form our main conclusions in section 4.
where R I is the position of the Ith nucleus and μ 0 is the vacuum permeability. 30 Similarly, the vector potential A B (r) of an external static magnetic field is where R O is the chosen magnetic gauge origin. The magnetic flux density B and the magnetic dipole moment m I are uniquely defined by the vector potentials A B (r) and A r ( ) m I , whereas the reverse does not hold since all the vector potentials of the form A′ = A + ∇f(r) generate the same magnetic field B, as ∇ × ∇f(r) = 0 for any smooth function f(r).
Even though exact solutions of the Schrodinger equation are gauge invariant, the use of finite one-particle basis sets introduces a gauge dependence in quantum chemical calculations of magnetic properties. The CDT can be made gauge origin independent by using gauge-including atomic orbitals (GIAOs), also called London atomic orbitals (LAOs). The GIAOs are defined as 8,31,32 where i is the imaginary unit and χ μ (0) (r) is a standard Gaussian-type basis function centered at R μ . The use of GIAOs eliminates the gauge origin from the expression we use for calculating the CDT (∂J α In eq 8, D is the density matrix in the atomic orbital basis, ∂D/∂B are the magnetically perturbed density matrices, ϵ αβγ is the Levi−Civita symbol, h(r) denotes the magnetic interaction operator without the |r − R I | −3 denominator with h r m r R p ( ) ( ) and h r m B r R r R 1 r R r R ( ) The Journal of Physical Chemistry A pubs.acs.org/JPCA Article and R I is the position of nucleus I. Finally, the nuclear magnetic shielding tensor of nucleus I, σ αβ I , can be calculated from eqs 3 and 5 as (11) It is important to note that all terms that contain the gauge origin R O cancel in eq 8, making the CDT calculation as well as eq 11 independent of the gauge origin. Analogously, all terms in eq 8 containing the nuclear position R I also cancel, eliminating explicit references to the coordinates of the nucleus I from the current density (the physical implicit dependence still remains). As a result, the integrated second-order magnetic properties have no reference to the gauge origin or the nuclear coordinates.
The Biot−Savart expression in eq 11 has advantages over the corresponding second-derivative expression. Contributions to the tensor elements can be visually interpreted by plotting the positive and negative parts of the integrand separately, yielding information about shielding and deshielding contributions to the elements of the magnetic shielding tensor. For example, in a system with a ring current, the σ zz I contribution given by will consist of both positive and negative shielding contributions due to the relative direction of the current density with respect to the investigated atom I. 11,33 The Biot−Savart expression in eq 11 can be calculated by quadrature when the CDT is known. Because established gradient-theory implementations of NMR shielding constants are typically used to compute the CDT, the shielding constants from eq 11 do not provide any new physical information; however, the numerically evaluated shielding constants can be compared to the analytically evaluated values to assess the accuracy of the numerical integration of the Biot−Savart expressions, which is useful for applications to other secondorder magnetic properties. For instance, a similar approach has recently been used to calculate and assess the accuracy of magnetizabilities from new density functional approximations, even though analytical methods to calculate the magnetizability tensor were not available in the used program. 19 2.2. Implementation. A numerical integration scheme for calculating spatial contributions to nuclear magnetic shieldings has been implemented into the freely available GIMIC program. 34 The atomic contributions to the magnetic shieldings are obtained by quadrature over atomic domains generated by the NUMGRID library, 35 which is based on the use of Becke's multicenter scheme. 36 The atomic domains were determined with the Becke partitioning scheme, 36 employing the iteration order k = 3 in the construction of the cutoff function as suggested by Becke. The radial integration points of the atom-centered grids are generated as suggested by Lindh et al., 37 and Lebedev's angular grids are used. 38 The CDT is constructed in GIMIC with eq 8 from the density matrix, the magnetically perturbed density matrices, and basis set information obtained from Turbomole 39 calculations of NMR shielding constants.
2.3. Computational Methods. The molecular structures of C 6 H 6 , C 4 H 4 , and B 3 N 3 H 6 were optimized with Turbomole 39 version 7.5 employing the B3LYP density functional, 40−42 the def2-TZVP basis set, 43 and the m5 quadrature grid; 44,45 the optimized molecular structures are given in the Supporting Information. Nuclear magnetic resonance (NMR) shielding constants were also calculated with Turbomole at the same level of theory by using GIAOs. 31,32,46,47 In the NUMGRID calculations, 21042 grid points were used for each carbon and 19234 grid points for each hydrogen.
The B3LYP/def2-TZVP level of theory has been found to yield good agreement compared to second-order Møller− Plesset (MP2) theory for the 1 H NMR magnetic shielding in tetramethylsilane (TMS, Si(CH 3 ) 4 ), as the 13 C shielding in TMS reproduced by the method deviates by only 7% (roughly 12 ppm) from the one obtained at the MP2/def2-TZVP and MP2/def2-TZVPP levels of theory. 48 Although we are aware that these results are not fully converged to the complete basis set limit, especially for the carbon shieldings, 49 the B3LYP/ def2-TZVP level of theory suffices for our present purposes of illustrating the spatial origins of magnetic shieldings: the accurate reproduction of 13 C shieldings is known to be challenging, 50 and the functional error is likely of the same order of magnitude as the basis set truncation error.
The methods presented in section 2.2 and their GIMIC implementation, however, can be applied in combination with any basis set or level of theory for which the density and perturbed density matrices are available. Basis set truncation errors for the def2-TZVP shieldings and their effects on the atomic contributions will be discussed in section 3.3, showing that the truncation errors in def2-TZVP only affect the contribution to the shielding of the same atom, whereas the contributions to the shieldings of the other atoms are reproduced accurately in the def2-TZVP basis set.

RESULTS AND DISCUSSION
3.1. Benzene. The magnetic shielding density for the 1 H NMR shielding in Figure 1a shows that the main shielding contribution in the molecular plane originates from the outer regions of the molecular electron density, where the diatropic ring current is strong. 11,33 Deshielding contributions arise close to the hydrogen nucleus and close to its adjacent (ipso) carbon. Shielding and deshielding contributions also arise from the valence electrons of the ipso and the nearest-neighbor (ortho) carbon atoms due to their local atomic current-density fluxes. All carbons have both shielding and deshielding contributions The Journal of Physical Chemistry A pubs.acs.org/JPCA Article for 1 H arising from the core electrons due to atomic current densities around the nucleus. The two core contributions cancel almost completely because the atomic current density has the same strength on both sides of the nucleus, and the relative distance to the positive (shielding, blue) and negative (deshielding, red) areas from the studied hydrogen nucleus is almost the same for the carbon atoms in the meta and para positions.
The zz contribution to the magnetic shielding density in the molecular plane for a 13 C nucleus in Figure 1b has an onionlike shell structure of shielding and deshielding contributions. The shielding contribution close to the nucleus arises from the core electrons, whereas the valence electrons deshield the nucleus. In the next shell, the shielding contribution originates from the diatropic ring current that flows on the outer side of the molecular ring near the hydrogen as well as from the paratropic ring current inside the C 6 H 6 ring. The atomic current density in the valence orbitals of the ortho carbon atoms also contributes to the 13 C shielding on closer side of the ortho carbon, while the contributions are deshielding on the remote side.
The ring-current contribution to the nuclear magnetic shielding constants can be analyzed by plotting the spatial distribution of the σ zz component to the nuclear magnetic shielding density. The zz component of the 1 H NMR shielding density calculated in a plane 1 a 0 above the molecular plane is shown in Figure 2a. The diatropic ring current flowing on the outside of the hydrogen shields the hydrogen nucleus. The ring current on the other side of the ring also shields it, while the diatropic ring current flowing on the inside of the hydrogen is deshielding. The paratropic ring current inside the C 6 H 6 ring deshields the hydrogen nucleus on the remote half of the ring, whereas inside the ipso carbon atom the paratropic ring current shields the hydrogen nucleus. The sign of the shielding contributions depends on the direction of the current density with respect to the studied nucleus according to the Biot− Savart expression in eq 12.
The ring-current contribution to the 13 C NMR shielding is seen in Figure 2b, where shielding contributions appear along the outer perimeter of the carbon ring. The paratropic ring current inside the C 6 H 6 ring leads to a shielding contribution near the studied carbon atom, whereas it is deshielding on the remote interior part of the ring. The deshielding contribution in the vicinity of the studied carbon originates from the diatropic ring current passing on the inside of the carbon atom.
The absolute value of the nuclear magnetic shielding density is illustrated by using a contour surface in Figure 3, where blue represents the shielding density of the hydrogen atom, while yellow is used to illustrate the shielding of the carbon atom. Figure 3 reveals that the shielding density near the ortho atoms contributes significantly, whereas the more distant atoms have negligible contributions, as expected due to the |r − R I | −3 denominator in the vector potential of m I .
Atomic contributions to the isotropic nuclear magnetic shielding constants can be analyzed by integration over atomic subdomains, yielding a compact representation of the spatial distribution of the shielding density. In this work, the atomic subdomains are defined by the Becke partitioning, 36 as discussed in section 2.2. Even though Becke partitioning was originally aimed for efficient numerical integration of density functionals, it has been shown to be useful for e.g. constructing mathematically well-based Pipek−Mezey orbital localization techniques, 51,52 and with a careful choice of the partitioning function it yields chemically sound atomic charges and bond orders. 53 The decomposition depends on the partitioning, i.e., the choice for the atomic weight functions. The original Becke partitioning yields a rough idea of the atomic decomposition of the shielding density; more sophisticated atomic decompositions are left to further work.
The resulting atomic contributions to the 1 H NMR and 13 C NMR magnetic shieldings of C 6 H 6 are given in Tables 1 and 2, respectively. These data suggest that the main contributions to the shielding originate from the vicinity of the studied atom   The contribution to the 1 H NMR shielding from the atomic domain of the studied hydrogen is 78.13% of the total shielding, while the contribution assigned to each ipso carbon is 6.11%. Contributions from all other atoms are in the interval of [0.61, 2.64]%. The contribution to the 13 C NMR shielding from the studied carbon is 70.49%. The ipso hydrogen and ortho carbons contribute with 5.61% and 7.00%, respectively, whereas the 13 C NMR contributions from the rest of the atoms are in the interval of [0.49, 2.08]%.
As a side note, although the molecular structure of borazine (B 3 N 3 H 6 ) is similar to that of benzene, a previous study of shielding densities suggested that borazine is nonaromatic. 15 However, a follow-up study showed that B 3 N 3 H 6 does sustain a diatropic ring current, although its strength is only 25% of that in C 6 H 6 . 54 A comparison of the zz contribution to the shielding densities of C 6 H 6 and B 3 N 3 H 6 (shown in the Supporting Information) reveals that B 3 N 3 H 6 has a similar but weaker ring-current contribution to the shielding density as for C 6 H 6 . Thus, B 3 N 3 H 6 cannot be considered to be nonaromatic.
3.2. Cyclobutadiene. The zz contribution to the 1 H NMR shielding density in the molecular plane of C 4 H 4 shown in Figure 4a is similar to the one for C 6 H 6 in Figure 1a. Even though C 4 H 4 is antiaromatic, it sustains a diatropic ring current along the outer edge of the molecule outside the hydrogen giving rise to a similar shielding contribution outside the hydrogen like in C 6 H 6 . 33 The ring current is paratropic inside the ring as in C 6 H 6 . Deshielding contributions appear at the hydrogen nucleus as well as at the ipso and ortho carbons due to local current densities. The atomic current density in the core of the carbon atoms leads to shielding and deshielding contributions that practically cancel, as for C 6 H 6 .
The contributions to the 13 C magnetic shielding density in the molecular plane of C 4 H 4 in Figure 4b also remind of those for C 6 H 6 . The onion structure of the alternating shielding and deshielding contributions around the studied carbon atom originates from current densities with different flux directions in the vicinity of the atom. The diatropic atomic current density in its core orbitals, the diatropic ring current flowing on the outside of hydrogen atom, and the paratropic ring current inside the C 4 H 4 ring shield the carbon nucleus, whereas the atomic current density of the valence orbitals deshields it.
In contrast, the magnetic shielding density in a plane 1 a 0 above (or below) the molecular plane of C 4 H 4 differs completely from the one for C 6 H 6 because C 6 H 6 sustains a diatropic ring current in the π orbitals, while the current density of C 4 H 4 is paratropic there. The 1 H and 13 C magnetic shielding densities of C 4 H 4 in Figures 5a and 5b show that the diatropic ring current along the outer edge of the molecule leads to a shielding contribution to 1 H NMR and 13 C NMR shieldings. The strong paratropic ring current which resides mainly inside the molecular ring leads to a shielding contribution to 1 H NMR in the closer half of the ring and a deshielding contribution from the remote part of the ring due to the different directions of the current-density fluxes relative to the studied hydrogen nucleus.
The deshielding contribution to the 13 C NMR shielding from the paratropic ring current dominates above the ring on the inside of it. A small shielding area is seen in Figure 5b, where the relative direction of the paratropic ring current leads to magnetic shielding. The paratropic ring current on the outside of the studied carbon deshields the carbon nucleus. The diatropic ring current along the outer edge of the molecule results in a weak shielding contribution in the vicinity of the hydrogen atom.
The absolute value of the nuclear magnetic shielding density is illustrated by using a contour surface in Figure 6, again showing that the most significant contributions arise from the ipso atoms, with some contributions from the ortho atoms.
The atomic contributions to the 1 H NMR and 13 C NMR magnetic shieldings of C 4 H 4 are given in Tables 3 and 4, respectively. Table 3 shows that 70.98% of the 1 H NMR shielding of C 4 H 4 originates from the atomic domain of the   The Journal of Physical Chemistry A pubs.acs.org/JPCA Article studied hydrogen. The contribution from the ipso carbon is 24.98%. The rest of the atoms contribute with less than 3.66%. The contributions from the para carbon and the ortho carbon with a formal single bond to the studied carbon are even negative. The contribution to the 13 C NMR shielding from the studied carbon is 87.50%. The ipso hydrogen contributes with 6.69%, and the ipso carbon with a formal double bond to the studied carbon contributes with 7.71%. Contributions to 13 C NMR from the rest of the atoms are small. The contributions from the ortho carbon with a formal single bond to the studied carbon and the carbon in the para position are also in this case negative.
3.3. Basis Set Dependence. We investigated the basis set truncation error in the def2-TZVP basis set with additional calculations using the fully uncontracted pc-n (unpc-n) polarization consistent basis sets series 55 and their augmented versions. 56 The basis set study was performed with Gaussian,57 and all basis sets were obtained from the Basis Set Exchange. 58 The full set of results is shown in the Supporting Information.
The resulting B3LYP complete basis set estimates from the quintuple-ζ unpc-4 set, which has a 11s6p3d2f1g and 18s11p6d3f2g1h composition for H and C, respectively, were found to be 42.35 and 24.03 ppm for the 13  Because of the noticeable basis set truncation error for the carbon shieldings, additional calculations were performed with the generally contracted pc-n basis sets, 55 their newer versions based on segmented contractions 59 (pcseg-n) and specializations thereof to the reproduction of nuclear magnetic shieldings 60 (pcSseg-n), as well as with the Karlsruhe def2 family of basis sets. 43 The pcseg-3 basis set 59 was found to yield excellent agreement with the unpc-4 values: the pcseg-3 basis set yields 13  Because of the good accuracy of the pcseg-3 basis set, spatial decompositions for C 4 H 4 and C 6 H 6 were recomputed in this basis; the decompositions are shown in the Supporting Information. Comparison of these data to the values in Tables  2 and 4 shows that basis set truncation error in def2-TZVP significantly affects only the shielding contribution from the ipso carbon, while the shielding contributions from the other atoms are strikingly similar, differing only up to 0.05 ppm for C 4 H 4 and 0.03 ppm for C 6 H 6 . This strongly suggests that the differences originate from orbitals localized to the ipso carbon, that is, an insufficient flexibility in the semicore region of the def2-TZVP basis set of carbon. Because the truncation error changes significantly the absolute nuclear magnetic shielding of the studied carbon, this also affects the relative percentages of the atomic contributions. Similar conclusions can also be made for the hydrogen shieldings by comparison of the data in the Supporting Information to Tables 1 and 3: the largest change (0.27 ppm) originates from the ipso hydrogen, while the contributions from all other atoms are negligible: less than 0.05 ppm for C 6 H 6 and less than 0.03 ppm for C 4 H 4 .

SUMMARY AND CONCLUSIONS
We have implemented methods for calculating and visualizing nuclear magnetic shielding densities in the GIMIC program. Studies of the shielding densities of benzene (C 6 H 6 ) and cyclobutadiene (C 4 H 4 ) show that the direction of the currentdensity flux relative to the studied nucleus determines whether the current density shields or deshields the nuclear magnetic moment. The paratropic ring current in the molecular plane within the C 6 H 6 and C 4 H 4 rings shields the studied nucleus when the current flows in the vicinity of the nucleus, while the current becomes deshielding on the remote side of the ring. The paratropic ring current inside the ring is much weaker in the aromatic benzene molecule than in the antiaromatic   The Journal of Physical Chemistry A pubs.acs.org/JPCA Article cyclobutadiene molecule. Benzene sustains a strong diatropic ring current in the π orbitals above and below the molecular ring, which results in shielding contributions to the 1 H NMR and 13 C NMR shieldings. However, the ring current passing the ipso carbon deshields the 1 H nuclear magnetic moment because it is a diatropic ring current near the studied 1 H nucleus that flows on the inside of it. The same holds for the 13 C NMR shielding. However, the diatropic ring current passing on the inside of the carbon is weaker and leads only to a small deshielding contribution.
The 1 H NMR and 13 C NMR shielding densities in the molecular plane of C 4 H 4 are similar to the ones of C 6 H 6 , whereas 1 a 0 from the molecular plane the shielding densities are completely different. C 4 H 4 sustains a strong paratropic ring current in the π orbitals inside the ring, whereas the ring current in C 6 H 6 is diatropic and flows mainly on the outside of the carbon ring.
Calculations of atomic contributions to the nuclear magnetic shielding constants using Becke's partitioning show that the largest contributions originate from the ipso atoms and its nearest neighbors. The ipso carbon contributes with 70.49% and 87.50% to the 13 C NMR shielding of C 6 H 6 and C 4 H 4 , respectively. The contribution from the ipso hydrogen to the 1 H NMR shielding is 78.13% and 70.98% for C 6 H 6 and C 4 H 4 , respectively. Even for small molecules like C 4 H 4 and C 6 H 6 , contributions from more distant atoms are only a few percent, which is utilized in local methods to calculate nuclear magnetic shielding constants.
Although the B3LYP/def2-TZVP level of theory was used for the most part of the present work, the methods presented herein can also be used with larger basis sets and post-Hartree−Fock levels of theory. We repeated the analysis in the pcseg-3 basis set, which we found to yield shielding constants in good agreement with our complete basis set estimates, which showed that most of the deficiencies in the def2-TZVP data originate from the atom under study, while the contributions from all other atomic domains are essentially already converged in def2-TZVP.
B3LYP/def2-TZVP optimized molecular structures for C 4 H 4 , C 6 H 6 , and B 3 N 3 H 6 , atomic shielding contributions calculated by using the pcseg-3 basis set, a basis set study, and the shielding density of borazine (PDF)