TThe Lorentz sphere visualised
S. Sturniolo a) and J.R. Yates b) Scientific Computing Department, UKRI Department of Materials, University of Oxford, Parks Road, Oxford, OX1 3PH,United Kingdom (Dated: 14 February 2019)
From the inception of Nuclear Magnetic Resonance as a spectroscopic technique, the local origin of chemicalshifts has been a topic of discussion. A useful concept employed to describe it has been that of the “Lorentzsphere”, the approximately spherical volume surrounding a given nucleus in which the electronic currentscontribute significantly to the chemical shift, whereas the outside can be considered as an uniformly magne-tised “bulk”.In this paper we use the output of the plane wave DFT code CASTEP to get a quantitative estimate of theLorentz sphere in periodic systems. We outline a mathematical description of a radial buildup function for themagnetic shielding starting from the electronic currents and the simple assumption of periodicity. We providean approximate upper bound for the Lorentz sphere’s size in any crystal, then compute buildup functions fora number of sites in two molecular crystals, showing how various chemical features such as hydrogen bondsinfluence to convergence to the final shielding value.PACS numbers: 74.25.nj, 71.15.MbKeywords: NMR, DFT, chemical shifts, Lorentz sphere, crystals
I. INTRODUCTION
Since the early days of NMR it has been vital to under-stand the way in which chemical structure gives rise tochemical shifts . Physical considerations lead one to sep-arate two main contributions: one originating from thebulk of the sample, an average effect of macroscopic mag-netisation and sample shape, and another coming fromlocal electronic currents, whose structure on the shortrange causes the observed differences between sites in thesame sample. In molecules tumbling around in a solventin the liquid state, there is a natural separation betweenintra- and intermolecular contributions for these local ef-fects, due to motional averaging of the latter. In thesolid state, however, this distinction is more blurred. Therange within which these local contributions are relevantis referred to as the “Lorentz sphere” in literature , butsuch considerations typically do not go beyond a quali-tative level. It would be very useful for the purpose ofNMR crystallography to have a better, quantitative un-derstanding of the Lorentz sphere, as it would allow usto better understand how much and in what ways eachspectral line is sensitive to the overall structure.The authors have already suggested one such approach making use of electronic currents computed by DensityFunctional Theory, in which the volumes of space con-tributing the most to the chemical shift at a given sitewere highlighted. This method was similar to the ap-proach of plotting current density maps that has beenused in the past . In this work however a differentapproach, specifically suited to crystalline solids, is pro- a) Electronic mail: [email protected] b) Electronic mail: [email protected] posed. The properties of the current in a periodic crystalare exploited to compute a radial buildup function. Aswe expand the volume of the crystal whose contributionwe account for in a sphere centred on the site, we seehow shielding tensors converge to their final value. Theanalytic form of this function gives an upper bound forthe size of the Lorentz sphere that depends only on thelattice parameters of the crystal itself. Finally, DensityFunctional Theory (DFT) and the Gauge-Including Pro-jector Augmented Wave (GIPAW) approach are usedto compute the currents induced by valence electrons intwo well known molecular crystals, ethanol and aspirin,and buildup functions are created for isotropic shielding,anisotropy and asymmetry at various sites. The shapeand speed of convergence of these curves highlights theway each site interacts with its chemical environment,confirming physical intuition and providing a strikinglysimple visual representation of the Lorentz sphere.
II. THEORY
Applying a static magnetic field to a diamagnetic sam-ple will induce microscopic electronic currents. Thesecurrents give rise to local magnetic fields which are expe-rienced by atomic spins and shift slightly their resonancefrequencies, thus generating the NMR spectrum. If weassume that the current density, J ( r (cid:48) ) is known, thenthe magnetic shielding tensor at a given position, σ ( r ),describing the relationship between external applied field B ext and local field B in , can be found simply by applyingthe Biot-Savart law : B in ( r ) = − σ ( r ) B ext = 1 c (cid:90) d r (cid:48) J ( r (cid:48) ) × r − r (cid:48) | r − r (cid:48) | (1) a r X i v : . [ phy s i c s . c o m p - ph ] F e b where the integral is carried out over the entirety ofspace. In the case of an infinite periodic system, how-ever, we can exploit the properties of the current densityto write it as a Fourier series over the reciprocal latticevectors G : B in ( r ) = 1 c (cid:90) d ρ ρ ρ × (cid:88) G J ( G ) exp( i G ρ ) exp( i Gr ) (2)where a change of coordinate ρ = r (cid:48) − r has been madefor convenience.Passing to polar coordinates, and using the plane waveexpansion : B in ( r ) = 1 c (cid:90) ρ dρd Ω ˆ ρ ρ × (3) (cid:88) G π J ( G ) exp( i Gr ) (cid:88) l,m i l j l ( Gρ ) Y l,m ( ˆ G ) Y ∗ l,m ( ˆ ρ )where the Y l,m and j l are respectively the sphericalharmonics and the spherical Bessel function with angularmomentum numbers l, m .Since the ˆ ρ versor can be entirely expressed as a functionof real spherical harmonics with l = 1:ˆ ρ = (cid:104) Y , − − Y , , i ( Y , − + Y , ) , √ Y , (cid:105) (4)(where dependency of the Ys on ( ˆ ρ ) is left implicit),when integrating over the solid angle Ω the sum over l, m can be removed as all other terms disappear, andthese are projected and replaced with the correspondingspherical harmonics dependent on ˆ G . In the end we have B in ( r ) = 1 c (cid:88) G πi ˆ G × J ( G ) exp( i Gr ) (5) (cid:90) dρj ( Gρ ) . Exploiting the properties of the Bessel functions: j ( x ) = − ∂∂x sin( x ) x (6)we find (cid:90) R dρj ( Gρ ) = − G (cid:20) sin( GR ) GR − (cid:21) (7)leading to B in ( r , R ) = 1 c (cid:88) G πi ˆ G × J ( G ) exp( i Gr ) (8)1 G (cid:20) − sin( GR ) GR (cid:21) . In the limit of R → ∞ this reduces to the known result B in ( r ) = 1 c (cid:88) G πi ˆ G G × J ( G ) exp( i Gr ) (9)which can be obtained from Eq. 1 also simply by ap-plying the convolution theorem to the individual Fouriertransforms of the two factors of the integral.Equation 8 gives us a way to compute a buildup functionfor the magnetic shielding tensor at a given site, assum-ing that the electronic current density is known. Theonly expensive operation is computing the current den-sity itself; given that, Eq. 8 is trivial to compute at anyvalue of R in just a few seconds. Full buildup functionswith good resolution can be obtained in minutes for aver-age sized systems. Even without considering the detailsof the current density for a specific system, though, wecan make some useful observations just from the form ofthe expression. The sum over all G vectors appearing inEq. 2 to 9 is a sum over all the points of the reciprocallattice, except for G = 0, since we assume that there isno net current in the crystal. This means that the con-vergence of the sum in Eq. 8 will be often dominated bythe shortest G in the lattice, for two reasons. The first isthat the G factor significantly enhances its relative mag-nitude. The second is that the Bessel function sin ( x ) /x converges quickly to zero as x goes to infinity; for x ∼ G , the faster the convergence of a specificterm. One can consider the simple case of a cubic lat-tice with parameter a and see that the argument will be x = 2 πR/a and therefore will be already within the 20%error for R ∼ a . In the more general case, Eq. 8 suggeststhat the size of the Lorentz sphere in a periodic systemwill tend to be on the scale of the inverse of the length ofthe shortest reciprocal lattice vector. This is by no meansa hard rule, as in special cases there could be construc-tive interference effects that cause the convergence to beslower. The exact size of the sphere will always dependon the level of precision for which we want to accountand the current density of the system of interest. III. COMPUTATION OF BUILDUP FUNCTIONS
During the last two decades, routine computation ofmagnetic shieldings by use of Density Functional The-ory (DFT) software has been made possible by the de-velopment of the Gauge-Including Projector AugmentedWave (GIPAW) method . This method solves the is-sue of computing shieldings when representing the coreelectrons of an atom with pseudopotentials. This is arequired step to only run the DFT calculation proper tocompute the wavefunction of the valence electrons, thusdecreasing the computational load significantly. Whenoperating within this approximation, the final computedmagnetic field can be split in three parts: B tot = B in + B susc = B loc + B lr + B susc . (10)Here B susc represents a contribution from the suscep-tibility of the sample, which is simply a function of itsshape, and B loc and B lr represent the local and longrange contributions respectively. The local contributionhere is used to mean an extremely short-ranged contribu-tion due to the core electrons and to part of the valenceelectons. This is a consequence of the way the GIPAWapproximation works, as the wavefunction near the nu-cleus is pseudised. Therefore, in the upcoming exampleswe will focus only on the long-range, which is the one thatcontains all the effects of valence electrons away from thenucleus. This means that the buildup functions shouldnot considered fully physical within the core radius foreach element, which are 0.32 ˚A for H, 0.74 ˚A for C and0.58 ˚A for O. A. Details of the calculations
Calculations were carried out using CASTEP ver-sion 16.1 . Structures for solid state ethanol andaspirin (acetylsalicylic acid) were retrieved from theCrystallography Open Database. Calculations were car-ried out using the PBE exchange-correlation functionalwith Tkatchenko-Scheffler dispersion corrections andthe default ultrasoft pseudopotentials. After convergencechecks, the cutoff energy was set to 800 eV for ethanoland 600 eV for aspirin; the k-point grid was set to 3x2x2for ethanol and 2x3x2 for aspirin. A full geometry op-timisation to a tolerance of 0.01 eV/˚A was carried outfor both systems, and then the valence electron cur-rent densities were calculated and written in a binaryfile using CASTEP’s Magres module with the keyword MAGRES WRITE RESPONSE = TRUE .Shielding buildup functions were constructed using aPython script making use of the NumPy, SciPy
ASE and Soprano libraries for computing and pars-ing of structure files. The code has been made availableon GitHub, currently still as a prototype, with the name pynics .Full output structures, the coordinates of the atoms forwhich buildup functions were calculated, and all otherresults are provided as supplementary information. B. Conventions and notation
We clarify here some of the notation used through-out this paper. First, while most NMR literature refer-ences chemical shifts, in this paper we will compute andplot magnetic shieldings. The difference between thesetwo quantities is simply a constant. Magnetic shield-ing tensors express the absolute magnetic response to anexternal field experienced at a site; chemical shifts are the same response referenced to that of another, knownsample. Therefore, we will focus on magnetic shieldings,which are the most natural output of a DFT calculation.The magnetic shielding is a rank 2 tensor with nine inde-pendent components, though for most NMR experimentsonly the symmetric part is of interest, reducing the in-dependent components to only six. In order to expressthem more conveniently, the so-called Haeberlen conven-tion is used . In this convention, the eigenvalues ofthe shielding tensors are computed. The isotropic shield-ing is defined as their average, or one third of the traceof the tensor: σ iso = σ + σ + σ . (11)The eigenvalues are ordered so that | σ − σ iso | > | σ − σ iso | > | σ − σ iso | (12)and two other quantities, anisotropy and asymmetry,are defined respectively as: σ aniso = σ − σ + σ σ asymm = σ − σ σ aniso . (14)Hydrogen atoms will be referred to by the site theyare bonded to. When discussing ethanol, this meansthey will be referred to as CH , CH and OH hydrogensrespectively. Despite being chemically equivalent, thesesites are not crystallographically equivalent, and have dif-ferent shieldings; however, they retain enough similaritythat treating them in groups still makes sense.For aspirin, the situation is more complex. In figure 1 amolecule of aspirin is represented; since the only site thatwas analysed in this case is the CH group, it is labelledappropriately. FIG. 1. Aspirin molecular structure. The CH group is la-belled with a star. C. Results
Shielding buildup functions were computed from theDFT valence electron current densities for all atomic sitesin one molecule of ethanol and aspirin inside the givencrystal structures. Table I shows the final lattice param-eters for the optimised structures as well as the inverse ofthe shortest reciprocal lattice vector, which should giveus a scale for convergence of the shielding.
Molecule a (˚A) b (˚A) c (˚A) 1 / min( | G | ) (˚A)Ethanol 5.33 6.81 8.06 1.25Aspirin 11.32 6.54 11.38 1.8TABLE I. Lattice parameters and inverse shortest reciprocallattice vector for ethanol and aspirin. -70-60-50-40-30-20-10 0 10 20 0.1 1 10 100(a) ∆ σ i s o ( pp m ) R (Ang) HCO -60-40-20 0 20 40 60 80 0.1 1 10 100(b) ∆ σ i s o ( pp m ) R (Ang) HCO
FIG. 2. Long range convergence of σ iso ( R ) − σ iso for all nucleiin one ethanol (a) and one aspirin (b) molecule, grouped byspecies, in their respective crystals. The value of 1 / min( | G | )for each system is marked with a vertical line In Figure 2, the buildup functions for all atoms inethanol and aspirin up to a distance of 100 ˚A are shown,expressed as deviation from the final value, ∆ σ iso ( R ) = σ iso ( R ) − σ iso . This plot makes the scale of the conver-gence apparent. For most nuclei, the largest part of theirshielding converges to a value close to the final one ona very short radius, comparable with the predicted con-vergence scale 1 / min( | G | ). Smaller oscillations continueup to a much larger range, but ultimately, they becomeimperceptible far before the 100 ˚A mark. This confirmsthe initial intuition that the large part of the shielding isjustified by local currents. More insights about the originof various contributions to the shielding can be gained byobserving the buildup functions on a smaller scale.Ethanol gives us a good insight in the way in whichinter-molecular interactions affect the size of the Lorentzsphere. In Figure 3 we can see the buildup functions forthe shielding isotropy on each of the six protons in anethanol molecule. The curves are grouped by carbon sitethat the protons are bonded to. One obvious observa-tion is that five of the curves converge much faster thanthe last one, the one for the OH proton. This is dueto the fact that this proton participates in a hydrogenbond with the oxygen from another molecule. The dis-tance between the proton and the acceptor oxygen is of1.67 ˚A, which corresponds roughly to the point where thecurve converges to a value very close to its final one. Thislong distance effect agrees with the observations made ina similar situation in a previous paper by the authors .The other five shieldings converge much quicker, on alength comparable to 1.1 ˚A, which is the length of the C-H covalent bond. It is also interesting to note that thesefive curves basically overlap on that distance, reflectingthe close similarity of the local electronic wavefunctionfor these five sites, and diverge in distinct groups for theCH and CH protons only afterwards, influenced differ-ently by the rest of the molecule. To the right side of thefigure, the final shieldings are marked by short segments,showing how the curves are almost effectively convergedto their final values.In Figure 4 we see the buildup functions for the twocarbons in ethanol. These curves are very similar inshape, but end up leading to two very distinct final val-ues because of two major differences. The first is that theCH carbon curve descends much quicker and lower ona very short range, within the core radius of the carbonpseudopotential. The second difference is a bump thatis visible around the 2 ˚A mark, in a range highlighted inthe figure and corresponding approximately to the rangeof distances comprising the OH group from the CH car-bon’s position. This suggests that the origin of this bumpis effectively the electronic cloud surrounding that group.The carbon curves also show much better than the pro-ton ones the effect of the rest of the crystal. A feature isvisible starting around 3.4 ˚A. This is approximately thedistance from either of the carbon atoms and the closestatom of a different molecule from the one they belongto. The range of distances representing the closest neigh-bouring molecule is highlighted. This shows very clearlythe importance of the direct effect of the crystal on theshielding, which appears to be similar for both sites. C - H b o n d O H .. O b o n d σ i s o ( pp m ) R (Ang) OHCH CH FIG. 3. Isotropic shielding buildup functions for hydrogenatoms in ethanol, by chemical group. On the right, the pre-dicted final values are shown as lines. Vertical lines are usedto mark the typical length of a C-H covalent bond and of thedistance between the hydrogen and the acceptor oxygen inthe hydrogen bond. -40-30-20-10 0 10 0 1 2 3 4 5 6 7 8 9OH Neighbour σ i s o ( pp m ) R (Ang) CH CH FIG. 4. Isotropic shielding buildup functions for carbon atomsin ethanol, by chemical group. On the right, the predictedfinal values are shown as lines. The range of distances atwhich the OH group appears with respect to the CH one, andthe general range of distances of the next closest neighbouringmolecule for both carbon atoms, are highlighted. Aspirin has been chosen as a test example because ofthe presence of one aromatic ring in its structure, whichis likely to cause strong intermolecular effects. This ismost visible in the protons belonging to the CH group.The arrangement of these protons is shown in Figure 5.All three protons are within a radius of 4 ˚A from the centre of the nearby molecule’s aromatic ring, but oneof them is much closer than the other two, and its C-Hbond is aligned in a direction almost orthogonal to theplane of the ring. To highlight the effect of the ring cur-rents, which is strongly dependent on direction, insteadof plotting the isotropy we plotted the absolute value ofthe anisotropy. The absolute value was picked becauseit avoids the issue of sudden changes in sign which canhappen when two eigenvalues cross each other due to theway the Haeberlen convention is defined. Figure 6 showsthe plot of | σ aniso | for all three protons, with their respec-tive distances from the centre of the nearby ring markedwith lines. It is easy to see that all three protons dis-play strong new features at this range, starting usually abit before the mark (since the mark labels the distancefrom the centre of the ring, the growing integration ra-dius intersects the ring’s electronic cloud a bit earlier).Predictably, the closest one displays the strongest effect. FIG. 5. Detail of aspirin crystalline structure with the CH protons highlighted. The neighbouring aromatic ring is visibleon the right. Finally, Figure 7 provides a clearer overview of the rateat which the shieldings converge. For all atoms in bothethanol and aspirin we can see how the size of the spherenecessary to achieve convergence grows as the toleranceis lowered. Save for a few outliers, the rough estimatesbased on 1 /min ( | G | ) account for convergence up to a 4ppm precision. The size tends then to grow much fasterfor smaller tolerances. IV. CONCLUSIONS
An analytic method to compute a radial buildup func-tion for magnetic shielding starting from the electronic | σ a n i s o | ( pp m ) R (Ang)
FIG. 6. Absolute value of the shielding anisotropy for theCH protons in aspirin. The distances of each proton fromthe centre of the nearby aromatic ring are marked by lines. r c o n v ( A n g ) | σ iso (r)- σ iso | (ppm)AspirinEthanol1/min(|G|) (Ethanol)1/min(|G|) (Aspirin) FIG. 7. Convergence radii for all nuclear shieldings in ethanoland aspirin at different tolerances. The values of 1 / min( | G | )for both are plotted for reference. The limits appear to holdwell for most nuclei for an accuracy down to 4 ppm, whileradii up to more than 7 ˚A are necessary to reach an accuracyof 1 ppm. current density calculated by electronic structure soft-ware in periodic systems has been provided. This methodallows us to have a rough estimate of the scale of what isknown as the Lorentz sphere for a given crystal, and tomake more accurate estimates for specific sites. Obser-vation of practical cases for two molecular crystals con-firms this technique’s potential for quantifying the ef-fects of chemical groups and molecules based on theirdistance, thus giving us a deeper understanding of how crystal structure contributes to NMR spectra. Regard-less of the specific examples used in this work, this tech-nique is applicable to periodic systems in general, includ-ing inorganic framework materials. A rigorous definitionof the size of the Lorentz sphere is proposed, and it isshown that this effectively corresponds to a rather shortrange for most chemical sites in the tested systems. Thisconfirms some long held beliefs in the field of NMR spec-troscopy that however have been historically hard to ex-press in a quantitative way. V. THANKS AND ACKNOWLEDGEMENTS
Computing resources provided by STFC ScientificComputing Department’s SCARF cluster. This work wasconducted within the framework of the CCP for NMRcrystallography, which is funded by the EPSRC grantsEP/J010510/1 and EP/M022501/1. J. R. Zimmerman and M. R. Foster. Standardization of n.m.r.high resolution spectra.
The Journal of Physical Chemistry ,61(3):282–289, 1957. C.J. Durrant, M.P. Hertzberg, and P.W. Kuchel. Magnetic sus-ceptibility: Further insights into macroscopic and microscopicfields and the sphere of lorentz.
Concepts in Magnetic Reso-nance Part A , 18A(1):72–95, 2005. Reinhard Ulrich, Ralf W. Glaser, and Anne S. Ulrich. Suscep-tibility corrections in solid state nmr experiments with orientedmembrane samples. part ii: Theory.
Journal of Magnetic Reso-nance , 164(1):115 – 127, 2003. Miri Zilka, Simone Sturniolo, Steven P. Brown, and Jonathan R.Yates. Visualising crystal packing interactions in solid-stateNMR: Concepts and applications.
The Journal of ChemicalPhysics , 147(14):144203, oct 2017. Cynthia J. Jameson and A. D. Buckingham. Nuclear magneticshielding density.
Journal of Physical Chemistry , 83(26):3366–3371, 1979. Cynthia J. Jameson and A. D. Buckingham. Molecular electronicproperty density functions: The nuclear magnetic shielding den-sity.
The Journal of Chemical Physics , 73(11):5684–5692, dec1980. Inmaculada Gar´ıca Cuesta, Alfredo S´anchez De Mer´as, StefanoPelloni, and Paolo Lazzeretti. Understanding the ring currenteffects on magnetic shielding of hydrogen and carbon nuclei innaphthalene and anthracene.
Journal of Computational Chem-istry , 30(4):551–564, mar 2009. Christopher James Pickard and F Mauri. All-electron mag-netic response with pseudopotentials: Nmr chemical shifts.
Physical Review. B, Condensed matter and materials physics ,63(24):2451011–24510113, 6 2001. Jonathan R. Yates, Chris J. Pickard, and Francesco Mauri. Cal-culation of nmr chemical shifts for extended systems using ultra-soft pseudopotentials.
Phys. Rev. B , 76:024401, Jul 2007. Christian Bonhomme, Christel Gervais, Florence Babonneau,Cristina Coelho, Fr´ed´erique Pourpoint, Thierry Aza¨ıs, Sharon EAshbrook, John M Griffin, Jonathan R Yates, Francesco Mauri,and Chris J Pickard. First-principles calculation of NMR pa-rameters using the gauge including projector augmented wavemethod: a chemist’s point of view.
Chem. Rev. , 112(11):5733–79, November 2012. John David Jackson.
Classical Electrodynamics Third Edition .Wiley, 1998. Stewart J. Clark, Matthew D. Segall, Chris J. Pickard, Phil J.Hasnip, Matt I. J. Probert, Keith Refson, and Mike C. Payne.First principles methods using CASTEP.
Zeitschrift f¨ur Kristal-lographie - Crystalline Materials , 220(5/6), jan 2005. P. G. J¨onsson. Hydrogen bond studies. CXIII. The crystal struc-ture of ethanol at 87 K.
Acta Crystallographica Section B ,32(1):232–235, Jan 1976. New Journalof Chemistry , 26(12):1733–1739, nov 2002. Erik R. McNellis, J¨org Meyer, and Karsten Reuter. Azobenzeneat coinage metal surfaces: Role of dispersive van der waals inter-actions.
Phys. Rev. B , 80:205414, Nov 2009. Travis E. Oliphant. Python for scientific computing.
Computingin Science & Engineering , 9(3):10–20, 2007. K. Jarrod Millman and Michael Aivazis. Python for scientistsand engineers.
Computing in Science & Engineering , 13(2):9–12, mar 2011. Ask Hjorth Larsen, Jens Jrgen Mortensen, Jakob Blomqvist,Ivano E Castelli, Rune Christensen, Marcin Duak, Jesper Friis,Michael N Groves, Bjrk Hammer, Cory Hargus, Eric D Her-mes, Paul C Jennings, Peter Bjerre Jensen, James Kermode,John R Kitchin, Esben Leonhard Kolsbjerg, Joseph Kubal, Kris-ten Kaasbjerg, Steen Lysgaard, Jn Bergmann Maronsson, Tris-tan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson,Carsten Rostgaard, Jakob Schitz, Ole Schtt, Mikkel Strange,Kristian S Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Wal-ter, Zhenhua Zeng, and Karsten W Jacobsen. The atomic simula-tion environmenta python library for working with atoms.
Jour-nal of Physics: Condensed Matter , 29(27):273002, 2017. Soprano - a library developed by the ccp for nmr crystallography. https://ccp-nc.github.io/soprano/ . Accessed: 2018-04-25. Simone Sturniolo. pynics. https://github.com/CCP-NC/pynics ,2018. Accessed: 2018-11-01. Ulrich Haeberlen.
Advances in Magnetic Resonance: Supple-ment . Academic Pr, 1976. Michael Mehring. Nuclear spin interactions in solids. In
Prin-ciples of High Resolution NMR in Solids , pages 8–62. Springer,1983.23