First-Principles Calculation of NMR Parameters Using the Gauge

Review pubs.acs.org/CR

First-Principles Calculation of NMR Parameters Using the Gauge Including Projector Augmented Wave Method: A Chemist’s Point of View Christian Bonhomme,*,† Christel Gervais,*,† Florence Babonneau,† Cristina Coelho,‡ Frédérique Pourpoint,† Thierry Azaïs,† Sharon E. Ashbrook,*,§ John M. Griffin,§ Jonathan R. Yates,*,∥ Francesco Mauri,⊥ and Chris J. Pickard# †

Laboratoire de Chimie de la Matière Condensée de Paris, Université Pierre et Marie Curie, Paris 06, CNRS UMR 7574, Collège de France, 75005 Paris, France ‡ IMPC, Institut des Matériaux de Paris Centre, FR2482, UPMC Université Pierre et Marie Curie Paris 06, Collège de France, 11 place Marcelin Berthelot, 75231 Paris Cedex 05, France § School of Chemistry and EaStCHEM, University of St. Andrews, North Haugh, St. Andrews KY16 9ST, United Kingdom ∥ Department of Materials, University of Oxford, Oxford OX1 3PH, United Kingdom ⊥ Laboratoire de Minéralogie Crystallographie, UMR CNRS 7590, Université Pierre et Marie Curie, UPMC, 75015 Paris, France # Department of Physics and Astronomy, University College London, London WC1E 6BT, United Kingdom 4. Applications of First-Principles Calculations in Chemistry 4.1. GIPAW: A Tool for Chemistry and Materials Science 4.2. Organic Structures and Assemblies 4.2.1. Amino Acids, Peptides, and Nucleic Acids 4.2.2. Carbohydrates and Oligosaccharides 4.2.3. Drugs and Polymorphism 4.2.4. Supramolecular Assemblies 4.2.5. Graphene and Carbon Nanotubes 4.2.6. Disorder and Local Dynamics 4.3. Inorganic Architectures 4.3.1. Si-Containing Structures 4.3.2. Al-Containing Structures 4.3.3. Borates and Borides 4.3.4. Niobiates 4.3.5. Phosphates 4.3.6. Li-Containing Materials 4.3.7. F-Containing Materials 4.3.8. Microporous Materials and Molecular Sieves 4.3.9. Local Dynamics 4.3.10. Disorder in Materials 4.4. Polymers and Hybrid Materials 4.5. Biomaterials and Nanomaterials 4.5.1. Biocompatible Materials and Biomaterials 4.5.2. Nanomaterials and Nanoparticles 4.6. Surfaces and Catalysis 4.7. From NMR and GIPAW to 3D Structures 4.7.1. Inorganic Structures 4.7.2. Organic Structures 4.8. Other Applications of GIPAW 4.8.1. Antisymmetric Parts of Tensors

CONTENTS 1. Introduction 2. Theory and Methods 2.1. Physical Basis of NMR 2.1.1. Magnetic Shielding 2.1.2. Spin−Spin Coupling 2.1.3. Electric Field Gradients 2.2. Electronic Structure Calculations 2.2.1. Hartree−Fock 2.2.2. DFT 2.2.3. Modeling of Solid Materials 2.2.4. Basis Sets 2.2.5. Pseudopotentials 2.3. Calculations of NMR Parameters 2.3.1. Shielding 2.3.2. Other NMR Parameters 2.3.3. Overview of GIPAW Method 2.3.4. Discussion of “Accuracy” 3. Experimental Solid-State NMR: Methods and Implementation 3.1. Basic Experimental Techniques (Spin I = 1/2) 3.2. Quadrupolar Nuclei (Spin I > 1/2) 3.3. Wideline NMR 3.4. Measurement of Interaction Parameters 3.5. Two-Dimensional Correlation Experiments 3.6. Combination of GIPAW Calculations and Magnetic Resonance

5734 5735 5735 5735 5736 5736 5737 5737 5737 5737 5738 5738 5739 5739 5741 5741 5742 5743 5743 5744 5745 5745 5746

5747 5747 5751 5751 5752 5753 5753 5754 5754 5756 5756 5758 5758 5758 5759 5759 5759 5760 5760 5761 5763 5763 5763 5764 5764 5765 5765 5766 5767 5767

5747 Received: March 13, 2012 Published: November 1, 2012

© 2012 American Chemical Society

5733

dx.doi.org/10.1021/cr300108a | Chem. Rev. 2012, 112, 5733−5779

Chemical Reviews 4.8.2. Stress in NMR 4.8.3. NQR Spectroscopy 4.8.4. g and Hyperfine Tensors 4.8.5. Mö ssbauer Spectroscopy 4.8.6. Other Spectroscopies 5. Future Work and Conclusions Author Information Corresponding Author Notes Biographies References

Review

can provide this necessary link and enable the full power of NMR spectroscopy to be exploited. The theoretical foundations for the calculation of NMR parameters were set by Ramsey in a series of pioneering papers in the early 1950s.11−13 Over subsequent decades, practical schemes to compute NMR parameters were developed within the quantum chemistry community and have been of great interest to experimentalists.14−26 Generally speaking, two major approaches are used for calculations. The first corresponds to wave function-based methods, which start with Hartree−Fock (HF) theory and encompass sophisticated correlated techniques such as configuration interaction. The second approach is density functional theory (DFT), which as the name suggests uses the electronic density as the fundamental variable. Wave function-based methods provide a systematic path to improve the accuracy of the calculations. However, this improvement increases significantly the computational cost of the calculations. DFT calculations provide a practical balance between accuracy and computational cost and are now suitable for medium to large molecules. For example, it was reported in 2004 that using linear scaling techniques, a calculation involving more than 1000 atoms was possible.27 At this stage, the reader is referred to several reviews dealing with HF theory and DFT for the calculation of chemical shifts and J coupling constants in molecular systems.28−33 A very detailed theoretical review is presented in ref 34. In the case of solid-state NMR, it is important to take account of the crystalline nature of the material.35 Traditional quantum chemical methods can be applied through use of a cluster of atoms constructed such that those in the center experience an environment similar to that in true extended solid. Convergence with respect to the size of the cluster must of course be tested carefully. There have been numerous applications of cluster models to interpret solid-state NMR experiments.36−40 Several techniques have been developed to include the longrange effects of the crystal lattice through a series of point dipoles.41,42 These electrostatic models have also been combined with cluster calculations.43 A third approach is to employ electronic structure techniques that model extended solids using periodic boundary conditions. The linear augmented plane-wave approach (LAPW) through the Wien series of codes has long been used to compute quadrupolar coupling constants.44 The first approach to compute magnetic shieldings under period boundary conditions was developed by Mauri et al.45 Later, Sebastiani et al.46 developed a method using the Carr−Parrinello Molecular Dynamics (CPMD) code47 to calculate chemical shifts under periodic boundary conditions.48−51 The same authors made also many important contributions in the field of hybrid quantum QM/MM molecular dynamics simulations.52,53 A major breakthrough in calculations using periodic boundary conditions occurred in 2001, with the pioneering work of Pickard and Mauri (Phys. Rev. B 2001, 63, 245101), leading to the now well-established GIPAW (Gauge Including Projector Augmented Wave) method. Based on DFT, a planewave basis set, pseudopotentials, and PAW,54 this new approach provided the ability to compute NMR tensors in extended systems with all-electron accuracy. It follows that periodic boundary conditions allow a truly solid-state calculation to be performed, which can be usefully applied for structural and chemical purposes. This highly innovative approach has led to a true revolution in the chemical

5767 5768 5768 5768 5769 5769 5770 5770 5770 5770 5773

1. INTRODUCTION First-principles calculations of NMR parameters in solids (including the chemical shift, the quadrupolar interaction for nuclei with spin quantum number I > 1/2, and the indirect J coupling constant) are of crucial importance in chemistry and materials science. During the past two decades, major methodological and instrumental advances have been achieved in solid-state NMR.1,2 High-resolution methods have been improved for I = 1/2 spins (1H, 19F, 31P, etc.):1 ultrafast magicangle spinning (MAS) up to 70 kHz at ultrahigh field (23.3 T), and efficient homonuclear decoupling in the case of strongly coupled nuclei such as 1H. In parallel, ingenious high-resolution techniques have been implemented for quadrupolar nuclei (I > 1/2: 17O, 23Na, 27Al, etc.)2 such as double rotation (DOR), dynamic angle spinning (DAS), multiple-quantum magic-angle spinning (MQMAS), and satellite-transition magic-angle spinning (STMAS). All of these methods allow the experimentalist to determine (at least in principle) the asymmetric unit of a given structure, or the distribution of isotropic chemical shifts and quadrupolar parameters in the case of amorphous materials. The acquisition of a high-resolution spectrum, however, does come with a cost, that is, the suppression (averaging) of all anisotropies. The interactions that affect NMR spectra are inherently anisotropic and can be represented mathematically by second- and fourth-rank tensors.3−5 It is well-known that anisotropies reflect detailed structural features and that their measurement and analysis can be valuable for an in-depth description of solid-state structures. The apparent contradiction between the need to acquire high-resolution spectra (without any remaining anisotropy) to resolve distinct resonances, and the necessity to measure simultaneously the anisotropic information was solved elegantly by the so-called “recoupling” techniques and formalized by Levitt.6 Generally, two-dimensional (2D) experiments are performed with one highly resolved dimension (obtained under MAS or MQMAS), while the second dimension contains information related to a selected anisotropic interaction. The use of recoupling experiments in combination with high-resolution schemes has led recently to major breakthroughs in the detailed description of through-bond and spatial connectivities for biosolids,7 inorganic/organic structures,8,9and glasses.10 It follows that solid-state NMR experiments offer new perspectives in terms of structural characterization when compared to powder X-ray diffraction (XRD) due to their sensitivity to the local, atomicscale environment. However, it is clear that an additional tool is required to fully exploit the information available, that is, to “bridge” between structural models and experimental NMR data. It is becoming apparent that first-principles calculations 5734


Chemical Reviews

Review

such an effective Hamiltonian can be obtained from the full crystal Hamiltonian by integrating over all degrees of freedom except for the nuclear spins and external fields. The effect of the electrons and positions of the nuclei is now incorporated into a small number of tensor properties, which define the key interactions in NMR. These tensors can be obtained from electronic structure calculations, and we now examine them in turn. 2.1.1. Magnetic Shielding. The Zeeman interaction between a magnetic field, B, and a set of spin 1/2 nuclei K is given by:

application of solid-state NMR, enabling the accurate calculations of NMR tensors in solids. In this Review, all applications to date (i.e., 2001−2012) of GIPAW in chemistry and materials science are presented. We have decided to ensure this Review is self-contained for the reader, with a particular emphasis on the theory and methods supporting the GIPAW approach (section 2), a brief introduction to the current methodological advances in solidstate NMR (section 3), and the various chemical topics relevant to GIPAW (section 4). Section 2 provides a comprehensive introduction to the physical basis of NMR, electronic structure calculations, and the subsequent calculation of NMR parameters in the context of DFT and periodic boundary conditions. These parameters include mainly the shielding of the nucleus, the electric field gradient, and the J coupling. Most importantly, a discussion of the “accuracy” of the approach is presented. Section 3 is devoted to the basics of solid-state NMR. Highresolution and two-dimensional (2D) correlation experiments (based on J and dipolar interactions) are discussed for spin I = 1/2 nuclei as well as for quadrupolar nuclei (I > 1/2). Recoupling experiments are introduced in the context of the measurement of interaction parameters, that is, internuclear distances and anisotropies. The main goal of this section is to highlight the fundamental concepts involved, and the methods required, for experimental measurement parameters. This will aid the reader in the understanding and appreciation of the examples discussed in section 4. Section 4 reviews the application of GIPAW in the NMR and chemistry communities. We note that, very recently, Charpentier55 published an interesting overview of the role of GIPAW calculations with particular focus on the characterization of crystals and glasses. However, in this work, we undertake an exhaustive review of the application of GIPAW to all chemical systems. Particular emphasis is placed on important classes of systems such as organic structures and assemblies, inorganic architectures, polymers and hybrids, biomaterials, and nanomaterials. The application of NMR to areas such as surfaces, interfaces, and catalysis is also covered. In our view, one of the major achievements of GIPAW in recent years has been to strengthen the link between NMR and crystallography, leading to the formalization of what is called now “NMR crystallography”.56−59 The effective combination of NMR and GIPAW to provide 3D structure determination is described in depth. Additional applications of GIPAW are presented as well, including less developed but original topics in NMR and NQR (nuclear quadrupolar resonance), ESR (electron spin resonance), and Mössbauer spectroscopies. Section 5 highlights the future of GIPAW in chemistry and materials science from both methodological and applications points of view. It is hoped that this work will not only provide an effective reference for current practitioners in the field, but will demonstrate the power that first-principles calculations have to unlock the potential of NMR spectroscopy for understanding solid-state structure and will inspire many new exciting investigations in the future.

