NMR Shielding in Metals Using the Augmented Plane Wave Method.

We present calculations of solid state NMR magnetic shielding in metals, which includes both the orbital and the complete spin response of the system in a consistent way. The latter contains an induced spin-polarization of the core states and needs an all-electron self-consistent treatment. In particular, for transition metals, the spin hyperfine field originates not only from the polarization of the valence s-electrons, but the induced magnetic moment of the d-electrons polarizes the core s-states in opposite direction. The method is based on DFT and the augmented plane wave approach as implemented in the WIEN2k code. A comparison between calculated and measured NMR shifts indicates that first-principle calculations can obtain converged results and are more reliable than initially concluded based on previous publications. Nevertheless large k-meshes (up to 2 000 000 k-points in the full Brillouin-zone) and some Fermi-broadening are necessary. Our results show that, in general, both spin and orbital components of the NMR shielding must be evaluated in order to reproduce experimental shifts, because the orbital part cancels the shift of the usually highly ionic reference compound only for simple sp-elements but not for transition metals. This development paves the way for routine NMR calculations of metallic systems.

Introduction

Nuclear magnetic resonance is a widely used experimental technique providing indirect information about the local structure around a selected nucleus for molecules and solids. NMR detects the response of a material to an external magnetic field by measuring the transition energies related to the reorientation of the nuclear magnetic moment.[1] The external magnetic field induces an electric current in the sample, which is the source of an induced magnetic field and partially screens or enhances the external field. The induced current and the resulting internal field are very sensitive to the details of the electronic and atomic structure.

The interpretation of the measured spectra is based on an assignment of NMR shifts to a particular atomic site, usually by comparing the measured shifts with that of similar, but simpler, reference compounds. This procedure may not be trivial when dealing with complicated compounds. Therefore, first-principles calculations of the corresponding magnetic shielding, in particular based on density functional theory (DFT), have been proven very useful. However, such calculations may provide much more information than just the final values of the shielding parameters and allow additional insight into the relation between electronic structure and NMR shielding.[2−4]

Ab initio calculations of NMR shielding of insulating solids are nowadays fairly common. There are several ab initio methods described in the literature for such calculations.[5−12] While the use of hybrid-DFT[13] is quite common in calculations for molecules using Gaussian basis sets, for solid state NMR, they do not seem to provide a systematic improvement[11] and usually NMR calculations are based on the standard DFT[14,15] framework. For metallic systems and the Knight shifts, the literature is far less extensive considering both methods and applications. In this case, besides the orbital motion, also the response from the electron spin has to be considered. Older works either ignored the orbital contribution to the shift[16,17] or use more or less severe approximations for it.[18−21] The modern gauge-invariant projector augmented wave (GIPAW) approach originally implemented for systems with band gaps[7,8] has been also formulated for metals by d’Avezac et al.[22] However, its application was limited to only Li, Al, and Cu, and the quality of the approach is limited (except for Li) by the difficulty to reach k-point convergence and the frozen core pseudopotential method, which cannot include the important effects of core polarization.

In this work, we apply and benchmark our DFT-based full-potential implementation for computing NMR shielding in metals using the all-electron augmented plane wave (APW) basis set as implemented into the WIEN2k package.[23] Following the usual approach, the orbital and spin contributions are computed separately. For the orbital part, we apply a gauge-invariant perturbation method.[11,12] Whereas, for the spin part, a direct self-consistent all-electron approach is employed. This allows one to include all possible contributions to the magnetic shielding. We establish the reliability of our implementation and demonstrate that it is possible to perform such calculations in a routine way.

The paper is organized as follows. In the next section, we present in more detail our theoretical approach for computing both the spin and orbital part of the response. We then demonstrate the nontrivial convergence of the NMR shifts with respect to k-mesh and Fermi smearing and benchmark our method by comparing our computed NMR shifts for several metals with available experimental data. Finally, we summarize and conclude our findings.

Theoretical Approach

Orbital Component

The NMR shielding tensor is defined as a proportionality constant between the induced magnetic field Bind at the nucleus at site R and the external uniform field B:Often only the isotropic shielding can be accessed experimentally. The actual measured quantity is the chemical shift δ, which is the NMR isotropic shielding σ(R) with respect to a reference compound, δ(R) = σref – σ(R).

