Large-Scale Computation of Nuclear Magnetic Resonance Shifts for

Article Cite This: J. Chem. Theory Comput. 2018, 14, 377−394

pubs.acs.org/JCTC

Large-Scale Computation of Nuclear Magnetic Resonance Shifts for Paramagnetic Solids Using CP2K Arobendo Mondal,† Michael W. Gaultois,‡ Andrew J. Pell,§ Marcella Iannuzzi,∥ Clare P. Grey,‡ Jürg Hutter,∥ and Martin Kaupp*,† †

Institut für Chemie, Theoretische Chemie/Quantenchemie, Technische Universität Berlin, Sekretariat C7, Strasse des 17 Juni 135, D-10623 Berlin, Germany ‡ Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom § Department of Materials and Environmental Chemistry, Arrhenius Laboratory, Stockholm University, SE-10691 Stockholm, Sweden ∥ Institut für Chemie, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland S Supporting Information *

ABSTRACT: Large-scale computations of nuclear magnetic resonance (NMR) shifts for extended paramagnetic solids (pNMR) are reported using the highly efficient Gaussian-augmented plane-wave implementation of the CP2K code. Combining hyperfine couplings obtained with hybrid functionals with g-tensors and orbital shieldings computed using gradient-corrected functionals, contact, pseudocontact, and orbital-shift contributions to pNMR shifts are accessible. Due to the efficient and highly parallel performance of CP2K, a wide variety of materials with large unit cells can be studied with extended Gaussian basis sets. Validation of various approaches for the different contributions to pNMR shifts is done first for molecules in a large supercell in comparison with typical quantum-chemical codes. This is then extended to a detailed study of g-tensors for extended solid transition-metal fluorides and for a series of complex lithium vanadium phosphates. Finally, lithium pNMR shifts are computed for Li3V2(PO4)3, for which detailed experimental data are available. This has allowed an in-depth study of different approaches (e.g., full periodic versus incremental cluster computations of g-tensors and different functionals and basis sets for hyperfine computations) as well as a thorough analysis of the different contributions to the pNMR shifts. This study paves the way for a more-widespread computational treatment of NMR shifts for paramagnetic materials.

1. INTRODUCTION

greater complexity of extended paramagnetic solids, progress in that field has been slower despite the importance of such materials and the complexity of their NMR spectra. However, substantial steps have been made here recently as well.28−36 The earliest solid-state DFT calculations that attempted to model the Fermi-contact (FC) contribution to the hyperfine shift used plane-wave pseudopotential calculations.28 The needed spin density at the NMR nuclei was approximated by an empirical correlation with the valence-shell spin-density distribution around the atom of interest. It has been shown that all-electron calculations are preferable for a direct computation

The underlying physics of NMR shifts in paramagnetic systems is invariably more complex than in diamagnetic ones due to interactions between electron and nuclear spin moments, which often dominate the observed quantities.1−9 These “hyperfine shifts”, however, encode important information about the molecular and electronic structure of a system.2,5,10 It is thus important to be able to compute such NMR shifts in paramagnetic systems (”pNMR shifts”) from first principles to complement experimental studies in interpretative, affirmative, or even in predictive mode. Consequently, the past 15 years have seen substantial efforts to improve theoretical formalisms and implementations in computer programs for pNMR shifts in molecular systems.2,11−27 Due to the much© 2017 American Chemical Society

Received: September 23, 2017 Published: November 28, 2017 377

DOI: 10.1021/acs.jctc.7b00991 J. Chem. Theory Comput. 2018, 14, 377−394

Article

Journal of Chemical Theory and Computation of the FC shifts.30 Consequently, the Gray and Carlier groups subsequently turned to all-electron calculations, either within the projector-augmented wave (PAW)37 reconstruction of the core density,38−40 the full-potential PAW implementation in WIEN2k,31,41 or the linear-combination of atomic orbital (LCAO) implementation in the CRYSTAL code42,43 (using Gaussian-type orbital basis sets, GTOs). Focusing on lithium transition-metal oxides, phosphates, and silicates, with emphasis on systems of relevance for battery materials, these computations allowed an analysis of the FC contributions in terms of spin-density transfer pathways from the transitionmetal centers to the NMR nucleus of interest (e.g., Li, P, and O).29−33 Notably, in contrast to isolated molecules, such calculations on extended paramagnetic solids need to be adjusted for the ferro- or antiferromagnetic couplings between the spin centers in the Curie−Weiss paramagnetic temperature region that frequently corresponds to the situation of the NMR measurement.29−31 For molecular systems, several studies 14−17,22−26,44−51 showed contributions beyond the FC shift to become important, particularly when the former is, for various reasons, small. These additional contributions consist of a) orbital shifts (analogous to the Ramsey-type expression52,53 of shifts for diamagnetic systems) and what has been termed “pseudocontact” or “dipolar” shifts.54,55 The latter arise from interactions that couple anisotropic hyperfine interactions at the nucleus of interest with the magnetic anisotropy at the paramagnetic metal center (encodable by the electronic gtensor and, for S > 1/2, by the zero-field splitting, ZFS). The first attempts to include pseudocontact (PC) shifts for extended solids using g-tensors obtained in DFT calculations have been published recently. Mali et al.34 included g-tensors obtained from GIPAW56,57 calculations in Quantum Espresso58 to compute PC shifts in a study of Li2FeSiO4 polymorphs.34 Pigliapochi et al.59 combined GIPAW g-tensors from Quantum Espresso within the generalized gradient approximation (GGA, here the PBE functional) with hyperfine couplings obtained at hybrid-DFT levels with CRYSTAL (GTO LCAO) for a number of olivine-type LiMPO4 materials. Here, we report our first systematic studies using a different type of implementation of pNMR shifts using the Gaussianaugmented-plane-wave (GAPW)60,61 approach available within the CP2K code.62 GAPW combines the advantages of GTOtype LCAO-basis sets for the orbitals with those of an augmented plane-wave basis for the density.63,64 A large advantage of CP2K is that it has fewer limitations regarding the permitted GTO basis sets compared to CRYSTAL, and it permits efficient computations with hybrid functionals and large basis sets for the hyperfine couplings. In fact, the extremely efficient parallelization of CP2K makes the code particularly suitable for highly parallel large-scale computations on supercomputers and allows us to address large unit cells.65 gTensor calculations with CP2K using GGA functionals are available and have been reported before.66,67 In addition, here we report inclusion of orbital shifts for extended solids at the same level. Notably, the present report illustrates the bridges between molecular and periodic pNMR shift computations. The setup of this paper is as follows: after describing the theoretical background of the computations and providing computational details, discussion of the results starts from validation of CP2K computations for various separate contributions and of full pNMR shifts for isolated molecules in a large supercell against typical computations using

molecular quantum-chemical codes. This is then extended to g-tensor computations for extended transition-metal fluorides and lithium vanadium phosphate materials of general composition LixV2(PO4)3 (x = 3, 2.75, 2.5, 2.25, 2, and 1). Finally, pNMR shifts for different Li sites are validated and analyzed at various computational levels for Li3V2(PO4)3. The choice of this material was partly motivated by the expectation that it exhibits small ZFS (the latter is not included in the present computations, calculations including ZFS will be reported elsewhere) while being closely related to the abovementioned battery materials. For g-tensor computations, we also have compared full computations using periodic boundary conditions against an incremental approach based on cluster models.

2. THEORY pNMR Shift Theory. Our pNMR shielding calculations follow the doublet-state framework proposed by Moon and Patchkovskii,11 found to be extendable beyond doublet states when the ZFS is negligible.12 In this formalism (note ref 22 for alternatives), the hyperfine shielding is expressed via a matrix product of the g- and A-tensors. Eq 1 formulates this for exchange-coupled solids in the Curie−Weiss temperature regime: I σ I = σorb −

μB S(S + 1) I 3kBg N μ N (T − Θ)

gAI (1)

where S is the total electronic spin quantum number, μB the Bohr magneton, kB the Boltzmann constant, gIN the g-factor of nucleus I, μN the nuclear magneton, T temperature, Θ the Weiss constant of the system, g the electronic g-tensor, and AI the hyperfine coupling tensor of nucleus I. In contrast to ref 11, which deals with isolated paramagnetic molecules, eq 1 replaces T by T − Θ in the denominator of the hyperfine shift term to account for the residual exchange couplings of the extended solid-state system in question, assuming Curie−Weiss behavior at the temperatures probed experimentally.68,69 Such a temperature dependence is also observed for the paramagnetic shifts in exchange-coupled systems and thus introduced into eq 1.29,30,33,70 For each material, the Weiss constant obtained experimentally by low-temperature extrapolation of the hightemperature susceptibility curves will be used. For the moment, this adds an empirical correction to the computations. In the long run, it will be desirable to obtain analogous information from a statistical first-principles treatment of the exchange interactions in the solid (see, e.g., refs 36, 71). Note, however, that at the simulation temperature (320 K) used, a change in the Weiss constant of 20 K in one or the other direction affects the hyperfine shifts only by about 5%. Both g and AI are general 3 × 3 matrices composed of a scalar (rank zero), a true anisotropic symmetrical second-rank tensor (g- and A-tensor anisotropies), and an asymmetric (rank 1) part (the latter is zero for the nonrelativistically computed AI used here). We may separate the g- and A-matrices into their individual contributions: g may be written as (ge + Δgiso)1 + Δg,̃ where ge is the isotropic free-electron g-value, Δgiso the isotropic deviation, and Δg̃ the traceless anisotropic part. Similarly, the nonrelativistic AI is separated into the isotropic Fermi-contact part AIFC1 and the dipolar tensor AIdip. Because AIdip and Δg̃ are traceless, eq 2 gives only four nonzero contributions to the isotropic shielding: 378


Article

Journal of Chemical Theory and Computation I I = σorb,iso − σiso

electrons locally, in the way done in eqs 5 and 6.2,43,74,77 We thus have to divide the spin density provided by CP2K by a factor of 2S, where S is the spin of the paramagnetic center.30 Aspects of g-tensor calculations. Because the formalism we use for pNMR shielding (cf. eq 1 above) requires the computation of g- and A-tensors, it is important to examine aspects that may affect their accuracy. The solid-state g-tensor implementation in CP2K62 is based on the second-order perturbation treatment by Pickard and Mauri57 for paramagnetic defects in solids using periodic boundary conditions (PBCs). The different contributing terms are:


I I I (geAFC )iso 1 + Δg isoAFC 1 + ΔgA ̃ dip

(2)

A total of four terms from the matrix product gAI contribute to the shielding anisotropy, giving: I I = σorb,aniso − σaniso


I I I I (ge A dip )aniso + Δg iso A dip + ΔgA + ΔgA ̃ FC ̃ dip

(3)

where the subscript “aniso” indicates the traceless anisotropic part of the sum in parentheses. Traditionally, the experimental isotropic pNMR shift is decomposed into three parts: orbital shift (δorb), Fermi-contact shift (δFC), and pseudocontact shift (δPC): exp exp exp δexp = δorb + δ FC + δ PC

g = ge1 + Δg ZKE + Δg SO + Δg SOO = ge1 + Δg