H = −∑ γK IK ·B

(1)

K

where ℏIK is the spin angular momentum and γK is the gyromagnetic ratio of nucleus K. If we consider B as the field at the nucleus due to the presence of an externally applied field Bext, we can express eq 1 as: H = −∑ γK IK (1 − σK⃡ )Bext

(2)

K

The first term is the interaction of the bare nucleus with the applied field, while the second accounts for the response of the electrons to the field. This electronic response is characterized by the absolute magnetic shielding tensor σ⃡ K, which relates the induced field to the applied field: Bin (R K ) = −σK⃡ Bext

(3)

In a diamagnet, the induced field arises solely from orbital currents j(r), induced by the applied field: Bin (r) =

1 c

∫

d3r′ j(r′) ×

r − r′ |r − r′|3

(4)

The shielding tensor can equivalently be written as a second derivative of the electronic energy of the system: σ⃡K =

∂ 2E ∂mK ∂Bext

(5)

In solution state NMR, or for powdered solids under magic angle spinning (MAS) conditions (see section 3), we are mainly concerned with the isotropic part of the shielding tensor σiso = 1/3Tr[σ⃡]. The magnetic shielding leads to nuclei in different chemical environments resonating at frequencies that are slightly different from the Larmor frequency of the bare nucleus. Rather than report directly the change in resonant frequency (which would depend on the operating magnetic field of the spectrometer), a normalized chemical shift is reported: νsample − νref δ= νref (6) where νref is the resonance frequency of a standard reference sample. The magnetic shielding and chemical shift are related by: σref − σsample δ= 1 − σref (7)

2. THEORY AND METHODS 2.1. Physical Basis of NMR

For all but very heavy elements, |σref| ≪ 1 and so:

The correct language to describe NMR is that of effective nuclear Hamiltonians (see, e.g., ref 60). In a conceptual sense,

δ ≈ σref − σsample 5735

(8)


Chemical Reviews

Review

Experimentally derived absolute shielding scales have been proposed for several nuclei, for example, 19F61 and 13C.62 Using such an absolute scale provides a strict test of the calculations; however, it is more common to determine σref indirectly by comparing calculated magnetic shielding and observed chemical shift in some way as this benefits from cancelation of errors. Strategies for doing this have been discussed in refs 56 and 63. An aspect of NMR that is particularly relevant for solid-state studies is that the key interactions are anisotropic, that is, orientation dependent. For example, the chemical shift tensor (or equivalently the shielding tensor) is a general second-rank tensor and can be divided into a symmetric and an antisymmetric contribution. The influence of the antisymmetric part of the tensor is discussed in section 4.8.1. The symmetric part of the tensor can be diagonalized to give three principal components, δii (δii = σref − σii), together with three Euler angles, which describe the orientation of the tensor in the crystal frame. From the three principal values, the chemical shift anisotropy (CSA) δaniso = δ33 − δiso and the asymmetry parameter ηδ = (δ22 − δ11)/δaniso) can be determined. δiso is the isotropic chemical shift (1/3 (δ11 + δ22 + δ33)). The following convention is used here: |δ33 − δiso| ≥ | δ11 − δiso| ≥ | δ22 − δiso|,64 although a number of other conventions can be found in the literature.65,66 2.1.2. Spin−Spin Coupling. In the previous section, we considered the effect of the magnetic field at a nucleus resulting from an externally applied field. However, there may also be a contribution to the magnetic field at a nucleus arising from the magnetic moments of the other nuclei in the system. In an effective spin Hamiltonian, we may associate this spin−spin coupling with a term of the form:

∑

H=