In metals, the external magnetic field modifies not only the orbital motion of the electrons, but also redistributes the spins. For a weak magnetic field, the two effects can be separated.[22] The orbital part of the induced field (Bindorb) is obtained according to the Biot–Savart law using the induced electric current jindorb(r) (in atomic units, with c as speed of light):

The formalism to calculate the induced current is based on a linear response approach[5,7,8] originally developed by Mauri, Pfrommer, and Louie (MPL)[5] within the pseudopotential or projector augmented wave method. This formalism has been adapted for the all-electron, full potential augmented plane wave method and implemented into our WIEN2k code.[23,24] The details of the implementation for insulators have been described in previous publications.[11,12] Formally our approach belongs to a set of gauge transformation methods, often referred as IGCV (individual gauge for core and valence) with a “d(r) = r” gauge choice for the valence electrons.[25]

The current density is evaluated as an expectation value of the current operator:The expression for the induced current involves only the first order terms with respect to the external field B:where Ψo(0) is an unperturbed Kohn–Sham (KS) occupied orbital, J0(r′) is the paramagnetic part of the current operator (the first term in eq ), J1(r′) is the diamagnetic component of the current operator (the second term in eq ). Ψ̃o(1) is the first order perturbation of Ψo(0) given by the standard formula:where the Greens function = Σ[|Ψe(0)⟩⟨Ψe(0)|]/(ϵ – ϵe) is calculated as a sum over all empty (e) states, and H(1) = (1/2c)r × p·B is the perturbation due to the external magnetic field in symmetric gauge. The last component in eq is proportional to the electron density which may lead to the gauge origin dependence of the total response. The generalized sum rule[7] allows to express this term using wave functions of the occupied bands:Finally, eq can be cast into a simple form:In the actual implementation, the position operator r is replaced by the limit r·û = lim(1/2q)(e – e–) to avoid the divergences for extended systems. The integral in eq is evaluated using a modified approach developed by Weinert[11,26] without imposing any approximation to the induced current.

The uniform component of the induced field (G = 0) is related to the macroscopic susceptibility:[7]Mauri and Pickard[7] calculate the macroscopic susceptibility using an expression derived by replacing r with r·û = lim(1/2q)(e – e–) in the integral evaluating the induced magnetic moment M = ∫d3rr × j.where F = (2−δ)Q, and i and j are the indices of the Cartesian coordinates. The tensor is calculated withwhere Ao are the matrix elements between the periodic part of the ground state (uo,) and perturbed (uo,) Bloch part of the wave functions at k and k + qαThis method works well for insulators, but for metals the expression in eq is very difficult to converge with respect to the Brillouin zone (BZ) sampling. This is most likely a consequence of the fact that eq has the form of a second derivative with respect to q, while the induced current itself is proportional to the first derivative only. To overcome these convergence problems we calculate by direct integration of the total induced orbital moment m = r × j. For this the unit cell is divided into nonoverlapping space-filling atomic basins and each basin is integrated individually.

In the APW method, the unit cell is decomposed into nonoverlapping atomic spheres and an interstitial region. The unperturbed wave functions as well as their first order perturbations are expressed using plane waves augmented with an atomic like angular momentum expansion inside the atomic spheres Sα:Inside the atomic spheres, APW uses numerical radial functions W(r) computed at predefined linearization energies,[24] which are chosen to match the energies of the corresponding occupied bands. This approach yields basically the exact radial wave functions for all occupied and low-energy conduction band states. However, the computation of the orbital part of the shielding tensor relies on the expansion of the perturbation due to an external magnetic field into a complete basis set (unoccupied states). Since our basis is not complete, such an expansion is not accurate, but we have solved this problem by supplying several additional local orbitals (NMR-LO) with radial wave functions evaluated at higher energies[11] and by augmenting the Green’s function in eq by radial functions proportional to r(∂/∂r)u(r) (DUDR).[12]

Spin Component

