Lionel Truflandier
email: lionel.truflandier@cnrs-imn.fr

Florent Michaël Paris

Florent Boucher
email: florent.boucher@cnrs-imn.fr

DFT investigation of 3d transition metal NMR shielding tensors in diamagnetic systems using the gauge-including projector augmented-wave method

come

DFT investigation of 3d transition metal NMR shielding

tensors in diamagnetic systems using the gauge-including projector augmented-wave method

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 nuclear-spin larger than ½, electric field gradient (EFG) tensors. The development of theoretical methods to calculate NMR properties is currently underway in several scientific communities. [2][3][4][5] 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 time-scale 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 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,[8][9][10][11][12][13] 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. [START_REF] Tossell | Computational Materials Chemistry[END_REF] To overcome such difficulties, Pickard and Mauri have developed the so-called "gaugeincluding projector augmented-wave" (GIPAW) pseudopotential approach in which the periodicity of the system is explicitly taken into account using a plane-wave 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 [START_REF] Mauri | [END_REF]18] is the possibility of keeping the nodal 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, [20][21][22][23] due to the difficulties involved in efficient pseudopotential development.

Nowadays, NMR spectroscopy applied to transition metals is widely used in the fields of coordination, bio-chemistry and solid state materials. Among the 3d transition metals, numerous NMR measurements on 51 V 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 49 Ti and 51 V 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 benchmark of titanium and vanadium complexes and validated by comparison to allelectron calculations. Application to 51 V containing extended systems will be discussed in Sec. IV. A first example of such an application has been published recently on the AlVO 4 system. [27] In this last part, we will finally concentrate on the relation between exchange-correlation 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,