Deviations of Δg from the free-electron g-factor (first term) are provided by the (isotropic) Zeeman kinetic energy (ZKE) contribution (second term), the “standard” spin−orbit (SO) contributions (third term, contains diamagnetic “gauge” and paramagnetic parts), and the spin-other-orbit (SOO) contribution (fourth term). Pickard and Mauri define these terms as:57,66 g ΔgijZKE = − 2e (T α − T β)δij (8) c

(4)

Conversion from nuclear shieldings to relative chemical shifts is done in the usual way by subtracting the shielding from that of a suitable reference compound (see the Computational Details I section). Thus, comparing eq 4 to eq 2, obviously σorb,iso exp corresponds to δorb. The sum of the two terms depending on geAIFC and ΔgisoAIFC corresponds to δexp FC , whereas in the doublet formalism (i.e., without ZFS), the term depending on the isotropic part of ΔgA ̃ Idip gives δexp PC . The shift tensor can be characterized by the shift anisotropy, Δ = δzz − (1/2)(δxx + δyy), and by the asymmetry parameter (deviation from axial symmetry), η = (δyy − δxx)/(δxx − δiso), where the components are arranged according to the Haeberlen convention to give |δzz − δiso| ≥ |δxx − δiso| ≥ |δyy − δiso|, where the orientations of the δii components correspond to the principal axis frame of the symmetric part of the shift tensor.72,73 Hyperfine Coupling Tensors. Here, we so far are restricted to the nonrelativistic formalism, where A iso corresponds to the Fermi-contact term and is thus proportional to the electronic spin density ρIs(0) at the position of nucleus I:74−77 I A iso

I 4π geμB gNμ N I I = ρs (0) = AFC 3 ⟨Sz⟩

ΔgijSO =

I 1 geμB gNμ N 2 ⟨Sz⟩

∫ d rρsI (r)

r5

c

ΔgijSOO = 2

∫Ω [jiα (r) × ∇Veffα (r) − jiβ (r) × ∇Veffβ (r)]j d3r c

∫Ω Bijcorr (r)ρs (r)d3r c

(10)

Here, ρs = ρα − ρβ is the spin density, α and β denote the spin channels σ, and Tσ, jσi , and ρσ are the unperturbed kinetic energy, induced current density, and electron density of the spin-σ channel, respectively. The integrals span the volume of the simulation cell, Ωc. Vσeff is the effective Kohn−Sham potential for the spin-σ electrons. Bcorr represents the magnetic ij field that originates from the corresponding total (selfinteraction-corrected) induced current density. The spin−orbit contributions dominate the g-tensor. Upon reduction of the full Dirac−Coulomb−Breit Hamiltonian for many-electron systems to the Breit−Pauli Hamiltonian, one obtains (a) the one electron spin−orbit nucleus term, (b) the two-electron spin−orbit term resulting from Coulomb and exchange screening of the nuclear charge, and (c) the spinother-orbit term arising from the Breit interaction (relativistic corrections to electron repulsion). Disregarding here the fact that relativistic exchange-correlation functionals would be required for a complete DFT computation of SO effects and assuming that we may use the Kohn−Sham determinant to compute SO matrix elements, we nevertheless end up with a large number of two-electron SO integrals. Several approximations have been suggested to simplify matters. Complete replacement of the two-body interactions by the effective oneelectron Kohn−Sham potential (Veff) is the simplest and most often used. However, this neglects contributions arising from antisymmetric (exchange-like) terms, including the SOO term.78,79 Mean-field approaches including the exchange and SOO terms include the spin−orbit mean-field (SOMF)78 method, additionally invoking a one-center approximation gives the so-called atomic mean-field integrals (AMFI) method.80 CP2K has the Pickard−Mauri57,66 Veff implementation augmented by an approximation to the SOO term. Veff is obtained by integrating products of the induced current

(5)

3rri j − δijr 2

(ge − 1)

(9)

where ⟨Sz⟩ is the expectation value of the z component of the total electronic spin. The components AIaniso,ij of the anisotropic part of the hyperfine matrix result from dipole−dipole interactions,74−77 We note in passing that in molecular computations of pNMR shifts, spin−orbit, and scalar relativistic corrections to the HFCs have sometimes been included.12,14,15,17,18,24,25 They are, however, currently not available in solid-state codes and, moreover, are not expected to be important for the lithium HFCs in the present work. I A aniso, ij =

(7)

I = A dip, ij

(6)

Subscripts i and j refer to Cartesian coordinates x, y, and z. While local spin and, thus, ⟨Sz⟩ is easily defined for molecules, it is less-straightforward for condensed-phase systems. In fact, the implementation of hyperfine tensors in solid-state codes does not normalize the spin density to the number of unpaired 379


Article

Journal of Chemical Theory and Computation densities and the gradient of the effective potential Vτeff (α, β)81 over the simulation cell. Vτeff is defined as: τ Veff (r ) = νext(r ) + νH(r ) + νxcτ (r )

where c is the speed of light in vacuum, i and j are the indices to represent x, y and z directions, r′ is the position of nucleus I, ji is the current density induced by a constant external magnetic field applied along the i axis, and ΩS is the volume of the entire integration domain, including the periodic replicas of the simulation cell. The other tensor components can be obtained by changing the indices accordingly. This corresponds to the open-shell generalization of the previously reported nuclear shielding implementation for diamagnetic systems,66 and we have used IGAIM gauges with CP2K (GIAOs with the Gaussian09 molecular code for comparison).

(11)

where νext is the Coulombic potential from the nuclei, νH the Hartree potential, and ντxc the exchange-correlation potential. For further elaboration and implementation, see refs 66 and 82. In eq 10, for the missing SOO term,57 Bcorr represents the i magnetic field that originates from the corresponding total (self-interaction-corrected) induced current density, ji,66 as experienced by the unpaired electron. That is, the SOO contribution to a g-shift tensor component becomes: ΔgxySOO

1 = S

∫

B1,y x (r )[ρα

− ρβ ]dr

3. COMPUTATIONAL AND EXPERIMENTAL DETAILS General Aspects. Calculations with the CP2K/Quickstep package62,89 used periodic boundary conditions for both molecules (large supercell) and solids. Optimizations of cell parameters used the hybrid Gaussian and plane waves (GPW) formalism together with the pseudopotential approximation, applying the PBE GGA exchange-correlation functional.90 Goedecker−Teter−Hutter (GTH) pseudopotentials91 and double-ζ MOLOPT basis sets (DZVP−MOLOPT−SR− GTH)92 have been used for all elements. For the expansion of the charge density in plane waves, an energy cutoff of 500 Ry has been used, and the convergence criterion over the maximum component of the wave function gradient was set to 1.0 × 10−7. Calculations have been carried out on a 2 × 2 × 2 supercell to minimize artifacts of the periodic boundary conditions. Note that, while the GAPW ansatz in CP2K uses minimal k-point sampling, the code enables the use of extended supercells to compensate for this and to thus improve accuracy. Keeping the cell parameters fixed, the atom positions were further optimized using the all-electron Gaussian-augmentedplane-waves (GAPW)60,61 formulation with PBE and def2TZVPD/def2-TZVP basis sets. 94,95 The wave function convergence criterion was set to 1.0 × 10−6, and the energy cutoff has been kept at 500 Ry. Magnetic-resonance parameter calculations employed the allelectron GAPW implementation of CP2K, which is particularly suitable for the properties needed here. Computation of orbital shieldings followed the existing CP2K implementation for diamagnetic systems66 (related to GIPAW) in its open-shell extension, using IGAIM. g-Tensor calculations used the related CP2K implementation,93 with spin−orbit treatments described in the Theory section, and IGAIM for the gauge treatment. Both g-tensors and orbital shieldings were obtained at the PBE level. Hyperfine coupling (HFC) tensor computations were based on the implementation of ref 77. In addition to PBE, PBE-based global hybrids were used, varying the exactexchange (EXX) admixture from 5% (PBE5) to 40% (PBE40).With minimal k-point sampling, the use of larger supercells is also known to help improve the accuracy for sensitive quantities such as HFCs.93 Basis sets to be used in the different comparisons of nuclear magnetic resonance and electron paramagnetic resonance (NMR and EPR) parameters vary and will be detailed further below. Experimental Determination of the Structure of Lithium Lactate. Solid Li lactate is an excellent standard for Li shift referencing (discussed later), so the crystal structure of Li lactate was determined by X-ray diffraction (XRD) to serve as the input structure for orbital shielding calculations (see below). Lithium L-lactate (Sigma-Aldrich, > 98%) was dissolved in ethanol, and recrystallization by isopropanol vapor transport produced needles and rectangular plates after

(12)

where B1, x is the magnetic field induced by the electronic currents when a magnetic field with unit magnitude is applied in the x direction (referred to as Veff,PM). This may miss the above-mentioned other exchange-like contributions, which is why one may expect a somewhat overestimated overall SO term compared to the more-complete treatments. The vector potential of the external magnetic field is defined up to an arbitrary translation. For computations in a finite basis set, this leads to the well-known gauge dependence. While this is thought to be less pronounced for g-tensors than, for example, for nuclear shieldings, it nevertheless needs to be accounted for. The most refined method for molecular calculations is the use of gauge-including atomic orbitals (GIAO).83 The individual gauges for atoms in molecules (IGAIM)84 or the closely related, simplified continuous set of gauge transformation (CSGT) approach85 have been found to perform almost as well.86 A common gauge origin (e.g., at the center of electronic charge, the center of nuclear charges or the position of a specific atom) is clearly the crudest possibility. The proper choice of gauge origin becomes important when the g-shifts are very small and/or spin orbit contributions from the atoms in different parts of a rather unsymmetrical molecule are of similar magnitude.87,88 CP2K offers extensions of IGAIM and CSGT to periodic boundary conditions. While the CSGT method has higher computational efficiency,66 in this work we have focused on using IGAIM. So far solid-state g-tensor codes offer only local or generalized gradient approximation (GGA) implementations but no hybrid functionals. The existing implementations of g-tensors in solid-state codes have furthermore been geared toward paramagnetic defects with only one paramagnetic center in the simulation cell. When using larger cells with more than one center, we found an incorrect size dependence, where the individual site gtensors simply added up (as also mentioned in ref 59). To calculate the g-tensor for periodic paramagnetic solids with multiple paramagnetic centers n, we thus have to normalize the resulting g-tensor of the simulation cell as: g = ge1 +

1 1 (Δg ZKE + Δg SO + Δg SOO) = ge1 + Δg n n (13)

σIorb,ij

Orbital Shielding Tensor. The components of the orbital shielding tensor of nucleus I for periodic solids, as implemented in CP2K, are defined as: I σorb,ij =

1 c

⎡

∫Ω ⎢⎣ |rr −− rr′|′3 S

⎤ × ji (r )⎥ d3r ⎦j

(14) 380


Article

Journal of Chemical Theory and Computation

Figure 1. Supercell (2 × 2 × 1) of nasicon-like lithium vanadium phosphate Li3V2(PO4)3 with 320 atoms (48 lithium atoms with three distinct sites) used for the HFC calculations.