The induced spin density, responsible for the spin part of the NMR shielding tensor, can be computed within a linear response formalism as proposed in ref (22). However, in this paper, we apply a more direct approach and perform self-consistent spin polarized calculations with a finite external magnetic field acting only on the electronic spin. This interaction with the external field can be cast into a spin-dependent potential leading to a shift of the effective exchange-correlation potentials for the two spins and a finite spin magnetization. It does not break the symmetry of the solid and therefore such calculations are straightforward and do not require large modifications of the existing WIEN2k code. The induced magnetic field at a given nucleus is computed using an expression for the magnetic hyperfine field (Bhf, already available in WIEN2k):[27]where the first term is the Fermi contact term (B = (8π/3) mav), and the second represents the dipole field Bdip. The Fermi contact term is related to the average spin density (mav) over a region near the nucleus with a diameter equal to the Thomson radius.[27] This is also the dominating contribution. The dipole term vanishes for high symmetry structures (as is the case for most metals considered in this study). But we have also noticed that, in cases when it is nonzero, the contribution is small and the value of the integral comes nearly entirely from within the atomic sphere, which further simplifies the calculations. The value of the dipole contribution is in this case related to the population matrix. The spin component of the magnetic susceptibility is proportional to the induced spin moment. This approach makes sense only when the response for both the hyperfine field and the induced spin moment depends linearly on the external magnetic field. Our tests show that this is the case for all systems considered in this work; as an example, Figure shows the results computed for Al. The dependence is linear even for fields as large as 200 T, with a negligible standard error for the slopes σs and ms/Bext. In order to obtain a strong enough response to evaluate the NMR shielding with a precision better than 1 ppm, we apply in our calculations a field equal to 100 T, which induces a spin-splitting of approximately 1 mRy. The spin susceptibility χs can easily be obtained from the total induced spin magnetic moment in the cell.

Figure 1

Dependence of (a) the induced hyperfine field Bhf and (b) the induced spin moment on the external magnetic field computed for fcc Al. σs and ms/Bext correspond to the slopes of a linear fit, and Δ denotes the corresponding standard deviations.

Computational Details

The orbital and spin components of the induced magnetic field are very sensitive to the details of the Fermi surface and require very dense samplings of the BZ. The convergence can be accelerated by an appropriate Fermi–Dirac “smearing” of the occupancy around the Fermi level. This issue has been discussed earlier by d’Avezac et al. in ref (22)., but while in this paper huge broadening parameters kbT = 10–80 mRy have been applied, we use a Fermi–Dirac function with equal temperatures for spin and orbital components ranging from kbT = 2–8 mRy (2 mRy corresponds to room temperature) and thus to a physically meaningful smearing. Figures and 3 present some convergence tests of σ and χ with respect to the k-mesh sampling in the BZ for three different smearing values. We have selected fcc Al as the example, since this was the most difficult case mentioned in ref (22). Naturally, broader “smearing” (larger temperature) leads to faster and smoother k-point convergence, but still, typical k-meshes above 60 × 60 × 60 are necessary. Since WIEN2k does not need to keep all wave functions in memory, there is no particular barrier. Unfortunately, the value of the “smearing” parameter affects the final converged values of the shielding. The actual dependency may vary from case to case, but Figures and 5 show the convergence of σ and χ with respect to the smearing for Al. The variation of the shielding and susceptibility is in this case fairly linear but also quite large. In the range between 2 and 8 mRy, the spin and orbital shielding varies by as much as 50 and 20 ppm, respectively. Therefore, it is necessary to estimate the effect of the smearing in all systems. In order to do that, we will quote the values of the spin and orbital shielding computed broadening parameters equal to 8 mRy and the coefficient describing the dependence on the smearing. In order to automatize the calculations we use a large number of about 2 × 106 (126 × 126 × 126) k-points in the full BZ for all systems. The errors due to finite BZ sampling, estimated by comparing results obtained with 2 × 106 and 1 × 106 k-points in the full BZ, do not exceed 1 ppm.

Figure 2

K-point convergence of the (a) spin and (b) orbital components of the isotropic shielding in fcc Al. kbT (in mRy) is the temperature in a Fermi–Dirac function used for smearing the band occupancies around the Fermi level.

Figure 3

k-Point convergence of the (a) spin and (b) orbital components of the macroscopic susceptibility of fcc Al.

Figure 4

Dependence of the (a) spin and (b) orbital components of the isotropic shielding in fcc Al on the smearing of the occupancy at the Fermi level (Fermi–Dirac temperature in mRy).

Figure 5

Dependence of the (a) spin and (b) orbital components of the macroscopic susceptibility of Al on the smearing of the occupancy at the Fermi level (Fermi–Dirac temperature in mRy).