⃡ )IL IK (D⃡KL + JKL

K 1/2, the NMR response will include an interaction between the quadrupole moment of the nucleus, eQ, and the electric field gradient (EFG) generated by the surrounding electronic structure.3 The EFG is a second rank, symmetric, traceless tensor G⃡ (r) given by

(11)

An equivalent expression arises from considering one nuclear spin (L) as perturbation, which creates a magnetic field at a second (receiving) nucleus (K): 2π ⃡ J ·mL ℏγK γL KL

4π

+

(9)

ℏγK γL

B(1) in (R K ) =

μ0 +

where rKL = RK − RL with RL the position of nucleus L and μ0 the vacuum permeability. D⃡ KL is a traceless tensor, and its effects will be averaged out under MAS when the spinning speed is much greater than D⃡ KL. J ⃡ KL is the indirect coupling and represents an interaction of nuclear spins mediated by the electrons. The J coupling is a small perturbation to the electronic ground state of the system, and it can be identified as a derivative of the total electronic energy E of the system: ⃡ = JKL

(14)

Here rL = r − RL with RL the position of nucleus L, μ0 the permeability of vacuum, δ the Dirac delta function, S the Pauli spin operator, g the Lande g-factor, and β the Bohr magneton. These perturbations give rise to a spin density m(1)(r) in the materials, which in turn creates a magnetic field through a second hyperfine interaction. By working to first order in these quantities, we can write the magnetic field at atom K induced by the magnetic moment of atom L as:

D⃡ KL is the direct dipolar coupling between the two nuclei and is a function of only the nuclear constants and the internuclear distance: 2 ℏ μ0 γK γL 3rKLrKL − 1 |rKL| D⃡KL = − 5 2π 4π |rKL|

⎛ 3r (μ ·r ) − |r |2 μ 1 ⎞ L L L L L ⎟ S·⎜⎜ ⎟ 4π ⎝ |rL|5 ⎠ μ0

Gαβ(r) =

∂Eγ (r) ∂Eα(r) 1 − δαβ ∑ ∂rβ 3 ∂rγ γ

(17)

where α,β,γ denote the Cartesian coordinates x,y,z, and Eα(r) is the local electric field at the position r, which can be calculated from the sum of nuclear and electron charge density n(r):

(12)

Equation 12 shows that the question of computing J ⃡ is essentially that of computing the magnetic field induced indirectly by a nuclear magnetic moment. The first complete analysis of this indirect coupling was provided by Ramsey.13,67

Eα(r) = 5736

∫

d3r

n(r) (rα − rα′) |r − r′|3

(18)


Chemical Reviews

Review

many-body problem can be exactly expressed in terms of a noninteracting set of fictitious particles whose charge density, and hence total energy, is identical to that of the true system. Using atomic units:

The EFG tensor is then equal to: Gαβ(r) =

∫

d3r

(rα − rα′)(rβ − rβ′) ⎤ n(r) ⎡ ⎢δαβ − 3 ⎥ 3 |r − r′| ⎣ |r − r′|2 ⎦ (19)

1 − ∇2 Ψ n(r) + veff (r)Ψ n(r) = εn Ψ n(r) 2

The quadrupolar coupling constant, CQ, and the asymmetry parameter, ηQ, can be obtained from the diagonalized electric field gradient tensor whose eigenvalues are labeled Vxx, Vyy, Vzz, such that |Vzz| > |Vyy| > |Vxx|:69

CQ =

eVzzQ h

with the charge density ρ obtained from a sum of occupied states, o, ρ = ∑o|Ψo(r)|2 and: veff (r ) = Vnuc + Vhartree[ρ] + Vxc[ρ]

(20)

Vxx − Vyy Vzz

(21)

2.2. Electronic Structure Calculations

First-principles calculations aim to provide a description of the properties of a material using only the fundamental assumptions of quantum mechanics.70 A common starting point is the nonrelativistic electronic Schrödinger equation in the Born−Oppenheimer approximation.70 For a system of electrons and nuclei: H(R)Ψ(r ; R) = E(R)Ψ(r ; R)

(24)

Here, −(1/2)∇2 is the operator that corresponds to the kinetic energy of the noninteracting electrons, Vnuc is the electrostatic interaction with the nuclei, and Vhartree[ρ] is the mean-field electrostatic interaction with the other particles in the system. All of these terms can be computed exactly. The final term Vxc[ρ] is known as the exchange-correlation contribution and represents the many-body interactions. The true functional form of Vxc[ρ] is not known, and to make progress simple forms have been developed that satisfy known physical constraints. The simplest exchange-correlation functional is the local density approximation (LDA).75 The LDA takes a piecewise approach to compute the exchange-correlation energy. Each point in space gives a contribution that is equal to the exact exchange-correlation energy of a uniform electron gas having the density of the given point. Given that the electron density varies widely throughout a material, this would appear to be a very crude approximation. However, it has been found to be a rather good approximation for the energy and structure of solids. For many properties (e.g., the binding energy of molecules), adding terms dependent on the gradient of the density (the generalized gradient approximation, GGA) provides an improvement.74 Numerous GGAs have been proposed; the one by Perdew, Burke, and Ernzerhof (PBE)76 has been most widely used in GIPAW calculations. 2.2.3. Modeling of Solid Materials. Macroscopic samples of solid materials contain uncountably large numbers of electrons, and a direct simulation of such systems would not be possible. Fortunately, many solids possess symmetries at the atomistic level. Translational symmetry means that it is sufficient to consider only the unit cell of material subject to periodic boundary conditions (see Figure 1). This typically reduces the number of atoms (and hence electrons) that must be considered to manageable numbers. However, while the charge density has the periodicity of the lattice (i.e., ρ(r + R) = ρ(r) for all lattice vectors R), the wave functions must be treated with care. For single-particle wave functions, Bloch’s theorem tells us that such objects are only quasi-periodic, that is, Ψnk(r + R) = eik·RΨnk(r), or equivalently:

and ηQ =

(23)

(22)

The electronic Hamiltonian H(R) depends parametrically on the nuclear positions R and describes the kinetic energy of the electrons together with their electrostatic interactions. The many-body wave function Ψ(r;R) is a function of the coordinates of the electrons r and nuclei R in the system. The total energy E provides a means to compare the relative stability of different phases of a material. However, derivatives of the total energy provide a rich variety of properties. For example, the force on an ion is given by the derivative of the total energy with respect to the position of the ion. As shown in section 2.1, the fundamental NMR interaction tensors can be expressed as derivatives of the total energy, and thus can, in principle, be determined through electronic structure calculations. 2.2.1. Hartree−Fock. The Schrödinger equation is far too complicated to solve exactly for more than a few electrons. One approach widely used in the quantum chemistry community is to treat the interaction between the electrons in an average (mean field) way. This leads to the Hartree−Fock set of equations. A range of post-Hartree−Fock methods have been developed to give a more accurate description of electron correlation. Methods such as MP2 add electron correlation as a perturbation to the Hartree−Fock result, while the configuration interaction approach builds the many-body wave function from the single-particle orbitals.71 However, the increased accuracy of these methods comes with a dramatic increase in computational cost. Wave function-based methods beyond Hartree−Fock are difficult and costly to apply to solid materials and have received only limited attention.72 2.2.2. DFT. An alternate approach is to use the electronic density as the fundamental variable. Hohenberg and Kohn73 proved that the total energy of a system can, in principle, be determined by the knowledge of the electronic density.70,74 As the density is a function of position, it is a much simpler quantity than a full many body wave function. A practical approach to so-called density functional theory (DFT) calculations is due to Kohn and Sham.75 They noted that the

Ψ nk (r) = ei k . rukn(r)

unk

(25)

unk(r)

where is a function periodic in the unit cell such that = unk(r + R). A consequence is that physical properties are calculated as an average over all values of the wavenumber k. The only unique values of k lie within the reciprocal unit cell or, equivalently and by convention, within the first Brillouin zone (BZ).77 As most properties in insulators vary smoothly across the BZ, it is possible to truly average with a summation over a set of regular spaced points in the BZ (i.e., a set of k-points). The point group symmetry of the crystal may reduce the number of unique points that must be considered. A common choice is the 5737


Chemical Reviews

Review

ditions. In contrast to GTO, planewaves do not look very much like atomic or molecular orbitals. However, calculations using planewaves will use significantly more basis functions than a calculation based on GTOs. A large calculation might use the order of 105 planewaves to represent the wave functions. To truncate the basis set, the sum is limited to a set of reciprocal lattice vectors contained within a sphere with a radius defined by the cutoff energy, Ecut: 1 k+G 2 Figure 1. (a) The unit cell of the lattice is shown in black. In (b) the atoms that must be included in the simulation are highlighted. (c) The reciprocal lattice of the crystal with the origin (gamma point) highlighted in red. Both the reciprocal unit cell (black) and the first Brillouin zone (BZ) (red) are highlighted. (d) First BZ divided by a uniform mesh of k-points. Because of the point group symmetry, only those within the gray triangle need to be explicitly considered in the calculation.

Figure 2. Examples of supercells for aperiodic structures: (a) defect in a lattice, and (b) surface in a vacuum supercell.

each case, the size of the supercell should be sufficiently large to minimize the interaction between the defect and its periodic images. Including areas of vacuum in the supercell enables the study of lower dimensional systems such as surfaces or nano objects. 2.2.4. Basis Sets. We now need to express the Kohn−Sham equations in terms of a practical algorithm that can be implemented in a computer program. A typical first step is to write the wave functions and charge density as a linear combination of simple mathematical functions. In practice, a finite number of such functions are used, meaning the wave functions are represented approximately, but as more functions are used the representation of the wave function should become more accurate. In quantum chemistry codes, a common choice is a set of localized atomic centered orbitals such as Gaussian-type orbitals (GTO) or Slater-type orbitals (STO).71,74 While localized orbitals can be, and are, used in periodic simulations (for example, in the Crystal80 and Siesta81 codes), the GIPAW technique has so far been implemented using a set of planewaves as basis functions, that is:

∑ ckn(G)ei(k+ G)·r G

≤ Ecut

(27)

Hence, the basis set is defined by the maximum kinetic energy of the waves it contains. This systematically controllable convergence is a key strength of the planewave basis. It is a particularly relevant point for the calculation of NMR parameters, which requires very precise description of several regions of space. For example, to accurately describe the Fermicontact contribution to the J coupling, a precise description of the wave function at the nucleus is needed, while the desciption of the induced orbital currents requires a more extended description. Jensen82 has developed a series of Gaussian basis sets especially optimized to provide a systematically improvable description of J. For planewave functions all that is needed for any property is to increase Ecut until the desired convergence is reached. The value of Ecut required for convergence in a particular system is primarily determined by the atomic species present in the system (or more precisely by the forms of the atomic pseudopotentials used; see section 2.2.5). This arises as the highest energy planewaves contribute to the regions of the wave function closest to the nucleus, which is not affected by the bonding. As an illustration, the Ecut required would increase across the series B to F due to the increasing localization of the 2p wave function. However, the Ecut required for two different systems containing F would be very similar. In any case, the degree of convergence can, and should, be tested to ensure the accuracy of any calculated parameters. An additional advantage to a planewave basis set is that they lend themselves to algorithms, which can exploit multicore and parallel computer architectures: for example, one can distribute planewave coefficients between computer elements. 2.2.5. Pseudopotentials. One potential drawback of a planewave basis is that representing the tightly bound core electrons would require a vast number of planewave coefficients extending to high cutoff energies. Exactly the same issue occurs in describing the oscillations of the valence wave functions close to the nucleus (see Figure 3). These two issues can be addressed by two distinct yet interrelated approximations.83 The frozen core approximation assumes that those electrons labeled as “core” do not take part in chemical bonding, and that they are essentially unchanged in different chemical environments. Gregor et al. have shown the validity of this approximation in the context of magnetic shielding.84 The core electrons can therefore be removed from a simple calculation on a free atom. However, the partitioning of states into core and valence is not always straightforward, and sometimes for accurate results it is necessary to also treat filled atomic shells as valence states (so-called semicore states). Specific examples include the 1s states in Li,85 and the 3p states in 3d transition metals.86 In a similar way, the oscillations of the valence wave functions close to the nucleus do not directly contribute to the bonding. Taking advantage of the frozen core approximation,

set of points determined according to the scheme of Monkhorst and Pack.78 It should be noted that metallic systems require more careful sampling of the BZ.79 In practice, many important materials do not have complete three-dimensional periodicity. The use of supercells allows calculations on systems containing defects (see Figure 2). In

Ψ nk(r) =

2

(26)

If G are a set of reciprocal lattice vectors, then the planewave functions automatically satisfy the periodic boundary con5738


Chemical Reviews

Review

In simple terms, the PAW (Projector Augmented Wave) transformation works by computing the component of a certain atomic-like state in a pseudowave function, and replacing the pseudized component by its all-electron form. This may seem like a rather approximate procedure, but because the atomic states form a good basis for the wave function in the region close to the nucleus, it can be made highly accurate by using multiple projectors. Within the PAW scheme, for an all-electron local or semilocal operator O, the corresponding pseudo-operator Õ is given by: Õ = O +

(29)

As constructed in eq 28, the pseudo-operator Õ acting on pseudowave functions will give the same matrix elements as the all-electron operator O acting on all-electron wave functions.

we can replace both the strong nuclear Coulomb potential and the interaction between the core and valence electrons by a smooth effective potential, which only acts within a defined region around the nucleus. Such an effective potential is known as a pseudopotential in the solid-state physics community and effective core-potential in the quantum chemistry literature. The exact form of the pseudopotential is not uniquely defined. For planewave calculations, the usual requirement is to choose the pseudopotential such that it leads to wave functions that are smooth in the core region, that is, that require only a minimum number of planewaves to represent the wave function. Numerous schemes have been devised to construct pseudopotentials, for example, the one due to Troullier and Martins.87 A more radical scheme for constructing pseudopotentials was devised by Vanderbilt,88 who noted that even smoother wave functions could be obtained if the condition that the norm of the wave functions should give the charge density is relaxed. In this “ultrasoft” scheme, atom-centered charges are included to maintain the total charge in the system. Ultrasoft pseudpotentials add significant complications for code developers, but for users they represent the state-of-art in terms of efficiency. They are of particular benefit when describing localized semicore states. The pseudopotential approximation is valid when dealing with properties that are weighted outside of the core region. However, for the NMR-related properties, the details of the electronic structure close to the nucleus are of prime importance. Fortunately, it turns out that all of this information can be recovered from a calculation based on the use of pseudopotentials. The basis for such an approach was provided by van de Walle and Blöchl,54 who introduced a linear transformation, T, which maps the valence pseudo wave functions |Ψ̃⟩ onto the corresponding all-electron wave functions, |Ψ̃⟩ = T|Ψ̃⟩:

∑ [|ϕR,n⟩ − |ϕR̃ ,n⟩]⟨pR̃ ,n | R, n

|pR̃ , n ⟩[⟨ϕR , n|O|ϕR , m⟩ − ⟨ϕR̃ , n|O|ϕR̃ , m⟩]⟨pR̃ , m |

R ,n,m

Figure 3. Atomic states in Na highlighting the difference between core (1s), semicore (2s), and valence (3s). States from a neighboring Na atom are superimposed (dashed lines).

T=1+

∑

2.3. Calculations of NMR Parameters

We now discuss how electronic structure techniques can be used to compute NMR tensors in solid materials. 2.3.1. Shielding. As noted in section 2.1, the question of computing the magnetic shielding in a diamagnet is essentially that of computing the electronic current induced by an applied field. The electronic changes due to experimentally realizable magnetic fields are small, and so perturbation theory is the appropriate tool. We adopt a notation in which the electronic wave function in the presence of an applied field is expanded as ψ(r) = ψ(0)(r) + ψ(1)(r) + O(B2), where ψ(0)(r) is the groundstate wave function in the absence of the field. The first-order change in the wave function ψ(1)(r) can be expressed as an admixture of unoccupied states, e:

ψ (1)(r) =

∑ aeψe(0)(r)

(30)

e

The current operator, J(r′), is obtained from the quantum mechanical probability current. It can be written as the sum of diamagnetic and paramagnetic terms: J(r′) = J d (r′) + J p (r′)

(31)

1 A(r′)|r′⟩⟨r′| c

(32)

J d (r′) =

Jp (r′) = −

p|r′⟩⟨r′| + |r′⟩⟨r′|p 2

(33)

To first order in B, the induced orbital current, j (r), is (1)

p (1) j(1)(r′) = 4 ∑ Re[⟨Ψ(0) o |J (r′)|Ψ o ⟩] o d (0) + 2 ∑ ⟨Ψ(0) o |J (r′)|Ψ o ⟩

(34)

o

(28)

The first-order change in the wave function by:

Here, |φR,n⟩ and |ϕ̃ R,n⟩ are, respectively, all-electron and pseudo partial waves derived from an isolated atomic calculation, and ⟨p̃R,n| are a set of projectors such that ⟨p̃R′,n|ϕ̃ R,m⟩ = δR,R′δn,m. Each projector and partial wave is an atomic-like function centered on an atomic site R, and the index n refers to the angular momentum quantum numbers and to an additional number, used if there is more than one projector per angular momentum channel. Both projectors and partial waves are obtained during the initial construction of the pseudopotential.

|Ψ(1) o ⟩ =

∑ e

|Ψ(1) o ⟩

is given

(0) |Ψ(0) e ⟩⟨Ψ e | (1) (0) H |Ψ o ⟩ = .(εo(0))H (1)|Ψ(0) o ⟩ εo − εe

(35)

where H(1) = (1/(2c))(p·A + A·p). Using the symmetric gauge for the vector potential, A(r) = (1/2)B × r, we arrive at the following expression for the induced current: 5739


Chemical Reviews j(1)(r ′) = 4

1 2c

Review

this problem45 is to consider the response to a magnetic field with a finite wavelength that is, B = sin(q·r)q̂. In the limit that q→0, the uniform field result is recovered. For a practical calculation, this enables one to work with cell-periodic functions, at the cost that a calculation at a point in the BZ k will require knowledge of the wave functions at k ± q (i.e., six extra calculations for the full tensors for all atomic sites). A complete derivation was presented in refs 45,92 leading to the final result for the current:

p (0) (0) ∑ Re[⟨Ψ(0) o |J (r ′).(εo )r × p|Ψ o ⟩] o

1 − ρ(r ′)B × r ′ 2c

(36)

(0) 2∑o|⟨Ψ(0) o |r′⟩⟨r′|⟨Ψo ⟩

where ρ(r′) = with the factor of 2 accounting for spin degeneracy. Equation 36 shows a division of the current into two terms, the first paramagnetic, the second diamagnetic. However, this derivation relies on a specific choice for A(r) which, unlike B, is not uniquely defined. In particular, one can choose A(r) = (1/2)B × (r − d) where the constant vector d is known as the gauge origin. Different choices for d and A(r) will alter the balance between the paramagnetic and diamagnetic currents. Only the total current is well-defined. The paramagnetic term converges more slowly than the diamagnetic term with respect to the size of an atomic orbital basis. This has an important practical consequence for calculations using finite basis sets, and in particular those based on localized orbitals. Two chemically identical sites at different distances from the gauge origin will have the same total shielding but a different decomposition into diamagnetic and paramagnetic terms. For an infinite basis, this would not be a problem, but for a finite basis, the two sites will have different computed shieldings. This is known as the “gauge origin problem”. This was a severe problem for early calculations of magnetic shielding, and several techniques have been developed to restore gauge invariance. Two common approaches are GIAO (Gauge-Including Atomic Orbitals)89,90 and IGLO (Individual Gauges for Local Orbitals).91 A planewave basis does not directly suffer from the gauge-origin problem; however, the PAW transformation introduces a set of localized functions and hence reintroduced the gauge dependence. To address this problem, Pickard and Mauri92 introduced a fielddependent transformation operator TB, which, by construction, imposes the translational invariance exactly: TB = 1 +

j(1)(r ′) = lim

q→0

2 cNk

S(r′, q) =

⎡

∑ ∑ Re⎢⎣ 1 ⟨uo(0),k |Jkp,k+ q (r′)

i=x ,y ,z o,k

i

i

⎤ .k + q (εo , k )B × u i ·(p + k)|uo̅ (0) , k ⟩⎥ i ⎦

(40)

qi = qui, Nk is the number of k-points included in the BZ summation and Jkp, k + q = − i

(p + k)|r ′⟩⟨r ′| + |r ′⟩⟨r ′|(p + k + q i) 2

(41)

Equivalent expressions valid when using separable normconserving pseudopotentials are given in ref 92 and for ultrasoft potentials in ref 93. We now highlight alternative methods to GIPAW for computing shielding tensors in crystals. Rather than use the long-wavelength approach to apply the position operator, Sebastiani et al. developed a method using localized Wannier orbitals.94This has been implemented in the planewave CPMD code.47 More information concerning the approach of Sebastiani can be found in refs 46,48−53. The CP2K program has a recent implementation using the Gaussian and Augmented Planewave method.95 In both cases, the shieldings are still calculated using perturbation theory (linear response). Very recently, calculation of shielding tensors has been implemented in the Linear Augmented Planewave + local orbitals (LAPW+lo) code Wien2k. The method uses perturbation theory and the long-wavelength approach to applying the position operator.96 The LAPW+lo method represents the “gold standard” in terms of numerical accuracy and will provide an important route to test the pseudopotentialbased GIPAW approach. Recently, a method that avoids linear response, the so-called “converse approach”, has been developed and implemented in the Quantum-Espresso package.97 In this approach, the shielding tensors are obtained from the macroscopic magnetization induced by magnetic point dipoles placed at the nuclear sites of interest. The “converse approach” is made possible by the recent developments that have led to the Modern Theory of Orbital Magnetization,98−100which provides an explicit quantummechanical expression for the orbital magnetization of periodic systems in terms of the Bloch wave functions and Hamiltonian, in the absence of any external magnetic field. The converse and linear response approaches should give the same shielding tensors if the same electronic structure method is used (e.g., the LDA). The main advantage of the converse approach is that it can be coupled easily to advanced electronic structure methods and situations where a linear-response

R, n

(37)

The resulting approach is known as the Gauge Including Projector Augmented Wave (GIPAW) method. We might hope to apply eq 36 directly to solid-state systems; however, for an extended system, there is an obvious problem with the second (diamagnetic) term of eq 36; the presence of the position operator r will generate a large contribution far away from r = 0, and the term will diverge in an infinite system. The situation is saved by recognizing that an equal but opposite divergence occurs in the first (paramagnetic) term of eq 36, and so only the sum of the two terms is well-defined. Through the use of a sum-rule,45,92,93 we arrive at an alternative expression for the current: 1 2c

(39)

where

∑ ei/2c r·R × B[|ϕR,n⟩ − |ϕR̃ ,n⟩]⟨pR̃ ,n |e−( i/2c)r·R × B

j(1)(r ′) = 4

1 [S(r′, q) − S(r′, −q)] 2q

p (0) (0) ∑ Re[⟨Ψ(0) o |J (r ′).(εo )(r − r ′)×p|Ψ o ⟩] o

(38)

. (ε(0) o )

In an insulator, the Green function is localized, and so j(1)(r′) remains finite at large values of (r − r′). At this point, there still remains the question of the practical computation of the current. For reasons of efficiency, it is desirable to work with just the cell-periodic part of the Bloch function. Equation 38 is not suitable for such a calculation as the position operator cannot be expressed as a cell-periodic function. One solution to 5740


Chemical Reviews

Review

of the relevant interactions in paramagnetic solids at a consistent computational level. In metallic systems, the electronic spin also plays an important role as the external field will create a net spin density (Pauli susceptibility), which will in turn create a magnetic field at the nucleus. This additional contribution to the magnetic shielding is known as the Knight shift. The linearresponse GIPAW approach has been extended115 to compute the magnetic shielding and Knight shift in metallic systems, although there have been few applications to date. In Mössbauer spectroscopy, two of the key observables are the quadrupole coupling and the isomer shift. The quadrupole coupling is given by the electric field gradient at the site, and the isomer shift arises from a hyperfine interaction. An implementation within the ABINIT code107 has been described by Zwanziger.116 The two key interaction tensors in ESR spectroscopy are the hyperfine tensor and the g tensor. The g tensor arises from the interaction of the electronic spin with external magnetic field. This term plays a somewhat similar role to the shielding in NMR; induced electronic currents in the sample modify the g tensor from its vacuum value. The GIPAW approach has been used to compute g tensors in several crystalline materials including defects in α-quartz and zirconia.117,118 More recently, the “converse approach” has been used to compute the g tensor.119 The specific advantage of this formulation is that it can be applied in cases where linear-response theory fails, for example, radicals and defects with an orbital-degenerate ground state or those containing heavy atoms. 2.3.3. Overview of GIPAW Method. At this point, we briefly pause to summarize how a typical calculation might be performed. Generally, one starts from an initial guess at the crystal structure of the material with a specification of the unit cell and the positions of atoms within this cell. This may come from a diffraction-based study or from previous computational prediction. A minimization of the ground-state energy will provide the forces on all of the atoms and the stress on the unit cell. If any of these are “large”, it might be desirable to perform a geometry optimization, that is, to allow the positions of the atoms (and potentially the unit cell dimensions) to vary so as to minimize the forces (stress) in the system. Having chosen the appropriate geometry, a GIPAW calculation provides the desired NMR tensors. It should be stressed that here, and in the remainder of this review, we use GIPAW to refer to the ensemble of techniques used to compute NMR tensors in the planewave pseudopotentials approach. There are many fully functional periodic DFT codes available. The GIPAW approach was originally implemented for norm-conserving pseudopotentials in the PARATEC code.120 This was a development code that is not generally widely available. An implementation for both norm-conserving and ultrasoft pseudopotentials was introduced into the CASTEP code,121 and this functionality is now also available in the GIPAW/PWSCF code within QUANTUM ESPRESSO.122 Both codes can run on desktop PCs and on massively parallel supercomputers. As was already mentioned above, calculations of EFG and Mössbauer isomer shifts can be performed in the ABINIT code.107 Generally speaking, the most intensive part of the calculation will be the optimization of the crystal structure. The calculation of the EFG tensors typically takes only a very small fraction of this time. Many factors influence the time taken for a calculation. As an indication, a large calculation on a molecular crystal with approximately 700 atoms took 1.5 h to find the ground-state

formulation is cumbersome or unfeasible. For example, in the case of DFT+U101 or hybrid functionals, the converse method should provide a convenient shortcut from the point of view of program coding. It should be noted, however, that a converseapproach calculation should be repeated for each inequivalent NMR site, while in the linear response approach a single calculation is sufficient to obtain the shielding tensors of all sites. Finally, a very recent approach implemented in the ADFBAND code uses gauge-including atomic orbitals and computes the shielding tensor as the second derivative of the total electronic energy with respect to a nuclear magnetic moment and an external magnetic field.102 2.3.2. Other NMR Parameters. 2.3.2.1. Electric Field Gradients. The computation of electric field gradient tensors is less demanding than either shielding or J coupling tensors as it requires only knowledge of the ground-state charge density. This is directly available in all-electron (i.e., non pseudopotential) codes. For example, the Wien series of codes44 employing the Linear Augmented Planewave method (LAPW) has been widely used and shown to reliably predict EFG tensors in solids.103 It should be noticed that in some studies GIPAW is combined with LAPW for EFG calculations, showing that both approaches are in very good agreement.104,105 The equivalent formalism for pseudopotential techniques using the PAW approach to reconstruct the all-electron charge is reported in ref 106. It should also be mentioned here that properties at the nuclei such as the EFG (but also the Mössbauer isomer shifts; see below) can be calculated in the PAW formalism within the ABINIT code.107 2.3.2.2. J Coupling. Several quantum chemistry packages provide the ability to compute J coupling tensors in molecular systems (see ref 31 for a review of methods). An approach to compute J tensors within the planewave pseudopotential approach has recently been developed.108 A comprehensive review of the formalism and its applications is provided in ref 109. As with Pickard and Mauri’s approach for computing the shielding tensor, density functional perturbation theory is used to compute the response quantities (here the induced magnetization density in addition to the current). PAW is needed to account for the use of pseudopotentials. While calculations of the shielding and EFG tensor can be performed using the primitive unit cell of the crystal, this is not always the case for the calculation of J in the above formulation. The perturbing nuclear dipole breaks the translational symmetry, and for crystals with small unit cells, it is necessary to construct a supercell of sufficient size to inhibit the interaction of the dipole with its periodic images. 2.3.2.3. Paramagnetic Coupling and Other Spectroscopic Observables. If a material contains an unpaired electron, then this net electronic spin can create an additional magnetic field at a nucleus via Fermi contact and spin-dipolar mechanisms. Paramagnetism introduces several difficulties from the solidstate NMR point of view: for example, the resonances can exhibit significant broadening. However, NMR has been used to analyze local magnetic interactions, for example, in manganites110 and lithium battery materials.111 There are several reports of calculations of NMR parameters in paramagnetic systems, for example, paramagnetic shifts of 6 112 Li and EFGs of layered vanadium phosphates.113 However, to the best of our knowledge, unlike the case of paramagnetic molecules,114 there is currently no methodology to predict all 5741


Chemical Reviews

Review

Burke, and Ernzerhof (PBE).76 This is a relatively simple functional whose form has been obtained using purely physical considerations. Usually, the agreement with experiment is reasonably good; a rough rule of thumb is that errors in the chemical shift are within 2−3% of the typical shift range for that element. There are, however, some notable exceptions: several groups have shown123,124 that while present functionals can predict the trends in 19F chemical shifts, a plot of experimental against calculated shifts has a slope significantly less than 1 (see section 4.1). Another example is the calculation of 17O chemical shifts125 in calcium oxide and calcium aluminosilicates. There are significant errors in the 17O shifts due to the failure of GGA PBE to treat the unoccupied Ca 3d states correctly. In ref 125, it was found that a simple empirical adjustment of the Ca 3d levels via the pseudopotential was sufficient to bring the 17O chemical shifts into good agreement with experiment. However, in both cases, it is clear that current GGAs do not describe all of the relevant physics. The “converse approach” to compute NMR parameters should provide an easy route to include exactexchange in the calculation of magnetic shielding, and it will be interesting to see if this can improve the treatment of these known difficult cases. 2.3.4.3. Relativistic Effects. The study of systems containing elements with large atomic numbers requires the treatment of the effect of special relativity as both low lying core-states and valence states in the vicinity of the nucleus will attain high values of kinetic energy. The behavior of such states is more accurately described by the Dirac equation. Scalar relativistic (i.e., spin−orbit free) effects can be accounted for in the construction of pseudopotentials, for example, by solving the radial Dirac equation and performing the appropriate average over spin−orbit split states. However, Yates et al.126 showed that, in this approach, it was also necessary to modify the GIPAW operators to account for the relativistic nature of the valence electrons close to the core region. The form of the operators was derived using the zeroth-order regular approximation (ZORA) method for expanding the Dirac Hamiltonian. The importance of relativity grows as one decends the periodic table. Quantum chemical studies have highlighted the importance of spin−orbit-induced effects (see, e.g., ref 127). However, to date there has not been a fully relativistic implementation of the GIPAW approach. Light atoms bonded to heavy atoms have been found to have a significant proportion of their magnetic shielding arising from spin− orbit-induced effect.127 The spin−orbit coupling on the heavy atom creates a spin-polarization, which in turn influences the magnetic field at the light atom nucleus via a hyperfine mechanism. 2.3.4.4. Geometry and Thermal Motion. NMR spectra are commonly obtained at room temperature. Given that firstprinciples calculations typically use a static configuration of atoms (e.g., obtained from diffraction), this raises questions about the influence of thermal motion on NMR spectra, even if it is thought that there are no specific motional processes (such as exchange) involved. Dumez and Pickard128 have examined two ways of including motional effects: by averaging NMR parameters over snapshots taken from molecular dynamics simulations, and by averaging over vibrational modes (as previously used by Rossano et al.129 to study the effects of temperature on 17O and 25Mg NMR parameters in MgO). They found the effects of zero-point motion to be significant as well as the influence of thermal effects on shielding anisotropies. Webber et al.130 measured the

energy, 75 h to optimize the geometry, and 20 h to obtain the EFG and shielding tensors (using 8 dual quad-core 2.27 GHz Intel Xeon E5520 CPUs, which corresponds to a total of 64 cores). 2.3.4. Discussion of “Accuracy”. When considering the accuracy of first-principles calculations, it is helpful to first make a distinction between computational and physical approximations. Physical approximations include the use of approximate exchange-correlation functionals, the frozen core approximation, neglect of thermal motion, etc. The quality of these approximations can typically only be judged by careful comparison to experimental data. On the other hand, computational approximations, such as the use of finite basis sets, can be systematically controlled. 2.3.4.1. Computational Approximations. The important concept here is that of convergence: any computational approximations should be tested such that their impact on the final calculated parameters can be determined. For example, when using a finite basis set to represent the electronic wave function and related properties, the number of basis functions should be increased until the point at which adding additional functions does not change the property of interest by a significant amount. Convergence can be achieved using a planewave basis by simply increasing the energy of the maximum planewave included in the basis (see Figure 4). In

Figure 4. Convergence with planewave cutoff energy of the 17O isotropic magnetic shielding in a O(SiH3)2 cluster derived from the αquartz structure. The all-electron result is taken from ref 106. Results for both a norm conserving and an ultrasoft pseudopotential are shown, demonstrating the faster convergence of the latter. Adapted from ref 93.

a similar sense, the spacing between the k-points used to take averages of the properties over the electron BZ should be decreased to obtain convergence of the final results. The use of pseudopotentials is the third important approximation relevant for the GIPAW approach. A pseudopotential for a given element is tested by careful comparison to all-electron results, for example, the comparison of magnetic shieldings with molecular results using quantum chemistry software92 and the same for EFG with the LAPW approach.105 2.3.4.2. Exchange-Correlation Functionals. An important physical approximation lies in the choice of the exchangecorrelation functional. The better is the functional, the closer the result will be to the exact solution of the Schrödinger equation. The majority of GIPAW calculations of NMR parameters have employed the GGA functional due to Perdew, 5742


Chemical Reviews

Review

changes in 1H and 13C chemical shifts (in the case of βmaltose) in the range from 348 to 248 K. By extrapolating the results to 0 K, the change in 1H shifts for the hydroxyl protons with respect to room temperature was estimated to be 0.5 ppm. NMR can be used to study motional processes in solids. One commonly used technique is the use of deuterium NMR. 2H has spin I = 1, and the magnitude of its quadrupolar coupling (typically 250 kHz) makes it suitable to study motional processes on the micro- and millisecond time scale. Griffin et al.131 have used 2H solid-state NMR to study the dynamic disorder of hydroxyl groups in hydroxyl-clinohumite. In this material, the deuteron can exchange between two crystallographic sites. By combining first-principles calculations, a simple model of the effect of motion on the NMR linebroadening, and experimental 2H NMR spectra, it was possible to obtain the activation energy for the exchange processes. All of the examples cited above will be described in depth in section 4.

3. EXPERIMENTAL SOLID-STATE NMR: METHODS AND IMPLEMENTATION Although the interactions that affect nuclear spins can provide a range of information on the local environment, for powdered solids their anisotropic nature limits the sensitivity and resolution of the NMR spectra observed, particularly when more than one nuclear species is present. A number of experimental approaches are commonly used, therefore, both to increase resolution through the removal of anisotropic broadening and to improve the sensitivity of the resulting NMR spectra. Once high-resolution spectra are obtained, it is also then possible to utilize experiments that “reintroduce” the anisotropic interactions (usually in a selective way), enabling (i) the accurate measurement of interactions and the extraction of relevant structural parameters and (ii) the transfer of magnetization between nuclear spins to probe structure and bonding. Here, we will briefly review some of the commonly used experimental approaches to aid the understanding of the examples described in the following section. However, more detailed descriptions of individual techniques can be found in the references given and in the following reviews and references therein.1,2,5,132−134

Figure 5. 31P NMR spectra of ammonium dihydrogen phosphate recorded under (a) static conditions and at MAS rates of (b) 1.25 kHz and (c) 10 kHz.

In some cases, it may not be possible to spin at sufficiently rapid rates to remove the anisotropic interactions. This can be a particular problem for dipolar couplings involving highly abundant nuclei with high γ, such as 1H or 19F. In this case, MAS can be combined with an approach termed dipolar “decoupling”. For heteronuclear dipolar couplings, e.g., 1H/13C, decoupling involves the irradiation of one spin (e.g., 1H) during the acquisition of the FID for the second (e.g., 13C).132 In its simplest form (termed continuous wave decoupling), high power radiofrequency (rf) radiation is applied throughout acquisition, although over the years significant effort has been invested in the implementation of ever more elaborate decoupling schemes. In many cases, modulating the duration and phase of the applied rf has been shown to enable more efficient removal of the heteronuclear dipolar coupling.136−140 The removal of homonuclear dipolar couplings (e.g., 1H−1H or 19F−19F) poses more of a practical challenge, as both spins within the dipolar-coupled pair are affected by the application of any rf pulses during the acquisition of an FID.132 To remove such interactions, it is necessary to combine rf irradiation with sample rotation, in an approach known as CRAMPS (combined rotation and multiple pulse spectroscopy). The high-resolution spectrum can be acquired either in the indirect dimension of a two-dimensional experiment or in the direct dimension through the use of “windows”, which allow point-by-point acquisition of the FID.141 Recent advances in this area include the introduction of DUMBO142,143 and PMLG144 techniques. It should be noted that many of these approaches introduce a scaling factor into the spectrum that has to be accounted for when determining chemical shifts. The technique of cross-polarization (CP) is routinely used in solid-state NMR of nuclei with low natural abundance, such as 13 C or 15N. In this approach, magnetization is transferred usually from a high-γ nucleus with high natural abundance, e.g.,

3.1. Basic Experimental Techniques (Spin I = 1/2)

The standard technique for removing anisotropic broadening in powdered solids is MAS;135 sample rotation about an axis inclined at 54.736° to the external magnetic field. Providing the rotation is sufficiently rapid, the resulting time-dependence enables the removal of anisotropic interactions such as dipolar couplings and CSA (i.e., those with an orientation dependence proportional to P2(cos θ) = 1/2(3 cos2 θ − 1)). For many spin I = 1/2 nuclei, this is able to produce truly high-resolution or “isotropic” spectra, containing sharp, narrow resonances at the isotropic or average chemical shift, as shown in Figure 5. The rate of sample rotation required is dependent upon the magnitude of the interactions present, with typical rotation rates of 5−70 kHz available using commercial probeheads. If the rotation is not sufficiently rapid, the anisotropic line shape will break up into a series of “spinning sidebands”, spaced equally at the rotation frequency from the isotropic resonance, as shown in Figure 5. Although this may result in spectral overlap if more than one resonance is present, the intensities of the spinning sideband manifold can also provide information relating to the interaction tensors. 5743


Chemical Reviews

Review

H, through the dipolar coupling between the two spins.145−148 The main aim of cross-polarization is to improve sensitivity, with a theoretical intensity gain of γI/γS, assuming that the abundance of I is much greater than that of S. CP transfer only occurs when a matching condition (termed the Hartmann− Hahn condition145) is satisfied, that is, ω1I = ω1S, by the choice of appropriate rf field strengths. If CP is performed under MAS, this condition is modified to account for sample rotation giving ω1I = ω1S + nωR, where n is typically 1 or 2.147,148 In addition to the improvement in sensitivity, it is typically found that the increased relaxation rate of the more abundant I spin also enables the experiment to be carried out more rapidly, producing an additional gain in signal per unit time. CP transfer can only occur when the two spins involved are close in space. Therefore, the S-spin spectrum is inherently “edited”, according to the proximity to spin I, providing structural information on the system under study. 1

3.2. Quadrupolar Nuclei (Spin I > 1/2)

Nuclei with spin quantum number I > 1/2 possess a nuclear electric quadrupole moment, which can interact with any electric field gradient generated at the nucleus by the surrounding charges within the system.2,133 The perturbation of the Zeeman energy levels by the quadrupolar interaction results in two types of transitions within the NMR spectrum: the central transition (CT) unaffected by the quadrupolar interaction to a first-order approximation, and the satellite transitions (ST) with a frequency that depends upon the magnitude of the quadrupolar coupling. In a powdered sample, these latter transitions can be broadened over many MHz, and, consequently, spectral acquisition is typically only concerned with the CT. In many cases, the magnitude of the quadrupolar interaction is so large that its effect upon the spectrum is best described by a second-order perturbation of the Zeeman energy levels; both the CT and the ST are affected by the second-order quadrupolar interaction, with the CT now both anisotropically broadened and shifted from the isotropic chemical shift. Because of the more complex angular dependence of the second-order quadrupolar interaction (dependent upon both P2(cos θ) and P4(cos θ) = 1/8(35 cos4 θ − 30 cos2 θ + 3)), it is not possible to remove the anisotropic broadening by MAS, and powder-pattern lineshapes remain even when spinning rapidly. Although the lineshapes can provide information on the magnitude (CQ) and asymmetry (ηQ) of the quadrupolar tensor, it can be difficult to extract this information when more than one spectral resonance is present. As shown in Figure 6a, a27Al (I = 5/2) MAS NMR spectrum of AlPO4-14 exhibits a complex resonance, resulting from the overlap of four quadrupolar broadened lineshapes.2,149 An enduring theme in NMR of quadrupolar nuclei is the quest for high-resolution spectra and the complete removal of the second-order quadrupolar broadening. The earliest approaches utilized composite sample rotation techniques to achieve this, that is, rotation around two different angles. In DOR (double rotation), the two rotations take place simultaneously, with an inner rotor spinning inside an outer rotor with a larger diameter.150,151 The axes of rotation are chosen such that the conditions P2(cos θ) = 0 and P4(cos θ) = 0 are satisfied simultaneously, that is, 54.736° and 30.56°. A significant advantage of DOR is that high resolution (i.e., isotropic spectrum) can be acquired in real time in a simple one-dimensional experiment. However, the low filling factor and low rf field strength available limit sensitivity, and the

Figure 6. 27Al (I = 5/2) (a) MAS and (b) two-dimensional triplequantum MAS NMR spectrum and corresponding isotropic projection of AlPO4-14, recorded at B0 = 14.1 T.2,149 Reproduced with permission from refs 2,149. Copyright 2009 and 2008 PCCP Owner Societies.

relatively slow spinning speed of the outer rotor (∼1−2 kHz) can result in complicated sideband manifolds. In dynamic angle spinning (DAS) experiments, rotation around the two different angles occurs sequentially, with two quadrupolar broadened powder pattern lineshapes, correlated in a two-dimensional experiment.152,153 The DAS experiment may be performed with a MAS detection period by means of a further switch in rotor angle to 54.736°. A severe limitation of the DAS experiment involves the storage of the magnetization in a population state during the periods where the rotor angle is switched, making long T1 relaxation times a requirement. In contrast to DAS and DOR where the quadrupolar broadening is removed by increasing the complexity of the physical manipulation of the sample, MQMAS (multiplequantum MAS)133,154 and STMAS (satellite-transition MAS)155,156experiments combine simple (MAS) sample rotation with manipulation of the nuclear spins. Each experiment exploits the similarity of the form, but not the exact magnitude, of the second-order quadrupolar broadening for different transitions within the spin system to remove the fourth-rank anisotropic components of the interaction within a twodimensional correlation experiment. For MQMAS, multiplequantum coherences (typically triple-quantum) are excited and evolve in t1, before conversion to observable CT coherences. A similar experiment is performed in STMAS, but with excitation and evolution of the satellite transitions in the t1 period. In the most recent experiments, a double-quantum filter (DQF) is employed to remove the uninformative autocorrelation diagonal.157 Fourier transformation yields a two-dimensional spectrum containing (after an appropriate shearing transformation) ridge-like lineshapes lying parallel to F2. The F1 dimension contains an isotropic, high-resolution spectrum, as shown in Figure 6b, a 27Al triple-quantum MAS spectrum of AlPO-14, where the four crystallographically distinct Al species 5744


Chemical Reviews

Review

can now be easily resolved.2,149 The advantage of MQMAS and STMAS experiments over DAS and DOR is that they can be implemented using a conventional probehead, and MQMAS in particular has seen widespread application. The main limitation of MQMAS is the relatively poor sensitivity, although a number of sensitivity enhancement techniques, notably FAM, DFS, and SPAM, have been employed to attempt to address this problem.158−161 The sensitivity of STMAS is typically significantly better (often up to a factor of ∼6−8), but it does require that the pulses are very accurately timed, the spinning speed is very stable, and that the spinning angle is precisely (i.e., to within ±0.002°) adjusted to remove the large first-order quadrupolar interaction that is also present in the ST.156,162

are determined by the inverse of the coupling present, that is, typically 100 ns to 10 μs for the first-order quadrupolar coupling, and 100 μs to 1 ms for lineshapes broadened by the second-order quadrupolar coupling. In particular, 2H (I = 1) NMR has been widely used as a sensitive probe of dynamics in solids, while the possibility of selective isotopic labeling provides good resolution and spectral simplification in many cases.2,169 3.4. Measurement of Interaction Parameters

Detailed information about structure and dynamics in the solid state can be extracted if the anisotropic interactions that affect the NMR spectrum can be accurately measured. For simple materials, this can be possible from a static spectrum, but in many cases the MAS required to provide site resolution will average or remove much of the desired information. In this case, “recoupling” sequences that use rf pulses, typically synchronized in some way with the rotation of the sample, can be employed to reintroduce the interactions in a selective and controlled way. Of particular interest is the dipolar coupling, which, due to its strong dependence upon the distance between two spins (proportional to r−3), encodes information on spatial proximity and internuclear distances. A number of techniques have been proposed to reintroduce this interaction into the NMR spectrum.1,132 For example, REDOR (rotational echo double resonance) experiments involve the use of a simple spin echo for one spin (I), with the application of a series of π pulses, synchronized with the sample rotation to a dipolar coupled spin, S.170 This reintroduces the heteronuclear IS dipolar coupling, causing those I spins close to an S neighbor to “dephase”, that is, decrease in intensity. The magnitude of the dephasing, measured as a function of echo duration, can then be used to provide information on the internuclear distance between I and S. Double-quantum coherence, that is, the creation of coherences between two spins, in techniques such as BABA (back-to-back pulses) has also provided information on internuclear distances, particularly for high abundant nuclei such as 1H, 19F, and 31P.171 Recently, a theoretical framework to describe recoupling has been introduced, with pulse sequences described in terms of their symmetry, and has enabled a more focused design of particular sequences for specific applications. These “symmetry-based recoupling sequences” have found widespread use in the twodimensional experiments discussed in a later section.6 In addition to the information provided by the dipolar coupling, measurement of the CSA can also provide a useful structural probe. As this interaction is removed by fast MAS, slower sample rotation is required, although this can become more difficult to use as the number of species in the spectrum increases, with potential overlap of centerbands and sidebands. A number of two-dimensional methods have been proposed to overcome this problem, resulting in spectra where the isotropic resonances in one dimension are correlated with their corresponding sideband manifold in the second. Many of these methods use several π pulses applied at timings synchronous with the MAS (as in the phase-adjusted spinning sideband (PASS) approach), allowing the manipulation of the phase and intensity of the MAS sidebands.172−174 Twodimensional PASS-type experiments offer a simple, efficient method for the site-specific determination of the CSA, provided a sufficient number of sidebands are observed. However, the measurement of small anisotropic interactions remains a significant challenge. In many cases, the MAS rates required

3.3. Wideline NMR

While the majority of solid-state NMR spectra are acquired using some form of sample rotation, for example, MAS, to improve resolution, in some cases spectra are acquired for static samples.2,163 This may be the case when information on an anisotropic interaction is required or if an interaction is so large that it is simply not possible to spin sufficiently rapidly to remove it. However, acquisition of the broad lineshapes encountered for static samples does pose practical challenges. Typically echo sequences are required to enable acquisition of undistorted lineshapes (i.e., to overcome “deadtime” problems), with the nature of the refocusing pulse dependent on the interaction that is to be refocused.164 For I = 1/2 nuclei, static spectra will provide information on the CSA interaction. For nuclei with spin quantum number I > 1/2, static spectra are typically used to extract information on the quadrupolar interaction, although it is often necessary to include contributions from the CSA tensor when fitting the spectral line shape. A range of techniques can also be used to overcome the sensitivity challenges associated with the acquisition of broad lineshapes. For example, Carr−Purcell− Meiboom−Gill (CPMG) echo trains can be used in acquisition, resulting in a manifold of “spikelets” increasing the peak height sensitivity.165 In addition, it is possible to increase the sensitivity of CT spectra for quadrupolar spins with halfinteger spin quantum number by manipulation of the ST using approaches such as DFS or FAM pulses, prior to the initial pulse.161 (Note it is also possible to use these approaches to improve the sensitivity of MAS and DOR spectra for halfinteger spins.) In addition to the sensitivity challenges for very wide lines, the limited excitation bandwidth of the rf pulses can also pose a problem. This is usually overcome by implementing the so-called variable offset cumulative spectroscopy (VOCS) approach,166,167 whereby the spectrum is acquired in a number of separate experiments, changing the transmitter frequency in each step. It has recently been shown that the use of WURST pulses alleviates this problem and enables much broader lines to be acquired in a single acquisition.168 The acquisition of broadened lineshapes is often used to investigate dynamics in the solid-state. In general, NMR spectra are sensitive to dynamics on a variety of time scales; fast motions are typically measured by relaxation methods, while much slower motion can be probed by exchange experiments. However, it is motion on time scales intermediate to these extremes that typically affects the width and shape of anisotropically broadened lines. Because of its magnitude, it is the quadrupolar interaction that is the focus of many studies. The dynamic time scales that are probed by line shape analysis 5745


Chemical Reviews

Review

Figure 7. (a) Rotor-synchronized 19F (14.1 T) through-space DQ MAS NMR spectrum of a fluorinated magnesium silicate mineral, recorded using one cycle of BABA recoupling, with a MAS rate of 30 kHz. Adapted with permission from ref 124. Copyright 2010 American Chemical Society. (b) Two-dimensional 27Al/31P MQMAS-J-INEPT NMR spectrum of calcined-dehydrated AlPO4-14, recorded (at B0 = 9.4 T) with τ = 2.5 ms. Adapted with permission from ref 149. Copyright 2008 PCCP Owner Societies.

interaction parameters can in some cases also be measured by rotation of the sample around an angle that is slightly offset from the magic-angle. In this approach, termed “off-MAS”, a small component of the interaction is reintroduced, while high resolution is largely maintained. This method has been used to measure both dipolar and J couplings between 13C−13C spin pairs in labeled amino acids, with the advantage that the relative signs of the dipolar and J couplings can also be determined.191,192 Off-MAS experiments have also been used to correlate chemical shift and quadrupolar interactions for 2H species in deuterated L-histidine.193

to achieve the necessary number of sidebands are either prohibitively low or are insufficient to remove the dipolar interactions. More recently, a number of experiments have been proposed that achieve an effect termed “CSA amplification”, where the conventional MAS spectrum (MAS rate, ωr) in the direct dimension is correlated with a sideband manifold in the indirect dimension with an effective MAS rate of ωr/N.175−177 These approaches can be thought of as amplifying the CSA by a factor N, enabling an accurate measurement to be obtained at faster MAS rates. To date, these approaches have been evaluated for, and applied to, 13C NMR of sugars, amino acids, and other small organic molecules, and 31P NMR of simple inorganic phosphates.176,177 The scalar or J coupling interaction can provide a highly sensitive probe of chemical bonding and geometry. This interaction differs from the dipolar coupling and CSA, as it has an isotropic component that survives MAS. However, J couplings can often be overlooked in solid-state NMR due to their weak intensities as compared to other larger anisotropic interactions. Indeed, direct measurement in the solid-state NMR spectrum is often not possible, because the observed MAS line width is usually larger than the magnitude of the splitting arising from the J coupling. However, a spin−echo experiment can be used to refocus a number of anisotropic interactions and other broadening mechanisms, to give transverse dephasing times that are often much longer than the decay constants of the directly acquired NMR signal. This results in so-called “refocused” line widths that can be sufficiently narrow as to enable measurement of splittings due to J couplings.178,179 J couplings can be measured either by analytical fitting of the modulated spin−echo intensity or by Fourier transformation of the spin−echo dimension to yield a two-dimensional “J-resolved” spectrum. The J-resolved approach has been used widely for the measurement of homonuclear J couplings between pairs of spin-1/2 nuclei.111,180−187 The J-resolved approach has also been used to measure heteronuclear J couplings both between spin-1/2 nuclei and those involving quadrupolar nuclei.188−190 The methods for the measurement of interaction parameters described above are based upon the assumption of perfect averaging of the anisotropic interactions by MAS. However,

3.5. Two-Dimensional Correlation Experiments

While high-resolution NMR spectra reveal information about the number of different species, their relative intensities, and respective chemical environments, in many cases more detailed structural information can be obtained using two-dimensional correlation experiments. In these experiments, magnetization is transferred between nuclear spins, through either the dipolar coupling (a through-space interaction) or the J coupling (mediated by chemical bonds), providing information on spatial proximity, connectivity, and topology.1,2 Although the J coupling is typically much smaller than the dipolar coupling, it does have an isotropic component that survives MAS. Furthermore, although observed linewidths in solids can be large, the refocused line width is generally smaller than the J coupling (see above), enabling a range of standard solution-state techniques to be employed for solid samples.1 In general, sequences are usually synchronized with the sample rotation, and heteronuclear decoupling is employed if strong dipolar interactions are present to ensure the selectivity of the transfer mechanism. In the solution state, spectral resonances are routinely assigned using two-dimensional correlation spectroscopy (COSY),194 with “cross peaks” indicating a scalar coupling between two homonuclear spins. A number of variants of this experiment have been developed for specific application to solids, such as UC2QF-COSY195 and, more recently, SARCOSY.196 As an alternative, homonuclear correlation spectra can also be obtained using experiments based on the solution-state INADEQUATE method,1,197 which uses doublequantum coherence to correlate pairs of spins that share a J 5746


Chemical Reviews

Review

14,149 correlating an isotropic 27Al spectrum with a 31P spectrum is shown in Figure 7b and enables the connectivity of the microporous framework to be deduced.

coupling. For solids, this experiment has been modified to include a refocusing step, ensuring undistorted lineshapes and improved sensitivity.197 The application of heteronuclear correlation experiments in the solid state has also been demonstrated, with a range of techniques adapted from solution-state HSQC and HMQC experiments.1,198,199 More recently, the INEPT experiment has also been applied to solids, again with the use of a refocusing step to ensure the acquisition of undistorted lineshapes.200 These experiments are usually denoted by the prefix MAS-J- to indicate application to solidstate samples. Many of the homonuclear and heteronuclear experiments described above can be adapted to allow magnetization to be transferred through the dipolar coupling, a much larger interaction but one with no isotropic component.1 This can be achieved using one of the many different types of recoupling approaches discussed above. These experiments are usually denoted using the prefix D- or MAS-D- to distinguish them from their J-based counterparts. Transfer via the through-space dipolar coupling is by definition less selective than that through covalent bond, but the additional correlation peaks observed aid not only in spectral assignment but also in determining spatial topology. Dipolar-recoupled versions of both homonuclear and heteronuclear experiments have been implemented. In the case of the former, these include both Double Quantum/ Single Quantum (DQ/SQ) correlations (most notably exploiting BABA171 and POST-C76 recoupling), and also SQ/SQ experiments, usually using some form of spin diffusion.1,201 An example of a 19F through-space DQ MAS NMR spectrum of a fluorinated magnesium silicate mineral is shown in Figure 7a, showing the spatial proximities of the various 19F species.124 For heteronuclear correlations, variants of the HSQC, HMQC, and INEPT have also been implemented (termed D-INEPT, D-HSQC, etc.), with the choice of recoupling technique tuned to the nuclei and system of interest.1 Not all heteronuclear experiments that utilize the dipolar coupling are simple adaptations of well-known solutionstate techniques. One of the simplest experiments, often termed a HETCOR experiment, can be viewed as an adaptation of the basic cross-polarization experiment described above, with the introduction of a t1 evolution period to encode chemical shift information for the I spin. A second example is the TEDOR (transferred echo double-resonance) experiment,202 which is a two-dimensional version of the REDOR experiment described above. Although the discussion above has focused primarily upon techniques introduced for spin I = 1/2 nuclei, many of the experiments can be easily adapted for application to quadrupolar nuclei.2,5,203 Much of the early work in this area focused on the implementation of heteronuclear correlation experiments, in particular those between a quadrupolar spin and a spin-1/2 nucleus, many exploiting the dipolar interaction.5,203 While the use of heteronuclear decoupling is often less of a necessity in many systems where quadrupolar nuclei are involved, the spin dynamics of recoupling and magnetization transfer are often more complex.204 However, in many cases, it is the lack of resolution that limits the information available from heteronuclear correlation spectra involving quadrupolar nuclei, unless the shift dispersion is relatively large. This problem can be overcome by combining the experiments described above with high-resolution approaches such as MQMAS or STMAS.203 An example of this approach, a MQMAS-J-INEPT spectrum of calcined AlPO-

3.6. Combination of GIPAW Calculations and Magnetic Resonance

Although NMR spectroscopy is sensitive to the local structural environment, the interpretation of the high-resolution spectra produced using the varied techniques described above can still pose a considerable challenge. This is particularly true for inorganic solids where many different nuclei with a wide variety of coordination environments are utilized, and there is relatively little information available in the literature. Many of the methods described provide information on NMR parameters, but it is necessary to utilize an additional tool, such as the GIPAW calculations discussed here, to help understand the implications for solid-state structure. In the following section, we will demonstrate how GIPAW calculations can be used to help understand the spectra acquired using the experimental methods discussed in this section to provide new and important information on structure, disorder, and dynamics in a range of solids.

4. APPLICATIONS OF FIRST-PRINCIPLES CALCULATIONS IN CHEMISTRY This section will review all applications of the GIPAW methodology to date (i.e., for the 2001−2011 time period). The year 2001 corresponds to the publication of the pioneering work by Pickard and Mauri92 establishing the concepts of GIPAW applied to NMR (see section 2). Since 2001, GIPAW has been applied to investigate a large variety of chemical topics. Here, we focus intentionally on the point of view of the chemist, including the advantages and drawbacks of the method. In particular, an explicit classification of molecules, materials, and chemical topics is proposed, which should be of broad interest for the majority of chemists. Section 4.1 shows the current impact of GIPAW in chemistry and materials science by reviewing the NMR active nuclei studied by this first-principles method. Section 4.2 is related to organic structures, and important classes of derivatives are discussed individually. The important roles of temperature and dynamics in GIPAW are explicitly reviewed. Section 4.3 is an exhaustive overview of the applications of GIPAW to inorganic architectures. Local dynamics and disorder are also considered. Sections 4.4 and 4.5 are devoted to the description of complex materials, such as polymers, hybrid materials, biomaterials, and nanomaterials. Such materials are usually characterized by interfaces and surface effects. These last points are reviewed in section 4.6. Section 4.7 is a fundamental one as it describes the impact of GIPAW in the development of the emerging field of NMR crystallography. In particular, it is shown how the combination of NMR experiments, first-principles calculations, and molecular dynamics can lead to the determination of crystal structures. Finally, section 4.8 describes other GIPAW applications related to important but less developed NMR topics and to other spectroscopies such as NQR, ESR, and Mössbauer. 4.1. GIPAW: A Tool for Chemistry and Materials Science

After the publication of the GIPAW methodology, the first studies demonstrated its applications to 1H, 11B, 13C, 29Si, and 17 O NMR.92,106,205−207 The main spectroscopic targets here were the isotropic chemical shifts, the CSA parameters, as well as the quadrupolar parameters (for 17O, spin I = 5/2). During 5747


Chemical Reviews

Review

Table 1 calculations nucleus 14

25

N

Mg

35,37

Cl

experimental

calculations characteristics

ref

experimental accuracy

ref

model organic and inorganic systems (silicon nitride, dicyanobenzene, hydrazine dihydrochloride, urea, imidazole, melamine) δiso(14N)exp = 1.05 δiso calc − 146.7 (ppm) (R = −0.997, chemical shift range: ∼350 ppm) CQ(14N): CQ exp = |0.978CQ calc| + 0.11 (ranging from 0 to ∼4 MHz, R = 0.997) ηQ(14N): ηQ exp = 0.978ηQ calc (R = 0.990)

211

CQ ± 0.05 MHz

211

organic and inorganic magnesium oxyanion compounds CQ calc = 1.225CQ exp − 0.13 (MHz) (R2 = 0.987, standard deviation: 0.784 MHz) σcalc = −1.049δexp. + 565.23 (ppm) (R2 = 0.927, standard deviation: 11.9 ppm, σref = 565 ppm with a standard deviation of 3 ppm)

212

Mg-containing compounds (including molybdates, tungstates, and vanadates) σcalc = −1.072(±0.036)δexp + 566.06 (ppm) (global chemical shift range of ∼200 ppm)

213

±0.05 < δiso < ±1.5 ppm ±0.01 < CQ < ±0.1 MHz ±0.02 < ηQ < ±0.1

213

alkaline earth chloride hydrates δiso calc = 1.4016δiso exp + 3.4923 (ppm) (with R2 = 0.921) (systematic overestimation of the calculated values) CQ calc = 1.3941CQ exp + 0.1562 (MHz) (with R2 = 0.995)

214

±1 < δiso < ±5 ppm ±0.05 < CQ(35Cl) < ±0.1 MHz ±0.02 < CQ(37Cl) < ±0.08 MHz ±0.03 < ηQ < ±0.1

214

recent additional observation of CQ calc overestimation in 1-butyl-3-methylimidazolium chloride complex of meso-octamethylcalix[4]pyrrole 43

45

Ca

Sc

69,71

Ga

ηQ ± 0.05

215

inorganic compounds exhibiting Ca−O bonds (including phosphates, silicates, borates, and carbonates) δcalc = 0.956δex. (ppm) (R2 = 0.981, standard deviation: 5.8 ppm, chemical shift range: ∼200 ppm) taking into account the error in the PBE DFT position of the Ca 3d levels with respect to the O 2p states125

216

δiso < ±5 ppm PQ < ±0.1 MHz

216

Ca-containing structures (including CaB6, CaBr2, CaCl2...) σcalc = −1.245δexp + 1122.5 (with R2 = 0.983)

217

δiso ± 2 ppm ±0.1 < CQ < ±0.5 MHz ±0.05 < ηQ < ±0.1

217

scandium sulfate pentahydrate neither CQ nor ηQ calculated data lead to the complete assignment of the 45Sc experimental spectrum (fair agreement) calculated δiso(45Sc) in agreement with experimental chemical shifts

218

δiso ± 0.5 ppm CQ ± 0.05 MHz ηQ ± 0.05

218

metal oxides and gallates different linear fits between calculated and experimental data lead to various rms: in the case of free slope fitting, rms reduced by a factor of 2

219

±1 < δiso < ±7 ppm ±0.01 < CQ < ±0.1 MHz ±0.01 < ηQ < ±0.1

220

73

Ge

germanium halides fair agreement between computed and experimental quadrupolar parameters reproduction of the chemical shifts not as precise

221

±1 < δiso < ±3 ppm ±0.01 < CQ < ±0.09 MHz

221

77

Se

Se-containing compounds (both organic and inorganic) δiso exp = 1.02δiso cal. + 1626 (R2 = 0.959, rms deviation: 92.4 ppm, chemical shift range ∼1600 ppm, σref = 1608 ppm) nonrelativistic approach sufficient in the systems studied

222

δiso ± 0.5 ppm ±10 < δii < ±40 ppm

222

5748


Chemical Reviews

Review

Table 1. continued calculations nucleus

experimental


ref


ref

bromine salts (CaBr2, SrBr2, BaBr2, and corresponding hydrates) CQ calc = 1.1235CQ exp (R2 = 0.982, rms deviation = 5.55 MHz) ηQ calc = 1.0737ηQ exp (R2 = 0.946, rms deviation: 0.0649) δiso calc = 1.1395δiso exp (R2 = 0.946, rms deviation: 57.7 ppm) excellent agreement between computed and experimental Euler angles between the CSA and EFG tensors

223

±5 < δiso < ±30 ppm ±0.02 < CQ < ±0.06 MHz ±0.01 < ηQ < ±0.04

223

CSA calculations: no systematic under- or overestimation of principal components δii 79,81

Br

89

Y

Y2O3, Y2Sn2O7, Y2Ti2O7, Y2O2S, YF3, YAlO3, α-Y2Si2O7, β-Y2Si2O7 good correlation between experimental and computed data (errors of only 15−40 ppm or 0.3−1% of the ∼4000 ppm chemical shift range)

224

δiso ± 0.5 ppm ΔCSA ± 3 ppm

224

91

Zr

ZrO2 polymorphs use of ZORA approach126 fair agreement for δiso (CSA contribution being usually negligible) and CQ

225

δiso ± 20 ppm CQ ± 0.2 MHz ηQ ± 0.1

225

93

Nb

niobate materials core electrons are treated in a scalar relativistic manner CQ: good linear correlation between experiment and calculation for values ranging from ∼20 to >80 MHz (R2 = 0.977) increased scatter for ηQ (R2 = 0.840) CSA: very good agreement (R2 = 0.986), δiso: modest correlation (R2 = 0.784)

105

δiso ± 5 ppm δii ± 5 ppm CQ ± 0.03 MHz ηQ ± 0.003

105

95

Mo

ternary molybdates (AMoO4, A′2MoO4) implementation of ZORA126 relativistic effects during pseudopotential generation mean difference between computed and experimental CQ values ∼0.51 MHz (for CQ values ranging from 1.1 to 3.1 MHz) σiso calc = −0.95δiso exp − 844 (ppm) (R2 = 0.960, standard deviation: 15.5 ppm, chemical shift range; ∼350 ppm)

226

±0.2 < δiso < ±0.6 ppm ±0.02 < CQ < ±0.05 MHz ±0.05 < ηQ < ±0.08

227

inorganic Sn-containing structures (Sn2+ and Sn4+) no additional relativistic effect (apart from the scalar relativity effect included in the pseudopotential approximation) δiso exp = 0.88δiso calc − 33 (ppm) (with R2 = 0.95, chemical shift range: ∼6500 ppm)

228

δiso ± 1 ppm

229

230

±10 < δiso < ±100 ppm ±0.1 < CQ < ±1 MHz ±0.002 < ηQ < ±0.02

230

119

Sn

127

I CQ calc

metal iodides (MgI2, CaI2, SrI2, BaI2), hydrates (BaI2·2H2O and SrI2·6H2O), CdI2 = 1.135CQ exp − 3.127 (R2 = 0.9755, rms: 13.8 MHz) with CQ values up to 2000 MHz ηQ calc = 0.758ηQ exp + 0.016 (R2 = 0.9457) δiso calc = 1.537δiso exp − 227.48 (R2 = 0.9581)

137

Ba

inorganic materials (including nitrates, carbonates, etc.) use of ZORA approach126 CQ: excellent linear correlation between experimental and calculated CQ (R2 = 0.997, overall range: from 5 to ∼30 MHz CSA parameters difficult to compare (minimal influence on the 137Ba spectra, and difficult to determine accurately experimentally) representation of some tensor orientations: see Figure 8

163

±0.1 < CQ < ±1 MHz

163

195

Pt

cisplatin and related square-planar compounds dependence of the NMR parameters on the local

231

±15 < δiso < ±50 ppm

231

5749


Chemical Reviews

Review

Table 1. continued calculations nucleus

experimental


ref


ref

232

±30 < δiso < ±200 ppm ±1 < CQ < ±10 MHz

232

hydrogen-bonding environments σiso calc = −2.41δiso exp − 1950, R2 = 0.7301, chemical shift range; ∼600 ppm) 209

Bi

BiOX (X = Cl, Br, I), nitrate, triflate, and acetate strong discrepancies for CQ due to: (i) nature of the used PP; (ii) uncertainty of 209Bi nuclear quadrupole moment (over a range from 78 to 256 MHz) experimental trend of increasing CQ in the series X = I, Br, Cl well described by calculations second-order perturbation theory not entirely valid in that particular case to fit the experimental spectra

absolute deviation with respect to experiment is estimated to be ∼15 ppm. It was noted that discrepancies between experimental and calculated values can be partially masked in the average isotropic chemical shifts. The computed Euler angles between CSA and EFG tensors were also extracted showing major discrepancies with experimental data: one cause could be related to the overall small CSA for the three sites in AlVO4. Here, it has to be emphasized that from the experimental point of view, Euler angles can be rather difficult to extract with high accuracy. Subsequent work by Truflandier et al.86 addressed the accuracy of the GIPAW approach for 3d transition metals (such as 49Ti and 51V) including: (i) the accuracy of the pseudopotential, (ii) the transferability of a given pseudopotential to a large series of Ti and V complexes, and (iii) the accuracy of the exchange-correlation functional (see section 2.2.2). Ultrasoft pseudopotential GIPAW results are compared to those obtained by the all-electron approach IGAIM (Individual Gauges for Atoms In Molecules)209 using geometry-optimized structures. The core−valence partition necessary for pseudopotential construction was carefully checked for the 3d elements, and as was observed for AlVO4, the nonrigid core state contribution of the 3p level was demonstrated (51V). It was also shown that the use of ultrasoft pseudopotential is necessary for convergence at a reasonable planewave cutoff energy. The convergence of the results was also checked with respect to the number of projectors used for each valence state (especially for the 3p and 3d states). For this purpose, VOCl3 (for 51V) and TiCl3CH3 (for 49Ti) were chosen as standard molecules. The application to periodic systems was also performed through the detailed study of 13 inorganic structures containing vanadium atoms (for a total of 18 V sites), including ortho- and pyro-vanadates (4-fold coordination), as well as 6-fold coordinated vanadium atoms. The correlation between σiso calc (ppm) and δiso exp (ppm) was clearly linear with a slope −1.047, correlation coefficient −0.988, and rms deviation 28 ppm. Moreover, σref = −1939 (59) ppm is in good agreement with the value obtained for VOCl3 from all-electron calculations. For assignment purposes, σref can be adjusted for a given small range of chemical shifts. Whereas δiso(51V) appears accurately reproduced by GIPAW, some discrepancies remain for the magnitude and asymmetry of the anisotropic shielding. Finally, using quantum chemical methods, the potential of hybrid functionals (including Hartree−Fock exchange) was examined using VOCl3. It was found to give a partial improvement in the calculation of the anisotropy. An investigation into the origin and nature of any

the past decade, numerous authors have applied GIPAW calculations to study NMR active nuclei across the periodic table (∼20 different nuclei to date). For a given isotope, a series of simple structures (involving one or more inequivalent crystallographic sites) is typically used as a benchmark for comparison to experimental data. As outlined in section 3, MAS is the method of choice for obtaining information on δiso. for spin I = 1/2 nuclei (such as 13C, 19F, 29Si, 31P, etc.). CSA parameters can be extracted from static, slow MAS, or CSArecoupled 2D experiments. The case of 1H (and to some extent 19 F) remains special as the determination of δiso usually requires much more involved multiple-pulses techniques. For quadrupolar nuclei (such as 17O, 35Cl, etc.), more complex experimental schemes are necessary for the precise determination of δiso, CQ, and ηQ (e.g., MQMAS, STMAS). The subsequent calculation of NMR parameters by GIPAW requires structural models (generally from XRD, neutron diffraction data, or in silico computation). Structures can be used as obtained or can be geometry-optimized depending upon the accuracy of the initial data. For a detailed discussion on this topic, see, for example, ref 56. It is generally accepted that for H-containing structures, prior optimization of the 1H positions is mandatory for accurate results. In most cases, experimental and calculated data are compared, and linear regression is used to obtain a value for σref (in ppm), enabling the calculations to be referenced correctly (more details on choice of reference shielding can be found in section 2.1.1). In terms of computation, a crucial step lies in the choice of the adequate pseudopotential. Early studies used norm-conserving pseudopotentials in the Troullier−Martins formulation, while more recent works typically employ the more efficient ultrasoft pseudopotentials (see section 2). It should be noted that some authors have performed detailed studies related to: (i) the improvement of the generated pseudopotential, (ii) the transferability of a given pseudopotential, (iii) the improvement of the correlation functional, (iv) the application of the pseudopotential approach to 3d and 4d elements, and (v) the relativistic effects for heavy nuclei. For illustrative purposes, we highlight here the case of 49Ti and 51V. The application of GIPAW to 51V has been carefully and systematically studied in its role as prototype for transition metal nuclei in diamagnetic solid-state systems.208 Using the PW91 functional, ultrasoft pseudopotentials, with the semicore 3s and 3p states treated as valence, GIPAW results were obtained for the three crystallographically inequivalent V sites in AlVO4. For the nine corresponding principal components of the CSA, the mean 5750


Chemical Reviews

Review

Figure 8. Calculated 137Ba quadrupolar tensor orientations: (a) barium nitrate, (b) barium carbonate, (c) barium chlorate monohydrate, (d) barium chloride dihydrate, and (e) anhydrous barium chloride. Adapted with permission from ref 163. Copyright 2010 American Chemical Society.

shifts. It is well established that the 1H isotropic chemical shifts are sensitive to long-range electrostatic effects. Several studies233−236 have quantified the effect on 1H shifts due to hydrogen bonding (both conventional and “weak”) as well as ring currents. Using a set of 15 organic compounds, Salager et al.237 estimated the rms error on 1H chemical shifts to be 0.33 ppm. For the dipeptide β-AspAla, 1H CSA calculated data have been used to simulate the build-up DQ coherences in a 1H 2D correlation experiment (see section 3.5).238 Deuterium NMR parameters (such as the quadrupolar constants, CQ) have also been computed by the GIPAW approach.131 13 C chemical shifts have been extensively studied with the GIPAW approach. In 2009, a systematic study of 13C chemical shifts was performed.239 13C chemical shift tensors of 14 molecules (aromatics and carbohydrates) whose crystal structures have been precisely determined by neutron diffraction were calculated using both a cluster approach and GIPAW. This study showed that the GIPAW 13C predictions are the most accurate (i.e., in best agreement with the experimental data) with an error on the isotropic shift of ∼1 ppm. As in many contributions, it was emphasized that the intermolecular shielding arising from the lattice (i.e., crystal packing) is significant (especially in the case of hydrogenbonded systems). 4.2.1. Amino Acids, Peptides, and Nucleic Acids. Amino acids were the first bio-organic molecules investigated by the GIPAW approach. Glutamic acid polymorphs240 as well as glycine, valine, tyrosine, and alanine and their hydrochlorides241 were studied with a particular emphasis on the calculation of the 17O NMR parameters. Indeed, 17O is a particularly powerful structural probe,242 as its chemical shift range covers almost 1000 ppm in organic molecules, as compared to ∼250 ppm for 13C. Moreover, its EFG is very sensitive to the local geometry. Generally speaking, the studies demonstrated the existence of correlations between the calculated 17O NMR parameters and the strength of the CO····H hydrogen bonding. Moreover, for the case of

systematic errors within the theoretical framework was also undertaken,210 with particular emphasis on (i) the quality of the structural models, (ii) the functional used, and (iii) the reliability and transferrability of the pseudopotentials. For the eigenvalues of the quadrupolar tensor, the following systematic errors were extracted: ±400 kHz for CQ; ±0.2 for ηQ; ±40 and ±17 ppm for δiso.(51V); and ±48 and ±58 ppm for CSA, obtained from XRD and DFT-optimized structures, respectively; the chemical shift anisotropies are systematically underestimated with a shift of 85 ppm. Finally, the quadrupolar moment Q of 51V was estimated to 53 ± 1 mbarn in very good agreement with the observed experimental value. In Table 1, the results related to the GIPAW approach are described in detail for different nuclei, summarizing calculated correlations versus experimental data as well as experimental typical accuracies. 4.2. Organic Structures and Assemblies

One of the first uses of the GIPAW/NMR combination for organic crystalline structures was dedicated to phenylphosphonic acid,233 a system with chelation properties useful for the synthesis of hybrid materials by surface modification. The accuracy of the calculated 1H, 13C, 17O, and 31P NMR parameters, including chemical shift tensors and quadrupolar parameters for 17O, enabled the unambiguous assignment of all of the sites present in the structure. Moreover, the effects of intermolecular interactions on the NMR parameters (in particular for 1H and 17O) were investigated by comparing the results of the calculations in the crystal and in an isolated molecule of phenylphosphonic acid. In the following sections, a detailed and exhaustive discussion of the application of GIPAW calculations to organic systems (including amino acids and peptides, carbohydrates and oligosaccharides, drugs and polymorphs, supramolecular assemblies, graphene, and nanotubes) is presented. Here, the main spectroscopic targets are typically 1H, 13C, 14,15N, and 17O. As might be expected for 1H, there have been numerous applications of the GIPAW method to determine chemical 5751


Chemical Reviews

Review

hydrochloride amino acids,241 the influence of the Cl····H hydrogen bonding on the 35Cl quadrupolar parameters was also investigated and showed a consistent increase in the chlorine quadrupolar coupling constants with the number of close Cl····H contacts. Monosodium L-glutamate monohydrate was also investigated more recently by 17O NMR243 combining DOR and MQMAS experiments (see section 3.2) with GIPAW calculations. This combination of techniques was able to resolve the eight distinct 17O sites despite them lying within a narrow chemical shift range of

First-Principles Calculation of NMR Parameters Using the Gauge

Recommend Documents