atoms) (see Figure S2). We note in passing that for the large orthorhombic cells of Li lactate, converged PBC-based computations of the orbital shieldings for large supercells are currently precluded. The shieldings have been computed at PBE/IGLO-II level in a supercell of sizes of 40 Å × 40 Å × 40 Å and 60 Å × 60 Å × 60 Å for the small and the large cluster, respectively. The resulting shieldings are 90.5 and 90.4 ppm for the small and large cluster, respectively (the shieldings of the inner six Li sites were averaged for the large cluster), indicating a converged cluster model. These data should be useful for further analyses of lithium shifts in solids. Comparisons of CP2K Supercell Results for Molecules against Quantum-Chemical Codes. To validate the computation of the different ingredients contributing to pNMR shifts obtained with CP2K against standard quantumchemical codes, orbital shieldings, HFC tensors, and g-tensors for two doublet molecules (cobaltocene and the TEMPO nitroxide radical) have been computed in a large cubic supercell of size 30 Å × 30 Å × 30 Å, using large def2-QZVP basis sets.94 This supercell size ensures negligible interactions with molecules in neighboring cells. These calculations used structures optimized at B3LYP/6-311++G(d,p)100−104 level with the Gaussian09105 code. Comparisons for HFCs and orbital shieldings were done against Gaussian09 results, g-tensors against ORCA106 and MAG107 data, using the same def2-QZVP basis sets. While ORCA provided only a common gauge for the external magnetic field (placed at the center of charge), the MAG calculations allowed additionally the use of gauge-including atomic orbitals (GIAOs).83 This turned out to be of importance for some comparisons (see below). Orbital shielding and gtensor comparisons were done with the BLYP functional,100,108 while HFCs were computed with both BLYP and PBE0109 functionals. Some additional explorations of gauge aspects and