Other computational parameters like atomic sphere radii, angular momentum expansions of charge densities, potentials and wave functions inside the atomic spheres as well as the linearization energies are kept as set by WIEN2k_15.1 defaults. The plane wave basis set size was determined using RKmax = 8, and we use in all cases the generalized gradient approximation by Perdew et al.[28]

Results and Discussion

The calculated isotropic shielding for selected simple metals are reported in Table . The spin (σs) and orbital (σo) contribution to the shielding are calculated using a Fermi–Dirac “smearing” of the occupancy around the Fermi level with a temperature set to 8 mRy. The coefficient describing the dependency of the shielding on the “smearing” (shown in parentheses) has been evaluated as [σ(T1) – σ(T2)]/(T1 – T2), where T1 and T2 are set to 8 and 4 mRy. This ensures k-point converged shielding values, which could not be achieved in all cases for smaller temperatures. For most cases the sensitivity to the “smearing” is of the order of a few tens of ppm/mRy. The exceptional case is Cs, where both the value for σs (−16 177 ppm) and its temperature dependency (587 ppm/mRy) are very large. In fact, Cs shows also in experiment a strong temperature dependency of the NMR shift[29,30] in the range of 1.57–1.49% for temperatures from 4 to 300 K.

Table 1

Isotropic NMR Shielding and Shift of Elemental Metalsa

	reference	σo	σs	σref	δth	δexp	Δth–ex
Li	Li2O	81 (0)	–264 (0)	95.5	279	260	19
Na	NaBr	518 (0)	–1021 (−1)	551.3	1054	1070	–16
K	KBr	1126 (8)	–2560 (−2)	1153.2	2589	2500	89
Rb	RbCl	3031 (3)	–6826 (−1)	3027.5	6822	6460	362
Cs	CsCl	5473 (4)	–16 177 (587)	5380.1	16 083	15 700	383
Mg	MgCl2	505 (−0)	–1078 (4)	552.1	1124	1120	4
Ba	BaCl2	5730 (−28)	–4160 (−80)	5661.0	4092	4030	62
Al	AlPO4	519 (−4)	–1591 (11)	511.6	1584	1595	–11
In	In2(SO4)3	2807 (10)	–8012 (−43)	3676.1	8881	8300	581
V	NaVO3	–5988 (62)	488 (6)	–1453.25	4046	5800	–1754
Cr	Na2CrO4	–9847 (−46)	461 (0)	–2567.3	6818	6900	–82
Mo	K2MoO4	–5795 (−44)	–27 (−30)	–824.9	4997	6100	–1103
Cu	CuBr	–330 (−8)	–1568 (2)	492.0	2390	2380	10
Ag	AgNO3	2219 (5)	–3670 (−1)	3771.5	5223	5210	13

σo and σs are the spin and orbital contributions to the shielding computed using temperature smearing of 8 mRy. The coefficient (in ppm/mRy) describing the dependence of the shielding on the “smearing” (calculated as[σ(T1) – σ(T2)]/(T1 – T2), where T1 and T2 are the temperatures in the Fermi–Dirac distribution and set to 8 and 4 mRy) is given in parentheses. σref is the isotropic shielding of the reference compound. The last three columns compare our calculated chemical shift (δ = σref – σo – σs) with the experimental values taken from refs (29) and (30) (when available high temperature data are cited) . Δth–ex gives the difference between calculated (δth) and experimental (δex) values.

Three of the elements in Table , Li, Al, and Cu, have been considered in a previous attempt to compute the NMR shielding in metals by d’Avezac et al. in ref (22). However, when comparing our results with the published results, we see good agreement only for Li, where the differences for all components of σ do not exceed 2 ppm. For Al and Cu, some shifts can be different by as much as 300–400 ppm. Reference (22) reports spin and orbital components of the shielding in Al equal to −1858 ± 70 and 548 ± 8 ppm, respectively. Comparing these numbers to ours (Table ), we see that most of the discrepancy comes from the spin component. The orbital contribution deviates only by 44 ppm. For Cu, ref (22) reports −2336 ± 20 and −874 ± 10 ppm for σs and σo, whereas our values are −1568 and −329 ppm, respectively. In this case, also the reported[22] value of 424 ppm for the reference compound (CuBr) is considerably different from our 492 ppm. On the other hand, our computed components of the macroscopic susceptibility, shown in Table , are in good agreement with ref (22) except for the orbital part of Cu. It is difficult to speculate about the source of those discrepancies (k-point convergence, too large smearing parameters, pseudopotentials); however, it is clear that the comparison between calculated and measured NMR shifts are in favor of our results.