{ } r r r A r r p r J ) ( 1 2 1 ) ( c , - - = . (1)
Here { } . A(r) is a vector potential connected to B through

) (r A B × ∇ = or ) ( ) ( 0 2 1 r r B r A - × =
, where r 0 is the gauge origin. The first and the 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 first-order-induced current j (1) (r). It yields a nonuniform induced magnetic field (1) in B which shields each nucleus N from B. The nuclear magnetic shielding tensor σ or the so-called chemical shift tensor defined as

3 N N (1) 3 N N (1) in ) ( 1 ) ( ) ( ∫ r r r r r j B r r B - - × = - = r d c σ , (2)
is a second-order magnetic response. The first-order induced current density ) ( (1) r j is obtained by means of perturbation theory applied to ) ( ' r J , [28] { }

) ( 1 Ψ Ψ ) ( 0 (1) (0) (1) ∑ r A r r p r j ρ c , o o o - - = ) ( ) ( ) 1 ( d ) 1 ( p r j r j + = . (3)
In this equation, the summation is over the occupied states o and 0 ρ is the unperturbed electron density. The ground state wave function

Ψ o is its corresponding first-order correction due to the magnetic field perturbation. ) ( , which depends on the first-order 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, ) ( (1) r j must be independent of the choice of gauge origin r 0 . Both ) ( is implicitly present in (1) Ψ o . Different approaches can be used to evaluate (1) Ψ o 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 (1) H . For a perturbation due to an external magnetic field,

B L H ⋅ = c 2 1 (1) , (4)
where

p r L × =
is the angular-momentum operator. Thus, the presence of (1) H in the calculation of (1) Ψ o is responsible for the implicit gauge dependence of ) (

) 1 ( p r j .
Due to incomplete basis set, the gauge origin independence on ) ( (1) r j is usually not completely verified, and it could in principle yield numerical divergence of the calculation of ) ( (1) r j

. 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 ) ( (1) r j is expected. The problem of different convergence rates is entirely solved when considering the simple case of an isolated closed-shell atom:

) ( ) 1 ( p r j
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 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 one-electron 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 ) ( (1) r j by performing a gauge transformation for each point of space. [36] The "continuous set of gauge transformations" method (CSGT) achieves gauge-invariance via a parametric function d(r) which is defined in real space and shifts continuously the gauge origin. The potential vector is redefined as,

)) ( ( 2
1 ) ( 0 r d r r B r A - - × = . ( 5
)
The type of CSGT method is determined by the choice of the ) (r d function. [2,19,36] If

) (r d
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 ) ( (1) r j is calculated. In other words, the function ) (r d 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. [37][38][39][40][41][42] 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 the difficulties associated with solid state systems, an alternative approach was proposed using the fully periodic GIPAW method. [5] B. The gauge-including projector augmented-wave approach In order to discuss the approximations introduced in the magnetic field dependant GIPAW approach, we need first to briefly describe the projector augmented-wave (PAW) electronic structure calculation method elaborated by Blöchl. [12] Within the frozen core approximation and the pseudopotential plane-wave formalism, the PAW method was developed by introducing an operator T that maps the true valence wave functions

Ψ onto pseudo-wave functions Ψ ~, Ψ Ψ Τ = .
The construction of T is carried out through the use of all-electron (AE) and pseudo (PS) atomic wave functions (so-called AE and PS partial waves), respectively i φ and i φ ~. As in other pseudopotential methods, a cutoff radius c r N, (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,

n n n n p ~N, N, N, N, ∑ ⎟ ⎠ ⎞ - ⎜ ⎝ ⎛ + = φ φ 1 Τ , (6)
Local projector functions The ability of the PAW method to reconstruct an AE wave function has allowed the use of the pseudopotential plane-wave 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. [44][45][46][47] Nevertheless, when considering a second-order 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,

∑ N, N, N, N, 2 2 n n n n c i c i e p ẽ B × r r B × r r B N N 1 Τ ⋅ - ⋅ ⎟ ⎟ ⎠ ⎞ - ⎜ ⎜ ⎝ ⎛ + = φ φ . ( 7
)
As a result, the GIPAW pseudo-eigenvector Ψ associated with the all-electroneigenvector Ψ is defined by Ψ Ψ is given by

+ + + ⎟ ⎟ ⎠ ⎞ - ⎜ ⎜ ⎝ ⎛ + = N N, , N, N, N N N, N, N N N, N, N ∑ ξ ξ O ξ ξ O ξ ξ O O m m n m n m n n p p ~φ φ φ φ . ( 8
)
If one applies the transformation given in Eq. ( 8) on the operator J(r) described in Eq.

(1), the GIPAW current density operator becomes

{ } ∑ N d N p N N ) ( ) ( Δ ) ( Δ ) ( ) ( 1 2 1 ) ( r ξ r J r J r ξ r r r A r r p r J N + ⎟ ⎠ ⎞ + ⎜ ⎝ ⎛ + - - = ' ' ' ' ' c ' ' , ' . (9)
The GIPAW nodal structure reconstruction leads to the introduction of the paramagnetic

) ( Δ p N ' r J and diamagnetic ) ( Δ d N ' r J
operators defined in the augmentation region N Ω ;

{ } { } 2 1 ) ( Δ N, N, N, N, N, , N, p N ∑ m m n m n m n n p ' ' , ' ' , p ' ⎟ ⎠ ⎞ - ⎜ ⎝ ⎛ = φ φ φ φ r r p r r p r J , (10) 2 )
( ) ( Δ N, N, N, N, N, , N, N d N ∑ m n n m n m n n p ' ' ' ' p c ' ' ⎟ ⎠ ⎞ - ⎜ ⎝ ⎛ - × = φ φ φ φ r r r r r r B r J . ( 11
)
If one develops J in powers of B and uses density-functional perturbation theory, [29] the GIPAW first-order current density is obtained and expressed in different contributions, [5] )

( ) ( ) ( ) ( (1) Δd (1) Δp (1) bare (1) ' ' ' ' r j r j r j r j + + = , (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 pseudo-Hamiltonian

B B HT T H + =
. Obviously, the expression for H depends entirely on the pseudopotential approach used: Either the norm-conserving [48,49] (NCPP) or the ultrasoft [50] (USPP) schemes. In this latter case, the relaxation of the norm-constraint imposes an additional generalised orthonormality constraint which must be solved via an overlap operator S,

m n, m , n , ~δ φ φ = N N S . ( 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 USPP-GIPAW approach. The work of Yates has permitted development and implementation of the USPP-GIPAW formalism. [51] Due to the introduction of the generalized orthonormality constraint, the first-order perturbed wave function

Ψ n given in Eq. ( 32) of Ref. [5] is redefined as,

(0) (1) (1) (1) Ψ ) )( ( Ψ n n n n S H G ε ε - = , (14)
with the Green function operator ( )

n ε G expressed through ( ) ∑ - = e e n ) ( e ) ( e n ε ε ε 0 0 Ψ Ψ G , (15)
with the sum running over empty states e. (1) H and (1) S are respectively the first-order perturbed GIPAW Hamiltonian and the first-order perturbed overlap matrix. The Green function involving virtual subspace is only used here for convenience in order to express the first-order 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 conjugate-gradient 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) (1) (0) (0) (0) ps (0) (1) (1) (0) (1) bare Ψ Ψ Ψ Ψ ) ' ( 2
1 ) Ψ ) )( ( Ψ ( Re 2 ) ( ∑ ∑ o o' o' oo' o o o o o o ' ' , ' c ' ' , ' S r r p r B r S H G r r p r j - × - - = ρ ε ε , ( 16
)
where o runs over the occupied states.

∑ = o o o ' ' (0) (0) ps Ψ Ψ 2 ) ' ( r r r ρ is the ground state pseudo-density. The paramagnetic augmentation current is given by [ ] (0) (1) (0) (0) p N (0) (0) p N N (0) (0) (1) (1) p N (0) N, (1) Δp Ψ Ψ Ψ ) ( Δ Ψ 2 Ψ ) ( Δ 2 1 Ψ 2 ) Ψ ) ( ) ( ) ( Δ Ψ ( 4Re ) ( ∑ ∑ o o' o' oo' o o o o o o o o ' ' , ' ic ' ' S r J r J r r B S H G r J r j - ⎟ ⎠ ⎞ ⋅ × + - ⎜ ⎝ ⎛ = ε ε
, [START_REF] Mauri | [END_REF] and the diamagnetic augmentation current is

Ψ ) ( Δ Ψ 2 ) ( ∑ N, (0) d N (0) (1) Δd o o o ' ' r J r j = . ( 18
)
The introduction of extra terms in the expression of ) (

(1) bare r j , resulting from the additional orthonormality constraint, yields more awkward calculations compared to the normconserving GIPAW method. The NCPP-GIPAW equations can be recovered by putting

1 = S
(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 r 0 = r N . [19] By reformulating Eqs. ( 36) and (37) of Ref. [5], it has been shown that the first-order 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 ) (

(1) bare ' r j and ) (

(1)

Δp ' r j
(the 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 first-order induced magnetic field is divided into four contributions which can be individually calculated, taking advantage of the linearity of Eq. ( 2),

) ( ) ( ) ( ) ( ) ( N (1) Δd N (1) Δp N (1) bare N core N (1) in r B r B r B r B r B + + + = . ( 19
) ) ( N core r B
, which depends only on the core electrons, i.e. of the isolated atom, is calculated once using the Lamb formula. [54] At this stage, several approximations are introduced to compute NMR chemical shift tensors from the GIPAW approach. Firstly, to evaluate the correction to ) ( N (1) in r B 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 on-site approximation neglects the effects of the augmentation currents of the neighbouring atoms to the shielding of the studied atom. Secondly, within periodic conditions,

(1) in B is formulated in reciprocal space using the Biot-Savart law.

Unfortunately, for a null vector of the reciprocal lattice (G = 0), ( )

0 (1) in = G B
becomes a macroscopic quantity. [START_REF] Mauri | [END_REF] 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 ) (

(1) bare ' r j contribution. Finally, the pseudopotential 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 plane-wave expansion, USPP-GIPAW calculations are able to reproduce NCPP-GIPAW 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: all-electron and USPP-GIPAW 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.

In order to validate the shielding tensor GIPAW calculations for the titanium and vanadium atoms, the USPP-GIPAW results have been compared to those obtained through the AE approach. The Gaussian 03 suite of programs [55] was used to compute all-electron magnetic response of molecules within the IGAIM approach, [30] combined with the "Perdew-Wang 91" exchange and correlation functional PW91. [56,57] Molecular geometries were optimized with symmetry constraints, using the B3LYP hybrid functional [58,59] with the 6-311+G(2d,p) basis set. [60][61][62][63] 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-ζ 6-311++G(3df,3pd) Pople's basis set developed by Watchers and Hay [62][63][64] for the first-row 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-ζ correlation-consistent basis set (cc-pCV5Z) 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 NMR-CASTEP codes, [5,20,67] respectively. The Brillouin zone was sampled using Monkhorst-Pack technique. [68] Relaxation of ionic positions were performed at an energy cutoff of 600 eV, using a k-point 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 k-point by the use of a super-cell approach, checking that the super-cell is large enough to avoid spurious interaction between periodic images. This condition was in general satisfied in a 12000

Bohr 3 (~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 r c , beyond which the pseudo-wave functions match the all-electron 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). r c 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×r c . 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:

( )

zz yy xx iso δ δ δ δ + + = 3 1 , ( 20
)
anisotropy component:

zz iso aniso δ δ δ - = , (21)
asymmetry component:

aniso yy xx δ δ δ η δ - = , (22)
with

iso yy iso xx iso zz δ δ δ δ δ δ - ≥ - ≥ - , (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

] [ ref ij ij a σ σ δ - ⋅ - = (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, first-principles 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 all-electron 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. [START_REF] Mauri | [END_REF]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 plane-wave basis sets. For second and third row elements, the core-valence 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 frozen-core PAW and fully relaxed calculations, previous studies have demonstrated that favorable choice of core-valence 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 first-order 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 53 Cr, 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 29 Si and 31 P 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 core-valence partition involved in the GIPAW chemical shielding tensors calculations for 3d elements. We present in Table I the shielding tensor calculated for 51 V in the well-known VOCl 3 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 non-rigid 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 all-electron 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 USPP-GIPAW 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 core-valence interaction was described by the Troullier-Martin [49] scheme, in the Kleinman-Bylander [75] form. To be consistent with the previous calculation, we used the same core-valence separation and projector allocation as for the (3s3p4s3d)-USPP. The cutoff 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 non-optimized (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 pseudo-atoms, 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 pseudo-wave functions. Keeping the constraints of normalization and continuity of two derivatives at r c , the pseudization function is optimized in order to minimize the kinetic energy beyond the cutoff wave vector q c . For the non-optimized USPP, using a default value of q c = 12.7 au ½ , the 51 V 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 q c 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 core-valence 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 plane-wave energy cutoff and on the AE and PS partial-wave function expansions. With respect to the "additive augmentation principle", [12] Blöchl has shown that the truncation of the partial-wave extension does not affect the completeness of the basis set, assuming the complementary participation of the plane-wave 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 plane-wave 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 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 super-cells. 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 plane-wave 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.

) ( ) ( ) ( ) ( ) ( ) ( N Δd N Δp N 0 G bare N 0 G bare N core N iso
The augmentation of the 4s state with two projectors has no effect on the isotropic shielding component. Indeed, since the paramagnetic augmentation current ) ' (

Δp r j is proportional to the angular momentum (Eqs. ( 4) and ( 17)), for a s angular momentum, only the bare ) (

Δd r j terms are dependent of the projector extension. Finally the 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 two-augmented 3p state, a shielding effect is obtained for a two-augmented 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 correction expressed in Eq. ( 25). 0 G bare ≠ σ is slightly affected by the pseudo-partial 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 VOCl 3 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 TiCl 3 CH 3 molecule are also discussed (Table III). The titanium USPP was built using the same core-valence separation and projector allocation as for the vanadium USPP 3s P 3p 2×P 4s 2×P 3d 2×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 correlation-consistent basis sets agree very well with the 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 USPP-GIPAW method in the calculation of shielding parameters for two molecules, namely VOCl 3 for the 51 V and TiCl 3 CH 3 for the 49 Ti, 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 all-electron calculations of the 51 V and 49 Ti isotropic shielding values have been reported in the literature for organometallic systems. [78][79][80][81][82][83][84] 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 3s P 3p 2×P 4s 2×P 3d 2×P and 3s P 3p 3×P 4s 2×P 3d 3×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 VF 5 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 TiCl 4 (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 %. 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 3s P 3p 2×P 4s 2×P 3d 2×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 13 C and 8.8 ppm for 31 P. The average deviation decreases from 17 to 13 ppm for the anisotropy parameters when we used the 3s P 3p 3×P 4s 3×P 3d 3×P USPP, but unfortunately the value related to shielding parameters increases to 10 ppm.

A global analysis of our results is given in Table

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 29 Si and 3.2 % for the 13 C. The same conclusions can be drawn for the 49 Ti 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 VOCl 3 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 VF 5 systems.

As a result, the 3s P 3p 2×P 4s 2×P 3d 2×P USPP is sufficient to achieve accurate 51 V isotropic chemical shift calculations with a reduced computational effort compared to a 3s P 3p 3×P 4s 2×P 3d 3×P USPP calculation. Furthermore, the calculation time using GIPAW method is of the order of IGAIM with the 6-311++G(3df,3pd) basis set, while it is considerably smaller when more extended basis sets such as cc-pCVXZ 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 plane-wave 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 USPP-GIPAW 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 Kohn-Sham (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 field-dependant Hamiltonian which yields additional terms in the shielding tensor expressions. [86][87][88] Moreover, these terms can be separated in scalar and spin-orbit coupling parts, depending on the level of approximation used. [87,89] Obviously, for a non-consistent use of methods, i.e. if two 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, all-electron 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 Koelling-Hammond equation. [90] Bouten and co-workers have studied NMR shielding predictions of 3d metal oxide (MO 4 n-with M = Cr, Mn, Fe) coupling zero-order 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 USPP-GIPAW and IGAIM results observed in Tables I andII.

Previous studies combining the ZORA and GIPAW methods [88,[91][92][93][94] have shown the influence of scalar relativity on 77 Se 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.

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 49 Ti and 51 V are in progress.

IV. APPLICATION TO EXTENDED SYSTEMS A. Results and discussion

Having validated, on various molecular systems, the NMR shielding calculations for 49 Ti and 51 V using the USPP-GIPAW 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 51 V nucleus, using a (3s P 3p 2×P 4s 2×P 3d 2×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] 3.

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 least-squares 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 VOCl 3 using an all-electron 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 NH 4 VO 3 , β-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 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 vanadium-based 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 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. [106][107][108][109][110] Linear response of crystalline or molecular orbitals to the magnetic field perturbation are strongly dependant on the occupied-virtual energy gap (E gap ) and the shape of the virtual orbitals, through the first-order corrected wave function (Eq. 14). Recent studies have shown that hybrid density functionals, which include a portion of Hartree-Fock (HF) exchange, partially overcome the "band gap error" problem in solid state systems. [111][112][113][114][115] 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 VOCl 3 inorganic system (bulk-optimized geometry have been kept, see Sec. II.C.). Shielding calculations were performed through the use of IGAIM method coupled with the 6-311++G(3df,3pd) basis set. Investigation of the influence of the exact exchange on shielding parameters has been performed using different exchange-correlation functionals. For GGAs, we have used the "Perdew-Wang 91" exchange and correlation functional PW91, [56,57] and the BLYP functional, which combined the "Becke's 1988" exchange and the "Lee-Yang-Parr" correlation functionals. [58,117] Hybrid XC functionals are defined by the following exchange-

correlation approximation, GGA C GGA X LDA X HF X hybrid XC E E )E 1 ( E E + Δ + - + = β α α . ( 26
)
Where HF X E is the "exact" HF exchange, LDA X E is the LDA exchange, Lee-Yang-Parr. [58] Results are collected in Table VI. Firstly, in order to probe the packing effect on the 51 V shielding parameters, we have used the cluster approximation using ten additional VOCl 3 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 IGAIM-cluster results, and validates both approaches.

By isolating a unique VOCl 3 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 half-and-half functional proposed by

Becke, [118] and defined with the following relation,

LYP C LDA X HF X HandH XC E )E 1 ( E E + - + = α α . ( 27
)
Evolution of the occupied-virtual 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 occupied-virtual energy splitting. E gap discrepancy between pure DFT exchange (α = 0, called HandH 0 ), and quasi-full HF exchange (α = 0.9, called HandH 0.9 ) is about 0.30 a.u.. Calculation using Hartree-Fock level of theory (results not shown)

gives a value of 0.50 a.u. compared to 0.19 and 0.41 for HandH 0 and HandH 0.9 . These results agree with the well-known LDA-GGA underestimation and Hartree-Fock overestimation of occupied-virtual 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 occupied-virtual energy gap induced by the use of hybrid XC functionals.

Nevertheless, according to an extensive study of the influence of pure exchange on 57 Fe isotropic shielding through GIAO-DFT calculations, [107] Schreckenbach has demonstrated that three factors are responsible for the improvement induced by the use of hybrid functionals: enhancement of the occupied-virtual 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 51 V anisotropy and asymmetry NMR parameters are probably linked to a fundamental DFT deficiency rather than GIPAW built-in 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 34 51 V and 49 Ti in diamagnetic molecular-like 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. , where x and n are the shielding parameters and the number of molecules respectively. The VO(OCH 2 CH 2 ) 3 N molecule was dismissed from the statistic calculation.

b Relative mean absolute deviation calculated using c Calculations were performed using norm-conserving pseudopotential with the LDA exchange-correlation functional. Δ and Δ r calculations related to the { 31 P, 29 Si, 13 C} nuclei were accomplished with n = {3, 7, 5}, from Ref. [5].

H

  associated with the eigenvalue o ε and

  the unperturbed charge density 0 ρ , is called the "diamagnetic" contribution and corresponds to the uniform circulation of the electrons.

  expand the pseudo-wave function locally onto the pseudo-atomic orbitals. The index n refers to the angular-momentum quantum numbers and to an additional number, which is used if there is more than one projector per angular-momentum channel. Constraints[12] are imposed by the PAW method: The evaluation of an observable quantity represented by an operator O can be expressed in terms of pseudowith 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.

  an isolated molecular system such as the VOCl 3 molecule, there are no surface currents (see Sec. II.B.) and the 0 G bare = σ component of Eq. (25) should tend to zero.

  the exchange correction and correlation parts of GGA functional. We use the three-parameter B3 exchange functional defined by Becke,[59] leading to a value of α = 0.2. The correlation GGA functionals GGA C E are taken as the Perdew-Wang 91,[56,57] and

Fig. 1 .(

 1 Fig. 1. GIPAW method convergence using different vanadium pseudopotentials (see the Sec. II.A. & B. for the pseudopotential setting details). Calculated 51 V isotropic shielding in VOCl 3 molecule is plotted versus the plane-wave energy cutoff E c .

Fig. 4 .

 4 Fig. 4. Plot of the 51 V GIPAW absolute isotropic shielding versus experimental chemical

Fig. 5 .

 5 Fig. 5. Experimental versus calculated 51 V chemical shift tensor eigenvalues for the

Fig. 6 .

 6 Fig.6. Evolutions of the occupied-virtual energy gap and51 V shielding tensor

  Table V for all the compounds having more than one vanadium site. The results are given between parentheses and allow straightforward assignments of the 51 V resonances in the AlVO 4 , 23 αand β-Mg 2 V 2 O 7 , and Ca 2 V 2 O 7 compounds. With a discrepancy of the order of a few ppm, we are able to discriminate inequivalent vanadium sites exhibiting close isotropic chemical shifts.

TABLE I .

 I Convergence of the51 V absolute isotropic and anisotropic shielding parameters as a function of the vanadium valence states involved in USPP-GIPAW calculations for the VOCl 3 molecule. The multi-projector USPP is defined by the notation nl k×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). Core radius r c and atomic reference energies (in parentheses) of the local atomic pseudopotential v loc are given in a.u.. TABLE II.51 V NMR shielding parameters in various molecular systems. The GIPAW calculations were performed using the (1) 3s P 3p 2×P 4s 2×P 3d 2×P and (2) 3s P 3p 3×P 4s 2×P 3d 3×P ultrasoft pseudopotentials and compared to the IGAIM calculations performed with the (1) 6-311++G(3df,3pd) and (2) ANO-3ζ LAO basis sets. Abbreviation for the VO(OCH 2 CH 2 ) 3 N complex. Computation of the NMR shielding parameters was not tractable for this molecule using the ANO-3ζ basis set.

	51 V				σ iso (ppm)			δ σ (ppm)	η δ
	Molecule	Valence State USPP/LAO	r c c GIPAW v loc c	Number of projectors σ iso (ppm) δ σ (ppm) IGAIM GIPAW IGAIM GIPAW	IGAIM
	VOCl 3	(4s3d) (1) (2)	2.4 -1910 p(-0.5) -1947	-1904 -1952	4s 2×P 3d 2×P -455 -463	-1806 -483 -464	0.00 -353 0.00	0.00 0.00
	[V(CO) 6 ] -	(3p4s3d) (1)	2.5	97	f(0.0)	91	3p 2×P 4s 2×P 3d 2×P 0	-1910 0	-434 n/a	n/a
		(3s3p4s3d) (2)	2.0	89	f(0.0)	3s P 3p 2×P 4s 2×P 3d 2×P 76 0	-1910 0	-455 n/a	n/a
	VF 5	(2p3s3p4s3d) a 0.8/2.0 f(0.0) 2p P 3s P 3p 2×P 4s 2×P 3d 2×P (1) -1220 -1233 -9 (2) -1280 -1258 1	-1920 -6 11	0.01 -461 0.14	0.00 0.00
	VOF 3	all-electron b (1)	-	-1177	-	-1177	-	336	-1904 317	-483 0.00	0.00
		(2)		-1212		-1214		335	348	0.00	0.00
	VOClF 2	(1)		-1415		-1418		293	290	0.37	0.27
		(2)		-1451		-1458		293	297	0.36	0.35
	VON a	(1)		-1546		-1548		-42	-61	0.02	0.01
		(2)		-1584		_			-46	_	0.02	_
	VOCl 2 F	(1)		-1663		-1663		-345	-358	0.46	0.47
		(2)		-1700		-1707		-349	-347	0.46	0.45
	VO(CH 3 ) 3	(1)		-3034		-3020		-1641	-1652	0.00	0.00
		(2)		-3074		-3057		-1647	-1615	0.00	0.00

a A core radius of 0.8 and 2.0 a.u. was used for the 2p and for the remaining states respectively. b IGAIM/6-311++G(3df,3pd). c a

TABLE IV .

 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 GIPAW51 V and49 Ti NMR results were computed using (1) 3s P 3p 2×P 4s 2×P 3d 2×P and (2) 3s P 3p 3×P 4s 2×P 3d 3×P ultrasoft pseudopotentials.

		GIPAW	IGAIM			Δ a (ppm)	Δ r b (%)
	Nucleus	PP-PW	LAO basis set		σ iso	δ σ	σ iso
	51 V	(1)	6-311++G(3df,3pd)	5.9	16.9	0.3
		(2)	ANO-3ζ		10.4	13.3	0.6
	49 Ti	(1)	6-311++G(3df,3pd)	9.4	2.9	0.6
		(2)	ANO-3ζ		9.5	2.7	0.7
	31 P c	Ref. [5]	cc-pCVQZ		8.8	-	2.6
	29 Si c	Ref. [5]	cc-pCVQZ		0.8	-	0.3
	13 C c	Ref. [5]	cc-pCVQZ		1.5	-	3.2
	a Mean absolute deviation calculated using	= x 1 n Δ	n i ∑	IGAIM i x	-	GIPAW i x

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.

Ti NMR shielding parameters in various molecular systems. The GIPAW calculations were performed using the (1) 3s P 3p 2×P 4s 2×P 3d 2×P and (2,3) 3s P 3p 3×P 4s 2×P 3d 3×P ultrasoft pseudopotentials and compared to the IGAIM calculations performed with the (1) 6-311++G(3df,3pd), (2) ANO-3ζ and (3) cc-pCV5Z LAO basis sets. 49 Ti a Predictive 51 V chemical shifts have been calculated with respect to the Eq. ( 24), using a


Thierry Dumas
email: thierry.dumas@inria.fr

Aline Roumy
email: aline.roumy@inria.fr

Christine Guillemot
email: christine.guillemot@inria.fr

IMAGE COMPRESSION WITH STOCHASTIC WINNER-TAKE-ALL AUTO-ENCODER

Keywords: Image compression, sparse representations, auto-encoders, Orthogonal Matching Pursuit

This paper addresses the problem of image compression using sparse representations. We propose a variant of autoencoder called Stochastic Winner-Take-All Auto-Encoder (SWTA AE). "Winner-Take-All" means that image patches compete with one another when computing their sparse representation and "Stochastic" indicates that a stochastic hyperparameter rules this competition during training. Unlike auto-encoders, SWTA AE performs variable rate image compression for images of any size after a single training, which is fundamental for compression. For comparison, we also propose a variant of Orthogonal Matching Pursuit (OMP) called Winner-Take-All Orthogonal Matching Pursuit (WTA OMP). In terms of rate-distortion trade-off, SWTA AE outperforms auto-encoders but it is worse than WTA OMP. Besides, SWTA AE can compete with JPEG in terms of ratedistortion.

INTRODUCTION

Auto-encoders are powerful tools for reducing the dimensionality of data. Deep fully-connected auto-encoders [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] are traditionally used for this task. However, two issues have so far prevented them from becoming efficient image compression algorithms: they can only be trained for one image size and one compression rate [START_REF] Krizhevsky | Using very deep autoencoders for content-based image retrieval[END_REF][START_REF] Deng | Binary coding of speech spectrograms using a deep auto-encoder[END_REF]. [START_REF] Toderici | Variable rate image compression with recurrent neural networks[END_REF] attempts to solve both issues. The authors train an auto-encoder on image patches so that images of various sizes can be compressed. Their auto-encoder is a recurrent [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF] residual auto-encoder that performs variable rate image compression after a single training. But all image patches have the same rate and therefore different distortions due to the texture complexity variety in image patches. In addition, recurrence, which is equivalent to scalability in image compression, is not optimal in terms of rate-distortion trade-off [START_REF] Schwarz | Overview of the scalability video coding extension of the H.264/AVC standard[END_REF][START_REF] Sullivan | Standardized extensions of high efficiency video coding (HEVC)[END_REF].

Instead, we propose to perform learning on whole images under a global rate-distortion constraint. This is done through This work has been supported by the French Defense Procurement Agency (DGA).

Winner-Take-All (WTA), which can be viewed as a competition between image patches when computing their representation. Furthermore, auto-encoders architecture must adapt to different rates. Therefore, during training, the WTA parameter that controls the rate is stochastically driven. These contributions give rise to Stochastic Winner-Take-All Auto-Encoder (SWTA AE).

Notation

Vectors are denoted by bold lower case letters and matrices by upper case ones. X j denotes the j th column of a matrix X. X F is the Frobenius norm of X. X 0 counts the number of non-zero elements in X. The support of a vector x is supp (x) = {i | x i = 0}.

STOCHASTIC WINNER-TAKE-ALL AUTO-ENCODER (SWTA AE)

We now present our Stochastic Winner-Take-All Auto-Encoder (SWTA AE) whose architecture is shown in Figure 1. SWTA AE is a type of auto-encoder. An auto-encoder is a neural network that takes an input and provides a reconstruction of this input. We justify below two of the most critical choices for the SWTA AE architecture.

Strided convolution

A compression algorithm must process images of various sizes. However, the most efficient neural networks [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF][START_REF] He | Delving deep into rectifiers: surpassing humanlevel performance on imagenet classification[END_REF] require that all images have the same size. Indeed, they include both convolutional layers and fully-connected layers, and the number of parameters of the latters directly depends on the image size. This imposes to train one architecture per image size. That is why our proposed SWTA AE only contains convolutional layers. Its encoder has two convolutional layers and its decoder has two deconvolutional layers [START_REF] Noh | Learning deconvolution network for semantic segmentation[END_REF]. Each layer i ∈ 1, 4 consists in convolving the layer input with the bank of filters W (i) , adding the biases b (i) and applying a mapping g (i) , producing the layer output. For the borders of the layer input, zero-padding of width p (i) is used. Max-pooling is a core component of neural networks [START_REF] Scherer | Evaluation of the pooling operations in convolutional architectures for object recognition[END_REF] that downsamples its input representation by appling a max filter to non-overlapping sub-regions. But max-pooling increases the rate. Indeed, if the encoder contains a maxpooling layer, the locations of maximum activations selected during pooling operations must be recorded and transmitted to the corresponding unpooling layer in the decoder [START_REF] Zeiler | Visualizing and understanding convolutional networks[END_REF][START_REF] Badrinarayanan | SegNet: a deep convolutional encoderdecoder architecture for image segmentation[END_REF]. Instead, for i ∈ 1, 2 , we downsample using a fixed stride s (i) > 1 for convolution, which does not need any signaling.

Semi-sparse bottleneck

The bottleneck is the stack of feature maps denoted Z ∈ R h×w×65 in Figure 1. Z is the representation of the input image that is processed in Section 2.3 to give the bitstream.

We propose to apply a global sparse constraint that povides control over the coding cost of Z. This is called Winner-Take-All (WTA). Let us define WTA via a mapping g α : R h×w×64 → R h×w×64 , where α ∈]0, 1[ is the WTA parameter. g α keeps the α × h × w × 64 most representative coefficients in its input tensor, i.e. those whose absolute values are the largest, and sets the rest to 0. g α only applies to the output of the convolution in the second layer involving the first 64 filters in W (2) , producing the first 64 sparse feature maps in Z. Figure 1 displays these sparse feature maps in orange. Varying α leads to various coding costs of Z. Note that [START_REF] Makhzani | Winner-take-all autoencoders[END_REF] uses WTA, but our WTA rule is different and g α does not apply to specific dimensions of its input tensor as this constraint is not relevant for image compression.

A patch of the input image might be represented by a portion of the first 64 sparse feature maps in Z that only contains zeros. We want to ensure that each image patch has a minimum code in Z to guarantee a sufficient quality of reconstruction per patch. That is why the last feature map in Z is not sparse. Figure 1 displays it in red. We have noticed that, during the training in Section 4.2, SWTA AE learns to store in the last feature map a subsampled version of its input image.

Bitstream generation

The coefficients of the non-sparse feature map in Z are uniformly quantized over 8-bits and coded with a Huffman code. The non-zero coefficients of the 64 sparse feature maps in Z are uniformly quantized over 8-bits and coded with a Huffman code while their position is coded as explained hereafter. Figure 1 defines a coordinate system (x, y, z) for Z. The non-zero coefficients in Z are scanned along (x, y, z) where z changes the fastest. The position along z is coded with a fixed-length code and, for each pair (x, y), the number of non-zero coefficients along z is coded with a Huffman code. This unequivocally characterizes the position of each non-zero coefficient in Z. We have observed that this processing is effective in encoding the position of the non-zero coefficients.

WINNER-TAKE-ALL ORTHOGONAL MATCHING PURSUIT (WTA OMP)

SWTA AE is similar to Orthogonal Matching Pursuit (OMP) [START_REF] Davenport | Analysis of orthogonal matching pursuit using the restricted isometry property[END_REF], a common algorithm for image compression using sparse representations [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF]. The difference is that SWTA AE computes the sparse representation of an image by alternating convolutions and mappings whereas OMP runs an iterative decomposition of the image patches over a dictionary. More precisely, let x ∈ R m be an image patch. Given x and a dictionary D ∈ R m×n , OMP finds a vector of coefficients y ∈ R n with k < m non-zero coefficients so that Dy equals to x approximatively.

For the sake of comparison, we build a variant of OMP called Winner-Take-All Orthogonal Matching Pursuit (WTA OMP). More precisely, let X ∈ R m×p be a matrix whose columns are formed by p image patches of dimension m and Y ∈ R n×p be a matrix whose columns are formed by p vec-tors of coefficients of dimension n. WTA OMP first decomposes each image patch over D, see [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. Then, it keeps the γ × n × p coefficients with largest absolute value for the nlength sparse representation of the p patches and sets the rest to 0, see [START_REF] Krizhevsky | Using very deep autoencoders for content-based image retrieval[END_REF]. The support of the sparse representation of each patch has therefore been changed. Hence the need for a final least-square minimization, see [START_REF] Deng | Binary coding of speech spectrograms using a deep auto-encoder[END_REF].

Algorithm 1 : WTA OMP Inputs: X ∈ R m×p , D ∈ R m×n , k < m and γ ∈]0, 1[. For each j ∈ 1, p , Y j = OMP (X j , D, k) (1)
I = f γ (Y) (2)
For each j ∈ 1, p ,

Z j = min z∈R n X j -Dz 2 2 st. supp (z) = supp (I j ) (3)
Output: Z ∈ R n×p .

TRAINING

Before moving on to the image compression experiment in Section 5, SWTA AE needs training. Similarly, a dictionary D ∈ R m×n must be learned for WTA OMP.

Training data extraction

We extract 1.0 × 10 5 RGB images from the ILSVRC2012 ImageNet dataset [START_REF] Deng | ImageNet: a large-scale hierarchical image database[END_REF]. The RGB color space is transformed into YCbCr and we only keep the luminance channel.

For SWTA AE, the luminance images are resized to 321 × 321. M ∈ R 321×321 denotes the mean of all luminance images. σ ∈ R * + is the mean of the standard deviation over all luminance images. Each luminance image is subtracted by M and divided by σ. These images are concatenated into a training set ∆ ∈ R 321×321×(1.0×10 5 ) .

For WTA OMP, η = 1.2 × 10 6 image patches of size √ m × √ m are randomly sampled from the luminance images. We remove the DC component from each patch. These patches are concatenated into a training set Γ ∈ R m×η .

SWTA AE training

As explained in Section 2.2, α tunes the coding cost of Z. If α is fixed during training, all the filters and the biases of SWTA AE are learned for one rate. That is why we turn α into a stochastic hyperparameter during training. This justifies the prefix "Stochastic" in SWTA AE. Since there is no reason to favor some rates during training, we sample α according to the uniform distribution U[µ -, µ + ], where µ -> 0 and µ + < 1. We select µ = 1.8 × 10 -1 and = 1.7 × 10 -1 to make the support of α large. At each training epoch, α is drawn for each training image of ∆.

As shown in Section 2.1, SWTA AE can process images of various sizes. During training, we feed SWTA AE with random crops of size 49×49 of the training images of ∆. This accelerates training considerably. The training objective is to minimize the mean squared error between these cropped images and their reconstruction plus l 2 -norm weights decay. We use stochastic gradient descent. The gradient descent learning rate is fixed to 2.0 × 10 -5 , the momentum is 0.9 and the size of mini-batches is 5. The weights decay coefficient is 5.0 × 10 -4 . Our implementation is based on Caffe [START_REF] Jia | Caffe: convolutional architecture for fast feature embedding[END_REF]. It adds to Caffe the tools introduced in Sections 2.2 and 2.3.

Dictionary learning for WTA OMP

Given Γ, k < m and γ ∈ ]0, 1[, the dictionary learning problem is formulated as [START_REF] Toderici | Variable rate image compression with recurrent neural networks[END_REF].

min D,Z1,...,Zη 1 η η j=1 Γ j -DZ j 2 2 st. ∀j ∈ 1, η , Z j 0 ≤ k st. η i=j Z j 0 ≤ γ × n × η (4) (4)
is solved by Algorithm 2 which alternates between sparse coding steps that involve WTA OMP and dictionary updates that use stochastic gradient descent. Given Γ and p ∈ N * + , let φ be a function that randomly partitions Γ into η p = η / p mini-batches X (1) , ..., X (ηp) , where, for i ∈ 1, η p , X (i) ∈ R m×p . Mini-batches make learning very fast [START_REF] Hinton | Lecture 6a, overview of mini-batch gradient descent[END_REF].

Algorithm 2 : dictionary learning for WTA OMP. Inputs:

Γ ∈ R m×η , k < m, γ ∈]0, 1[, p ∈ N * + and ε ∈ R * + . D ∈ R m×n is randomly initialized. ∀j ∈ [|1, n|], D j ← D j / D j 2
For several epochs do:

X (1) , ..., X (ηp) = φ (Γ, p) ∀i ∈ [|1, η p |], Z (i) = WTA OMP X (i) , D, k, γ D ← D -ε ∂ X (i) -DZ (i) 2 F ∂D ∀j ∈ [|1, n|], D j ← D j / D j 2 Output: D ∈ R m×n .
For OMP, given Γ, a dictionary D ∈ R m×n is learned using K-SVD [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF]

IMAGE COMPRESSION EXPERIMENT

After training in Section 4, we compare the rate-distortion curves of OMP, WTA OMP, SWTA AE, JPEG and JPEG2000 on test luminance images.

Image CODEC for SWTA AE

Each input test luminance image is pre-processed similarly to the training in Section 4.1. The mean learned image M is interpolated to match the size of the input image. Then, the input image is subtracted by this interpolated mean image and divided by the learned σ. The encoder of SWTA AE computes Z. The bitstream is obtained by processing Z as detailed in Section 2.3.

Image CODEC for OMP and WTA OMP

A luminance image is split into 8×8 non-overlapping patches. The DC component is removed from each patch. The DC components are uniformly quantized over 8-bits and coded with a fixed-length code. OMP (or WTA OMP) finds the coefficients of the sparse decompositions of the image patches over D (or D). The non-zero coefficients are uniformly quantized over 8-bits and coded with a Huffman code while their position is coded with a fixed-length code. Then, for WTA OMP only, the number of non-zero coefficients of the sparse decomposition of each patch over D is coded with a Huffman code.

Comparison of rate-distortion curves

In the literature, there is no reference rate-distortion curve for auto-encoders. We compare SWTA AE with JPEG and JPEG2000 2 even though the image CODEC of SWTA AE is less optimized. Furthermore, we compare SWTA AE with its 2 JPEG and JPEG2000 code: http://www.imagemagick.org/script/index.php non-sparse Auto-Encoder counterpart (AE). AE has the same architecture as SWTA AE but its Z only contains non-sparse feature maps. Note that, to draw a new point in the AE ratedistortion curve, AE must be first re-trained with a different number of feature maps in Z.

Figure 2 shows the rate-distortion curves of OMP, WTA OMP, AE, SWTA AE, JPEG and JPEG2000 for two of the most common images: LENA and BARBARA. In terms of rate-distortion trade-off, SWTA AE outperforms AE and WTA OMP is better than OMP. This highlights the value of WTA for image compression. When we compare SWTA AE with WTA OMP, we see that iterative decomposition is more efficient for image compression using sparse representations. Moreover, SWTA AE can compete with JPEG. We also ran this image compression experiment on several crops of LENA and BARBARA and observed that the relative position of the six rate-distortion curves was comparable to the relative positioning in figure 2. The size of the test image does not affect the performance of SWTA AE. More simulation results and a complexity analysis for OMP, WTA OMP and SWTA AE can be found on the web page https://www.irisa.fr/temics/demos/ NeuralNets/AutoEncoders/swtaAE.htm.

CONCLUSIONS AND FUTURE WORK

We have shown that, SWTA AE is more adaptated to image compression than auto-encoders as it performs variable rate image compression for any size of image after a single training and provides better rate-distortion trade-offs.

So far, our work has focused on the layer of auto-encoders which is dedicated to coding. Yet, many avenues of research are still to be explored to improve auto-encoders for image compression. For instance, [START_REF] Tobias Springenberg | Striving for simplicity: the all convolutional net[END_REF] proves that removing a maxpooling layer and increasing the stride of the previous convolution, as we do, harms neural networks. This has to be addressed.

Fig. 1 :

 1 Fig. 1: SWTA AE architecture.

  1 , and the parameters m and n are optimized with an exhaustive search. This leads to m = 64 and n = 1024. For SWTA AE, the same values for m and n are used for training D via Algorithm 2. Moreover, k = 15, γ = 4.5 × 10 -3 , p = 10 and ε = 2.0 × 10 -2 .

Fig. 2 :

 2 Fig. 2: Evolution of PNSR with the rate.

K-SVD code: http://www.cs.technion.ac.il/ elad/software/