DFT investigation of 3d transition metal NMR shielding tensors in diamagnetic systems using the gaugeincluding projector augmentedwave method
Lionel Truflandier,* Micha?l Paris and Florent Boucher Institut des Matériaux Jean Rouxel, UMR 6502, Université de Nantes  CNRS, 2, rue de la Houssinière BP 32229, 44322 Nantes Cedex, France
We present a density functional theory based method for calculating NMR shielding tensors for 3d transition metal nuclei using periodic boundary conditions. Calculations employ the gaugeincluding projector augmentedwave pseudopotentials method. The effects of ultrasoft pseudopotential and induced approximations on the secondorder magnetic response are intensively examined. The reliability and the strength of the approach for
49
Ti and
51
V nuclei is shown by comparison with traditional quantum
chemical methods, using benchmarks of finite organometallic systems. Application to infinite systems is validated through comparison to experimental data for the 51V nucleus in various vanadium oxide based compounds. The successful agreement obtained for isotropic chemical shifts contrasts with full estimation of the shielding tensor eigenvalues, revealing the limitation of pure exchangecorrelation functionals compared to their exactexchange corrected analogues.
*
Electronic addresses: Lionel.Truflandier@cnrsimn.fr or Florent.Boucher@cnrsimn.fr
1
I. INTRODUCTION Nuclear magnetic resonance (NMR) spectroscopy is a powerful technique to investigate the structures of molecules, solids or biomolecular systems. For extended systems, the interpretation of spectra provides useful information with regard to the chemical local environment, the number of sites, the coordination number, the internuclear distances or the degree of distortion of polyhedra. In some cases, high resolution NMR measurements can even be used to determine crystallographic space groups.[1] However assignment and interpretation of the resonance lines often remains delicate. This problem can be partially overcome by performing first principles calculations of NMR parameters, i.e. shielding tensors and, for nuclearspin larger than ?, electric field gradient (EFG) tensors. The development of theoretical methods to calculate NMR properties is currently underway in several scientific communities.[25] To perform tractable NMR calculations, one has to deal with the size of the systems under investigation and with the high dependence of the methods with respect to the various levels of approximation, which can significantly affect the computational resources needed. Furthermore, the timescale for NMR spectroscopy is slow compared with the rovibrational effects of a chemical system. Thus, in order to get quantitative agreement between experimental and calculated results, we have to look beyond static calculations and internal motion contributions to NMR parameters have to be evaluated. Excluding dynamic disorder, those effects can usually be neglected in solid state NMR due to the restricted atomic motion compared with liquid measurement.[1] The reader may find discussions about the state of the art in NMR calculations in several reviews. The review by Helgaker et al.,[4] for instance, gives a broad description of the various quantum chemical methods developed in computational
2
chemistry. The primary effects involved in NMR calculations are described in the de Dios and Facelli reviews.[6,7] Concerning EFG tensors, it is now well established that they can obtained, at a high level of precision, by performing accurate ground state density calculations. The EFG is directly related to the asphericity of the electron density in the vicinity of the nucleus probe. Various approaches can be used to obtain the full tensor components, the choice being specifically dependant on the type of system under study.[3,813] For shielding tensors, the problem is much more complicated. Until recently, the common calculation methods have been based on a molecular approach using localized atomic orbitals (LAOs), the cluster approach being used to mimic infinite periodic systems. However, two important problems remain. Firstly, investigations of molecular materials are carried out by isolating a molecule from the bulk. As a consequence, the chemical environment is neglected in the calculations even though intermolecular interactions may contribute to the shielding and quadrupolar parameters.[14,15] Secondly, in the case of a nonmolecular material, the most common compounds in solid state chemistry, strong difficulties of calculations and convergence problem usually occur when using a finite size model. [16] To overcome such difficulties, Pickard and Mauri have developed the socalled "gaugeincluding projector augmentedwave" (GIPAW) pseudopotential approach in which the periodicity of the system is explicitly taken into account using a planewave basis set to expand the wave functions.[5] This approach was proposed within the framework of density functional perturbation theory (DFPT). The advantage of the GIPAW approach over other pseudopotential methods [17,18] is the possibility of keeping the nodal
3
properties of the wave functions in the neighbourhood of the core in the presence of a magnetic field. Considering the rigid contribution of core electrons with respect to NMR parameters,[19] accuracy comparable to all electron calculations can be achieved.[5] Nevertheless the application to extended systems was, to date, limited to elements belonging to the first three rows of the periodic table,[2023] due to the difficulties involved in efficient pseudopotential development. Nowadays, NMR spectroscopy applied to transition metals is widely used in the fields of coordination, biochemistry and solid state materials. Among the 3d transition metals, numerous NMR measurements on 51V nuclei have been performed in order to probe the vanadium(+V) sites in homogeneous and heterogeneous catalysis,[24] battery materials or metalloproteins.[25,26] In this paper we will investigate the calculation of 49Ti and 51V NMR shielding tensors in organometallic and diamagnetic inorganic systems, using complexes of titanium and vanadium, and vanadium(+V)based compounds as representative cases. We will explore for the first time the accuracy of the pseudopotential GIPAW approach on 3d transition metal referring to all electron calculations obtained from traditional quantum chemical methods, the purpose being to apply the computational methodology on extended systems. In Sec. II, we will briefly explore the theoretical methods commonly used in computational chemistry, in order, first, to outline the context in which the GIPAW method was developed, and second to underline approximations and difficulties inherent in the use of a pseudopotential planewave method and its application to 3d elements. In Sec. III, we will present the sensitivity of the shielding tensor components accuracies with respect to the level of improvement of the pseudopotential generation. Afterwards, transferability will be checked by means of a
4
benchmark of titanium and vanadium complexes and validated by comparison to allelectron calculations. Application to 51V containing extended systems will be discussed in Sec. IV. A first example of such an application has been published recently on the AlVO4 system.[27] In this last part, we will finally concentrate on the relation between exchangecorrelation functional improvement and reliability of the results. II. THEORETICAL METHODS A. The electronic current density and the gauge problem The response of matter to a uniform external magnetic field B can be represented by an induced electronic current density j(r) which is associated with the operator J(r), through the following relation given in atomic units,
1 1 J (r ) = ? { p , r r } ? A(r ) r r . 2 c
(1)
Here { p, r r } denotes the anticommutator of the momentum p and projection r r operators: { p , r r } = p r r + r r p . A(r) is a vector potential connected to B through
B = ? × A(r ) or A(r ) = 1 B × (r ? r0 ) , where r0 is the gauge origin. The first and the 2
second parts of the right hand side of the Eq. (1) are the paramagnetic and the diamagnetic current operators, respectively. In a closed shell molecule or insulating nonmagnetic material and within the field strengths typically used in NMR experiments, the induced electronic current density is calculated through the firstorderinduced current j(1)(r). It yields a nonuniform induced magnetic field B (1) which shields each nucleus N in from B. The nuclear magnetic shielding tensor σ or the socalled chemical shift tensor defined as
5
(1) B in (rN ) = ? σ (rN )B =
rN ? r 1 3 (1) ∫d r j (r) × r ? r 3 , c N
(2)
is a secondorder magnetic response. The firstorder induced current density j(1) (r ) is obtained by means of perturbation theory applied to J (r' ) ,[28]
1 j(1) (r ) = ?∑ Ψ (0) { p , r r } Ψ (1) ? ρ 0 A(r ) o o c o
= j(p1) (r ) + j(d1) (r ) .
(3)
In this equation, the summation is over the occupied states o and ρ 0 is the unperturbed
(0) electron density. The ground state wave function Ψ o is the eigenvector of the field
(1) independant Hamiltonian H (0) associated with the eigenvalue ε o and Ψ o
is its
corresponding firstorder correction due to the magnetic field perturbation. j(d1) (r ) , which
depends only on the unperturbed charge density ρ 0 , is called the "diamagnetic" contribution and corresponds to the uniform circulation of the electrons. j(p1) (r ) , which depends on the firstorder perturbed wave function, is called the "paramagnetic" contribution to the total current and is assumed to be a correction due to the molecular environment. The chemical shift tensor being an observable quantity, j(1) (r ) must be independent of the choice of gauge origin r0. Both j(p1) (r ) and j(d1) (r ) are separately gauge dependent,
nevertheless only their sum must satisfy the gauge invariance property. The gauge dependence of j(d1) (r ) is explicit through the presence of A, while the gauge dependence of j(p1) (r ) is implicitly present in Ψ (1) . Different approaches can be used to evaluate Ψ (1) o o
6
such as the Sternheimer equation, the Green's function method, the sum over states approach or the Hylleraas variational principle.[29] All these methods use the first order perturbed Hamiltonian H (1) . For a perturbation due to an external magnetic field, H (1) = 1 L?B , 2c (4)
where L = r × p is the angularmomentum operator. Thus, the presence of H (1) in the calculation of Ψ (1) is responsible for the implicit gauge dependence of j(p1) (r ) . o Due to incomplete basis set, the gauge origin independence on j(1) (r ) is usually not completely verified, and it could in principle yield numerical divergence of the calculation of j(1) (r ) . Actually, the diamagnetic term converges faster than the paramagnetic part with respect to the basis size. In fact, the diamagnetic term converges quite easily, since only an accurate determination of the ground state density is needed. Considering the paramagnetic contribution, careful choice of gauge origin[30] can lead to a decrease in its magnitude over a particular region of space. As a consequence, a smaller error in the calculated value of j(1) (r ) is expected. The problem of different convergence rates is entirely solved when considering the simple case of an isolated closedshell atom:
j(p1) (r ) vanishes when the intuitive choice of gauge origin is taken at the nucleus.
Several methods have been developed to solve the gauge problem for molecular systems, using localized atomic orbitals (LAOs). In the limit of complete basis sets, without dependence on the magnetic field, the calculated magnetic shielding tensor should be gauge invariant.[31] Nevertheless, only small molecules have been studied in such a way because of the prohibitive computational effort required.[32,33] An alternative and practical method has been developed through the use of LAOs including explicit field 7
dependence. This well known approach called “gauge invariant atomic orbital” (GIAO) was introduced first by London and generalized for molecular systems by Ditchfield over 30 years ago.[34,35] Each oneelectron function has its own local gauge origin represented by a multiplicative complex factor. Latter, Keith and Bader have presented new methods based on the calculation of j(1) (r ) by performing a gauge transformation for each point of space.[36] The “continuous set of gauge transformations” method (CSGT) achieves gaugeinvariance via a parametric function d(r) which is defined in real space and shifts continuously the gauge origin. The potential vector is redefined as, 1 A(r ) = B × (r ? r0 ? d(r )) . 2 (5)
The type of CSGT method is determined by the choice of the d(r ) function.[2,19,36] If
d(r ) is a constant, the single gauge origin method is obtained. In their first work, Keith and Bader proposed a partition of the induced current density into contributions of atoms in a molecule.[30] This method called “individual gauges for atoms in molecules” (IGAIM) is based on the displacement of the gauge origin to the position of the nearest nucleus to the point r at which j(1) (r ) is calculated. In other words, the function d(r ) takes discrete values equal to the atomic center positions present in the molecule. For chemical shift calculations, CSGT and IGAIM methods give similar results.[36] GIAO, CSGT and IGAIM methods have been developed for molecular NMR calculations using localized basis sets. The difficulty associated with application of localized methods to extended systems was circumvented by the use of a cluster approximation.[3742] The accuracy of the results is closely related to the basis quality and the cluster size, and limited convergence was reached despite extensive computational effort. To overcome
8
the difficulties associated with solid state systems, an alternative approach was proposed using the fully periodic GIPAW method.[5] B. The gaugeincluding projector augmentedwave approach In order to discuss the approximations introduced in the magnetic field dependant GIPAW approach, we need first to briefly describe the projector augmentedwave (PAW) electronic structure calculation method elaborated by Bl?chl.[12] Within the frozen core approximation and the pseudopotential planewave formalism, the PAW method was developed by introducing an operator T that maps the true valence wave functions Ψ
~ ~ onto pseudowave functions Ψ , Ψ = Τ Ψ . The construction of T is carried out
through the use of allelectron (AE) and pseudo (PS) atomic wave functions (socalled AE and PS partial waves), respectively φi
~ and φi . As in other pseudopotential
methods, a cutoff radius rN,c (for each nucleus N) is used to define the augmentation region ΩN where the operator T must restore the complete nodal structure of the AE wave functions,
~ ? Τ = 1 + ∑ ? φ N,n ? φ N,n N,n ? Local projector functions ~N, n p
? ~ ? p N,n , ?
(6)
are introduced to expand the pseudowave function
locally onto the pseudoatomic orbitals. The index n refers to the angularmomentum quantum numbers and to an additional number, which is used if there is more than one projector per angularmomentum channel. Constraints[12] are imposed by the PAW
p method: ~ N,n
~ and φN ,n
have to be orthogonal inward Ω N and vanish beyond this ~ are identical to φN ,n . The evaluation of an observable
9
region ~ N,n , whereas φN ,n p
quantity represented by an operator O can be expressed in terms of pseudowave functions by ~ ~ Ψ T + OT Ψ = Ψ O Ψ , with an accuracy comparable to an AE
calculation. However, within the framework of practical PAW calculations, completeness conditions can not be achieved. The results are dependent on the PS wave function planewave basis set expansion and on the AE and PS atomic wave function number. The ability of the PAW method to reconstruct an AE wave function has allowed the use of the pseudopotential planewave formalism for calculations of hyperfine and EFG parameters.[3,43] The efficiency of the EFG calculations has been demonstrated for a large series of nuclei.[4447] Nevertheless, when considering a secondorder magnetic response as the shielding tensor, intricacies appears. It was demonstrated that the PAW approach does not preserve the translational invariance of eigenvectors in the presence of a uniform magnetic field.[5] The solution proposed by Pickard and Mauri, similar to the GIAO method, is to introduce a field dependant phase factor to the GIPAW method. Here, the multiplicative complex factor is carried out by the operator,
ΤB = 1 + ∑ e 2c
i
r?rN ×B
N,n
? ~ ? φ N,n ? φ N,n ? ?
? ~ ? i r?r ×B ? p N,n e 2 c N . ? ?
(7)
As a result, the GIPAW pseudoeigenvector Ψ eigenvector Ψ is defined by
i
associated with the allelectron
Ψ = TB Ψ . For a local or semilocal operator,
introducing ξ N = e 2 c
r ?rN × B
+ , the GIPAW pseudooperator O = TB OTB is given by
O =O+
N,n ,m
∑ξ
N
~ ~ ~ ? φ ξ+ O ξ φ p N,n ? N,n N N N,m ? φN,n ξ + O ξ N φN,m N ? ?
? ~ ? p N,m ξ + . N ? ?
(8)
10
If one applies the transformation given in Eq. (8) on the operator J(r) described in Eq. (1), the GIPAW current density operator becomes 1 1 ? ? J (r' ) = ? { p , r' r' } ? A(r' ) r' r' + ∑ ξ N (r )? ΔJ p (r' ) + ΔJ d (r' ) ? ξ + (r ) . N N N 2 c ? ? N (9)
The GIPAW nodal structure reconstruction leads to the introduction of the paramagnetic ΔJ p (r' ) and diamagnetic ΔJ d (r' ) operators defined in the augmentation region Ω N ; N N ΔJ p (r' ) = N 1 p ? ∑ ~N,n ? 2 n ,m ?
~ ~ ? φ N,n { p , r' r' } φ N,m ? φ N,n { p , r' r' } φ N,m ? ~ N,m , (10) p
? ? ~ ? p N,m . ?
ΔJ d (r' ) = N
B × (r' ?rN ) ~ ~ p ? ∑ ~N,n ? φ N,n r' r' φ N,m ? φ N,n r' r' φ N,n 2c ? n ,m
(11)
If one develops J in powers of B and uses densityfunctional perturbation theory,[29] the GIPAW firstorder current density is obtained and expressed in different contributions,[5]
j(1) (r' ) = j(1) (r' ) + j(1) (r' ) + j(1) (r' ) , bare Δp Δd
(12)
As in Eq. (4) the first order perturbed Hamiltonian is required and expressed thanks to an
+ expansion in powers of B of the GIPAW pseudoHamiltonian H = TB HTB . Obviously,
the expression for H depends entirely on the pseudopotential approach used: Either the normconserving[48,49] (NCPP) or the ultrasoft[50] (USPP) schemes. In this latter case, the relaxation of the normconstraint imposes an additional generalised orthonormality constraint which must be solved via an overlap operator S,
φN ,n S φN ,m = δ n,m .
~
~
(13)
Due to this additional degree of freedom, the simplifications (see Eqs. (11) and (12) with the following discussion of Ref. [5]) which are valid for a NCPP are no longer valid within the USPPGIPAW approach. The work of Yates has permitted development and 11
implementation of the USPPGIPAW formalism.[51] Due to the introduction of the
(1) generalized orthonormality constraint, the firstorder perturbed wave function Ψ n
given in Eq. (32) of Ref. [5] is redefined as,
(1) (0) Ψ n = G (ε n )( H (1) ? ε n S (1) ) Ψ n ,
(14)
with the Green function operator G(ε n ) expressed through
G (ε n ) = ∑
e ( ( Ψe 0 ) Ψe 0 )
εn ? εe
,
(15)
with the sum running over empty states e. H (1) and S (1) are respectively the firstorder perturbed GIPAW Hamiltonian and the firstorder perturbed overlap matrix. The Green function involving virtual subspace is only used here for convenience in order to express
(1) the firstorder perturbed wave function Ψ n of Eq. (14). Practically, [51] the closure
relation based on the summation of the occupied and virtual subspaces, coupled with a conjugategradient minimization scheme leads to a simple linear system of equations, involving solely the occupied ground state wave functions.[52,53] This advantageous scheme, which reduces considerably the computational time, succeeds to express the three different contributions of Eq. (12) as
(0) (0) j(1) (r' ) = 2∑ Re ( Ψ o { p , r' r' }G (ε o )( H (1) ? ε o S (1) ) Ψ o ) ? bare
o
(0) (0) ? ∑ Ψ o { p , r' r' } Ψ o'
oo'
1 ps ρ (r ' )B × r' 2c , (16) (0) (0) Ψ o' S (1) Ψ o
(0) (0) where o runs over the occupied states. ρ ps (r ' ) = 2∑ Ψo r' r' Ψo is the ground state
o
pseudodensity. The paramagnetic augmentation current is given by
12
? (0) (0) j(1) (r' ) = ∑ ? 4Re ( Ψ o ΔJ p (r' )G (ε o ) ( H (1) ? ε o S (1) ) Ψ o ) Δp N N,o ? , (17) 1 ? (0) p (0) (0) p (0) (0) (1) (0) [B × rN ? r' , ΔJ N (r' )] Ψo ? ? 2∑ Ψo ΔJ N (r' ) Ψo' Ψ o' S Ψo + 2 Ψo 2ic ? oo' and the diamagnetic augmentation current is
(0) (0) j(1) (r' ) = 2∑ Ψ o ΔJ d (r' ) Ψ o . Δd N N,o
(18)
The introduction of extra terms in the expression of j(1) (r ) , resulting from the additional bare orthonormality constraint, yields more awkward calculations compared to the normconserving GIPAW method. The NCPPGIPAW equations can be recovered by putting
S = 1 (Eqs. (36) and (37) of Ref. [5]). In order to increase tractability and accuracy of
calculations, the gauge origin in the GIPAW approach is put at the nucleus center setting r0 = rN.[19] By reformulating Eqs. (36) and (37) of Ref. [5], it has been shown that the firstorder induced current expressed in Eq. (12) is invariant upon a rigid translation through the individual invariance of its three contributions. Then, for a sufficient basis set expansion, the same rate of convergence is observed for j(1) (r' ) and j(1) (r' ) (the Δp bare convergence is governed by the first terms of the right hand sides of Eqs. (16) and (17)). Finally, in order to reduce the computational resources required for the chemical shielding tensor calculations, the firstorder induced magnetic field is divided into four contributions which can be individually calculated, taking advantage of the linearity of Eq. (2),
(1) B in (rN ) = B core (rN ) + B (1) (rN ) + B (1) (rN ) + B (1) (rN ) . bare Δp Δd
(19)
B core (rN ) , which depends only on the core electrons, i.e. of the isolated atom, is calculated once using the Lamb formula.[54]
13
At this stage, several approximations are introduced to compute NMR chemical shift
(1) tensors from the GIPAW approach. Firstly, to evaluate the correction to B in (rN ) due to
the paramagnetic and diamagnetic augmentation terms, only the augmentation region ΩN of the nucleus N is considered, i.e. the sum on N in Eqs. (16) and (17) is no longer carried out. This onsite approximation neglects the effects of the augmentation currents of the neighbouring atoms to the shielding of the studied atom. Secondly, within periodic
(1) conditions, B in
is formulated in reciprocal space using the BiotSavart law.
Unfortunately, for a null vector of the reciprocal lattice (G = 0), B (1) (G = 0) becomes a in macroscopic quantity.[17] The induced field depends on the surface currents and, as a result, on the shape of the sample. Therefore, the macroscopic magnetic susceptibility χ has to be evaluated and no full GIPAW approach is available at the moment. Thus, this quantity is calculated using only the j(1) (r' ) contribution. Finally, the pseudopotential bare used for GIPAW calculation must be chosen with caution. Earlier studies show good agreement between all electron (IGAIM) and pseudopotential GIPAW (NCPP) calculations.[5] For a noteworthy reduction of planewave expansion, USPPGIPAW calculations are able to reproduce NCPPGIPAW results.[51] Without neglecting the intrinsic pseudopotential generation parameters, and especially for 3d elements, the choice of the valence states as well as the number of projectors must be precisely examined in order to reach converged NMR shielding parameters. This issue will be investigated in the next section. C. Computational details: allelectron and USPPGIPAW calculations In this part we review the default computational parameters employed for the study. If different settings are used, then the calculation details will be explicitly given in the text. 14
In order to validate the shielding tensor GIPAW calculations for the titanium and vanadium atoms, the USPPGIPAW results have been compared to those obtained through the AE approach. The Gaussian 03 suite of programs[55] was used to compute allelectron magnetic response of molecules within the IGAIM approach,[30] combined with the "PerdewWang 91" exchange and correlation functional PW91.[56,57] Molecular geometries were optimized with symmetry constraints, using the B3LYP hybrid functional[58,59] with the 6311+G(2d,p) basis set.[6063] The default force tolerance parameter of 0.02 eV/? was kept. We considered different kinds of LAOs in order to check the basis set dependence on the shielding tensor calculations of vanadium and titanium atoms. The tripleζ 6311++G(3df,3pd) Pople's basis set developed by Watchers and Hay[6264] for the firstrow transition elements, the augmented tripleζ atomic natural orbital (ANO) of Roos and coworkers, tabulated from Sc to Cu atoms,[65] as well as Dunning's quintuple ζ correlationconsistent basis set (ccpCV5Z) developed for the Ti atom by Bauschlicher[66] were used. The basis sets for elements in the first three rows were adapted in order to be consistent with those used for 3d transition metals. For extended systems, all the calculations were carried out using the PW91 functional. The geometry optimization and GIPAW investigations were performed using the CASTEP and NMRCASTEP codes,[5,20,67] respectively. The Brillouin zone was sampled using MonkhorstPack technique.[68] Relaxation of ionic positions were performed at an energy cutoff of 600 eV, using a kpoint spacing always smaller than 0.05 ?1 and keeping experimental unit cells. The residual forces on atom positions were converged within 0.05 eV/?. Molecules were studied with 1 kpoint by the use of a supercell approach, checking that the supercell is large enough to avoid spurious
15
interaction between periodic images. This condition was in general satisfied in a 12000 Bohr3 (~12×12×12 ?3) simulation cell. Shielding tensor calculations for molecular and extended systems were carried out through the crystal approach.[5] The interaction of nuclei and core states with the valence electrons was taken into account by the use of USPPs.[50,69] The selection of core levels were the common ones: 1s, [He]2s2p and [Kr]4d for the elements of the second row, third row and for calcium, and lanthanum, respectively. Two projectors were introduced for each remaining ns and np valence states and for the specific case of the 1s valence state of hydrogen two projectors were also used. The core radii rc, beyond which the pseudowave functions match the allelectron ones, are given in parentheses (a.u.) for the various atoms: H(0.8), C(1.4), N(1.5), O(1.3), F(1.4), Mg(2.8), P(1.8), S(1.7), Cl(1.7), Ca(1.8) and La(2.3). rc was set to the same value for all angular momentum channels of a given atom. Moreover, nonlinear core corrections were employed,[70] with a cutoff radius equal to 0.7×rc. Finally, the same USPP settings were used for DFT geometry optimization as well as shielding parameters calculations, apart from the 3d elements, where the USPP settings for the GIPAW calculations are given explicitly in the text. D. Conventions The conventions used to calculate the chemical shift parameters {δiso, δaniso, ηδ}, from chemical shift tensor eigenvalues {δxx, δyy, δzz}, are defined as follows, isotropic component: anisotropy component: asymmetry component:
δ iso =
1 (δ xx + δ yy + δ zz ), 3
(20) (21) (22) 16
δ aniso = δ iso ? δ zz , ηδ = δ xx ? δ yy , δ aniso
with
δ zz ? δ iso ≥ δ xx ? δ iso ≥ δ yy ? δ iso ,
(23)
The shielding parameters {σiso, σaniso, ησ} are deduced from the calculated eigenvalues using relations similar to (18), (19) and (20). One obtains σiso = ? (σxx+σyy+σzz) and ησ =
ηδ while σaniso = δaniso according to the relation
δ ij = ?a ? [σ ij ? σ ref ]
(24)
where δij and σij are the chemical shift and absolute shielding tensor components respectively, a is a slope (equal to unity in experiments) and σref is the isotropic shielding of a reference compound. Unfortunately, firstprinciples calculations of σref involve the consideration of rovibrational and intermolecular effects. In order to circumvent such tricky calculations, σref was evaluated assuming a linear regression between computed
σiso and experimental δiso values.
III. GIPAW: APPLICATION TO 3d TRANSITION METALS A. Validation of the frozen core approximation Within the framework of the pseudopotential approximation, the GIPAW method is able to converge towards allelectron magnetic response calculations. One contributing factor of this success is the assumption of a rigid contribution to the shielding NMR parameters of core electrons, i.e. the validity of the frozen core approximation.[17,19,71] The main concept is that the core electrons are not involved in the chemical reactivity, i.e. the core wave functions of an atom remain unmodified whatever its chemical environment is. Therefore, the AE atomic potential can be replaced by a pseudopotential which mimics the potential created by the nucleus surrounded by its inner electrons. The orthogonality condition between the valence and the core states being relaxed, the valence wave functions become smoother and easier to calculate using planewave basis sets. For 17
second and third row elements, the corevalence states separation is quite obvious and usual selections of core states are employed by the community for first principle PP calculations. Difficulties appear for the fourth row elements, especially for the 3d transition metals.[72] Comparing atomic total energies, using the frozencore PAW and fully relaxed calculations, previous studies have demonstrated that favorable choice of corevalence separation, in terms of computational cost, leads to less accurate results.[73] In the case of the vanadium atom, inaccurate results were found when keeping the {1s2s2p3s3p} states as core states (in the following discussion core and valence shells will be distinguished by the use of braces and parentheses, respectively) while including the 3p states into the valence improved the precision. Consequently, for the firstorder magnetic response calculation applied to 3d transition metals through DFPT calculations, one must carefully check the gap between core and valence states. Within the frozen core approximation and GIAO approach, Schreckenbach and Ziegler have concluded that,[74] for the third period nuclei, the 2p state must be included explicitly in the valence to get accurate results. They also mentioned that for a 3d transition metal like 53Cr, the 3s and 3p valence shells are necessary. More recently, using the IGAIM approach and choosing the gauge origin at the nucleus center (see Eq. (5)), investigations for 29Si and 31P atoms have demonstrated that the core contribution to the chemical shielding is purely diamagnetic,[19] corresponding to a rigid participation of the {1s2s2p} core shells to the shielding tensor. Those contradictory conclusions led us to study the influence of the corevalence partition involved in the GIPAW chemical shielding tensors calculations for 3d elements. We present in Table I the shielding tensor
18
calculated for
51
V in the wellknown VOCl3 molecule, using different vanadium
pseudopotentials, going from a large {1s2s2p3s3p} to a small {1s2s} core. As previously suggested in the literature,[5] two projectors per channel were used for each angularmomentum, except in the case of the (3s3p4s3d)GIPAW calculation and for (2p3s3p4s3d)GIPAW where only one projector is used for the both the 2p and 3s channels. In all cases, the energy cutoff was set large enough to reach convergence for the calculated shielding values with respect to the basis size. A dramatic discrepancy, compared to the AE calculation, is observed in Table I when only the (4s3d) shells are used for the valence. The nonrigid core state contribution of the 3p level is obvious when one compares the (4s3d) and (3p4s3d)GIPAW calculations. Furthermore, considering the anisotropy parameter, a better agreement between GIPAW and allelectron IGAIM calculations is obtained for an extension of the valence states up to the 3s and even 2p atomic functions. Unfortunately, comparing the
51
V isotropic shielding convergence for
the (3s3p4s3d) and (2p3s3p4s3d)GIPAW calculations (Figure 1), with respect to the cutoff energy, application to solid state systems is not tractable when including the 2p functions in the valence states. B. Pseudopotential optimization and convergence In order to demonstrate the computational efficiency of the USPPGIPAW approach applied to 3d element shielding tensor calculations, we have plotted in Figure 1 the convergence evolution of a NCPP and USPP. For the NCPP case, the corevalence interaction was described by the TroullierMartin[49] scheme, in the KleinmanBylander[75] form. To be consistent with the previous calculation, we used the same corevalence separation and projector allocation as for the (3s3p4s3d)USPP. The cutoff
19
radii were obviously reduced to a reasonable value of 0.9 au. Moreover, to also demonstrate the interest of using optimized USPP for the vanadium atom,[76,77] we present the convergence results obtained for a nonoptimized (3s3p4s3d)USPP. The method for generating optimized pseudopotentials was introduced by Rappe, Rabe, Kaxiras and Joannopoulos (RRKJ).[76] The RRKJ scheme is based on the statement that, for isolated pseudoatoms, the total energy convergence is mainly dependant of its kinetic part, which governs the total energy of extended systems. Therefore, to achieve optimal convergence, the authors have proposed a direct method to minimize the high Fourier components of the pseudowave functions. Keeping the constraints of normalization and continuity of two derivatives at rc, the pseudization function is optimized in order to minimize the kinetic energy beyond the cutoff wave vector qc. For the nonoptimized USPP, using a default value of qc = 12.7 au?, the 51V isotropic shielding is converged to within 0.5 ppm at a cutoff energy of 750 eV (Fig. 1). The optimized USPP obtained by setting the qc parameter to 5.3 au?, allows reduction of the energy cutoff by about 200 eV. For the same level of accuracy, using NCPP, the cutoff must be dramatically augmented to 3000 eV, which forbids definitively its use for 3d metal shielding calculations involving the (3s3p4s3d) valence states. Finally, whatever the selected GIPAW corevalence separation or pseudopotential scheme are, one should carefully check the convergence using extended basis sets. The same remarks stands for the IGAIM method. C. Completeness of the basis set Within the framework of the PAW method, the completeness of the basis set depends on both the planewave energy cutoff and on the AE and PS partialwave function
20
expansions. With respect to the “additive augmentation principle”,[12] Bl?chl has shown that the truncation of the partialwave extension does not affect the completeness of the basis set, assuming the complementary participation of the planewave expansion. In order to have a tractable implementation of the PAW formalism for electronic structure calculations, this author has demonstrated that the use of a finite number of partial wave functions yields negligible discrepancy by comparison to AE calculations. To check the transferability of those properties beyond the GIPAW method and to compare shielding parameters with fully converged IGAIM values, we have investigated the convergence of the method with regard to the number of projectors used for each valence state. The validation of the shielding convergence with respect to the planewave energy cutoff is quite obvious and has been shown previously in Figure 1. If we rewrite Eq. (19) in terms of the isotropic shielding components, we find
G≠ G= σ iso (rN ) = σ core (rN ) + σ bare0 (rN ) + σ bare0 (rN ) + σ Δp (rN ) + σ Δd (rN ) .
(25)
Clearly, for an isolated molecular system such as the VOCl3 molecule, there are no
G= surface currents (see Sec. II.B.) and the σ bare0 component of Eq. (25) should tend to zero.
Thus, the value of this component is a useful tool to check the absence of interactions between periodic images of the molecular system, in the limit of very large supercells. In our calculations the value was always smaller than 0.5 ppm. Figure 2 shows the projector dependence of the various components of Eq. (25). For all the tested configurations, the planewave energy cutoff was set to 700 eV and we used a {1s2s2p}(3s3p4s3d) state configuration for the USPP. As expected, the sensitivity of the paramagnetic correction term is larger than the diamagnetic one, with respect to the number of projectors used. The augmentation of the 4s state with two projectors has no effect on the isotropic 21
shielding component. Indeed, since the paramagnetic augmentation current j(1) (r ' ) is Δp proportional to the angular momentum (Eqs. (4) and (17)), for a s angular momentum, only the bare j(1) (r ) (which contains a diamagnetic part) and the diamagnetic bare augmentation j(1) (r ' ) terms are dependent of the projector extension. Finally the Δd scattering property of the 4s state is well reproduced with at least one projector. On the other hand, augmentation of the 3p and 3d states leads to strong variations of the isotropic shielding components, especially for the paramagnetic augmentation term. While a deshielding effect is observed for a twoaugmented 3p state, a shielding effect is obtained for a twoaugmented 3d state. Therefore, this antagonistic effect must be countered by a balanced choice of the number of projectors allocated to the 3p and 3d states. Opposite variations are observed (Figure 2) for the bare term and the diamagnetic augmentation
G≠ correction expressed in Eq. (25). σ bare0 is slightly affected by the pseudopartial wave
expansion, which yields variations within 2 ppm, against 30 ppm for σ Δd . Furthermore, three projectors are needed to achieve convergence of the paramagnetic augmentation term with respect to the 3p and 3d states. Now, if we compare the fully converged IGAIM and GIPAW results (Table II) for the VOCl3 molecule, fairly good agreement is observed between both series of shielding parameters. In order to improve the reliability of the method for 3d transition metals, the shielding parameters of
49
Ti in the simple
TiCl3CH3 molecule are also discussed (Table III). The titanium USPP was built using the same corevalence separation and projector allocation as for the vanadium USPP 3sP3p2×P4s2×P3d2×P (see caption of Table I for details). The cutoff radius was set to 1.8 a.u. for all the angular momentum channels. Concerning isotropic shielding, AE calculations performed with ANO as well as correlationconsistent basis sets agree very well with the 22
GIPAW results, whereas a weak discrepancy of 2 % is observed for the anisotropy parameter. D. Pseudopotential transferability: application to organometallic systems After having demonstrated the accuracy of the USPPGIPAW method in the calculation of shielding parameters for two molecules, namely VOCl3 for the 51V and TiCl3CH3 for the 49Ti, it is important to test the transferability of our approach in various electronic and geometric environments. Thus, we have worked with benchmarks of eight V and six Ti based molecular diamagnetic systems. Several allelectron calculations of the
49 51
V and
Ti isotropic shielding values have been reported in the literature for organometallic
systems.[7884] Here, we have focused our investigations on the complexes presented in Tables II and III, which have been studied in recent works by Bühl et al..[83] Computation of the NMR shielding parameters within the GIPAW approach was investigated through the use of 3sP3p2×P4s2×P3d2×P and 3sP3p3×P4s2×P3d3×P ultrasoftpseudopotentials (see caption of Table I) which leads to different convergence levels. As pointed out in Tables II and III and keeping in mind the extended range of the absolute shielding components observed for 3d transition metals, excellent agreement is found between the GIPAW and IGAIM approaches, whatever the level of chosen accuracy. For vanadium isotropic values (Table II), the most important relative discrepancies are observed for the [V(CO)6] and VF5 complexes (6 % and 1 %, respectively, for the first level of accuracy), which may be attributed to the singular electronic environment of the vanadium nucleus. This statement is also true for TiCl4 (Table III), which exhibits a discrepancy of 2 % for the second level of convergence, whereas the isotropic value of [Ti(CO)6]2 compared to AE calculation remains inferior to 1 %.
23
A global analysis of our results is given in Table IV which also gather previously published calculations on
31
P,
29
Si and
13
C. [5] Regarding the mean absolute deviations
between GIPAW and AE, the differences for the anisotropy parameters are larger than for the isotropic shieldings. In the case of the 3sP3p2×P4s2×P3d2×P USPP, GIPAW and AE isotropic shielding values differ by only 6 ppm which is acceptable for the
51
V atom
compared to 1.5 ppm for 13C and 8.8 ppm for 31P. The average deviation decreases from 17 to 13 ppm for the anisotropy parameters when we used the 3sP3p3×P4s3×P3d3×P USPP, but unfortunately the value related to shielding parameters increases to 10 ppm. Eventually, if we now assess the percentage of deviation of the
51
V isotropic shielding
parameters with respect to the calculated value, a comforting mean value of 0.3 % is found (0.6 % for the second level of convergence), against 0.3 % for 29Si and 3.2 % for the 13C. The same conclusions can be drawn for the 49Ti results, and we remark that the average deviation of the anisotropy parameter is divided by a factor 4 compared to the vanadium value. In an NMR experiment, we are not directly interested in absolute shielding values but rather in chemical shift parameters with regard to a reference. If we now choose VOCl3 as the reference system, then, using Eq. (24) with a = 1, we can calculate GIPAW and IGAIM
51
V chemical shifts. From the values reported in Table III, we found a mean
relative discrepancy of 1.6 % and 1.3 % between GIPAW and IGAIM calculations for both the levels of convergence and only 0.8 % between the two GIPAW calculations. This last value drops to 0.2 % when excluding the singular [V(CO)6] and VF5 systems. As a result, the 3sP3p2×P4s2×P3d2×P USPP is sufficient to achieve accurate
51
V isotropic
chemical shift calculations with a reduced computational effort compared to a
24
3sP3p3×P4s2×P3d3×P USPP calculation. Furthermore, the calculation time using GIPAW method is of the order of IGAIM with the 6311++G(3df,3pd) basis set, while it is considerably smaller when more extended basis sets such as ccpCVXZ or ANO are used. Fast and stable convergence of GIPAW calculations could be a promising alternative compared to time consuming LAO methods in the case of 3d elements. This leads us to consider the planewave DFT method as an accurate and efficient approach for the calculation of NMR chemical shift in finite organometallic systems E. Relativistic effects A complete investigation of the relativistic effects on vanadium and titanium shielding tensor calculations is beyond the scope of this paper, but some comments have to be given in order to keep in mind the level of approximation used in the GIPAW method. It will also give some hints to clarify the origin of the differences found between the allelectron IGAIM and the USPPGIPAW methods. Calculation of the NMR shielding tensor can be separated into two steps: The self consistent field (SCF) procedure which at least leads to the unperturbed KohnSham (KS) eigenvalues and orbitals, and the linear response of these orbitals due to the presence of the magnetic field. Thus, two kinds of relativistic effects are distinguished when calculating the shielding parameters:[85] The indirect term which is associated to the energy and shape modifications of the unperturbed KS orbitals induced by a relativistic SCF procedure,[13] and the direct relativistic effects associated to the use of a relativistic fielddependant Hamiltonian which yields additional terms in the shielding tensor expressions.[8688] Moreover, these terms can be separated in scalar and spinorbit coupling parts, depending on the level of approximation used.[87,89] Obviously, for a nonconsistent use of methods, i.e. if two
25
different levels of relativistic approximations are used for the SCF and shielding calculations, the analysis and comparison of results should be undertaken with caution. In our investigations, allelectron calculations are performed with no relativistic approximation, whereas in GIPAW method, introduction of indirect relativistic effects is performed through the pseudopotential approximation. Indeed, the atomic
pseudopotentials and wave functions are generated by resolving the scalar relativistic KoellingHammond equation.[90] Bouten and coworkers have studied NMR shielding predictions of 3d metal oxide (MO4n with M = Cr, Mn, Fe) coupling zeroorder regular approximation (ZORA) and GIAO methods.[85,87] They have shown that indirect relativistic effects are from three to four times larger than the direct ones with, on isotropic shieldings, an average magnitude of 63 ppm and 17 ppm for the indirect and direct effects, respectively. However, the indirect contribution does not seem to be rigid with respect to the 3d metal and the considered electronic environment. Therefore, this incomplete insertion of indirect effects could explain the small discrepancies toward the USPPGIPAW and IGAIM results observed in Tables I and II. Previous studies combining the ZORA and GIPAW methods[88,9194] have shown the influence of scalar relativity on 77Se molecular systems. By taking into account both the direct and indirect effects, an average increase of 69 ppm of the selenium isotropic shielding is observed. However, when calculating a relative chemical shift and comparing to experiments, either using a reference system, or better, by applying a linear regression (Eq. 24), no difference is then found between these two calculations. Similar conclusions can be drawn for the
125
Te, where the relativistic effect is even larger and increases the
chemical shielding by about 255 ppm.
26
As a consequence further work on the influence of GIPAW indirect relativistic effects is necessary, in particular to define the magnitude of the indirect contributions on the shielding parameters, but we are confident that reasonably good results can be obtained for the chemical shift when using the current implementation of the USPP GIPAW method. Investigations on third row elements, especially for 49Ti and 51V are in progress. IV. APPLICATION TO EXTENDED SYSTEMS A. Results and discussion Having validated, on various molecular systems, the NMR shielding calculations for 49Ti and
51
V using the USPPGIPAW method, we will explore now for the first time the
accuracy of the pseudopotential approach in calculating the shielding parameters of 3d transition metals in extended systems. We will only focus here on the 51V nucleus, using a (3sP3p2×P4s2×P3d2×P) USPP for the vanadium (see Sec. III.D.) and an energy cutoff of 700 eV. NMR shielding tensors were calculated for thirteen inorganic vanadium systems, chosen to span a large range of chemical shift for the
51
V. Consequently, a total of
eighteen distinct vanadium sites have been investigated. The list of compounds is collected in Table V. Considering previous experimental studies,[95] five different types of vanadium species have been established: orthovanadate[9698] with almost regular tetrahedral units, pyrovanadate[99] with slightly distorted tetrahedra,
metavanadates[100,101] with distorted tetrahedra, vanadates[95,101103] with distorted octahedra and crystal embedded complexes containing distorted vanadium polyhedra with different surrounded atoms[104] (O and N for VO2[acpyinh]; O, N and S for VO(OEt)(ONS). A schematic representation of the different structural types and local vanadium environments is shown in Figure 3.
27
For the eighteen vanadium sites, the correlation between calculated isotropic shielding coefficients and experimental isotropic chemical shifts is shown in Figure 4 and evaluated by a linear leastsquares fit according to Eq. (24). This regression displays the good accuracy of the GIPAW method considering the value of the slope 1.047(41) (the ideal value being 1.0) and the correlation coefficient of 0.988. The root mean square deviation of 28 ppm is an indication of the attainable precision for a predictive calculation of isotropic chemical shift in inorganic vanadium based systems. It is also important to note that the fitted σref value (1939(59) ppm) is in perfect agreement with the isotropic shielding parameter obtained for VOCl3 using an allelectron calculation (Table I). From this linear regression, the theoretical chemical shift parameters have been calculated for the eighteen vanadium sites and compared to the experimental values (Table V). The larger discrepancies between experimental and theoretical isotropic components, observed for NH4VO3, βVOPO4 and VO(OEt)(ONS), can be explained by the strong distortion of the first coordination sphere for the vanadium atom. Moreover, for the special case of VO(OEt)(ONS), the metal atom is located in a quite unusual distorted square pyramid environment formed by one sulfur, one nitrogen and three oxygen atoms. When many inequivalent sites are present in the same structure, the primary interest is not to predict the isotropic chemical shifts, but instead to assign NMR resonances to the different environments of the probe nucleus. As emphasized in Figure 4, when we focus on a short range of chemical shift (between 1450 to 1350 ppm, for instance), the agreement between calculated and experimental values can be improved by a small adjustment of the σref value. This has been done in Table V for all the compounds having 28
more than one vanadium site. The results are given between parentheses and allow straightforward assignments of the
51
V resonances in the AlVO4,23 α and βMg2V2O7,
and Ca2V2O7 compounds. With a discrepancy of the order of a few ppm, we are able to discriminate inequivalent vanadium sites exhibiting close isotropic chemical shifts. Unfortunately, the previous conclusions are not transferable to anisotropy and asymmetry parameters. Despite the quite reasonable agreement between experimental and theoretical anisotropy parameters obtained for ortho and pyrovanadates, huge differences are observed for the other families of vanadiumbased compounds. Moreover, the asymmetry parameters are generally poorly reproduced (large experimental deviation could be observed in TABLE V). These disagreements suggest the existence of an indirect relation between the degree of distortion of polyhedra and the theoretical δaniso reliability. Finally, especially for high anisotropy values, a significant trend of underestimation of the calculated parameter is revealed. In order to check the overall quality of the correlation between experimental and calculated shielding parameters, and to understand the lack of reliability observed for the calculated δaniso and ηδ, the eigenvalues of the chemical shift tensor have to be considered.[27] Experimental eigenvalues have been obtained from chemical shift parameters using Eqs. (20) to (22), whereas theoretical values have been deduced from absolute shielding eigenvalues, using Eq. (24) and the linear regression previously fitted. We have shown that the classification of chemical shift eigenvalues according to the relation (22) can lead to inversions of calculated components with regard to the experimental values.[27] In order to have a consistent comparison, incorrect assignments have been corrected when needed. The correlation is plotted in Figure 5. When all the eighteen vanadium sites are considered, poor agreement is observed
29
between the experimental and theoretical eigenvalues, which contrast with the very good correlation observed for δiso values (Figure 4). This contrasting behavior may be related to an error compensation induced by the averaging process bound to the isotropic values. Looking more carefully at the different classes of compounds, orthovanadates, which reveals low δaniso and high ηδ values, are characterized by a low dispersion of the eigenvalues and exhibit certainly the best agreement (see discussion Ref. [27]). This is in contrast to the pyrovanadates, where experimental and theoretical δaniso display quite good agreement, yet strong disparities are graphically observed due to the poor reproduction of the ηδ values (Table V). For the other families, the correlation in Figure 4 is even worse, in relation with the increasing polyhedron distortion. B. Improvement of DFT calculations The previous results demonstrate the difficulties to quantitatively reproduce the chemical shift components using DFT. Observed discrepancies between experiment and theory are masked by the average isotropic chemical shift and magnified by the anisotropy and asymmetry parameters. Moreover, difficulties in calculating, in specific cases, isotropic chemical shift for
17
O have been recently reported and discussed for solid state NMR.
[105] The authors have invoked the "band gap error" and the inaccuracy of the local density approximation (LDA) and generalized gradient approximation (GGA) exchangecorrelation (XC) functionals to properly describe excited state spectra. Other investigations on molecular systems have shown that calculated shielding parameters are highly dependant on the type of XC functionals.[106110] Linear response of crystalline or molecular orbitals to the magnetic field perturbation are strongly dependant on the occupiedvirtual energy gap (Egap) and the shape of the virtual orbitals, through the first30
order corrected wave function (Eq. 14). Recent studies have shown that hybrid density functionals, which include a portion of HartreeFock (HF) exchange, partially overcome the "band gap error" problem in solid state systems.[111115] In the case of quantum chemical NMR calculations, it was established that implementation of exact exchange in functionals leads to a huge improvement of calculated transition metal isotropic shieldings in organometallic systems.[116] To our knowledge, apart from an isolated computational investigation of the effect of the XC functionals on anisotropy for nuclei in organic molecules,[108] theoretical investigations have mainly been carried out considering the average isotropic component obtained from the three eigenvalues of the second rank shielding tensor. To discuss the influence of the HF exchange on anisotropy and asymmetry parameters, we now focus our attention on the VOCl3 inorganic system (bulkoptimized geometry have been kept, see Sec. II.C.). Shielding calculations were performed through the use of IGAIM method coupled with the 6311++G(3df,3pd) basis set. Investigation of the influence of the exact exchange on shielding parameters has been performed using different exchangecorrelation functionals. For GGAs, we have used the "PerdewWang 91" exchange and correlation functional PW91,[56,57] and the BLYP functional, which combined the "Becke’s 1988" exchange and the "LeeYangParr" correlation functionals.[58,117] Hybrid XC functionals are defined by the following exchangecorrelation approximation, E hybrid = αE HF + (1 ? α )E LDA + βΔE GGA + E GGA . XC X X X C (26)
Where E HF is the "exact" HF exchange, E LDA is the LDA exchange, ΔE GGA and E GGA are X X X C respectively the exchange correction and correlation parts of GGA functional. We use the 31
threeparameter B3 exchange functional defined by Becke,[59] leading to a value of α = 0.2. The correlation GGA functionals E GGA are taken as the PerdewWang 91,[56,57] and C LeeYangParr.[58] Results are collected in Table VI. Firstly, in order to probe the packing effect on the 51V shielding parameters, we have used the cluster approximation using ten additional VOCl3 entities which mimic the bulk environment on a central molecule (Table VI). This procedure works pretty well in the present case if we compare the GIPAW calculations and the IGAIMcluster results, and validates both approaches. By isolating a unique VOCl3 molecule and comparing to the cluster results, we conclude that the influence of the Van der Waals interactions on calculated shielding parameters are negligible. Thus, calculations carried out with an isolated molecule should be reliable enough to be extrapolated to the fully periodic GIPAW calculations. Inspection of Table VI reveals that the two GGAs as well as the two hybrid functionals give similar results. The differences between both sets of pure and hybrid functionals are around 35 ppm for
σiso and 20 ppm for δaniso. Considering a GGA and the corresponding hybrid functional,
we observe a fairly good improvement of δaniso with regard to experiment (Table V) when exact exchange is introduced. Afterwards, we have studied the dependence of the calculated shielding eigenvalues on the amount of exact exchange involved in the hybrid functional. This has been done using the halfandhalf functional proposed by Becke,[118] and defined with the following relation,
LYP E HandH = αE HF + (1 ? α )E LDA + E C . XC X X
(27)
Evolution of the occupiedvirtual gap and shielding eigenvalues with regard to the mixing coefficient α are displayed in Figure 6. Increase of the exact exchange leads to a linear widening of the occupiedvirtual energy splitting. Egap discrepancy between pure DFT 32
exchange (α = 0, called HandH0), and quasifull HF exchange (α = 0.9, called HandH0.9) is about 0.30 a.u.. Calculation using HartreeFock level of theory (results not shown) gives a value of 0.50 a.u. compared to 0.19 and 0.41 for HandH0 and HandH0.9. These results agree with the wellknown LDAGGA underestimation and HartreeFock overestimation of occupiedvirtual energy gap. Considering the shielding components results, we observed that σiso and δaniso are strongly dependent on the exact exchange, and the anisotropy parameter is the more affected. Following the above observations, we could suspect that the anisotropy improvement is closely bound to the correction of the occupiedvirtual energy gap induced by the use of hybrid XC functionals. Nevertheless, according to an extensive study of the influence of pure exchange on 57Fe isotropic shielding through GIAODFT calculations,[107] Schreckenbach has
demonstrated that three factors are responsible for the improvement induced by the use of hybrid functionals: enhancement of the occupiedvirtual gap, increase of the diffuse character of virtual molecular orbitals and the coupling contribution due to the HF exchange (Eq. (21) from Ref. [107]). All these contributions, and especially the last two, have an important effect on the paramagnetic part of the shielding tensor. As a result further work is in progress to understand quantitatively the influences of the exact exchange on the shielding tensor eigenvalues. At least we can deduce that the discrepancies found for the 51V anisotropy and asymmetry NMR parameters are probably linked to a fundamental DFT deficiency rather than GIPAW builtin approximations. V. CONCLUSION We have shown that extension of the GIPAW method to 3d nuclei in finite and infinite systems is reliable and reproduces with high accuracy the NMR isotropic shieldings of
33
51
V and
49
Ti in diamagnetic molecularlike and extended inorganic systems. The stable
and fast convergence of the pseudopotential method is able to overcome difficulties due to the incomplete expansion of the localized basis, reducing considerably the computational cost associated with traditional quantum chemical methods. Moreover the use of scalar relativistic pseudopotentials leads to the introduction of indirect relativistic corrections without increasing calculation time, which are the dominant contribution in 3d transition metals compared to fully relativistic calculations. Furthermore, direct assignment of
51
V solid state NMR resonances is allowed. We have demonstrated that
principal components of the shielding tensors should be considered in order to avoid erroneous conclusions on the quality of the theoretical model, when looking for correlation between calculated and experimental results. Despite a lack of reliability observed for anisotropy and asymmetry parameters, we are hopeful that future investigations will correct these limitations of DFT. Finally, we believe that this new approach will be a complementary and useful tool for experimental NMR research applied to organometallic and solid state chemistry. ACKNOWLEDGEMENTS The calculations presented in this work have been carried out at the Centre Régional de Calcul Intensif des Pays de la Loire financed by the French Research Ministry, the Région Pays de la Loire, and Nantes University. L.T. gratefully acknowledges C.J. Pickard for useful discussions and J.R. Yates for providing his PhD thesis manuscript. We also wish to thank C. Payen, N. Dupré and C. Ewels for careful reading of the manuscript.
34
References [1]R. K. Harris, Solid State Sci. 6, 1025 (2004). [2]J. R. Cheeseman, G. W. Trucks, T. A. Keith, and M. J. Frisch, J. Chem. Phys. 104, 5497 (1996). [3]H. M. Petrilli, P. E. Bl?chl, P. Blaha, and K. Schwarz, Phys. Rev. B 57, 14690 (1998). [4]T. Helgaker, M. Jaszunski, and K. Ruud, Chem. Rev. 99, 293 (1999). [5]C. J. Pickard and F. Mauri, Phys. Rev. B 63, 245101 (2001). [6]A. C. d. Dios, Prog. Nucl. Magn. Reson. Spectrosc. 29, 229 (1996). [7]J. C. Facelli, Concepts in Magn. Reson. Part A 20A, 42 (2003). [8]E. N. Kaufmann and R. J. Vianden, Rev. Mod. Phys. 51, 161 (1979). [9]E. Wimmer, H. Krakauer, M. Weinert, and A. J. Freeman, Phys. Rev. B 24, 864 (1981). [10]P. Blaha, K. Schwarz, and P. Herzig, Phys. Rev. Lett. 54, 1192 (1985). [11]P. Blaha, K. Schwarz, and P. H. Dederichs, Phys. Rev. B 37, 2792 (1988). [12]P. E. Bl?chl, Phys. Rev. B 50, 17953 (1994). [13]P. Pyykk?, Chem. Rev. 88, 563 (1988). [14]M. Strohmeier, D. Stueber, and D. M. Grant, J. Phys. Chem. A 107, 7629 (2003). [15]J. Czernek, J. Phys. Chem. A 107, 3952 (2003). [16]J. A. Tossell, in Computational Materials Chemistry, edited by L. A. Curtis and M. S. Gordon (Kluwer, Dordrecht, 2004). [17]F. Mauri, B. G. Pfrommer, and S. G. Louie, Phys. Rev. Lett. 77, 5300 (1996). [18]D. Sebastiani and M. Parrinello, J. Phys. Chem. A 105, 1951 (2001). [19]T. Gregor, F. Mauri, and R. Car, J. Chem. Phys. 111, 1815 (1999).
35
[20]M. Profeta, C. J. Pickard, and F. Mauri, J. Am. Chem. Soc. 125, 541 (2003). [21]T. Charpentier, S. Ispas, M. Profeta, F. Mauri, and C. J. Pickard, J. Phys. Chem. B 108, 4147 (2004). [22]C. Gervais, M. Profeta, F. Babonneau, C. J. Pickard, and F. Mauri, J. Phys. Chem. B 108, 13249 (2004). [23]J. R. Yates, S. E. Dobbins, C. J. Pickard, F. Mauri, P. Y. Ghi, and R. K. Harris, Phys. Chem. Chem. Phys. 7, 1402 (2005). [24]O. W. Howarth, Prog. Nucl. Magn. Reson. Spectrosc. 22, 453 (2001). [25]C. P. Gray and N. Dupré, Chem. Rev. 104, 4493 (2004). [26]D. C. Crans, J. J. Smee, E. Gaidamauskas, and L. Yang, Chem. Rev. 104, 849 (2004). [27]L. Truflandier, M. Paris, C. Payen, and F. Boucher, J. Phys. Chem. B 110, 21403 (2006). [28]C. J. Jameson and A. D. Buckingham, J. Chem. Phys. 73, 5684 (1980). [29]X. Gonze, Phys. Rev. A 52, 1096 (1995). [30]T. A. Keith and R. F. W. Bader, Chem. Phys. Lett. 194, 1 (1992). [31]W. N. Lipscomb, Advan. Magn. Res. 2, 137 (1966). [32]P. Lazzeretti and R. Zanasi, J. Chem. Phys. 68, 832 (1978). [33]F. Keith and R. Ahlrichs, J. Chem. Phys. 71, 2671 (1979). [34]F. London, J. Phys. Radium, Paris 8, 397 (1937). [35]R. Ditchfield, Mol. Phys. 27, 789 (1974). [36]T. A. Keith and R. F. W. Bader, Chem. Phys. Lett. 210, 223 (1993). [37]X. Xue and M. Kanzaki, J. Phys. Chem. B 103, 10816 (1999).
36
[38]G. Valerio, A. Goursot, R. Vitrivel, V. Malkina, and D. R. Salahub, J. Am. Chem. Soc. 120, 11426 (1998). [39]G. Valerio and A. Goursot, J. Phys. Chem. B 103, 51 (1999). [40]J. A. Tossell, Chem. Phys. Lett. 303, 435 (1999). [41]J. A. Tossell and R. E. Cohen, J. NonCryst. Solids 120, 13 (2001). [42]J. A. Tossell and H. J., J. Phys. Chem. B 109, 1794 (2005). [43]C. G. Van de Walle and P. E. Bl?chl, Phys. Rev. B 47, 4244 (1993). [44]K. Schwarz, C. AmbroschDraxl, and P. Blaha, Phys. Rev. B 42, 2051 (1990). [45]P. Dufek, P. Blaha, and K. Schwarz, Phys. Rev. Lett. 75, 3545 (1995). [46]T. J. Bastow, M. I. Burgar, and C. Maunders, Solid State Commun. 122, 629 (2002). [47]T. Bredow, P. Heitjans, and M. Wilkening, Phys. Rev. B 70, 115111 (2004). [48]L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982). [49]N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). [50]D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). [51]J. R. Yates, PhD Thesis, Cambridge University, 2003. [52]S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987). [53]S. Baroni, S. d. Gironcoli, A. D. Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001). [54]W. E. Lamb, Phys. Rev. 60, 817 (1941). [55]Gaussian 03, Revision C.02, Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.;
37
Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; 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.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; AlLaham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; and Pople, J. A.; Gaussian, Inc., Wallingford CT, 2004. [56]J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992). [57]J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). [58]C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988). [59]A. D. Becke, J. Chem. Phys. 98, 5648 (1993). [60]A. D. McLean and G. S. Chandler, J. Chem. Phys. 72, 5639 (1980). [61]R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. 72, 650 (1980). [62]A. J. H. Wachters, J. Chem. Phys. 52, 1033 (1970). [63]K. Raghavachari and G. W. Trucks, J. Chem. Phys. 91, 1062 (1989). [64]P. J. Hay, J. Chem. Phys. 66, 4377 (1977). [65]R. PouAmerigo, M. Merchan, I. NebotGil, P. O. Widmark, and B. Roos, Theor. Chim. Acta 92, 149 (1995).
38
[66]C. W. Bauschlicher, Theor. Chim. Acta 10, 141 (1999). [67]M. D. Segall, P. J. D. Lindan, M. J. Probert, C. J. Pickard, P. J. Hasnip, S. J. Clark, and M. C. Payne, Phys.: Cond. Matt. 14, 2717 (2002). [68]H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). [69]K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt, Phys. Rev. B 47, 10142 (1993). [70]S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982). [71]J. Richard, B. Levy, and P. Millie, Mol. Phys. 36, 1025 (1978). [72]G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). [73]N. A. W. Holzwarth, G. E. Matthews, R. B. Dunning, A. R. Tackett, and Y. Zeng, Phys. Rev. B 55, 2005 (1997). [74]G. Schreckenbach and T. Ziegler, Int. J. Quantum Chem. 60, 753 (1996). [75]L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982). [76]A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990). [77]J. S. Lin, A. Qteish, M. C. Payne, and V. Heine, Phys. Rev. B 47, 4174 (1993). [78]M. Bühl and F. A. Hamprecht, J. Comput. Chem. 19, 113 (1998). [79]M. Bühl, M. Kaupp, O. L. Malkina, and V. G. Malkin, J. Comput. Chem. 20, 91 (1998). [80]M. Bühl and M. Parrinello, Chem. Eur. J. 7, 4487 (2001). [81]M. Bühl, P. Imhof, and M. Repisky, Chem. Phys. Chem. 5, 410 (2004). [82]M. Bühl, R. Schurhammer, and P. Imhof, J. Am. Chem. Soc. 126, 3310 (2004). [83]S. Grigoleit and M. Bühl, Chem. Eur. J. 10, 5541 (2004).
39
[84]W. von Philipsborn, Chem. Soc. Rev. 28, 95 (1999). [85]R. Bouten, E. J. Baerends, E. v. Lenthe, L. Visscher, G. Schreckenbach, and T. Ziegler, J. Phys. Chem. A 104, 5600 (2000). [86]G. Schreckenbach and T. Ziegler, Int. J. Quantum Chem. 61, 899 (1997). [87]S. K. Wolff, T. Ziegler, E. van Lenthe, and E. J. Baerends, J. Chem. Phys. 110, 7689 (1999). [88]J. R. Yates, C. J. Pickard, M. C. Payne, and F. Mauri, J. Chem. Phys. 113, 5746 (2003). [89]S. K. Wolff and T. Ziegler, J. Chem. Phys. 109, 895 (1998). [90]D. D. Koelling and B. N. Harmon, J. Phys. C: Solid State Phys. 10, 3107 (1977). [91]E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 99, 4597 (1993). [92]E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 101, 9783 (1994). [93]R. van Leeuwen, E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 101, 1272 (1994). [94]A. J. Sadlej, J. G. Snijders, E. van Lenthe, and E. J. Baerends, J. Chem. Phys. 102, 1758 (1995). [95]O. B. Lapina, V. M. Masthikhin, A. A. Shubin, V. N. Krasilnikov, and K. I. Zamaraev, Prog. Nucl. Magn. Reson. Spectrosc. 24, 457 (1992). [96]J. Skibsted, C. J. H. Jacobsen, and H. J. Jakobsen, Inorg. Chem. 37, 3083 (1998). [97]U. G. Nielsen, A. Boisen, M. Brorson, C. J. H. Jacobsen, H. J. Jakobsen, and J. Skibsted, Inorg. Chem. 41, 6432 (2002). [98]U. G. Nielsen, H. J. Jakobsen, and J. Skibsted, Solid State Nucl. Magn. Reson. 23, 107 (2003).
40
[99]U. G. Nielsen, H. J. Jakobsen, and J. Skibsted, J. Phys. Chem. B 105, 420 (2001). [100]J. Skibsted, N. C. Nielsen, H. Bilds?e, and H. J. Jakobsen, J. Am. Chem. Soc. 115, 7351 (1993). [101]U. G. Nielsen, H. J. Jakobsen, and J. Skibsted, Inorg. Chem. 39, 2135 (2000). [102]J. Skibsted, N. C. Nielsen, H. Bilds?e, and H. J. Jakobsen, Chem. Phys. Lett. 188, 405 (1992). [103]O. B. Lapina, D. F. Khabibulin, A. A. Shubin, and V. M. Bondareva, J. Mol. Catal. A 162, 381 (2000). [104]N. Pooransingh, E. Pomerantseva, M. Ebel, S. Jantzen, D. Rehder, and T. Polenova, Inorg. Chem. 42, 1256 (2003). [105]M. Profeta, M. Benoit, F. Mauri, and C. J. Pickard, J. Am. Chem. Soc. 126, 12628 (2004). [106]Y. RuizMorales, G. Schreckenbach, and T. Ziegler, J. Phys. Chem. A 101, 4121 (1997). [107]G. Schreckenbach, J. Chem. Phys. 110, 11936 (1999). [108]P. J. Wilson, R. D. Amos, and N. C. Handy, Chem. Phys. Lett. 312, 475 (1999). [109]V. G. Malkin, O. L. Malkina, M. E. Casida, and D. R. Salahub, J. Am. Chem. Soc. 116, 5898 (1994). [110]T. Helgaker, P. J. Wilson, R. D. Amos, and N. C. Handy, J. Chem. Phys. 113, 2983 (2000). [111]R. L. Martin and F. Illas, Phys. Rev. Lett. 79, 1539 (1997). [112]T. Bredow and A. R. Gerson, Phys. Rev. B 61, 5194 (2000).
41
[113]J. K. Perry, J. TahirKheli, and W. A. Goddard, III, Phys. Rev. B 63, 144510 (2001). [114]J. Muscat, A. Wander, and N. M. Harrison, Chem. Phys. Lett. 342, 397 (2001). [115]J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, J. Chem. Phys. 123, 174101 (2005). [116]M. Bühl, Chem. Phys. Lett. 267, 251 (1997). [117]A. D. Becke, Phys. Rev. A 38, 3098 (1988). [118]A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
42
TABLE I. Convergence of the 51V absolute isotropic and anisotropic shielding parameters as a function of the vanadium valence states involved in USPPGIPAW calculations for the VOCl3 molecule. The multiprojector USPP is defined by the notation nlk×P where an integer k is associated to each nl atomic state and displays the number of projectors allocated (one projector is allocated to the 3s channel). Valence State (4s3d) (3p4s3d) (3s3p4s3d) rcc 2.4 2.5 2.0 vlocc p(0.5) f(0.0) f(0.0) f(0.0) Number of projectors 4s2×P3d2×P 3p2×P4s2×P3d2×P 3sP3p2×P4s2×P3d2×P 2pP3sP3p2×P4s2×P3d2×P 
σiso (ppm) δσ (ppm)
1806 1910 1910 1920 1904 353 434 455 461 483
(2p3s3p4s3d)a 0.8/2.0 allelectronb
a