Table 2

Calculated Macroscopic Susceptibilities (in 10–6 cm3/mol) of Simple Metalsa

	χo	χs	χth	χexp	Δth–ex
Li	0.0 (0.0)	29.3 (28.4)	29.3 (28.4)	24.5	4.8
Na	–4.6	24.6	20.0	16	4.0
K	–13.8	44.9	31.1	20.8	10.3
Rb	–23.0	54.7	31.7	17	14.7
Cs	–35.9	88.6	52.7	29	23.7
Mg	–5.3	37.6	32.3	13.1	19.2
Ba	–4.3	85.8	81.5	20.6	60.9
Al	–1.9 (−1.1)	17.4 (17.7)	15.5 (16.6)	16.5	–1.0
In	–19.1	17.8	–1.3	–10.2	8.9
V	100.6	189.5	290.1	285	5.1
Cr	124.6	35.0	159.6	167	–7.4
Mo	56.2	76.7	132.9	72	60.9
Cu	–2.4 (−17.6	11.2 (10.8)	8.8 (−6.8)	–5.6	14.4
Ag	–18.2	9.8	–8.4	–19.5	11.1

χo and χs correspond to the orbital and spin components of the total χth. The measured values of χexp have been taken from refs (32) and (33). The GIPAW values for Li, Al, and Cu taken from ref (22) are given in parentheses. The last column (Δth–exp) gives the difference between calculated and experimental values.

The results included in Table show that, in principle, both the orbital and the spin components of the shielding have considerable values. However, for most of the sp-metals (alkali and alkali-earth metals, Al, Ga), the orbital part σo of the metal and the corresponding (usually) highly ionic reference compound cancel out to a large extent, so that the final shift is fairly close to the absolute value given by σs. This can be understood due to the fact that the diamagnetic core contribution is more or less identical and also the valence contribution of the ion (formally no occupied valence electron) and the metal (only s-electrons) are always very small and nearly identical. This is, however, not true for transition metals, where large differences between σo of the metal and the reference compound exist. For instance, in Cu, σo is negative (paramagnetic) since the 3d electrons are not completely occupied and their contribution dominates over the diamagnetic core contribution, whereas the orbital response in CuBr is negative (diamagnetic), since in the Cu+ ion the 3d shell is completely occupied and behaves like a closed shell.

Interestingly, the transition metals with partly occupied d states (V, Cr, and Mo) show a quite exotic behavior which is very different from that of the other elements in this study. The orbital shielding is large and negative (dominated by the large paramagnetic shift of the d electrons), whereas the spin part is either rather small compared to other metals or even positive (for V and Cr) . Of course, also these metals have large paramagnetic Fermi contact contributions due to the spin-polarization of the valence 4s(5s for Mo) electrons, but in addition the magnetic field also introduces a sizable 3d(4d in Mo) spin-magnetic moment, which introduces a huge core polarization of opposite sign (see Table ) . A similar cancellation has been observed by Ebert and co-workers[18,19] applying multiple scattering theory and KKR-type band structure approach. We can capture this effect because in our all-electron self-consistent calculations we allow also a relaxation of the core electrons. The dominance of the core polarization over the valence 4s hyperfine field is well-known from hyperfine field calculations/measurements in ferromagnets (e.g., bcc Fe),[27] and thus, in Fe Mössbauer spectroscopy, one often assumes that the measured hyperfine field is directly related to the Fe 3d magnetic moment.

Table 3

Angular Momentum Decomposition (s,p,d) of the Induced Spin Magnetic Moment within the Atomic Spheres (for a Field of 200 T) and Valence and Core Components of the Spin Component of the Shielding (only the Contact Term Is Listed and “Valence” Includes All Atomic States with Energies Higher than −6 Ry, for Example, Also the ”Semicore” States (1s and 2s in Li, 2s and 3s in Na, ...))

	ms (10–3μB)			σs (ppm)
	s	p	d	valence	core
Li	0.4	1.1	0.0	264	0
Na	0.6	0.5	0.1	1033	–12
K	0.4	0.2	0.0	2556	4
Rb	0.4	0.2	0.2	6795	31
Cs	0.5	0.2	0.7	16 084	93
Mg	0.5	0.7	0.3	1097	–16
Ba	0.0	0.3	2.1	4078	82
Al	0.3	0.6	0.2	1584	7
In	0.4	0.9	0.0	7956	56
V	0.4	1.9	25.1	3439	–3927
Cr	0.0	0.2	5.3	613	–1074
Mo	0.0	0.2	3.1	150	–123
Cu	0.1	0.2	0.7	1677	–109
Ag	0.2	0.7	0.3	3708	–39

The overall comparison between calculated and experimental[29,30] NMR shifts is very good as displayed in Figure . In an ideal situation, a linear least-square fit of theory vs experiment should have a slope of one with zero offset. Our fit δexp = A0 + A1δth has A0 = 280 ± 258 ppm and A1 = 0.964 ± 0.042, which fulfills this constrain within one standard deviation. In particular, when considering the scale, the measured shifts vary from 260 ppm to nearly 16 000 ppm, the agreement is definitely very good, and only V and Mo visually deviate from the linear least-square fit. We attribute these exceptions to possible problems of the generalized gradient DFT approximation with the dramatically different bonding character between metal and the reference compound or the (huge) core response, which is also known to lead to some errors in hyperfine field calculations for Mössbauer spectroscopy.[31] Still, our theoretical approach captures all the main features with reasonable precision.

Figure 6

Correlation between calculated and measured NMR shifts for simple metals. The line shows least-square fit, δexp = A0 + A1δth, with A0 = 280 ± 258 ppm and A1 = 0.964 ± 0.042.

While our calculated isotropic shieldings/shifts reproduce the measured values very well, the agreement for the macroscopic susceptibility (χ), see Table , is in general not of the same quality. The susceptibilities for Li, Al, and Cu have been computed before by d’Avezac et al.,[22] and their spin components (χs) are very similar to our results for all three elements. Our calculations also agree with ref (22) on the small orbital contribution (χo) in Li and Al, but there is a quite large discrepancy for the orbital component in Cu. They reported a large negative value for χo resulting in a total diamagnetic χ in agreement with experiment, while our value of χo is much smaller resulting in a positive total χ. Interestingly, the computed and measured magnetic susceptibilities are particularly different for heavier sp-metals, while the errors for 3d elements are much smaller. In general, the calculated susceptibilities for metals seem to be in much poorer agreement with experiment than for insulators.[12] In any case, the effect of the magnetic susceptibility on the NMR shielding is in most cases only few ppm and the potentially introduced error for the NMR shielding is therefore relatively small.

Conclusions

NMR shielding of insulating solids can nowadays be routinely calculated by various ab initio packages. This was, however, not the case for metals, where the k-point convergence and the use of pseudopotential methods is a severe issue. A first approach using the GIPAW method has been published by d’Avezac et al.,[22] but the discussion was limited to Li, Cu, and Al and the results were not very promising showing large discrepancies between experiment and theory of 300–400 ppm for Al and Cu. In the current paper we revisited the subject in a more systematic way using our recent full potential APW implementation for NMR shielding. We have computed NMR shieldings and shifts for a set of metallic elements including alkali, alkali-earth, transition and noble metals. The overall agreement between measured and calculated shifts is very good, including the elements used in ref (22). The calculations are, however, more demanding than for insulators or semiconductors, especially in terms of k-point convergence, but together with a reasonable Fermi smearing converged results can be obtained. We find in particular for transition metals that in general both orbital and spin contributions to the NMR shielding are needed in order to properly reproduce experimental data. Only for alkali and alkali-earth metals, which are referenced to ionic insulators, most of the orbital contribution cancels out and only the spin contribution (Knight shift) determines the relative NMR shift of the metal against the insulator. An all-electron implementation which allows for core polarization effects is necessary in general, since the induced hyperfine field originating from the core electrons can dominate the direct valence electron polarization. We conclude that first-principles calculations of the NMR shielding in metals is now possible, and our scheme paves the way for routine calculations also in more complicated systems.