several days. X-ray diffraction (XRD) on a rectangular single crystal with dimensions 0.30 mm × 0.08 mm × 0.02 mm was performed by collecting ω and φ scans with a Bruker D8QUEST using Cu Kα radiation (λ = 1.5406 Å), and structure solution and refinement was performed with use of the SHELX program package.96 Material as purchased was lightly ground and powder XRD was performed to confirm the single crystal was representative of the bulk (PANalytical Empyrean, Cu Kα radiation). The structure was solved using a unit cell in the P21 space group, with dimensions a = 4.5412(4) Å, b = 4.8838(4) Å, c = 10.1982(10) Å, and β = 98.226(7)°; Li+ ions sit in a LiO4 tetrahedral environment (see Figure S1). H atoms were treated by a mixture of independent and constrained refinement, The absolute structure was assigned from known L-lactate material, with (S)-configuration at C2. Crystallographic details are listed in Table S1, atomic coordinates are listed in Table, S2, and anisotropic displacement parameters are listed in Table S3. Our structure determination differs from a previous report by Clough and Poldy, who observed a similar crystal habit of a thin rectangular plate but assigned a P21/c space group in their preliminary analysis.97 Due to symmetry considerations, the P21/c space group is unable to support the periodic packing of a single enantiomorph, whereas the P21 space group determined in the structure solution presented here is a nonenantiogenic Sohncke space group capable of supporting the packing of the L-lactate enantiomorph. Shift Referencing. 7Li shifts of the solid in question are obtained by subtracting its computed absolute shielding from that of lithium lactate (LiC3H5O3), i.e., δI = σIref − σI. The orbital shieldings for solid lithium lactate were computed for two cluster models cut out from the XRD supercell structure (see above), to evaluate convergence with cluster size. The smaller cluster has eight units of LiC3H5O3 (96 atoms). The second, much-larger cluster has 18 units of LiC3H5O3 (216 381


Article


Figure 2. Molecular models of the two distinct vanadium centers of Li3V2(PO4)3 solids, where the terminal oxygens have been saturated with hydrogen atoms. Model A and B have four and three lithium atoms, respectively.

ture. These computations used basis sets well-validated in molecular HFC calculations. A (14s11p6d)/[8s7p4d] basis was used for the transition metal.110 Only the most-diffuse sfunction from that work had to be removed to accommodate limitations of diffuse functions for compact solids in CP2K, reducing the size from [9s7p4d] to [8s7p4d]. Unmodified extended Huzinaga−Kutzelnigg-type IGLO-II99 basis sets were used for the main-group atoms. This possibility to use extended GTO basis sets is a distinct advantage of CP2K over previous FC shift calculations with CRYSTAL.30 To finally compute 7Li nuclear shieldings and, thus, NMR shifts of Li3V2(PO4)3, the computed EPR parameters and orbital shieldings were inserted into eq 1, using HFC tensors at various hybrid levels, as well as orbital shieldings and g-tensors at the PBE level, as described previously. The value of the Weiss constant (Θ = −37 K) has been taken from experiment.70 Incremental Cluster-Model Approach to g-Tensors for Periodic Solids. As an alternative to the periodic g-tensor calculations, we have evaluated also the use of an incremental scheme in which the g-tensor of the unit cell is built from clusters around a given metal center. This exploits the expected locality of g-tensors in insulators or semiconductors and opens the principal possibility to apply more easily the higher levels of theory for this property. The approach is analogous to the incremental scheme for electron correlation in solids.111 The cluster models were extracted from the optimized Li3V2(PO4)3 structures as shown in Figure 2. They include the core structure of the metal coordination, i.e., vanadium octahedrally coordinated by six tetrahedral phosphates (PO4). The oxygen valences of the phosphate groups were saturated with additional hydrogen atoms. The hydrogen-atom positions were optimized at PBE/def2-TZVP level using Turbomole.112 Because Li3V2(PO4)3 has two distinct vanadium environments, two sets of cluster models were constructed. For these, the number of lithium counterions was varied, including ions further and further removed from the phosphates. For comparability with the periodic calculations, g-tensors for each of the resulting clusters were computed with CP2K at PBE/def2-TZVP/IGLOII level in a supercell of dimension 40 Å × 40 Å × 40 Å. The resulting g-tensors were reinserted into the solid-state structures in appropriate orientation (using the corresponding

SO operators for g-tensors were done for a number of simple trigonal-planar transition-metal trifluoride complexes MF3 (M = Ti, V, Cr, Mn, and Fe), using the PBE functional and def2TZVP98 and IGLO-II basis sets for transition metals and fluorine atoms, respectively, based on B3LYP/6-311++G(d,p)optimized structures. Evaluation of g-Tensor Normalization for Periodic solids. To validate the normalization of g-tensors for multiple spin-centers in solids according to eq 13, two sets of test calculations were done. The first set has an increasing number of isolated trifluoride complexes in a supercell of dimension 40 Å × 40 Å × 40 Å using the B3LYP/6-311++g(d,p) structures for the individual molecules (see above). The second set consists of variation of the supercell dimension (and thus of the number of spin centers in the cell) for simple periodic solids, based on structures from X-ray crystallography. Here, we chose the four extended trifluorides MF3 (M = Ti, V, Mn, and Fe), with the transition metals octahedrally coordinated. The gtensor calculations were done at PBE level with IGLO-II and def2-TZVP basis sets for main-group and transition elements, respectively. Calculations on LixV2(PO4)3 (x = 3, 2.75, 2.5, 2.25, 2, and 1). The initial cell optimization (GPW level with PBE, GTH pseudopotentials, and DZVP-MOLOPT-SR basis sets) used a 2 × 2 × 2 supercell, leading to 512−640 atoms per cell, thus minimizing artifacts of the periodic boundary conditions (Figure 1). Subsequent GAPW optimization of atom positions at the PBE/def2-TZVPD/def2-TZVP level (the first basis designation holds for vanadium) used either a simple 1 × 1 × 1 unit cell with 72−80 atoms for the comparatively moreexpensive subsequent g-tensor and orbital shielding calculations or 2 × 2 × 1 supercells with 288−320 atoms for the moresensitive HFC calculations to ensure convergence with cell size. HFC values provided by CP2K have been normalized properly to the local spin state (see the Theory section). EPR and NMR parameters were computed for both experimental and optimized structures. g-Tensor calculations were done at the PBE/def2-TZVP/IGLO-II level. HFC tensors were calculated both at the PBE GGA level and using the PBEbased global hybrids, varying exact-exchange admixture from 5% (PBE5) to 40% (PBE40). This includes the well-known PBE0 functional (25%),90,109 termed PBE25 in this nomencla382


Article


Figure 3. Molecular test cases, (a) cobaltocene and (b) TEMPO radical, used to validate CP2K supercell results.

Figure 4. Relative (top) and absolute (bottom) deviations of isotropic HFCs for various main-group nuclei obtained with CP2K compared to Gaussian09 for cobaltocene (CoCp2) (left) and TEMPO (right). Results were collected with def2-QZVP basis set and BLYP or PBE0 functionals.

known computational package65 and used widely in the solidstate physics and quantum-chemistry community, it has so far not been used to compute NMR shifts for paramagnetic solidstate systems, although some experience for g-tensors66,67 and NMR shifts of diamagnetic systems113−115 exists. Before delving into explicit solid-state calculations, we thus need to establish the accuracy of the CP2K computations of HFC and g-tensors as well as orbital shieldings. This is best done by computations on isolated molecules in a sufficiently large supercell, in comparison with standard molecular quantum-chemical calculations using identical density functionals and basis sets. Let us thus focus on results for the doublet spin molecules TEMPO and cobaltocene (in staggered conformation; see Figure 3 and the Computational Details section), starting with isotropic HFCs, known to be very

rotational matrices), making use of point-group symmetry to replicate identical sites and to thus minimize the number of gtensor calculations needed. For example, the unit cell of Li3V2(PO4)3 has eight vanadium centers but only two distinct nonequivalent sites. Only two g-tensors thus need to be computed to generate the eight site g-tensors in the cell, which then lead size-consistently to the unit-cell g-tensor by suitable averaging. We note that the method is somewhat related to the approach in ref 59, in which site g-tensors were obtained from periodic calculations by replacing all but one paramagnetic metal ion in the cell by diamagnetic ions.

4. RESULTS AND DISCUSSION Validation of CP2K Results for Isolated Molecules against Quantum-Chemical Codes. While CP2K is a well383


Article

Journal of Chemical Theory and Computation Table 1. Comparison of g-Tensors between CP2K and ORCA for TEMPO and Cobaltocenea package CP2K ORCA

CP2K ORCA

functional

SO operator

BLYP PBE BLYP PBE PBE

Veff,PM Veff,PM SOMF SOMF AMFI

BLYP PBE BLYP PBE PBE

Veff,PM Veff,PM SOMF SOMF AMFI

gauge TEMPO IGAIM IGAIM common common common cobaltocene IGAIM IGAIM common common common

g11

g22

g33

giso

2.002 2.002 2.002 2.002 2.002

2.007 2.006 2.006 2.006 2.006

2.010 2.009 2.009 2.009 2.009

2.006 2.006 2.006 2.006 2.006

1.833 1.831 1.858 1.848 1.848

2.018 2.017 2.017 2.016 2.016

2.065 2.063 2.058 2.057 2.057

1.972 1.970 1.977 1.974 1.974

a

BLYP and PBE results with def2-QZVP basis sets. SOMF operators were used with ORCA and eq 8 with CP2K. The common gauge origin was placed at the center of charge.

sensitive to various computational parameters. Figure 4 compares relative and absolute deviations of isotropic HFCs (Fermi-contact contribution) for various main-group nuclei in both test molecules with CP2K against Gaussian09 as reference, using either the BLYP GGA or the PBE0 hybrid functional. Relative deviations are generally below 1%, with absolute deviations below 0.05 MHz (ca. 0.1% MHz for the nitroxide nitrogen atom in TEMPO). This confirms that (a) the use of a supercell in CP2K does not cause any notable deviations, (b) the GAPW ansatz in CP2K faithfully reproduces the Gaussian09 results, and (c) any deviations due to the different integration grids used in the two codes also remain minor. Moving on to g-tensors (Table 1), in addition to those computational parameters discussed above, the nature of the SO operators comes into play. A comparison is made to calculated g-tensors obtained with ORCA using the accurate spin−orbit mean-field (SOMF) approximation for the SO matrix elements. In contrast, CP2K utilizes the Mauri−Pickard SO approximations57 from eq 8. Agreement is good in both cases, but the CP2K results generally give slightly larger deviations from the free-electron g-value than the ORCA data. This is not unexpected, as the treatment of the SOO term is somewhat less complete in CP2K, which is thought to cause a slight overestimate of the SO matrix elements. It is important to point out, however, that previous tests suggest that this overestimate is much smaller than for other solid-state g-tensor implementations.67 CP2K is thus well-suited to also give us the needed g-tensor, with the current caveat that we cannot yet use hybrid functionals. Table 2 shows the results for triplet-state VF3, a metal complex reasonably close to a local environment of the solidstate systems in focus in this work (see below). MAG allows the use of a number of different SO operators, ranging from a full treatment (FULL) of the Breit−Pauli SO integrals via the atomic meanfield approximation (AMFI; ORCA’s SOMF may be viewed as intermediate between these two treatments) to Koseki’s semiempirical effective one-electron SO operators (1el.-eff.)116−118 to the bare one-electron SO term (1el.). Focusing on the GIAO-based results, we see that AMFI deviates relatively little from the full treatment (in spite of some known incompatibility between AMFI and GIAO),87 whereas 1el.-eff. and 1el. overshoot deviations from the free-electron gvalue drastically for this system (as two-electron SO contributions diminish the overall SO matrix elements, their incomplete treatment causes an overestimate). The IGAIM-

Table 2. Comparison of Different Codes, Gauges, and SO Operators for g-Tensor Calculationsa on VF3 package

gauge

CP2K MAG

IGAIM COMMON

GIAO

ORCA

COMMON

SO operator

g11

g22

g33

giso

Veff,PM FULL AMFI 1el.-eff. 1el. FULL AMFI 1el.-eff. 1-el. SOMF AMFI 1el.-eff.

1.874 1.975 1.968 1.926 1.931 1.911 1.888 1.751 1.767 1.954 1.954 1.940

1.899 1.986 1.981 1.958 1.960 1.927 1.909 1.804 1.815 1.963 1.963 1.953

1.976 1.996 1.994 1.984 1.985 1.982 1.977 1.946 1.950 1.992 1.992 1.990

1.916 1.986 1.981 1.956 1.959 1.940 1.924 1.834 1.844 1.970 1.970 1.961

a

PBE results with def2-TZVP and IGLO-II basis sets for vanadium and fluorine, respectively. See text for the abbreviations for SO operators and choices of gauge origin. Supercell calculations for CP2K, see text.

based CP2K results with the eq 8 SO operator are close to the GIAO-based MAG results with AMFI, confirming both the reasonable SO operator of CP2K and the good accuracy of the IGAIM gauge compared to GIAO. Much larger deviations are apparent for calculations that have a common gauge at the metal center: here the match between ORCA and MAG data is negligible, within the expected differences between codes, SO-operators, grids, and so on. However, deviations of common-gauge calculations from both GIAO- and IGAIM-based data are notable. This limits the possibility to obtain accurate solid-state g-tensors by inserting molecular model-cluster calculations with ORCA in an incremental scheme, at least for DFT treatments (see below), and points to a need to have a careful eye on the gauge-origin question. Notably, the common-gauge calculations tend to provide smaller deviations from the free-electron value than the GIAO- or IGAIM-based results, an observation we also make for related larger model vanadium complexes (see below). Table S4 has detailed computed values for other MF3 (M = Ti, Cr, Mn, and Fe) molecules, comparing different codes, gauges, and SO operators for g-tensor calculations. Turning finally to the orbital shieldings, Figure 5 shows for both cobaltocene and TEMPO that IGAIM supercell calculations with CP2K agree very well with GIAO-based 384


Article


Figure 5. Relative and absolute deviations of CP2K supercell IGAIM calculations of orbital shieldings from Gaussian09 GIAO data for cobaltocene (left) and TEMPO radical (right), at BLYP/def2-QZVP level.

calculations using Gaussian09: relative deviations are generally below 1% (slightly above for the small oxygen shielding in TEMPO), and absolute deviations are below 1 ppm. We can thus conclude from these molecular test calculations that all contributions we need for the doublet-state formalism of pNMR shift calculations (eq 1) are obtainable with CP2K in good accuracy compared with standard molecular quantumchemical codes. One possible limitation is the use of only GGAtype functionals for the g-tensor calculations. It is known that such functionals tend to delocalize too much spin density onto the ligands, leading to underestimated SO contributions from the metal center and thus too small g-tensors.119 This is why we will also explore the use of an incremental scheme based on molecular cluster g-tensor calculations for the different metal sites in the unit cell, as this opens more-facile access to the use of hybrid functionals or even of post-Hartree−Fock methods. Validation of g-Tensor Normalization for More than One Spin Site in a Unit Cell. Here, we validate the correctness of the normalization to the number of spin centers in eq 13. We first do this by placing an increasing number of molecular metal trifluoride complexes MF3 (Ti, V, Cr, Mn, and Fe) into a large supercell of 40 Å × 40 Å × 40 Å, so as to minimize interactions between the molecules. Calculations for up to four molecules per unit cell at PBE/def2-TZVP/IGLO-II level are given in Table S5 (see structural arrangements in Tables S15−S20). Indeed, in spite of the different spin states of the various metal fluorides, normalization renders the g-tensors for the different numbers of spin centers exactly equal, thus confirming the validity of eq 13 for well-separated molecules. Table 3 extends this type of test to extended solids, as exemplified by the solid trifluoride materials, using the same PBE/def2-TZVP/IGLO-II level as for the molecular test. Cell parameters and structures are taken from experimental XRD data,120−123 which were used to generate supercells of increasing sizes, having between 24 and 54 spin centers per supercell. Interaction between the spin centers is now of course more pronounced than for the more or less isolated molecules in the previous test case, and we cannot expect the normalized g-tensors (eq 13) to stay constant as we extend the supercell size. However, the values vary over only a very small range (supercell g-tensors exhibit relatively small anisotropy), thus confirming the procedure and providing the necessary setup for g-tensor calculations in supercells containing more than one spin center. Validation for and Application to Solid Lithium Vanadium Phosphates. Here we focus in particular on g-

Table 3. Comparison of Normalized Solid-State g-Tensors for Different Supercell Sizes of Extended MF3 (M = Ti, V, Mn, or Fe) Periodic Solids Using CP2K solid TiF3 VF3

MnF3

FeF3

unita

nb

g11

g22

g33

gisoc

321 231 331 322 232 121 131 221 221 321 231 331

36 36 54 24 24 24 36 48 24 36 36 54

1.941 1.955 1.945 1.924 1.918 1.975 1.973 1.977 2.026 2.025 2.026 2.024

1.956 1.954 1.957 1.900 1.903 1.987 1.988 1.987 2.026 2.026 2.026 2.024

1.951 1.954 1.954 1.906 1.894 1.987 1.987 1.989 2.025 2.026 2.025 2.023

1.949 1.954 1.952 1.910 1.905 1.983 1.983 1.984 2.025 2.026 2.025 2.024

a

Digits represent the number of unit-cell replica into the direction of the unit-cell vectors. bNormalization factor reflects the number of paramagnetic spin centers in the super cell. cEq 13 has been used for calculating the normalized g-tensor for periodic solids at the PBE/ def2-TZVP/IGLO-II level.

tensor computations for the series of lithium vanadium phosphates, LixV2(PO4)3 (x = 3, 2.75, 2.5, 2.25, 2, and 1), comparing full periodic calculations and the incremental cluster model. Subsequently, the full 7Li pNMR shifts of Li3V2(PO4)3 are analyzed. Table 4 provides results of the PBE/DZVP-MOLOPT-SRGTH GPW optimizations (cf. the Computational Details section) of unit-cell parameters in comparison to highresolution single-crystal XRD data of these phosphate materials.124 The optimizations were done for a 2 × 2 × 2 supercell, with dimensions of more than 15 Å in each direction, to minimize artifacts from the periodic-boundary conditions. Agreement to within ca. 1−2% in the cell parameters is generally achieved. Either these optimized parameters or the experimental cell parameters were used subsequently to optimize the atom positions at all-electron GAPW PBE/def2TZVPD/def2-TZVP level (cf. the Computational Details section). The same procedure has been used to optimize the structures for the more complicated Li2.75V2(PO4)3 and Li2.25V2(PO4)3 systems, for which experimental structural data are lacking. We note in passing the rather large unit-cell sizes and number of atoms (GAPW for atom positions: >300 atoms for 2 × 2 × 1 supercell and GPW for cell parameters: 2 × 2 × 2 385


Article


Table 4. Comparison of Optimized and XRD Unit-Cell Parameters and Volume for LixV2(PO4)3 (x = 3, 2.75, 2.5, 2.25, 2, and 1)a material Li3V2(PO4)3

Li2V2(PO4)3

LiV2(PO4)3 Li2.75V2(PO4)3 Li2.5V2(PO4)3 Li2.25V2(PO4)3 a

structure

a (Å)

b (Å)

c (Å)

β (deg)

volume (Å3)

124

8.597 8.611 0.16 8.461 8.529 0.80 8.301 8.382 0.98 8.624 8.604 8.614

8.592 8.606 0.16 8.619 8.645 0.30 8.518 8.560 0.49 8.617 8.636 8.654

12.001 12.098 0.80 11.892 11.965 0.61 11.653 11.745 0.79 12.015 11.956 11.905

90.57 90.57 0.00 90.00 90.00 0.00 89.61 89.61 0.00 90.57 90.57 90.57

886.415 896.448 1.13 867.228 882.129 1.72 823.833 842.702 2.29 892.834 888.240 887.288

XRD OPT percent difference XRD124 OPT percent difference XRD125 OPT percent difference OPT OPT OPT

Cell parameter optimizations at GPW-DZVP-MOLOPT-SR-GTH level for a 2 × 2 × 2 supercell (monoclinic cell, P21/n space group).

molecular units in the unit cell for these materials, giving compositions Li9V8(PO4)12, Li10V8(PO4)12, and Li11V8(PO4)12. This has subtle structural consequences that need to be taken into account in the NMR and EPR parameter calculations. Based on these structures, we may now compute the various contributions to the pNMR shifts, starting with the g-tensors obtained using the normalization to the number of spin centers in the cell, as discussed above. The PBE/def2-TZVP/IGLO-IIbased results for all lithium vanadium phosphate materials studied here are provided in Table 5, using both optimized and XRD structures where available and only the optimized ones in the remaining cases. While only scarce experimental g-tensor data are available for these materials (from Figure 4 of ref 128 we infer g|| = 1.9260, g⊥ = 2.0003 for Li2V2(PO4)3 and Li2.5V2(PO4)3), the overall magnitudes compare well to related VIII, and VIV sites with octahedral coordination by oxygen atoms (see, e.g., ref 126.), even though we have to assume that the PBE functional underestimates the deviations from the freeelectron value (slightly overestimated SO matrix elements may provide some error compensation). Interestingly, the predicted deviations from the free-electron value diminish with an increase of the number of VIV sites (Table 5). In view of these limitations regarding the functional in PBCbased g-tensor calculations, we evaluate additionally in Table 6 the performance of an incremental cluster model for obtaining the g-tensor in Li3V2(PO4)3. That is, for both structurally distinct vanadium centers in the unit cell, we have cut clusters with one vanadium site and the six coordinated phosphates

supercell optimizations, >620 atoms), made possible by the efficient CP2K implementation. Figure 6 shows for the example

Figure 6. Comparison of XRD (in yellow) and optimized (in blue) structures of Li3V2(PO4)3.

of Li3V2(PO4)3 the excellent agreement between XRD (in yellow) and fully optimized (in blue) structures (a similar agreement is found for LiV2(PO4)3 and Li2V2(PO4)3; see Figure S3). This gives us confidence in the optimized structures for those systems, where XRD data are not available, i.e., Li2.25V2(PO4)3 and Li2.75V2(PO4)3. Note the presence of four

Table 5. Computed g-Tensors for LixV2(PO4)3 (x = 3, 2.75, 2.5, 2.25, 2, and 1)a,b material Li3V2(PO4)3 Li2.75V2(PO4)3 Li2.5V2(PO4)3 Li2.25V2(PO4)3 Li2V2(PO4)3 LiV2(PO4)3

structure

V+III

V+IV

g11

g22

g33

giso

Δc

ηd

XRD OPT OPT OPT OPT XRD OPT XRD OPT

8 8 7 6 5 4 4 0 0

0 0 1 2 3 4 4 8 8

1.834 1.864 1.875 1.871 1.880 1.897 1.895 1.877 1.906

1.870 1.875 1.891 1.885 1.901 1.909 1.913 1.911 1.919

1.902 1.915 1.914 1.917 1.911 1.935 1.935 1.913 1.919

1.869 1.885 1.894 1.891 1.897 1.914 1.914 1.901 1.915

−0.05 0.05 0.03 0.04 −0.03 0.03 0.03 −0.03 −0.01

0.89 0.34 0.74 0.56 0.58 0.55 0.87 0.09 0.00

a

Unit-cell g-tensor obtained after normalization (eq 13). Columns V+III and V+IV show the number of vanadium atoms present in the unit cell having oxidation states of + III and + IV, respectively (with an overall total of eight vanadium sites present). bComputations on 1 × 1 × 1 cells at the PBE/ def2-TZVP/IGLO-II level. cg-Tensor anisotropy, Δ = gzz − (1/2)(gxx + gyy) the components ordered as | gzz − giso | ≥ | gyy − giso | ≥ | gxx − giso |, giso = (g11 + g22 + g33)/3 and g11 ≤ g22 ≤ g33. dAsymmetry parameter, η= (gxx − gyy)/(gxx − giso). 386


Article

Journal of Chemical Theory and Computation Table 6. Comparison of Calculateda Site g-Tensors for Clusters,b Unit-Cell g-Tensors Obtained from Cluster Modelsc, and PBC Unit-Cell g-Tensorsd for Li3V2(PO4)3

model VA

model VB

unit-cell model

PBC results

no. of Li atomse

Li−O distancef (Å)

g11

g22

g33

giso

0 1 2 3 4 0 1 2 3 0 (0−0)

− 1.94 1.96 1.97 2.09 − 1.97 1.94 2.02 −

1.792 1.788 1.828 1.821 1.820 1.839 1.834 1.843 1.838 1.822

1.832 1.847 1.856 1.844 1.852 1.839 1.857 1.856 1.860 1.836

1.929 1.931 1.933 1.926 1.936 1.945 1.962 1.956 1.943 1.929

1.851 1.855 1.873 1.864 1.870 1.874 1.884 1.885 1.881 1.862

3 (1−2) 5 (3−2) 7 (4−3) unit cell

− − − −

1.825 1.858 1.848 1.864

1.848 1.868 1.874 1.875

1.936 1.896 1.901 1.915

1.870 1.874 1.875 1.885

Computational Details section), together with hybrid functionals. Using IGLO-II and 8s7p4d basis sets for main-group elements and vanadium, respectively (see the Computational Details section), we thus evaluate the influence of EXX admixture for PBE-based hybrids with EXX admixtures between 5% and 40% for the 2 × 2 × 1 supercell from both XRD- and DFT-optimized structures. As shown in Figure 7, there are three distinct lithium sites in the structure (Lia, Lib, and Lic) and two distinct vanadium sites (V1 and V2) contributing to the hyperfine shifts. Lia is in a tetrahedral site sharing edges with two V1 octahedra, Lib is in a distorted trigonal bipyramidal site sharing a face with V2 octahedra and two corners with V1 and V2 octahedra, and Lic is also in a distorted trigonal bipyramidal site sharing a face with V2 octahedra and an edge with V1 octahedra.127,128 Isotropic (AFC) and dipolar (Adip) HFCs for all three lithium sites are provided in Tables S6 and S7. The AFC values decrease as AFC(Lia) > AFC(Lic) > AFC(Lib) for both XRD and optimized structures. In all cases, the computed AFC decreases almost linearly from PBE5 to PBE40, consistent with more-localized spin density around the metal center. We will see below from the pNMR shift analyses that significant EXX admixtures are needed for “best” agreement with experiment. Thus, while the computation of HFCs with hybrid functionals is appreciably more expensive than with the PBE GGA, this is clearly warranted by the results (see below). The AFC for the Lia and Lib sites is larger when computed at the XRD rather than at the optimized structure, while matters are reversed for Lic. The Adip tensor is not only less-sensitive to the EXX admixture, as expected, but also depends less on structure (XRD versus optimized). The asymmetry parameter η (see the Computational Details section) remains almost constant with the functional, and structural effects are also small. The orbital shifts are directly taken relative to the lithium lactate reference (PBE/Ahlrichs-VTZ/IGLO-II level with IGAIM; cf. Table S8). They are less affected by the nature of the spin centers but are characteristic of the overall structure of the solid. While they are usually small compared to the FC shifts, they may play an important role when the latter are small. Remembering the small differences between XRD and optimized structures (cf. Figure 6), we nevertheless see effects of these changes on the orbital shifts: while for the XRD orb orb structure δ(Liorb a ) < δ(Lic ) < δ(Lib ), for the optimized orb orb structure, we see δ(Lia ) < δ(Lib ) < δ(Liorb c ), i.e., a switch between Lib and Lic. The fact that the Lia orbital shift is less affected reflects its more rigid environment, while the other two Li positions vary more. Indeed, in the delithiation process of Li3V2(PO4)3, the first lithium atoms are removed from the Lic site followed by the Lib sites.128 Now having everything needed for the 7Li pNMR shieldings and shifts within the doublet formalism (lacking ZFS and currently also SO contributions to HFCs, which, however, are expected to be small), we may now proceed to discuss the total lithium shifts for Li3V2(PO4)3. This then also allows an evaluation of the relative importance of the different terms, of the importance of the functional on the HFC contributions and a comparison of PBC-based and cluster-model-based g-tensor calculations. For comparison with experimental values,128 we report shifts at 320 K temperature (a value slightly higher than room temperature being used to account for frictional heating due to magic-angle spinning). The vanadium ions exhibit the +III oxidation state (d2 configuration with local triplet spin) and an approximately octahedral coordination. We obtain seven

