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 fullpotential 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 B ind 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, 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 (B ind orb ) is obtained according to the Biot−Savart law using the induced electric current j ind orb (r) (in atomic units, with c as speed of light): (2) 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, J 0 (r′) is the paramagnetic part of the current operator (the first term in eq 3), J 1 (r′) is the diamagnetic component of the current operator (the second term in eq 3). Ψ̃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 4 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 4 can be cast into a simple form: In the actual implementation, the position operator r is replaced by the limit r·ûi = lim q→0 (1/2q)(e iquî·r − e −iqu i ·r ) to avoid the divergences for extended systems. The integral in eq 2 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 χ ←→ o using an expression derived by replacing r with r·ûi where F i,j = (2−δ i,j )Q i,j , and i and j are the indices of the Cartesian coordinates. The tensor ⃡ Q is calculated with where A k,q i o are the matrix elements between the periodic part of the ground state (u o,k ) and perturbed (u o,k (1) ) Bloch part of the wave functions at k and k + q α This method works well for insulators, but for metals the expression in eq 9 is very difficult to converge with respect to the Brillouin zone (BZ) sampling. This is most likely a consequence of the fact that eq 9 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 χ ←→ o 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 lm n,α,k (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 5 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 spindependent potential leading to a shift of the effective exchangecorrelation 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 (B hf , already available in WIEN2k): 27 hf av 3 (13) where the first term is the Fermi contact term (B c = (8π/3) m av ), and the second represents the dipole field B dip . The Fermi contact term is related to the average spin density (m av ) 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 1 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 m s /B ext . 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.
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 k b T = 10−80 mRy have been applied, we use a Fermi−Dirac function with equal temperatures for spin and orbital components ranging from k b T = 2−8 mRy (2 mRy corresponds to room temperature) and thus to a physically meaningful smearing. Figures 2 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  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 × 10 6 (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 × 10 6 and 1 × 10 6 k-points in the full BZ, do not exceed 1 ppm.
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 RK max = 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 1. 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 [σ(T 1 ) − σ(T 2 )]/(T 1 − T 2 ), where T 1 and T 2 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.
Three of the elements in Table 1 Table 1 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 spinpolarization 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 3) . 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 Mossbauer spectroscopy,   The Journal of Physical Chemistry C Article one often assumes that the measured hyperfine field is directly related to the Fe 3d magnetic moment. The overall comparison between calculated and experimental 29,30 NMR shifts is very good as displayed in Figure 6. 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 = A 0 + A 1 δ th has A 0 = 280 ± 258 ppm and A 1 = 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 Mossbauer spectroscopy. 31 Still, our theoretical approach captures all the main features with reasonable precision.
While our calculated isotropic shieldings/shifts reproduce the measured values very well, the agreement for the macroscopic susceptibility (χ), see Table 2, 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.

Notes
The authors declare no competing financial interest.