A core radius of 0.8 and 2.0 a.u. was used for the 2p and for the remaining states respectively.
b c
IGAIM/6311++G(3df,3pd).
Core radius rc and atomic reference energies (in parentheses) of the local atomic pseudopotential vloc are given in a.u..
43
TABLE II.
51
V NMR shielding parameters in various molecular systems. The GIPAW
calculations were performed using the (1) 3sP3p2×P4s2×P3d2×P and (2) 3sP3p3×P4s2×P3d3×P ultrasoft pseudopotentials and compared to the IGAIM calculations performed with the (1) 6311++G(3df,3pd) and (2) ANO3ζ? LAO basis sets.
51
V USPP/LAO (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) (1) (2)
σiso (ppm)
GIPAW 1910 1947 97 89 1220 1280 1177 1212 1415 1451 1546 1584 1663 1700 3034 3074 IGAIM 1904 1952 91 76 1233 1258 1177 1214 1418 1458 1548 _ 1663 1707 3020 3057
δσ (ppm)
GIPAW 455 463 0 0 9 1 336 335 293 293 42 46 345 349 1641 1647 IGAIM 483 464 0 0 6 11 317 348 290 297 61 _ 358 347 1652 1615 GIPAW 0.00 0.00 n/a n/a 0.01 0.14 0.00 0.00 0.37 0.36 0.02 0.02 0.46 0.46 0.00 0.00
ηδ
IGAIM 0.00 0.00 n/a n/a 0.00 0.00 0.00 0.00 0.27 0.35 0.01 _ 0.47 0.45 0.00 0.00
Molecule VOCl3
[V(CO)6]VF5 VOF3 VOClF2 VONa VOCl2F VO(CH3)3
a
Abbreviation for the VO(OCH2CH2)3N complex. Computation of the NMR shielding parameters was not tractable for this molecule using the ANO3ζ basis set.
44
TABLE III. 49Ti NMR shielding parameters in various molecular systems. The GIPAW calculations were performed using the (1) 3sP3p2×P4s2×P3d2×P and (2,3) 3sP3p3×P4s2×P3d3×P ultrasoft pseudopotentials and compared to the IGAIM calculations performed with the (1) 6311++G(3df,3pd), (2) ANO3ζ? and (3) ccpCV5Z LAO basis sets.
49
Ti USPP/LAO (1) (2) (3)
σiso (ppm)
GIPAW 1459 1491 623 626 2171 2206 1845 1879 2434 2473 778 796 IGAIM 1471 1494 1489 621 622 622 2183 2191 2208 1859 1885 1876 2448 2468 2451 780 780 781
δσ (ppm)
GIPAW 465 471 0 0 445 451 483 489 0 0 0 0 IGAIM 462 479 476 0 0 0 441 448 451 477 491 487 0 0 0 0 0 0 GIPAW 0.00 0.00 n/a n/a 0.00 0.00 0.78 0.78 n/a n/a n/a n/a
ηδ
IGAIM 0.00 0.00 0.00 n/a n/a n/a 0.00 0.00 0.80 0.80 0.80 n/a n/a n/a n/a n/a n/a
Molecule TiCl3CH3
[Ti(CO)6]2
(1) (2) (3)
TiCl(CH3)3
(1) (2) (3)
TiCl2(CH3)2
(1) (2) (3)
Ti(CH3)4
(1) (2) (3)
TiCl4
(1) (2) (3)
45
TABLE IV. Comparison between GIPAW and IGAIM methods for various nuclei, using benchmarks of molecules, through the consideration of the deviation Δ and relative mean absolute deviation Δr. The GIPAW
51
V and
49
Ti NMR results were computed using (1)
3sP3p2×P4s2×P3d2×P and (2) 3sP3p3×P4s2×P3d3×P ultrasoft pseudopotentials. GIPAW Nucleus
51
IGAIM LAO basis set 6311++G(3df,3pd) ANO3ζ 6311++G(3df,3pd) ANO3ζ ccpCVQZ ccpCVQZ ccpCVQZ
1 n
Δa (ppm)
Δrb (%)
PPPW (1) (2)
σiso
5.9 10.4 9.4 9.5 8.8 0.8 1.5
δσ
16.9 13.3 2.9 2.7 
σiso
0.3 0.6 0.6 0.7 2.6 0.3 3.2
V
49
Ti
(1) (2)
31 c
P
Ref. [5] Ref. [5] Ref. [5]
29
Sic Cc
13
a
Mean absolute deviation calculated using Δx =
∑
i
n
GIPAW xiIGAIM ? xi
, where x and n are the
shielding parameters and the number of molecules respectively. The VO(OCH2CH2)3N molecule was dismissed from the statistic calculation.
b
Relative mean absolute deviation calculated using Δ r x =
1 n
∑
i
n
xiIGAIM ? xiGIPAW xiIGAIM
× 100 .
c
Calculations were performed using normconserving pseudopotential with the LDA exchangecorrelation functional. Δ and Δr calculations related to the {31P, 29Si, 13C} nuclei were accomplished with n = {3, 7, 5}, from Ref. [5].
46
TABLE V. Experimental and calculated 51V shielding parameters (δ iso, δ aniso, ηδ) using the USPPGIPAW method, in various vanadate compounds.
Theoretical (ppm) compound orthovanadate AlVO4 V(1) V(2) V(3) LaVO4 pyrovanadate 705a(738b) 633(670) 742(773) 616 V(1) V(2) V(1) V(2) Ca2V2O7 metavanadate NH4VO3 Mg(VO3)2 Ca(VO3)2 vanadate V2O5 622 718 5 310 519 468 484 429 271 371 0.07 0.01 0.03 0.90 0.45 612 ± 1 755 7 369 ± 1 504 ± 2 645 ± 1 818 323 336 ± 68 485 ± 29 0.11 ± 0.05 0.00 0.03 0.35 ± 0.10 0.25 ± 0.25 [102] [103] [95] [104] [104] 601 544 567 156 263 414 0.37 0.21 0.39 570 ± 1 534 ± 1 563 ± 1 240 ± 5 310 ± 3 517 ± 5 0.70 ± 0.03 0.30 ± 0.03 0.18 ± 0.03 [100] [101] [101] V(1) V(2)
a
Experimental (ppm)
site
δiso
δaniso
96 77 62 49 73 73 113 264 72 473
ηδ
0.55 0.86 0.50 0.65 0.89 0.53 0.49 0.26 0.36 0.62
δiso
744 ± 1 661 ± 1 776 ± 1 605 ± 1 604 ± 1 549 ± 1 639 ± 1 494 ± 1 575 ± 1 534 ± 1
δaniso
120 ± 6 87 ± 8 82 ± 7 50 ± 5 103 ± 2 57 ± 3 113 ± 7 262 ± 3 71 ± 3 530 ± 10
ηδ
0.72 ± 0.10 0.74 ± 0.17 0.88 ± 0.11 0.71 ± 0.05 0.34 ± 0.16 0.91 ± 0.10 0.90 ± 0.10 0.10 ± 0.10 0.54 ± 0.35 0.50 ± 0.03
Ref.
[97]
[98] [99] [99] [99]
αMg2V2O7 βMg2V2O7
628(603) 570(549) 669(639) 517(495) 576(570) 543(539)
βVOPO4
VOCl3 (103 K) complexe VO(OEt)(ONS) VO2[acpyinh]
a
Predictive 51V chemical shifts have been calculated with respect to the Eq. (24), using a = 1.047 and σref = 1939.
b
Relative 51V chemical shifts for AlVO4, α and βMg2V2O7, and Ca2V2O7 are reported relative to the reference values –2004, 1943, 1940 and –1959 respectively (a = 1 in Eq. (24)).
47
TABLE VI. Influence of the XC functional on the calculated 51V shielding parameters in VOCl3. method GIPAW clustera XC functional PW91 PW91 PW91 BLYP molecule B3PW91 B3LYP
a
σiso (ppm)
1944 1941 1924 1959 2185 2226
δaniso (ppm)
429 445 415 434 366 390
ηδ
0.03 0.03 0.00 0.00 0.00 0.00
A cluster of eleven VOCl3 entities have been used, keeping the geometry used for the GIPAW periodic NMR calculations. Shielding parameters are related to the central molecule.
48
Fig. 1. GIPAW method convergence using different vanadium pseudopotentials (see the Sec. II.A. & B. for the pseudopotential setting details). Calculated 51V isotropic shielding in VOCl3 molecule is plotted versus the planewave energy cutoff Ec.
(Color online) Fig. 2. Evolution of the 51V isotropic shielding components as a function of the number of projectors used in the USPPGIPAW calculation, for the VOCl3
G≠ molecule. The scale of σ Δp was reduced by a factor of 15 compared to σ bare0 and σ Δd .
(Color online) Fig. 3. Polyhedral projection of the various classes of vanadiumbased inorganic systems using representative compounds of Table V: (a) LaVO4 for orthovanadate, (b) βMg2V2O7 for pyrovanadate, (c) NH4VO3 for metavanadate, (d) and (e) represent the vanadate class with βVOPO4 and CaVO3 and (f) an organometallic complexe with VO(OEt)(ONS). Structural distortions are shown in terms of distance (given in ?) with their first coordination sphere.
Fig. 4. Plot of the 51V GIPAW absolute isotropic shielding versus experimental chemical shifts for the 18 vanadium sites referenced in the TABLE V. The solid line represents the linear correlation. All the fitted parameters are given in the upper right panel.
49
Fig. 5. Experimental versus calculated
51
V chemical shift tensor eigenvalues for the
various vanadate compounds of TABLE V. The solid line represents perfect agreement between calculation and experiment.
Fig. 6. Evolutions of the occupiedvirtual energy gap and
51
V shielding tensor
eigenvalues as a function of the HartreeFock mixing coefficient involved in the HandH hybrid exchangecorrelation functional for VOCl3.
50
Figure 1
51
Figure 2
52
Figure 3
53
Figure 4
54
Figure 5
55
Figure 6
56