a

PBE/IGLO-II/def2-TZVP results. bCluster model with variable number of Li counterions; see the Computational Details section. c g-Tensor for unit cells obtained by replacing the cluster-model gtensors into the solid-state structure after averaging. dFull calculation of unit-cell g-tensor from periodic calculations (using proper normalization). eNumber of Li counterions included in cluster model. fShortest Li−O distance from the VO6 octahedron.

(terminated by hydroxy functions) from the solid (Figure 2). To include the effect of the lithium counterions, a variable number of ions were included (within 2.2 Å from the oxygen atoms in the direct vanadium coordination sphere). Up to four lithium ions had to be considered for model VA and up to three lithium ions for model VB (Figure 2). To be able to compare directly the g-tensor results for the unit cell obtained by placing the molecular g-tensors for these models into the solid-state structure against the full periodic calculations, we used the same PBE functional and the same basis sets and did the cluster computations in a large super cell in CP2K. Starting with the individual site g-tensors (Table 6), we see that for model VA, the overall deviations from the free-electron value (cf. giso) diminish slightly upon the inclusion of more Li counterions, while the opposite effect holds true for model VB. The generation of the overall unit-cell g-tensor shows a clearer effect of the counterions: as more of them are added, g11 and g22 increase, and g33 decreases somewhat. Overall, agreement with the full PBC result is good, suggesting the cluster modeling to be an adequate alternative that opens avenues toward improved electronic-structure levels from molecular computations. Likely, protonation of phosphate oxygen atoms and addition of Li ions as a means of saturating valencies in the clusters is a sufficiently small perturbation to the central part of the cluster to still provide accurate local g-tensors. The second and arguably most-important individual aspect contributing to the pNMR shifts in the lithium vanadium phosphates are the 7Li HFCs. The isotropic HFCs enter the often dominant FC term, whereas the HFC anisotropies influence not only the shift anisotropies but also the isotropic shift via the PC term. It is the isotropic HFCs that tend to be most sensitive to the computational level, however. Here, the strength of CP2K is the possibility to use extended Gaussiantype all-electron basis sets within the GAPW scheme (see the 387


Article


Figure 7. Local environment of three Li and two V sites in Li3V2(PO4)3. (a) Lia, (b) Lib, (c) Lic, (d) V1, and (e) V2; see the text for more details.

still part of the FC-shift. Obviously, this second term should not be neglected for materials with significant g-shifts. In the absence of ZFS, the PC-shift is represented by the third term, δΔgA ̃ dip. It reflects the trace of the interaction between the anisotropic parts of the HFC and g-tensors. Because Adip depends little on the functional (see above) and the g-tensor has been obtained at PBE level, this term remains almost constant for a given site. While it may seem small at first sight, it may become relevant for cases of small FC shifts. The difference between PC shifts for Lib and Lic is 7.1 and 4.9 ppm for XRD and optimized structures, respectively. Considering more complex systems with lower symmetry like Li2.75V2(PO4)3 (see the 7Li MAS NMR in Figure 7 of ref 128), shift separations may be less than 20 ppm, and the inclusion of PC shifts may be crucial for assigning the NMR signals. As we saw already earlier, the orbital shifts (fourth term) may differ between sites and for different structures. They can be negative or positive (for the given reference). While in the present case they are also smaller than the FC shifts, it is clear that their inclusion may become important when FC shifts are small. The question of which functional used for AFC provides the best agreement of total

terms contributing to the shift tensor and four to the isotropic part (eq 2). δge AFC and δΔgisoAFC make up the Fermi-contact shift, and δΔgA ̃ dip (the pseudocontact shift), and we have δorb as the orbital shift. Figure 8 illustrates the contributions to the isotropic 7Lishifts at various levels, comparing the total shifts with experimental values128 for the three lithium sites (Tables S9 and S10 in Supporting Information provide detailed numerical comparisons of all terms). For better contrast, Figure 9 (bottom) compares the total isotropic shifts for both XRD and optimized structures. Assignment of the shifts to the corresponding sites is clearly possible, giving the order δLia > δLic > δLib. Just by visual inspection of the structure, it would appear difficult to make this assignment, pointing to the usefulness of explicit computations. The almost-linear dependence of AFC on EXX admixture in the functional (see above) translates into a similar dependence of the total isotropic shifts, dominated by the δgeAFC term. Deviations of the isotropic gvalue from the free-electron value combine with AFC to give in most cases the second-largest contribution (only when δgeAFC is clearly the largest; see later discussion), δΔgisoAFC, which is 388


Article


Figure 8. Comparison of 7Li chemical shifts of three distinct lithium sites, Lia (top), Lib (middle), and Lic (bottom) for Li3V2(PO4)3 computed for both XRD and optimized (OPT) structures. Variations with EXX admixture to PBE-based functionals for the HFC tensors are shown. g-Tensor and orbital shielding obtained at PBE level. Shieldings converted to shifts relative to solid lithium lactate (LiC3H5O3). Experimental values come from ref 128 (see Tables S9 and S10 for numerical values).

pNMR shifts with experiment depends somewhat on the input structure, and results also differ between the three Li sites (Figures 8 and 9). For example, for Lia use of the XRD structure would suggest best performance of PBE25 (PBE0, 25% EXX admixture), while 15% EXX admixture might be

indicated upon using the optimized structure. A similar dependence holds for Lib but with a shift to lower optimum EXX admixtures (Figure 8). Finally, matters appear reversed for Lic, in which the use of the XRD structure would suggest 15% to 20% EXX admixture as optimal, while use of the optimized 389


Article


depends also much less on EXX admixture than the isotropic HFC, it is not surprising that the shift anisotropies, which tend to be dominated by the HFC anisotropies, also show a small dependence on the functional. For all three lithium sites, use of 15% to 40% EXX admixture provides good agreement with experiment. We keep in mind that the g-tensor was computed at PBE level (it appears in the anisotropic terms and is also included in the determination of the Weiss constant; see earlier discussion). Other choices might affect shift anisotropies to some extent.

5. CONCLUSIONS This work shows that the spin Hamiltonian terms needed to compute NMR shifts for paramagnetic solids based on density functional theory can be obtained in an efficient way from the CP2K program package. In fact, the extensive parallelization, general computational efficiency, the possibility to use extended Gaussian basis sets including relatively diffuse functions, and the possibility to use hybrid functionals efficiently for the computation of hyperfine tensors as well as the availability of gtensors (with reasonable spin−orbit operators) and orbital shieldings with distributed gauges (IGAIM and CSGT) for paramagnetic solids make CP2K a very valuable tool for future work in the field. Limitations apply currently to the use of only semilocal functionals for g-tensors and orbital shifts and the absence of zero-field splittings and spin−orbit corrections to hyperfine couplings. We also found that current standard solid-state implementations of g-tensors are geared toward one isolated defect per unit cell. Application to extended solids with more than one openshell site per cell requires thus a suitable normalization of gtensors. We could furthermore show that an incremental cluster model that constructs the g-tensor of the simulation cell from individual molecular g-tensors is a viable and accurate alternative to computations using periodic boundary conditions. This opens the possibility of using advanced molecular electronic-structure methods to construct a suitable g-tensor for a solid (insulator or semiconductor), making use of its relative locality. Similar approaches seem attractive to obtain zero-field splittings. We note that the gauge of the magnetic vector potential for g-tensors had a larger impact on the present results than anticipated. This also has to be kept in mind for future work. This paper has established the basis for reliable computations using CP2K by first comparing the various contributions to pNMR shifts (hyperfine and g-tensors and orbital shieldings) for molecules (computed in a supercell) against wellestablished quantum-chemical codes. Second, the various computational parameters have been evaluated for simple extended solid-state fluorides and, third, for more-complex ternary lithium vanadium phosphates. We have furthermore provided an accurate 7Li NMR shift standard for solid-state computations by modeling the absolute shieldings in solid lithium lactate. Finally, all relevant contributions to the 7Li pNMR shifts of the three lithium sites in Li3V2(PO4)3 have been examined in detail. While the Fermi-contact shifts dominate, it is clear that pseudocontact and orbital shifts may not be negligible for accurate calculations (in particular if contact shifts are small), and deviations of the isotropic g-value from ge may affect even the contact terms. The interrelation between input structure (XRD versus DFT-optimized) and computation of the pNMR parameters (exact-exchange admixture for hybrid functionals used for the hyperfine

Figure 9. Comparison of computed total 7Li chemical shifts for Li3V2(PO4)3 as a function of the EXX admixture to the functional for the HFC tensors. Top: comparison of results with g-tensors computed within PBC (solid line) or obtained from cluster models (dashed line). Bottom: the effect of the input structure (XRD vs optimized) is probed. g-Tensor and orbital shielding at the PBE level. Shifts relative to solid lithium lactate (LiC3H5O3). Experimental values come from ref 128 (see Tables S9−S11 for numerical values).

structure requires 35% to 40% for best agreement. This shows the subtle interdependence between structure, functional and further aspects of the computations (e.g., the level used for gtensor and orbital shift, SO operators, etc.). Taking, for the sake of argument, PBE25 (PBE0, 25% EXX admixture) as functional used for HFC computations, we arrive at the conclusion that with the XRD structure the shift for Lia is described well, whereas those for Lib and Lic are underestimated by 5−10 ppm. Using instead the optimized structure, the shifts for Lia and Lib are somewhat low (by approximately 10 ppm), and that for Lic 10 ppm too large. In any case, it can be stated clearly that with essentially all hybrid functionals, the order of shifts between the three sites is reproduced faithfully. A comparison of 7Li-shifts with PBC-based and clustermodel-based g-tensors shows essentially negligible differences (Figure 9, top), suggesting the incremental cluster scheme as a viable alternative (also for our ongoing efforts toward including ZFS contributions). Further information on isotropic shifts and shift anisotropies are found in Tables S11 and S12. A total of five terms contribute to the shift-tensor anisotropies (see eq 3). Detailed numerical results are provided in Tables S12−S14. Interestingly, the shift anisotropy (Δ) and asymmetry parameter (η) depend much less on the choice of XRD versus the optimized structures, likely because they are dominated by the less-sensitive Adip. In general, the clear trends ΔLia < ΔLib < ΔLic and ηLia ≈ ηLib > ηLic are found. Because Adip 390


Article


(2) Bertini, I.; Luchinat, C.; Parigi, G. Solution NMR of paramagnetic molecules; Elsevier: Amsterdam, The Netherlands, 2001. (3) Helm, L. Relaxivity in paramagnetic systems: Theory and mechanisms. Prog. Nucl. Magn. Reson. Spectrosc. 2006, 49, 45−64. (4) Cole, K. E.; Paik, Y.; Reeder, R. J.; Schoonen, M.; Grey, C. P. 2H MAS NMR studies of deuterated goethite (α-FeOOD). J. Phys. Chem. B 2004, 108, 6938−6940. (5) Grey, C. P.; Dupré, N. NMR studies of cathode materials for lithium-ion rechargeable batteries. Chem. Rev. 2004, 104, 4493−4512. (6) Nielsen, U. G.; Paik, Y.; Julmis, K.; Schoonen, M. A. A.; Reeder, R. J.; Grey, C. P. Investigating sorption on iron-oxyhydroxide soil minerals by solid-state NMR spectroscopy: A 6Li MAS NMR study of adsorption and absorption on goethite. J. Phys. Chem. B 2005, 109, 18310−18315. (7) Kowalewski, J.; Kruk, D.; Parigi, G. NMR relaxation in solution of paramagnetic complexes: recent theoretical progress for S ≥ 1. Adv. Inorg. Chem. 2005, 57, 41−104. (8) Kim, J.; Nielsen, U. G.; Grey, C. P. Local environments and lithium adsorption on the iron oxyhydroxides lepidocrocite (γFeOOH) and goethite (α-FeOOH): A 2H and 7Li solid-state MAS NMR study. J. Am. Chem. Soc. 2008, 130, 1285−1295. (9) Kaupp, M.; Köhler, F. H. Combining NMR spectroscopy and quantum chemistry as tools to quantify spin density distributions in molecular magnetic compounds. Coord. Chem. Rev. 2009, 253, 2376− 2386. (10) Mackenzie, K. J. D.; Smith, M. E. Multinuclear solid-state NMR of inorganic materials; Elsevier Science: Oxford, U.K., 2002. (11) Moon, S.; Patchkovskii, S. First-principles calculations of paramagnetic NMR shifts. In Calculation of NMR and EPR parameters. Theory and applications; Kaupp, M., Bühl, M., Malkin, V. G., Eds.; Wiley-VCH, Weinheim, Germany, 2004. (12) Hrobárik, P.; Reviakine, R.; Arbuznikov, A. V.; Malkina, O. L.; Malkin, V. G.; Köhler, F. H.; Kaupp, M. Density functional calculations of NMR shielding tensors for paramagnetic systems with arbitrary spin multiplicity: Validation on 3d metallocenes. J. Chem. Phys. 2007, 126, 024107. (13) Vaara, J. High-Resolution NMR spectroscopy; Contreras, R. H., Ed.; Elsevier: Amsterdam, The Netherlands, 2013. (14) Pennanen, T. O.; Vaara, J. Nuclear magnetic resonance chemical shift in an arbitrary electronic spin state. Phys. Rev. Lett. 2008, 100, 133002. (15) Autschbach, J.; Patchkovskii, S.; Pritchard, B. Calculation of hyperfine tensors and paramagnetic NMR shifts using the relativistic zeroth-order regular approximation and density functional theory. J. Chem. Theory Comput. 2011, 7, 2175−2188. (16) Van den Heuvel, W.; Soncini, A. NMR chemical shift in an electronic state with arbitrary degeneracy. Phys. Rev. Lett. 2012, 109, 073001. (17) Komorovský, S.; Repiský, M.; Ruud, K.; Malkina, O. L.; Malkin, V. G. Four-component relativistic density functional theory calculations of NMR shielding tensors for paramagnetic systems. J. Phys. Chem. A 2013, 117, 14209−14219. (18) Pennanen, T. O.; Vaara, J. Density-functional calculations of relativistic spin-orbit effects on nuclear magnetic shielding in paramagnetic molecules. J. Chem. Phys. 2005, 123, 174102. (19) Schmitt, S.; Jost, P.; van Wüllen, C. Zero-field splittings from density functional calculations: Analysis and improvement of known methods. J. Chem. Phys. 2011, 134, 194113. (20) Soncini, A.; Van den Heuvel, W. Communication: Paramagnetic NMR chemical shift in a spin state subject to zero-field splitting. J. Chem. Phys. 2013, 138, 021103. (21) Kubica, A.; Kowalewski, J.; Kruk, D.; Odelius, M. Zero-field splitting in nickel(II) complexes: A comparison of DFT and multiconfigurational wavefunction calculations. J. Chem. Phys. 2013, 138, 064304. (22) Van den Heuvel, W.; Soncini, A. NMR chemical shift as analytical derivative of the Helmholtz free energy. J. Chem. Phys. 2013, 138, 054113.

couplings) has been clearly exposed. Differences between the shifts of the three individual lithium sites are sufficient to get a clear computational assignment. The computations furthermore provide the entire shift tensor for comparison with suitable experimental data, at levels that go beyond previous approaches by including more of the relevant contributions. Our ongoing work focuses on the inclusion of zero-field splitting, which is expected to appreciably affect the pNMR shifts of systems with high-spin ions, particularly for later 3d elements such as Co or Ni.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00991. Tables and figures providing additional information on the comparison of different codes, gauges and SO operators, Cartesian coordinates, a comparison of XRD and optimized structures of lithium vanadium phosphate, all structures of lithium vanadium phosphate solids, detailed lithium site HFC tensor data, unit-cell and cluster model structures of lithium lactate, complete results of lithium orbital shift tensors, numerical data on total lithium isotropic shifts and contributions, and detailed lithium chemical shift tensor data. (PDF) Crystallographic information file of lithium L-lactate (CIF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Michael W. Gaultois: 0000-0003-2172-2507 Marcella Iannuzzi: 0000-0001-9717-2527 Clare P. Grey: 0000-0001-5572-192X Martin Kaupp: 0000-0003-1582-2819 Funding

This work has been funded by the EU 7th Framework Programme within the Marie Curie Actions Initial Training Network scheme under grant agreement no. 317127, the “pNMR project”. M.W.G. is grateful for support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska−Curie grant agreement no. 659764. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to Dr. Ladislav Benda, Sebastian Gohr, Ieuan D. Seymour, Roberta Pigliapochi, and the members of ITN− pNMR network for useful discussions and to Dr. Andrew Bond at the University of Cambridge for assistance with single crystal data collection, reduction, and solution. The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC resources that have contributed to the research results reported in this paper.

■

REFERENCES

(1) Berliner, L. NMR of paramagnetic molecules; Reuben, J., Ed.; Plenum Press: New York, 1993. 391


Article

Journal of Chemical Theory and Computation (23) Martin, B.; Autschbach, J. Temperature dependence of contact and dipolar NMR chemical shifts in paramagnetic molecules. J. Chem. Phys. 2015, 142, 054108. (24) Gendron, F.; Sharkas, K.; Autschbach, J. Calculating NMR chemical shifts for paramagnetic metal complexes from first-principles. J. Phys. Chem. Lett. 2015, 6, 2183−2188. (25) Gendron, F.; Autschbach, J. Ligand NMR Chemical Shift Calculations for Paramagnetic Metal Complexes: 5f1 vs 5f2 Actinides. J. Chem. Theory Comput. 2016, 12, 5309−5321. (26) Vaara, J.; Rouf, S. A.; Mareš, J. Magnetic couplings in the chemical shift of paramagnetic NMR. J. Chem. Theory Comput. 2015, 11, 4840−4849. (27) Martin, B.; Autschbach, J. Kohn-Sham calculations of NMR shifts for paramagnetic 3d metal complexes: Protocols, delocalization error, and the curious amide proton shifts of a high-spin iron(II) macrocycle complex. Phys. Chem. Chem. Phys. 2016, 18, 21051−21068. (28) Carlier, D.; Ménétrier, M.; Grey, C. P.; Delmas, C.; Ceder, G. Understanding the NMR shifts in paramagnetic transition metal oxides using density functional theory calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 174103. (29) Middlemiss, D. S.; Ilott, A. J.; Clément, R. J.; Strobridge, F. C.; Grey, C. P. Density functional theory-based bond pathway decompositions of hyperfine shifts: equipping solid-state NMR to characterize atomic environments in paramagnetic materials. Chem. Mater. 2013, 25, 1723−1734. (30) Kim, J.; Middlemiss, D. S.; Chernova, N. A.; Zhu, B. Y. X.; Masquelier, C.; Grey, C. P. Linking local environments and hyperfine shifts: A combined experimental and theoretical 31P and 7Li solid-state NMR study of paramagnetic Fe(III) phosphates. J. Am. Chem. Soc. 2010, 132, 16825−16840. (31) Zhang, Y.; Castets, A.; Carlier, D.; Ménétrier, M.; Boucher, F. Simulation of NMR fermi contact shifts for lithium battery materials: The need for an efficient hybrid functional approach. J. Phys. Chem. C 2012, 116, 17393−17402. (32) Strobridge, F. C.; Middlemiss, D. S.; Pell, A. J.; Leskes, M.; Clément, R. J.; Pourpoint, F.; Lu, Z.; Hanna, J. V.; Pintacuda, G.; Emsley, L.; Samoson, A.; Grey, C. P. Characterising local environments in high energy density Li-ion battery cathodes: A combined NMR and first principles study of LiFexCo1−xPO4. J. Mater. Chem. A 2014, 2, 11948−11957. (33) Clément, R. J.; Pell, A. J.; Middlemiss, D. S.; Strobridge, F. C.; Miller, J. K.; Whittingham, M. S.; Emsley, L.; Grey, C. P.; Pintacuda, G. Spin-transfer pathways in paramagnetic lithium transition-metal phosphates from combined broadband isotropic solid-state MAS NMR spectroscopy and DFT calculations. J. Am. Chem. Soc. 2012, 134, 17178−17185. (34) Mali, G.; Sirisopanaporn, C.; Masquelier, C.; Hanzel, D.; Dominko, R. Li2FeSiO4 polymorphs probed by 6Li MAS NMR and 57 Fe Mössbauer spectroscopy. Chem. Mater. 2011, 23, 2735−2744. The square of the g-tensor in eq. 2 has been used in that work. While this conforms with earlier formalisms, it is not fully consistent with the formalisms used for the quantum-chemical calculations, neither in that work nor here. (35) Seymour, I. D.; Middlemiss, D. S.; Halat, D. M.; Trease, N. M.; Pell, A. J.; Grey, C. P. Characterizing oxygen local environments in paramagnetic battery materials via 17O NMR and DFT calculations. J. Am. Chem. Soc. 2016, 138, 9405−9408. (36) Lee, J.; Seymour, I. D.; Pell, A. J.; Dutton, S. E.; Grey, C. P. A systematic study of 25Mg NMR in paramagnetic transition metal oxides: Applications to Mg-ion battery materials. Phys. Chem. Chem. Phys. 2017, 19, 613−625. (37) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (38) Yazyev, O. V.; Tavernelli, I.; Helm, L.; Röthlisberger, U. Core spin-polarization correction in pseudopotential-based electronic structure calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 115110.

(39) Bahramy, M. S.; Sluiter, M. H. F.; Kawazoe, Y. Pseudopotential hyperfine calculations through perturbative core-level polarization. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 035124. (40) Marsman, M.; Kresse, G. Relaxed core projector-augmentedwave method. J. Chem. Phys. 2006, 125, 104101. (41) Blaha, P.; Schwarz, K.; Madsen, G. K. H.; Kvasnicka, D.; Luitz, J. WIEN2k, an augmented plane wave + local orbitals program for calculating crystal properties; Technische Universität Wien: Vienna, Austria, 2001. (42) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Crystal: A computational tool for the ab initio study of the electronic properties of crystals. Z. Kristallogr. - Cryst. Mater. 2005, 220, 571−573. (43) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL09 User’s Manual; University of Torino: Torino, Italy, 2009. (44) Pennanen, T. O.; Machacek, J.; Taubert, S.; Vaara, J.; Hnyk, D. Ferrocene-like iron bis(dicarbollide), [3-FeIII-(1,2-C2B9H11)2]−. The first experimental and theoretical refinement of a paramagnetic 11B NMR spectrum. Phys. Chem. Chem. Phys. 2010, 12, 7018−7025. (45) Kruck, M.; Sauer, D. C.; Enders, M.; Wadepohl, H.; Gade, L. H. Bis(2-pyridylimino)isoindolato iron(II) and cobalt(II) complexes: Structural chemistry and paramagnetic NMR spectroscopy. Dalton Trans. 2011, 40, 10406−10415. (46) Mareš, J.; Liimatainen, H.; Pennanen, T. O.; Vaara, J. Magnetic properties of Ni2+(aq) from first principles. J. Chem. Theory Comput. 2011, 7, 3248−3260. (47) Bühl, M.; van Mourik, T. NMR spectroscopy: Quantumchemical calculations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 634−647. (48) Pritchard, B.; Autschbach, J. Theoretical investigation of paramagnetic NMR shifts in transition metal acetylacetonato complexes: Analysis of signs, magnitudes, and the role of the covalency of ligand-metal bonding. Inorg. Chem. 2012, 51, 8340−8351. (49) Rouf, S. A.; Mareš, J.; Vaara, J. 1H chemical shifts in paramagnetic Co(II) pyrazolylborate complexes: A first-principles study. J. Chem. Theory Comput. 2015, 11, 1683−1691. (50) Gohr, S.; Hrobárik, P.; Repiský, M.; Komorovský, S.; Ruud, K.; Kaupp, M. Four-component relativistic density functional theory calculations of EPR g- and hyperfine-coupling tensors using hybrid functionals: Validation on transition-metal complexes with large tensor anisotropies and higher-order spin-orbit effects. J. Phys. Chem. A 2015, 119, 12892−12905. (51) Benda, L.; Mareš, J.; Ravera, E.; Parigi, G.; Luchinat, C.; Kaupp, M.; Vaara, J. Pseudo-contact NMR shifts over the paramagnetic metalloprotein CoMMP-12 from first principles. Angew. Chem., Int. Ed. 2016, 55, 14713−14717. (52) Ramsey, N. F. A molecular beam resonance method with separated oscillating fields. Phys. Rev. 1950, 78, 695−699. (53) Ramsey, N. F. Electron coupled interactions between nuclear spins in molecules. Phys. Rev. 1953, 91, 303−307. (54) McConnell, H. M.; Robertson, R. E. Isotropic nuclear resonance shifts. J. Chem. Phys. 1958, 29, 1361−1365. (55) Kurland, R. J.; McGarvey, B. R. Isotropic NMR shifts in transition metal complexes: The calculation of the Fermi contact and pseudocontact terms. J. Magn. Reson. 1970, 2, 286−301. (56) Pickard, C. J.; Mauri, F. All-electron magnetic response with pseudopotentials: NMR chemical shifts. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 245101. (57) Pickard, C. J.; Mauri, F. First-principles theory of the EPR g tensor in solids: Defects in quartz. Phys. Rev. Lett. 2002, 88, 086403. (58) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM 392


Article


(80) Schimmelpfenning, B. AMFI - an atomic mean field integral program; University of Stockholm: Stockholm, Sweden, 1996. (81) Schreckenbach, G.; Ziegler, T. Calculation of the g-tensor of electron paramagnetic resonance spectroscopy using gauge-including atomic orbitals and density functional theory. J. Phys. Chem. A 1997, 101, 3388−3399. (82) Slater, J. C. A simplification of the Hartree-Fock method. Phys. Rev. 1951, 81, 385−390. (83) Ditchfield, R. Self-consistent perturbation theory of diamagnetism. Mol. Phys. 1974, 27, 789−807. (84) Keith, T.; Bader, R. Calculation of magnetic response properties using atoms in molecules. Chem. Phys. Lett. 1992, 194, 1−8. (85) Keith, T. A.; Bader, R. F. Calculation of magnetic response properties using a continuous set of gauge transformations. Chem. Phys. Lett. 1993, 210, 223−231. (86) Cheeseman, J. R.; Trucks, G. W.; Keith, T. A.; Frisch, M. J. A comparison of models for calculating nuclear magnetic resonance shielding tensors. J. Chem. Phys. 1996, 104, 5497−5509. (87) Kacprzak, S.; Reviakine, R.; Kaupp, M. Understanding the EPR parameters of Glycine-derived radicals: The case of N-acetylglycyl in the N-acetylglycine single-crystal environment. J. Phys. Chem. B 2007, 111, 811−819. (88) Kaupp, M.; Gress, T.; Reviakine, R.; Malkina, O. L.; Malkin, V. G. g-Tensor and spin density of the modified Tyrosyl radical in Galactose oxidase: A density functional study. J. Phys. Chem. B 2003, 107, 331−337. (89) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (90) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (91) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (92) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (93) Van Yperen-De Deyne, A.; Pauwels, E.; Van Speybroeck, V.; Waroquier, M. Accurate spin-orbit and spin-other-orbit contributions to the g-tensor for transition metal containing systems. Phys. Chem. Chem. Phys. 2012, 14, 10690−10704. (94) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (95) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (96) Sheldrick, G. M. Crystal structure refinement with SHELXL. Acta Crystallogr., Sect. C: Struct. Chem. 2015, 71, 3−8. (97) Clough, S.; Poldy, F. Tunnelling rotation rate of a weakly hindered methyl group (E.P.R. and ENDOR measurements). J. Phys. C: Solid State Phys. 1973, 6, 1953. (98) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (99) Kutzelnigg, W.; Fleischer, U.; Schindler, M. NMR basic principles and progress; Diehl, P., Fluck, E., Günther, H., Kosfeld, R., Seelig, J., Eds.; Springer: Berlin, Germany, 1991. (100) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (101) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (102) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: A critical analysis. Can. J. Phys. 1980, 58, 1200−1211.

ESPRESSO: A modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (59) Pigliapochi, R.; Pell, A. J.; Seymour, I. D.; Grey, C. P.; Ceresoli, D.; Kaupp, M. DFT investigation of the effect of spin-orbit coupling on the NMR shifts in paramagnetic solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2017, 95, 054412. (60) Lippert, G.; Hutter, J.; Parrinello, M. The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations. Theor. Chem. Acc. 1999, 103, 124− 140. (61) Iannuzzi, M.; Chassaing, T.; Wallman, T.; Hutter, J. Ground and excited state density functional calculations with the Gaussian and augmented-plane-wave method. Chimia 2005, 59, 499−503. (62) CP2K Developer Group. CP2K home page. http://www.cp2k. org/ (accessed). (63) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (64) Del Ben, M.; Hutter, J.; VandeVondele, J. Second-order MøllerPlesset perturbation theory in the condensed phase: An efficient and massively parallel Gaussian and plane waves approach. J. Chem. Theory Comput. 2012, 8, 4177−4188. (65) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: Atomistic simulations of condensed matter systems. WIREs Comput. Mol. Sci. 2014, 4, 15−25. (66) Weber, V.; Iannuzzi, M.; Giani, S.; Hutter, J.; Declerck, R.; Waroquier, M. Magnetic linear response properties calculations with the Gaussian and augmented-plane-wave method. J. Chem. Phys. 2009, 131, 014106. (67) Pauwels, E.; Asher, J.; Kaupp, M.; Waroquier, M. Cluster or periodic, static or dynamic-the challenge of calculating the g-tensor of the solid-state glycine radical. Phys. Chem. Chem. Phys. 2011, 13, 18638−18646. (68) Blundell, S. Magnetism in condensed matter (Oxford master series in physics), 1st ed.; Oxford University Press: New York, 2001. (69) Chernova, N. A.; Nolis, G. M.; Omenya, F. O.; Zhou, H.; Li, Z.; Whittingham, M. S. What can we learn about battery materials from their magnetic properties? J. Mater. Chem. 2011, 21, 9865−9875. (70) Cahill, S. L.; Chapman, P. R.; Kirby, W. C.; Goward, R. G. The challenge of paramagnetism in two-dimensional 6,7Li exchange NMR. Appl. Magn. Reson. 2007, 32, 565−581. (71) Clément, R. J.; Middlemiss, D. S.; Seymour, I. D.; Ilott, A. J.; Grey, C. P. Insights into the Nature and Evolution upon Electrochemical Cycling of Planar Defects in the β-NaMnO2 Na-Ion Battery Cathode: An NMR and First-Principles Density Functional Theory Approach. Chem. Mater. 2016, 28, 8228−8239. (72) Mehring, M. Principles of high resolution NMR in solids; SpringerVerlag: Berlin, Germany, 1983. (73) Haeberlen, U. High resolution NMR in solids; Academic Press: New York, 1976. (74) Blöchl, P. E. First-principles calculations of defects in oxygendeficient silica exposed to hydrogen. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 62, 6158−6179. (75) Blügel, S.; Akai, H.; Zeller, R.; Dederichs, P. H. Hyperfine fields of 3d and 4d impurities in nickel. Phys. Rev. B: Condens. Matter Mater. Phys. 1987, 35, 3271−3283. (76) Harriman, J. E. Theoretical foundations of electron spin resonance; Academic Press: New York, 1978. (77) Declerck, R.; Pauwels, E.; Van Speybroeck, V.; Waroquier, M. First-principles calculations of hyperfine parameters with the Gaussian and augmented-plane-wave method: Application to radicals embedded in a crystalline environment. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 245103. (78) Neese, F. Efficient and accurate approximations to the molecular spin-orbit coupling operator and their use in molecular g-tensor calculations. J. Chem. Phys. 2005, 122, 034107. (79) Malkina, O. L.; Vaara, J.; Schimmelpfennig, B.; Munzarova, M.; Malkin, V. G.; Kaupp, M. Density functional calculations of electronic g-tensors using spin orbit pseudopotentials and mean-field all-electron spin orbit operators. J. Am. Chem. Soc. 2000, 122, 9206−9218. 393


Article


(121) Hepworth, M. A.; Jack, K. H. The crystal structure of manganese trifluoride, MnF3. Acta Crystallogr. 1957, 10, 345−351. (122) Hoppe, R.; Becker, S. Neues über TiF3 und Fluorotitanate(III): Cs2K[TiF6], Rb2K[TiF6] und TiF3 Mit einer Bemerkung über Rb3[TiF6]. Z. Anorg. Allg. Chem. 1989, 568, 126−135. (123) Smith, R. I.; Jørgensen, J.-E. On the compression mechanism of FeF3. Acta Crystallogr., Sect. B: Struct. Sci. 2006, 62, 987−992. (124) Onoda, M.; Hirose, H. Crystal structures and 3d electron configurations for the LixV2(PO4)3 insertion electrode system with a semi-double-electron reaction. J. Phys. Soc. Jpn. 2012, 81, 094801. (125) Yin, S.-C.; Grondey, H.; Strobel, P.; Anne, M.; Nazar, L. F. Electrochemical property: Structure relationships in monoclinic Li3−yV2(PO4)3. J. Am. Chem. Soc. 2003, 125, 10402−10411. (126) Krzystek, J.; Ozarowski, A.; Telser, J.; Crans, D. C. Highfrequency and -field electron paramagnetic resonance of vanadium(IV, III, and II) complexes. Coord. Chem. Rev. 2015, 301, 123−133. (127) Patoux, S.; Wurm, C.; Morcrette, M.; Rousse, G.; Masquelier, C. A comparative structural and electrochemical study of monoclinic Li3Fe2(PO4)3 and Li3V2(PO4)3. J. Power Sources 2003, 119−121, 278− 284. (128) Yin, S.-C.; Strobel, P. S.; Grondey, H.; Nazar, L. F. Li2.5 V2(PO4)3 : A room-temperature analogue to the fast-ion conducting high-temperature γ-phase of Li3V2(PO4)3. Chem. Mater. 2004, 16, 1456−1465.

(103) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Selfconsistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72, 650−654. (104) Blaudeau, J.-P.; McGrath, M. P.; Curtiss, L. A.; Radom, L. Extension of Gaussian-2 (G2) theory to molecules containing thirdrow atoms K and Ca. J. Chem. Phys. 1997, 107, 5016−5021. (105) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B., Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision E.01; Gaussian Inc.: Wallingford CT, 2009. (106) Neese, F. ORCA, An ab Initio, density functional and semiempirical program package, version 3.0.1; 2012. (107) Malkin, V. G.; Malkina, O. L.; Reviakine, R.; Arbuznikov, A. V.; Kaupp, M.; Schimmelpfennig, B.; Malkin, I.; Repiský, M.; Komorovský, S.; Hrobárik, P.; Malkin, E.; Helgaker, T.; Ruud, K. MAG-ReSpect, version 2.3, 2010. (108) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (109) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (110) Munzarová, M.; Kaupp, M. A critical validation of density functional and coupled-cluster approaches for the calculation of EPR hyperfine coupling constants in transition metal complexes. J. Phys. Chem. A 1999, 103, 9966−9983. (111) Stoll, H.; Paulus, B.; Fulde, P. On the accuracy of correlationenergy expansions in terms of local increments. J. Chem. Phys. 2005, 123, 144108. (112) TURBOMOLE GmbH. TURBOMOLE version 7.0; University of Karlsruhe and Forschungszentrum Karlsruhe GmbH: Karlsruhe, Germany, 2015; available from http://www.turbomole.com. (113) Bonnet, M.-L.; Iannuzzi, M.; Sebastiani, D.; Hutter, J. Local disorder in lithium imide from density functional simulation and NMR spectroscopy. J. Phys. Chem. C 2012, 116, 18577−18583. (114) Allolio, C.; Salas-Illanes, N.; Desmukh, Y. S.; Hansen, M. R.; Sebastiani, D. H-bonding competition and clustering in aqueous LiI. J. Phys. Chem. B 2013, 117, 9939−9946. (115) Wang, X.; Hou, Z.; Ikeda, T.; Terakura, K. NMR chemical shifts of 15N-bearing graphene. J. Phys. Chem. C 2014, 118, 13929− 13935. (116) Koseki, S.; Schmidt, M. W.; Gordon, M. S. Effective nuclear charges for the first- through third-row transition metal elements in spin-orbit calculations. J. Phys. Chem. A 1998, 102, 10430−10435. (117) Koseki, S.; Schmidt, M. W.; Gordon, M. S. MCSCF/631G(d,p) calculations of one-electron spin-orbit coupling constants in diatomic molecules. J. Phys. Chem. 1992, 96, 10768−10772. (118) Koseki, S.; Gordon, M. S.; Schmidt, M. W.; Matsunaga, N. Main group effective nuclear charges for spin-orbit calculations. J. Phys. Chem. 1995, 99, 12764−12772. (119) Kaupp, M.; Reviakine, R.; Malkina, O. L.; Arbuznikov, A.; Schimmelpfennig, B.; Malkin, V. G. Calculation of electronic g-tensors for transition metal complexes using hybrid density functionals and atomic meanfield spin-orbit operators. J. Comput. Chem. 2002, 23, 794−803. (120) Jack, K. H.; Gutmann, V. The crystal structure of vanadium trifluoride VF3. Acta Crystallogr. 1951, 4, 246−249. 394


Large-Scale Computation of Nuclear Magnetic Resonance Shifts for

Recommend Documents