David Lauvergnat

André Nauts

Smolyak scheme for solving Schrödinger equation: application to malonaldehyde in full dimensionality

In 1963, Smolyak introduced an elegant approach to overcome the exponential scaling, in terms of the number of variables (i.e. degrees of freedom), of the size of direct products [S. A. Smolyak Soviet Mathematics Doklady, 4, 240 (1963)]. The main idea is to replace a single and large direct product by a sum of selected small direct products. This approach has been applied in various scientific domains and, in particular, in quantum dynamics, where Avila and Carrington used it for the first time in 2009 [G. Avila and T. Carrington, J. Chem. Phys., 131, 174103 (2009)]. Since then, several calculations on systems with 12 degrees of freedom have been published by Avila and Carrington, by other groups and by ourselves.

In the present study and to push the limit to larger and complex systems, the Smolyak scheme is combined with the use of an on-the-fly calculations of the kinetic energy operator [A. Nauts and D. Lauvergnat, Mol. Phys. 116, 3701 (2018)] in our home-made Fortran code, ElVibRot-Tnum-Tana. This procedure is applied to compute the tunneling splitting of malonaldehyde in full dimensionality (21D) using a recent potential of Mizukami et al. [W. Mizukami, S. Habershon, and D.P. Tew, J. Chem. Phys. 141, 144310 (2014)]. Our tunneling splitting calculations, 21.7±0.3 cm -1 and 2.9±0.1 cm -1 , show an excellent agreement with experimental values, 21.6 cm -1 and 2.9 cm -1 for the normal isotopologue and the mono-deuterated one, respectively.

I. INTRODUCTION:

In quantum dynamics, one of the major difficulties is to obtain a compact representation of the delocalized wave function or wave packet for a system with d degrees of freedom. Formally, those quantum states are delocalized over the entire space and to overcome this difficulty several strategies can be used. For instance, with a stochastic approach such as quantum Diffusion Monte Carlo (DMC), [START_REF] Suhm | Quantum Monte Carlo studies of vibrational states in molecules and clusters[END_REF][START_REF] Clary | Torsional diffusion Monte Carlo: A method for quantum simulations of proteins[END_REF][START_REF] Wang | Full-dimensional quantum calculations of ground-state tunneling splitting of malonaldehyde using an accurate ab initio potential energy surface[END_REF][START_REF] Mizukami | A compact and accurate semi-global potential energy surface for malonaldehyde from constrained least squares regression[END_REF][START_REF] Benoit | Quaternion formulation of diffusion quantum Monte Carlo for the rotation of rigid molecules in clusters[END_REF] the ground state of a quantum state can be obtained without the need of expanding the wave function on a basis set. However, the excited states are more difficult to compute and usually, the approximate fixed node approach is used. [START_REF] Suhm | Quantum Monte Carlo studies of vibrational states in molecules and clusters[END_REF] The Feynman Path integral Molecular Dynamics approaches [START_REF] Marx | Ab initio path integral molecular dynamics: Basic ideas[END_REF][START_REF] Beutier | Computing thermal Wigner densities with the phase integration method[END_REF] or related approaches, Ring Polymer MD [START_REF] Habershon | Ring-Polymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space[END_REF] or Centroid MD [START_REF] Cao | The formulation of quantum statistical mechanics based on the Feynman path centroid density. I. Equilibrium properties[END_REF] , enable to represent the delocalized wave packet by several coupled trajectories. Those methods are intensively used for large systems such as liquids. Finally, the wave function or wave packet can be expanded on a multidimensional basis set or grid with the usual standard approaches. In this case, those methods have to be adapted in order to avoid the exponential scaling, as a function of d, of the basis sets and/or the grid sizes. There are several ways to overcome this difficulty, without being exhaustive, one can use: (i) a highly contracted scheme, such as, VSCF [START_REF] Horn | Vibrational states of very floppy clusters : Approximate separability and the choice of good curvilinear coordinates for XeHe2 , I2He[END_REF][START_REF] Bowman | MULTIMODE: A code to calculate rovibrational energies of polyatomic molecules[END_REF][START_REF] Benoit | No Title[END_REF][START_REF] Christiansen | Vibrational structure theory: new vibrational wave function methods for calculation of anharmonic vibrational energies and vibrational contributions to molecular[END_REF] (Vibrational Self Consistent Field), MCTDH [START_REF] Beck | THE MULTICONFIGURATION TIME-DEPENDENT HARTREE (MCTDH) METHOD: A HIGHLY EFFICIENT ALGORITHM FOR PROPAGATING WAVEPACKETS[END_REF] (Multi-configuration time-dependent Hartree) or its multilayer version [START_REF] Wang | Nonperturbative quantum simulation of time-resolved nonlinear spectra: Methodology and application to electron transfer reactions in the condensed phase[END_REF][START_REF] Hammer | Intramolecular proton transfer in malonaldehyde: Accurate multilayer multi-configurational time-dependent Hartree calculations[END_REF][START_REF] Vendrell | Multilayer multiconfiguration time-dependent Hartree method: implementation and applications to a Henon-Heiles hamiltonian and to pyrazine[END_REF] (ii) a nondirect-product basis-set where the basis functions are selected in terms of excitations (like in most of the VSCF implementations) or when the basis set is pruned. [START_REF] Dawes | A multidimensional discrete variable representation basis obtained by simultaneous diagonalization[END_REF] However, with these schemes, the exponential scaling remains when a direct-product grid is required. (iii) a Smolyak scheme, [START_REF] Smolyak | Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions[END_REF] which avoids the exponential scaling as function of d and can be viewed as a combination of a non-direct product basis set and a sparse grid adapted to the basis set. Its use on realistic applications is recent. In 2009 Avila and Carrington [START_REF] Avila | Nonproduct quadrature grids for solving the vibrational Schrödinger equation[END_REF] showed the efficiency of this scheme for the vibrational levels of water and the six stretches of SF6. Later on, we used this combination (with a crude implementation) to study the torsional motion of nitric acid in a 9D-ab initio model. [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF] Since then, several studies have been published by Avila and Carrington, [START_REF] Avila | Using a pruned basis, a non-product quadrature grid, and the exact Watson normal-coordinate kinetic energy operator to solve the vibrational Schrödinger equation for C2H4[END_REF][START_REF] Avila | Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D[END_REF][START_REF] Avila | Reducing the cost of using collocation to compute vibrational energy levels: Results for CH 2 NH[END_REF] by Matyus group [START_REF] Avila | Full-dimensional (12D) variational vibrational states of CH 4 •F - : Interplay of anharmonicity and tunneling[END_REF][START_REF] Avila | Toward breaking the curse of dimensionality in (ro)vibrational computations of molecular systems with multiple large-amplitude motions[END_REF] and by ourselves. [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF][START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF][START_REF] Chen | Smolyak Algorithm Adapted to a System-Bath Separation: Application to an Encapsulated Molecule with Large-Amplitude Motions[END_REF] The main goal of the present study is to show that, with the Smolyak scheme, large realistic molecular systems (larger than 12 degrees of freedom) with complex coordinates can be studied.

However, to study those kinds of systems with the standard approaches, one needs a good set of coordinates to reduce the coupling between the modes of the Hamiltonian and therefore to use a compact basis set. In other words, with this "good" set of coordinates, the full Hamiltonian must be well approximated by a sum of uncoupled low dimensional (such as 1D or 2D) Hamiltonians. For instance, when rectilinear normal modes are well adapted to the system or to the chemical or physical processes, the full Hamiltonian is expressed as a sum of uncoupled 1D-quadratic Hamiltonians plus some anharmonic contributions. However, for some processes with large amplitude motions or for systems with several minima or with flat potential energy surfaces, the curvilinear coordinates are usually better adapted than the rectilinear usual ones. Unfortunately, with curvilinear coordinates, Q, the expression of the kinetic energy operator [START_REF] Podolsky | Quantum-mechanically correct form of Hamiltonian function for conservative systems[END_REF][START_REF] Wilson | Molecular Vibrations : The Theory of Infrared and Raman Vibrational Spectra[END_REF][START_REF] Nauts | Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians[END_REF] , 𝑇 " ! , with a Euclidean normalization, is not a simple expression anymore and is given as follows:

𝑇 " ! #𝑸, 𝜕 𝑸 ' = - ℏ # 2 - 1 𝐽(𝑸) 𝜕. 𝜕𝑄 $ 𝐽(𝑸) 𝐺 $% (𝑸) 𝜕. 𝜕𝑄 % &,& $,%() (1)
In the previous expression, 𝐺 $% (𝑸) are the contravariant components of the metric tensor. The volume element, 𝑑𝜏 ! , is Euclidean and given by 𝐽(𝑸)𝑑𝑄 ) ⋯ 𝑑𝑄 & , where 𝐽(𝑸) is the Jacobian. However, this Euclidean volume element is usually not adapted to the multidimensional basis set. Therefore, a new volume element, 𝑑𝜏 * = 𝜌(𝑸)𝑑𝑄 ) ⋯ 𝑑𝑄 & , must be imposed for the basis set to be orthonormal. This implies that the wave function and the operators must be transformed accordingly as follows: [START_REF] Podolsky | Quantum-mechanically correct form of Hamiltonian function for conservative systems[END_REF][START_REF] Wilson | Molecular Vibrations : The Theory of Infrared and Raman Vibrational Spectra[END_REF][START_REF] Nauts | Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians[END_REF][START_REF] Fábri | Numerically constructed internal-coordinate Hamiltonian with Eckart embedding and its application for the inversion tunneling of ammonia[END_REF]

Where Ψ ! and Ψ * are, respectively, the wave function normalized with the Euclidean volume element and the non-Euclidean one. Furthermore, 𝐎 < ! and 𝐎 < * are « operationally identical » operators, adapted to the volume elements 𝑑𝜏 ! and 𝑑𝜏 * , respectively, in order to preserve the values of the matrix elements of these operators. [START_REF] Nauts | Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians[END_REF][START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF] In the remainder of the text the index r of the operators and the wave functions will often not be mentioned when the volume element used is 𝑑𝜏 * . Nowadays, the KEO can be handled relatively easily, either numerically [START_REF] Harthcock | Calculation of Kinetic Energy Terms for the Vibrational Hamiltonian: Application to Large-Amplitude Vibrations Using One-, Two-, and Three-Dimensional Models[END_REF][START_REF] Senent | Determination of the kinetic energy parameters of non-rigid molecules[END_REF][START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Mátyus | Toward black-box-type full-and reduced-dimensional variational (ro)vibrational computations[END_REF] with possibly very complex coordinates or analytically with some specific kinds of coordinates [START_REF] Chapuisat | Vector parametrization of the N-body problem in quantum mechanics: Polyspherical coordinates[END_REF][START_REF] Gatti | A general expression of the exact kinetic energy operator in polyspherical coordinates[END_REF][START_REF] Ndong | Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach[END_REF][START_REF] Ndong | Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach: General formulation and removal of singularities[END_REF] (Jacobi, polyspherical, ...).

In the present study, the combination of the Smolyak scheme and a set of curvilinear coordinates for which the KEO is obtained numerically with TNUM [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF] , will be applied to the calculation of the tunneling splitting of malonaldehyde in full dimensionality i.e. with 21 degrees of freedom. Indeed, Malonaldehyde is a molecule which presents two equivalent minima associated to the proton transfer (see Figure 1). Therefore, its vibrational ground state energy level is split due to the tunneling through the barrier. The experimental tunneling splitting is known experimentally with great precision: 21.583 cm -1 for the normal istopologue [START_REF] Turner | Microwave spectroscopic study of malonaldehyde. 4. Vibration-rotation interaction in parent species[END_REF][START_REF] Varella | Real-time observation of intramolecular proton transfer in the electronic ground state of chloromalonaldehyde: an ab initio study of time-resolved photoelectron spectra[END_REF] and 2.915 cm -1 for the mono-deuterated istopologue (the transfer atom is the deuterium). [START_REF] Baughcum | Microwave spectroscopic study of malonaldehyde. 3. Vibration-rotation interaction and one-dimensional model for proton tunneling[END_REF] From a theoretical point of view this molecule has been intensively studied. In particular several simulations have been done in full dimensionality (21D): Reaction Surface Hamiltonian (RSH), [START_REF] Tew | A reaction surface Hamiltonian study of malonaldehyde[END_REF] instanton, 47-49 DMC, 3,4 MCTDH [START_REF] Hammer | Intramolecular proton transfer in malonaldehyde: Accurate multilayer multi-configurational time-dependent Hartree calculations[END_REF][START_REF] Hammer | Multiconfigurational time-dependent Hartree calculations for tunneling splittings of vibrational states: Theoretical considerations and application to malonaldehyde[END_REF][START_REF] Schröder | Theoretical studies of the tunneling splitting of malonaldehyde using the multiconfiguration time-dependent Hartree approach[END_REF] approaches. In particular and due to the availability of recent full dimensionality potential energy surfaces, [START_REF] Wang | Full-dimensional quantum calculations of ground-state tunneling splitting of malonaldehyde using an accurate ab initio potential energy surface[END_REF][START_REF] Mizukami | A compact and accurate semi-global potential energy surface for malonaldehyde from constrained least squares regression[END_REF] the DMC and MCTDH lead to tunneling splitting values in good agreement with the experimental ones (see Table I). [START_REF] Turner | Microwave spectroscopic study of malonaldehyde. 4. Vibration-rotation interaction in parent species[END_REF][START_REF] Varella | Real-time observation of intramolecular proton transfer in the electronic ground state of chloromalonaldehyde: an ab initio study of time-resolved photoelectron spectra[END_REF] , b) ref [START_REF] Baughcum | Microwave spectroscopic study of malonaldehyde. 3. Vibration-rotation interaction and one-dimensional model for proton tunneling[END_REF] In the second section, we will present some wave function notations needed in the next sections. In section III, we will describe the first Smolyak scheme, which is needed for understanding the second Smolyak scheme presented in section IV. In section V, the second Smolyak scheme will be applied to malonaldehyde and finally, a conclusion and perspectives will be given in section VI.

II. Wave function representation and notations

Let the quantum state of a system, often loosely called the wave function of a system, be described by the ket |𝜓⟩, which can be expanded on a set of basis functions expressed as a product of n @𝑏 + ! , B:

|𝜓⟩ = -𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ) • @𝑏 + " ) B ⋯ @𝑏 + ! , B ⋯ @𝑏 + # / B + " …+ # (3)
where 𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ) are the coefficients of the wave function expansion. To each @𝑏 + ! , B corresponds a primitive basis function, 𝑏 + ! , , belonging to a primitive basis set, 𝑩 , = H𝑏 ) , , 𝑏 # , ⋯ I.

Furthermore, each 𝑩 , is associated to one or several coordinates denoted collectively as 𝑸 , , and the value of a basis function, 𝑏 + ! , , at 𝑸 , , is expressed as 𝑏 + ! , (𝑸 , ).

We are now in a position to switch to a coordinate representation, 𝜓(𝑸 ) , ⋯ 𝑸 , , ⋯ 𝑸 / ), which can be obtained as follows:

𝜓(𝑸 ) , ⋯ 𝑸 , , ⋯ 𝑸 / ) = -𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ) • 𝑏 + " ) (𝑸 ) ) ⋯ 𝑏 + ! , (𝑸 , ) ⋯ 𝑏 + # / (𝑸 / ) + " …+ # (4)
Equation ( 4) provides a way to compute the value of the wave function, 𝜓(𝑸 ) , ⋯ 𝑸 , , ⋯ 𝑸 / ), at a given point from its basis set coefficients, 𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ). We can also perform the reverse operation, i.e. compute these coefficients from 𝜓(𝑸 ) , ⋯ 𝑸 , , ⋯ 𝑸 / ) using an integration:

𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ) = J𝑏 + " ) , … 𝑏 + ! , , … 𝑏 + # / @𝜓B = K 𝑏 + " ) (𝑸 ) ) ⋯ 𝑏 + ! , (𝑸 , ) ⋯ 𝑏 + # / (𝑸 / )𝜓(𝑸 ) , ⋯ 𝑸 , , ⋯ 𝑸 / )𝑑𝜏 (5)
Furthermore, to perform the numerical integration, we need a grid representation of the wave function. Therefore, each basis set, 𝑩 , , is associated to a grid with 𝑛𝑞 , points, 𝑮 , = O𝑸 ) , , 𝑸 # , ⋯ 𝑸 , are, respectively, the ui th grid point and its corresponding weight. The coefficients 𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ) are thus given by:

𝝍 𝑩 𝟏 ⋯𝑩 𝒏 (𝑘 ) , … , 𝑘 , , … 𝑘 / ) = -𝜔 1 " ) 𝑏 + " ) #𝑸 1 " ) ' ⋯ 𝜔 1 ! , 𝑏 + ! , #𝑸 1 ! , ' ⋯ 𝜔 1 # / 𝑏 + # / #𝑸 1 # / ' 𝜓(𝑸 1 " ) , ⋯ 𝑸 1 ! , , ⋯ 𝑸 1 # / ) 1 " …1 # (6)
Moreover, at these grid points, the wave function is simply given as:

𝝍 4 " ⋯4 # (𝑢 ) , … , 𝑢 , , … 𝑢 / ) = 𝜓(𝑸 1 " ) , ⋯ 𝑸 1 ! , , ⋯ 𝑸 1 # / ) (7)
Several remarks are in order:

(i) Each primitive basis set, 𝑩 , can be a 1D-basis set (such as a harmonic oscillator basis set, a Fourier basis ... ) or 2D or 3D-basis sets (such as the spherical harmonics or a Wigner basis set…). So, for instance, in the case of a 2D-basis set, the corresponding grid, 𝑮 , , is a 2D-grid (not necessarily a rectangular one). In this context, the physical dimension, i.e. the number of coordinates, d, is larger or equal to n, the number of basis sets, 𝑩 , . (ii) Usually, the grid points, 𝑸 1 ! , , associated to a 1D-grid are the gaussian quadrature points and the weights are the corresponding gaussian weights.

(iii) The transformations of a wave function from a basis representation to a coordinate representation, Eq (4), or the reverse, Eq (5), are not numerically efficient so that sequential transformations are used instead (see section IV for more details).

(iv) As already mentioned in the introduction and in order to avoid the exponential scaling as a function of n or d, it is essential to select some basis functions (or to prune them with respect to the full direct-product). Several ways to perform this selection exist. For instance, in the case of one of the simplest selection, ∑ (𝑘 , -1)

/ ,()

≤ 𝐿, the total number of basis functions is

(67/)! 6!•/!
instead of (𝐿 + 1) / . This selection is equivalent to limiting the polynomial degree of the basis function to L or to limiting the excitation to L. Unfortunately, in the case of this simple scheme, it is not trivial to select or prune the point grid so that, for instance, one should use the Smolyak scheme instead.

III. FIRST SMOLYAK SCHEME:

In Smolyak's paper of 1963, [START_REF] Smolyak | Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions[END_REF] the main expression of the first scheme is given by Eq (8), in which 𝑺 $;< ) is usually a multidimensional grid or a basis set. Its main advantage is, by performing a selection (see below), to circumvent the exponential scaling of the number of grid points or number of basis sets when the number of degrees of freedom, d, or the number of basis sets, n, increases. In our new implementation, 𝑺 $;< ) can represent either a multidimensional grid or a basis set.

𝑺 $;< ) = -∆𝑺 ℓ " ) ⨂ ⋯ ∆𝑺 ℓ ! , ⋯ ⨂ ∆𝑺 ℓ # / >?|𝓵|?6 (8) 0 ≤ ℓ ) ⋯ ℓ , + ⋯ ℓ / ≤ 𝐿 and ℓ , ∈ [0, 𝐿] (9)
In Eq (8), the index 1 of 𝑺 $;< ) stands for the first Smolyak scheme and the ∆𝑺 ℓ ! , are defined as follows:

_ ∆𝑺 > , = 𝑺 > , ∆𝑺 ℓ ! , = 𝑺 ℓ ! , -𝑺 ℓ ! B) , (10)
From a simple point of view, Eq. ( 8) is a sum (or, more precisely, a union) of direct-products

(∆𝑺 ℓ " ) ⨂ ⋯ ∆𝑺 ℓ ! , ⋯ ⨂ ∆𝑺 ℓ # / )
where the ∆𝑺 ℓ ! , are the "differences" between restricted primitive basis sets or grids, 𝑺 ℓ ! , (see below). Furthermore, |𝓵|, is the sum of the ℓ , , and the parameter 𝐿 adds a constraint on the ℓ , (Eq. ( 9)). This constraint is fundamental in the Smolyak scheme since it allows the reduction of the number of direct-products, ∆𝑺 ℓ " ) ⨂ ⋯ ∆𝑺 ℓ ! , ⋯ ⨂ ∆𝑺 ℓ # / . Indeed, this number grows as a polynomial in n, with degree 𝐿 and its expression is given by

(67/)! 6!•/! .
First, let us focus on the Smolyak scheme in the case of a basis expansion for which each 𝑺 ℓ ! , is a restricted finite basis set, with 𝑛𝑏 ℓ ! , basis functions,

𝑩 ℓ ! , = a𝑏 ) , , 𝑏 # , ⋯ 𝑏 /C ℓ ! ! , b.
In order to avoid any confusion, an index ℓ , is added to the basis set notation, 𝑩 , , to label the different possible restricted sets. Thus, using Eq (10), the ∆𝑩 ℓ ! , (ℓ , ∈ [0, 𝐿]) are given by:

⎩ ⎪ ⎨ ⎪ ⎧ ∆𝑩 > , = O𝑏 ) , ⋯ 𝑏 /C ( ! , P ∆𝑩 ) , = O𝑏 /C ( ! 7) , ⋯ 𝑏 /C " ! , P ⋮ ∆𝑩 ℓ ! , = a𝑏 /C ℓ ! )" ! 7) , ⋯ 𝑏 /C ℓ ! ! , b (11)
The relation between ℓ , and 𝑛𝑏 ℓ ! , is essential in the Smolyak schemes. Indeed, one can use any increasing integer sequences [START_REF] Avila | Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D[END_REF][START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF][START_REF] Avila | Pruned bases that are compatible with iterative eigensolvers and general potentials: New results for CH3CN[END_REF][START_REF] Powers | The effect of the condensed-phase environment on the vibrational frequency shift of a hydrogen molecule inside clathrate hydrates[END_REF] and, the sequences associated to two different restricted basis sets, 𝑩 ℓ ! , and 𝑩 ℓ * D (𝑖 ≠ 𝑗), can be different. Some particular sequences are worth mentioning:

(i) For 1D-basis sets (orthonormal polynomials, Fourier basis set ...), one of the simplest sequences is: 𝑛𝑏 ℓ ! , = 1 + ℓ , .

(ii) In the present study, we will also use arithmetic sequences, 𝑛𝑏 ℓ ! , = 𝐴 , + 𝐵 , • ℓ , , where 𝐴 , > 0 and 𝐵 , > 0.

(iii) For spherical harmonic basis sets (with 𝒀 D ! E ! as basis functions), we can define any increasing sequence between ℓ , and 𝑗 , EFG , where 𝑗 , EFG is the largest value for 𝑗 , . [START_REF] Powers | The effect of the condensed-phase environment on the vibrational frequency shift of a hydrogen molecule inside clathrate hydrates[END_REF][START_REF] Benoit | Does cage quantum delocalisation influence the translation-rotational bound states of molecular hydrogen in clathrate hydrate?[END_REF] Given the definitions of the 𝑛𝑏 ℓ ! , sequences, the ∆𝑩 ℓ ! , (Eqs (10) and ( 11)) and the constraint on the ℓ , (Eq. ( 9)), the Smolyak scheme applied to the corresponding multidimensional basis set (Eq (8)) allows the selection of an n-dimensional basis set, reduced in size with respect to the full directproduct basis set. In particular, when all the 𝑛𝑏 ℓ ! , sequences are ≤ 𝐿) at the end of section II, providing that 𝑘 , = 1 + ℓ , and the basis functions are 𝑏 + ! , (𝑸 , ).

𝑛𝑏 ℓ ! , = 1 + ℓ , ,
For example, with 𝑛 = 2, 𝐿 = 3, and 𝑛𝑏 ℓ ! , = 1 + ℓ , , the number of 2D basis functions i.e. the reduced number of direct-products, is 10 (see Table II), whereas the number of 2D basis functions would be 16 without the Smolyak scheme.

Table II. List of the direct-products of the first scheme in terms of basis functions (fourth column) in Eq (8) with 𝒏 = 𝟐 and 𝑳 = 𝟑, such that |𝓵| = 𝓵 𝟏 + 𝓵 𝟐 ≤ 𝟑. For both basis sets, 𝒏𝒃 𝓵 𝒊 𝒊 = 𝟏 + 𝓵 𝒊 (i=1, 2).

ℓ ) ℓ # |𝓵| = ℓ ) + ℓ # ∆𝑩 ℓ " ) ⨂∆𝑩 ℓ $ # 0 1 2 3 0 0 1 2 3
{𝑏 ) ) • 𝑏 ) # } {𝑏 # ) • 𝑏 ) # } {𝑏 H ) • 𝑏 ) # } {𝑏 I ) • 𝑏 ) # } 0 1 2 1 1 2 3
{𝑏 ) ) • 𝑏 # # } {𝑏 # ) • 𝑏 # # } {𝑏 H ) • 𝑏 # # } 0 1 2 2 3 {𝑏 ) ) • 𝑏 H # } {𝑏 # ) • 𝑏 H # } 0 3 3 {𝑏 ) ) • 𝑏 I # }
The selection by means of the Smolyak scheme seems intricate. However, it allows circumventing the exponential scaling of the number of grid points as well. More precisely, to adapt the Smolyak scheme to a grid in Eq (8), each 𝑺 ℓ ! , has to be viewed as a primitive restricted grid, 𝑮 ℓ ! , , with 𝑛𝑞 ℓ ! , grid points. When one uses Gaussian quadratures, the points of the grid 𝑮 ℓ ! , are different from those of another grid 𝑮 ℓ J ! , (ℓ r , ≠ ℓ , ), although some points might be identical. It is important to note, that the grid difference, ∆𝑮 ℓ ! , , (Eq (10)) has to be understood as the union of the grid points of 𝑮 ℓ ! , and 𝑮 ℓ ! B)

,

where the weights of 𝑮 ℓ ! B) , are negatives. However, if one can use nested grids, where all the points of one grid are included in a larger grid, [START_REF] Avila | Nonproduct quadrature grids for solving the vibrational Schrödinger equation[END_REF][START_REF] Avila | Using a pruned basis, a non-product quadrature grid, and the exact Watson normal-coordinate kinetic energy operator to solve the vibrational Schrödinger equation for C2H4[END_REF][START_REF] Avila | Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D[END_REF] the common points appear only once in ∆𝑮 ℓ ! , and their corresponding weights are the differences between the weights of 𝑮 ℓ ! , and 𝑮 ℓ ! B) ,

.

Finally, with this first Smolyak scheme, a wave function can be represented as a sum of Smolyak contributions and each contribution, @𝜓 K,ℓ " ⋯ℓ ! ⋯ℓ # B, is expanded on a small direct-product (given by the ℓ , ) of either basis sets or grids.

|𝜓⟩ = -@𝜓 K,ℓ " ⋯ℓ ! ⋯ℓ # B >?|𝓵|?6 (12)
IV. SECOND SMOLYAK SCHEME:

Smolyak 19 also proposed a second scheme, Eq (13), expressed directly in terms of the 𝑺 ℓ ! , , (

without the differences ∆𝑺 ℓ ! , ) and strictly equivalent to Eq (8). This new expression can be easily obtained by expanding the ∆𝑺 ℓ ! , in terms of the 𝑺 ℓ ! , . When the constraint on the ℓ , is defined as in Eq. ( 9), the coefficients, 𝐷 ℓ " ,ℓ $ ,…ℓ # , are (-1) 6B|𝓵| 𝐶 /B) 6B|𝓵| , [START_REF] Wasilkowski | Explicit Cost Bounds of Algorithms for Multivariate Tensor Product Problems[END_REF] where 𝐶 /B) 6B|𝓵| are binomial coefficients. When other constraints are used, the 𝐷 ℓ " ,ℓ $ ,…ℓ # can be computed numerically by directly transforming the first Smolyak scheme into the second one.

𝑺 $;< # = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • 𝑺 ℓ " ) ⨂ ⋯ 𝑺 ℓ ! , ⋯ ⨂𝑺 ℓ # / 6B/7)?|𝓵|?6 (13)
In Eq (13), the index 2 of 𝑺 $;< # stands for the second scheme and the boundaries on the ℓ , are slightly different from those of Eq. ( 9). Indeed, the lower bound of |𝓵| is 𝐿 -𝑛 + 1 for Eq (13) instead of 0 for Eq. ( 9). This feature results straightforwardly from transforming the first Smolyak scheme into the second one. Like for the first scheme, the previous expression can be used either with primitive basis sets or primitive grids. Some useful features of this second Smolyak scheme are the following:

(i) The Smolyak terms, 𝑩 ℓ " ) ⨂ ⋯ 𝑩 ℓ ! , ⋯ ⨂𝑩 ℓ # / , are small direct-products associated to (ℓ ) ⋯ ℓ , ⋯ ℓ / )
in which the size of the basis sets, 𝑩 ℓ ! , , are 𝑛𝑏 ℓ ! , . In a way, each small directproduct is an approximation of the full direct-product.

(ii) In each Smolyak term, 𝑩 ℓ " ) ⨂ ⋯ 𝑩 ℓ ! , ⋯ ⨂𝑩 ℓ # / , a n-dimensional basis function, 𝑏 + " ) (𝑸 ) ). 𝑏 + $ # (𝑸 # ) ⋯ 𝑏 + ! , (𝑸 , ) ⋯ 𝑏 + # / (𝑸 / )
, can be present only once (when 𝑘 , ≤ 𝑛𝑏 ℓ ! , ) or not present at all (when 𝑘 , > 𝑛𝑏 ℓ ! , ). Therefore, it is convenient to associate to this basis function and this Smolyak term, a parameter, 𝜀 + " ⋯+ ! ⋯+ # ℓ " ⋯ℓ ! ⋯ℓ # , equal to 1 or 0 according to whether the basis function is, respectively present or not in the Smolyak term.

(iii) In Eq. ( 13) , some n-dimensional basis functions,

𝑏 + " ) (𝑸 ) ). 𝑏 + $ # (𝑸 # ) ⋯ 𝑏 + ! , (𝑸 , ) ⋯ 𝑏 + # / (𝑸 / ),
are present in several Smolyak terms. For instance, the first basis function

(𝑏 ) ) (𝑸 ) ). 𝑏 ) # (𝑸 , ) ⋯ 𝑏 ) , (𝑸 , ) ⋯ 𝑏 ) / (𝑸 / )
) is present in all the terms of Eq (13) (see Table III).

However, since Eq. ( 13) is strictly equivalent to Eq (8), where a given n-dimensional basis function appears only once, if one regroups all identical basis functions from the different Smolyak terms including their coefficients, 𝐷 ℓ " ,ℓ $ ,…ℓ # , the basis function has to appear only once.

In other words, the following relation is fulfilled for all the n-dimensional basis functions:

- 𝐷 ℓ " ,ℓ $ ,…ℓ # • 𝜀 + " ⋯+ ! ⋯+ # ℓ " ⋯ℓ ! ⋯ℓ # 6B/7)?|𝓵|?6 = 1 (14)
This relation results automatically from the transformation of the first Smolyak scheme, where the coefficients, 𝐷 ℓ " ,ℓ $ ,…ℓ # , are absent, to the second one where the coefficients, 𝐷 ℓ " ,ℓ $ ,…ℓ # ,which can be negative, are present. Thus, as required, the n-dimensional basis functions selected in the second Smolyak scheme are strictly equivalent to those selected in the first Smolyak scheme.

𝒊 = 𝟏 + 𝓵 𝒊 (i=1, 2). 𝑘 ) , 𝑘 # 1,1 2,1 3,1 4,1 1,2 2,2 3,2 1,3 2,3 1,4 ℓ ) ℓ # 𝑩 ℓ " ) ⨂𝑩 ℓ $ # 𝐷 ℓ " ,ℓ $ 𝜀 + " ,+ $ ℓ " ,ℓ $ 2 3 0 0 {𝑏 ) ) , 𝑏 # ) , 𝑏 H ) }⨂{𝑏 ) # } {𝑏 ) ) , 𝑏 # ) , 𝑏 H ) , 𝑏 I ) }⨂{𝑏 ) # } -1 1 1 1 1 1 1 1 1 1 2 1 1 {𝑏 ) ) , 𝑏 # ) }⨂{𝑏 ) # , 𝑏 # # } {𝑏 ) ) , 𝑏 # ) , 𝑏 H ) }⨂{𝑏 ) # , 𝑏 # # } -1 1 1 1 1 1 1 1 1 1 1 1 0 1 2 2 {𝑏 ) ) }⨂{𝑏 ) # , 𝑏 # # , 𝑏 H # } {𝑏 ) ) , 𝑏 # ) }⨂{𝑏 ) # , 𝑏 # # , 𝑏 H # } -1 1 1 1 1 1 1 1 1 1 1 0 3 {𝑏 ) ) }⨂{𝑏 ) # , 𝑏 # # , 𝑏 H # , 𝑏 I # } 1 1 1 1 1 -𝐷 ℓ " ,ℓ $ • 𝜀 + " ,+ $ ℓ " ,ℓ $ #?|𝓵|?H 1 1 1 1 1 1 1 1 1 1
To illustrate these features, let us take the 2D example used for the first Smolyak scheme with n=2, L=3 and with the 𝑛𝑏 ℓ ! , = 1 + ℓ , (see Table III). With the constraints on |𝓵| ,(𝐿 -𝑛 + 1 ≤ |𝓵| ≤ 𝐿, see Eq.( 13)), the second Smolyak scheme has 7 direct-products 𝑩 ℓ " ) ⨂𝑩 ℓ $ # made up of the ten 2D-basis functions given in Table II of the first Smolyak scheme. Furthermore, Eq ( 14) is fulfilled for all the 10 basis functions (last line of Table III) and in particular, one can see that the first basis function, 𝑏 ) ) • 𝑏 ) # , is present in all the terms (𝜀 ),) ℓ " ,ℓ $ = 1).

(iv) Like for the first scheme, a wave function can be represented as a "weighted" sum of Smolyak contributions, @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B, and each contribution is expanded on a small direct-product (parameterized by the ℓ , ) of either restricted basis sets or restricted grids.

|𝜓⟩ = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B 6B/7)?|𝓵|?6 (15)
The @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B on the basis set is noted as, 𝜓

𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # #
(: ) , while on the grids it is noted,

𝜓 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (:
). Moreover, several mixed representations of |𝜓⟩ can be used: (a) for some indices i (for instance from 1 to i), the wave function is expanded on the basis and for the other indices, the wave function is expanded on the grid.

𝝍 𝑩 𝟏 ⋯𝑩 𝒊 𝑮 𝒊,𝟏 ⋯𝑮 𝒏 (… ) = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ ! ! 𝑮 ℓ !," !," ⋯𝑮 ℓ # # (⋯ ) 6B/7)?|𝓵|?6 (16)
(b) Some terms of the Smolyak scheme are expanded on the grids while other are on the basis sets. In the following relation, all Smolyak terms are on the basis sets except one term (𝓵 r ) which is on the grid.

𝝍 𝒎𝒊𝒙𝒆𝒅 (… ) = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (⋯ ) 6B/7)?|𝓵|?6 𝓵R𝓵 J + 𝐷 ℓ J " ,ℓ J $ ,…ℓ J # • 𝜓 𝑮 ℓ -" " ⋯𝑮 ℓ -! ! ⋯𝑮 ℓ -# # (⋯ ) (17)
(c) a mixture of (a) and (b). All, these features are very important for the computation of the action of an operator on the wave function (see next section).

(v) The expression of the scalar product, between two wave functions, |𝜓⟩ and |𝜒⟩, is rather simple:

⟨𝜒|𝜓⟩ = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • J𝜒 ℓ " ⋯ℓ ! ⋯ℓ # @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B 6B/7)?|𝓵|?6 (18)
In other words, the scalar product is a "weighted" sum of scalar product contributions on the Smolyak terms (see the proof in Appendix).

(vi) The number of elements (basis functions, 𝑁𝐸 -, or grid points, 𝑁𝐸 S ) of the Smolyak scheme are simply:

𝑁𝐸 -= - { 𝑛𝑏 ℓ ! , / ,() 6B/7)?|𝓵|?6 𝑁𝐸 S = - { 𝑛𝑞 ℓ ! , /
,() 6B/7)?|𝓵|?6 (19) 𝑁𝐸 S can be viewed as the total number of grid points whereas 𝑁𝐸 -, can be much larger than the total number of basis functions, NB, without the duplicated ones (see iii). Therefore, it is important to store a wave function in a compact basis representation (without duplicates), 𝝍 -" ⋯-# (: ). In this case, we need a mapping between a basis function in the small direct-products of the second Smolyak scheme and an index of the compact basis representation. Actually, this mapping is stored as a vector of 𝑁𝐸 -integers, which gives the indices of the nD-basis functions of the compact representation as a function of the indices of the basis functions of the second Smolyak scheme. This vector can be huge (see the malonaldehyde section) and is one of the memory bottlenecks of our implementation.

For instance, for n=21, with 𝑛𝑏 ℓ ! , = 𝑛𝑞 ℓ ! , = ℓ , + 1 and L=8 (see Table IV), the number of basis function elements, 𝑁𝐸 -, is half a billion, while the number of basis functions is only 4 million. (vi) When the grid is needed and in order to improve the numerical exactness of the action of an operator on a wave function, one needs more grid points than basis functions. Therefore, in the standard pseudo-spectral approach, for each index i, the number of grid points of 𝑮 , is larger than the number of basis function of 𝑩 , (𝑛𝑞 ℓ ! , > 𝑛𝑏 ℓ ! , ).

In the context of the Smolyak scheme, although it is possible to use the same feature for all small direct-products by using two different sequences in ℓ , for 𝑛𝑞 ℓ ! , and 𝑛𝑏 ℓ ! , . However, in this way, the total number of grid points, 𝑁𝐸 S , will increase dramatically. For example, with n=21, the sequence IV).

To overcome this difficulty, the same sequences (or almost the same sequences) can be used, but the L associated to the grid, LG, has to be larger than the L associated to the basis set, LB. In our implementation, LG and the sequences in ℓ , for the 𝑛𝑞 ℓ ! , enable to define the number of terms in the second Smolyak scheme, hence, 𝑁𝐸 S and 𝑁𝐸 -are of the same order (they are equal, when all

𝑛𝑏 ℓ ! , = 𝑛𝑞 ℓ ! ,
). Furthermore, the parameter, LB, and the sequences in ℓ , for the 𝑛𝑏 ℓ ! , enable to define the number of the total basis functions NB. Then, when the wave function needs to be stored in the compact basis representation, some basis functions from the Smolyak terms,

𝝍 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (:
), which are not present in the compact basis representation, 𝝍 𝑩 " ⋯𝑩 # (: ), are not

transferred from 𝝍 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ) to 𝝍 𝑩 " ⋯𝑩 # (:
). This procedure is certainly not optimal; however, it is extremely simple.

IV.a) Basis-grid transformation with the second Smolyak scheme:

The transformation of a quantum state, |𝜓⟩, from a basis representation to a grid (or the reverse) is fundamental in quantum dynamics. Indeed, the computational efficiency of the operator action on |𝜓⟩ depends strongly on the nature of the operator (local on the grid such as the potential or nonlocal such as the KEO) and the representation of the quantum state. Formally, the basis to grid transformation and reverse can be done using, respectively, Eq (4) and Eq (6). However, they are not efficient when n is large so that, for a direct-product expansion of the quantum state, one uses a sequential transformation. [START_REF] Light | Discrete-variable representations and their utilization[END_REF] Since within the second Smolyak scheme, the quantum state is expressed as a sum of small directproduct contributions on each Smolyak term (Eq (15)), one can use the usual sequential transformation for each term. For a given Smolyak term (parameterized by (ℓ ) , … ℓ , , … ℓ / )) and assuming that the grid to basis transformation is already performed for the index 1 to i-1, the transformation for the index i (𝑮 ℓ ! , to 𝑩 ℓ ! , ) is obtained as follows:

𝜓 𝑩 ℓ " " ⋯𝑩 ℓ !)" !)" 𝑩 ℓ ! ! 𝑮 ℓ !," !," ⋯𝑮 ℓ # # (… 𝑘 ,B) , 𝑘 , , 𝑢 ,7) … ) = -𝜔 1 ! , 𝑏 + ! , #𝑸 1 ! , ' 1 ! (),/2 ℓ ! ! • 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ !)" !)" 𝑮 ℓ ! ! 𝑮 ℓ !," !," ⋯𝑮 ℓ # # (… 𝑘 ,B) , 𝑢 , , 𝑢 ,7) … ) (20)
The reverse transformation for the index i (𝑩 ℓ ! , to 𝑮 ℓ ! , ) is obtained as follows:

𝜓 𝑮 ℓ " " ⋯𝑮 ℓ !)" !)" 𝑮 ℓ ! ! 𝑩 ℓ !," !," ⋯𝑩 ℓ # # (… 𝑢 ,B) , 𝑢 , , 𝑘 ,7) … ) = -𝑏 + ! , #𝑸 1 ! , ' • 𝜓 𝑮 ℓ " " ⋯𝑮 ℓ !)" !)" 𝑩 ℓ ! ! 𝑩 ℓ !," !," ⋯𝑩 ℓ # # (… 𝑢 ,B) , 𝑘 , , 𝑘 ,7) … ) + ! (),/C ℓ ! ! (21)
Each index uj on the grid, 𝑮 ℓ * D , runs from 1 to 𝑛𝑞 ℓ * D and the index kj on the basis, 𝑩 ℓ * D , runs from 1 to

𝑛𝑏 ℓ * D .
Remark: for both equations, ( 20) and (21), and for given i and ui, the values of the grid point 𝑸 1 ! , are usually different for different values of ℓ , . Formally, an index ℓ , is missing on 𝑸 1 ! , .

To summarize, a full transformation needs three steps: (i) a loop on the Smolyak terms i.e. the selected (ℓ ) ⋯ ℓ , ⋯ ℓ / ) sets. (ii) for each term, a loop on the index i. (iii) for a given i, the partial transformations (Eqs (20) ( 21)) which involve a multidimensional loop on the indices uj and kj.

IV.b) Operator action on a quantum state:

As for the calculation of the overlap (or scalar product, Eq (18)), the action of an operator, 𝑶 < , on a quantum state, |𝜓⟩, can be split and distributed on each Smolyak term:

𝑶 < |𝜓⟩ = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • 𝑶 < @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B 6B/7)?|𝓵|?6 (22)
Furthermore, from Eqs ( 18) and ( 22), the computation of J𝜑@𝑶 < @𝜓B, can be obtained as follows:

J𝜑@𝑶 < @𝜓B = - 𝐷 ℓ " ,ℓ $ ,…ℓ # • J𝜑 ℓ " ⋯ℓ ! ⋯ℓ # @𝑶 < @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B 6B/7)?|𝓵|?6 (23)
Eqs ( 22) and ( 23) need further explanation. Indeed, for a given Smolyak term, (ℓ ) ⋯ ℓ , ⋯ ℓ / ), the operator action, 𝑶 < @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B, is performed only with the nD-basis functions of the small direct product, 𝑩 ℓ " ) ⨂ ⋯ 𝑩 ℓ ! , ⋯ ⨂𝑩 ℓ # / and the result is not projected on the full nD-basis set (NB basis functions) but only on the basis functions of the considered Smolyak term (where

∏ 𝑛𝑏 ℓ ! , / ,()
is the number of basis functions). In a way, it is equivalent to computing a small operator matrix (∏ 𝑛𝑏 ℓ ! , / ,()

× ∏ 𝑛𝑏 ℓ ! , / ,() ) associated to the Smolyak term. Each element of this matrix is €𝑏 + " ) … 𝑏 + ! , … 𝑏 + # / •𝑶 < •𝑏 + " . ) … 𝑏 + ! . , … 𝑏 + # .
/ ', where

𝑏 + " ) … 𝑏 + ! , … 𝑏 + # / and 𝑏 + " . ) … 𝑏 + ! . , … 𝑏 + # . / are both present in 𝑩 ℓ " ) ⨂ ⋯ 𝑩 ℓ ! , ⋯ ⨂𝑩 ℓ # / .
Furthermore, when there is no Smolyak term which include both basis functions (𝑏

+ " ) … 𝑏 + ! , … 𝑏 + # / and 𝑏 + " . ) … 𝑏 + ! . , … 𝑏 + # . / ), the coupling element, €𝑏 + " ) … 𝑏 + ! , … 𝑏 + # / •𝑶 < •𝑏 + " . ) … 𝑏 + ! .
, … 𝑏 + # . / ', is not calculated at all, hence its value is zero.

In other words, some coupling elements are missing. However, when the Smolyak parameter, L, increases, some Smolyak terms with will be added and the overall operator matrix will include more coupling terms. Therefore, when L is sufficiently large all necessary coupling elements will be present and it will give converged results.

For instance, taking the 2D-example given in Table III:

-The basis functions, 𝑏 # ) • 𝑏 ) # and 𝑏 H ) • 𝑏 ) # are both present in the following Smolyak terms, (ℓ ) , ℓ # ): (2,0), (3,0) and (2,1). Therefore, J𝑏 # ) • 𝑏 ) # @𝑶 < @𝑏 H ) • 𝑏 ) # B is calculated as follows:

𝐷 #,> J𝑏 # ) • 𝑏 ) # @𝑶 < @𝑏 H ) • 𝑏 ) # B + 𝐷 H,> J𝑏 # ) • 𝑏 ) # @𝑶 < @𝑏 H ) • 𝑏 ) # B + 𝐷 #,) J𝑏 # ) • 𝑏 ) # @𝑶 < @𝑏 H ) • 𝑏 ) # B with 𝐷 #,> +𝐷 H,> +𝐷 #,) =1, so that the sum is equal to J𝑏 # ) • 𝑏 ) # @𝑶 < @𝑏 H ) • 𝑏 ) # B.
-For the basis functions, 𝑏 H ) • 𝑏 ) # and 𝑏 ) )

• 𝑏 H # , there is no Smolyak term where both basis functions are present. Therefore, J𝑏 H ) • 𝑏 ) # @𝑶 < @𝑏 ) )

• 𝑏 H # B is not calculated and never used for the given L (L=2).

However, when L increases, this element will be calculated.

This approximation is essential in our implementation. Indeed, it gives a faster numerical implementation and it enables a trivial parallelization.

IV.c) Numerical aspects of Hamiltonian action:

Indeed, an important point when computing the action of an operator for large systems, is to perform the numerical operations as fast as possible and to avoid large intermediate vectors or matrices.

In particular, in the case of the second scheme, it is important to avoid a full Smolyak representation because the number of elements, 𝑁𝐸 -(basis functions) or 𝑁𝐸 S (grid points) can be very large (see Table IV). Therefore, a quantum state is better stored on the n-dimensional basis set with the compact representation (NB basis functions). The numerical procedure to compute |𝜑⟩ = 𝑶 < |𝜓⟩ is split on each Smolyak term (for a given (ℓ ) ⋯ ℓ , ⋯ ℓ / )) as follows:

(i) From the |𝜓⟩ on the compact basis set representation (𝝍 -" ⋯-# (: ), Eq. ( 3)), extract all basis functions of one term of the Smolyak sum on the basis representation, 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ). The size of the vector, (∏ 𝑛𝑏 ℓ ! , / ,() ), is small compared to NB or 𝑁𝐸 -.

(ii) For the element, 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ), compute the operator action, 𝑶 < @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B. For this step to be efficient, it is important to use the right representation (grid or basis). For instance, when 𝑶 < is known in terms of the coordinates, Q, or when it is a local operator in Q (a function of Q), like the potential, 𝜓

𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ) must be transformed into 𝜓 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (: ) (Eq ( 21
)
) to calculate the action of 𝑶 < . Then, the result obtained in Q (𝜑 (20).

𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (: )) must be transformed into the basis set representation to get 𝜑 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ) using Eq
(iii) The result, 𝜑

𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (:
), is added to the compact basis set representation of |𝜑⟩. In a way, it is the reverse operation of step (i).

Remarks: These operations on the loop are independent from each other. However, the mapping used in steps (i) and (iii) leads to the coupling of the different 𝑶 < @𝜓 ℓ " ⋯ℓ ! ⋯ℓ # B calculations. In each of these three steps, |𝜓⟩ and |𝜑⟩ are never stored on the full Smolyak representation, but only on the compact basis representation. During the calculation, only manipulations and operations on small direct products are performed.

The action of the kinetic energy operator, KEO, on a quantum state requires special care, in particular when a numerical approach is used, the memory requirement may become huge. Indeed, the number of terms associated to the second order differential operator grows as d 2 /2 and each term is associated to a function in Q (the functions 𝑓 # $% (𝑸), 𝑓 ) $ (𝑸) and 𝑉 ;< (𝑸) in Eq. ( 25) or the 𝐺 $% (𝑸)

and 𝐽(𝑸) in Eq. ( 24)), which have to be stored on the full Smolyak grid (𝑁𝐸 S elements).

𝑇 " * #𝑸, 𝜕 𝑸 ' = - ℏ # 2 : 𝐽(𝑸) 𝜌(𝑸) - 1 𝐽(𝑸) 𝜕. 𝜕𝑄 $ 𝐽(𝑸) -𝐺 $% (𝑸) 𝜕. 𝜕𝑄 % : 𝜌(𝑸) 𝐽(𝑸) & %() & $() = -𝑓 # $% (𝑸) 𝜕 # . 𝜕𝑄 $ 𝜕𝑄 % & $T% + -𝑓 ) $ (𝑸) 𝜕. 𝜕𝑄 $ & $ + 𝑉 ;< (𝑸) and 𝑑𝜏 * = 𝜌(𝑸)𝑑𝑄 ) ⋯ 𝑑𝑄 & (24) (25) (26)
For instance, using the grid size from Table IV with L=8, the number of grid points is about half a billion, therefore one needs about one terabytes of memory, since with 21 degrees of freedom, 253 terms are present in Eq. ( 25). To avoid the requirement a large amount of memory, we recently proposed a strategy in which each term is recomputed when needed. [START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF] However, to be efficient, we use equation Eq. ( 24)) instead of Eq. ( 25) and only one Smolyak term (for a given (ℓ ) ⋯ ℓ , ⋯ ℓ / )) is considered in this scheme.

The numerical procedure to compute |𝜑⟩ = 𝑯 < |𝜓⟩ for a given Smolyak term defined with (ℓ ) ⋯ ℓ , ⋯ ℓ / ) can be summarized as follows:

(i) From the |𝜓⟩ on the compact basis set representation (𝝍 -" ⋯-# (𝑘 ) , … , 𝑘 , , … 𝑘 / ), Eq. ( 3)), extracting one term of the Smolyak sum in the basis representation, 𝜓 (21).

𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ). (ii) Transform 𝜓 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ) into 𝜓 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (: ) using Eq
(iii) Compute the action of the potential, V(Q), on 𝜓

𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (: ) to get 𝜑 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (:
). In this case, for all grid points (𝑢 ) , 𝑢 # … 𝑢 / ) of the small direct product, the action of the potential is: US / , on the grid. [START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF] Here, the operator US 0 on the grid and sum on r. (f) multiplication by -ℏ # / ‡4𝐽(𝑸)𝜌(𝑸) and update the vector 𝜑

𝜑 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (𝑢 ) , 𝑢 # … 𝑢 / ) = 𝑉(𝑸 1 " ) , 𝑸 1 $ # … 𝑸 1 # / ) • 𝜓 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (𝑢 ) , 𝑢 # … 𝑢 / ) (
𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # .
(vi) Then, the result obtained on the grid, 𝜑 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # (: ), must be transformed on the basis set to get 𝜑 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ) using Eq (20).

(vii) The result, 𝜑 𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # # (: ), is added to the compact basis set representation of |𝜑⟩.

In our previous publication, [START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF] this scheme was used on the full grid (with 𝑁𝐸 S elements) and 2d+4 intermediate vectors were required on the full grid. However, in the present implementation specific to the second Smolyak scheme, we need one and only one vector on the full grid i.e. the potential (this can be avoided if the potential is recalculated). Indeed, the sizes of intermediate vectors required in the Hamiltonian action scheme described above are of small size (∏ 𝑛𝑞 ℓ ! , / ,() ) compared to 𝑁𝐸 S , because the Smolyak terms are treated one after each other. Furthermore, one also needs to store a vector of integers (size 𝑁𝐸 -) for the mapping between the basis representation of one Smolyak term, 𝜓

𝑩 ℓ " " ⋯𝑩 ℓ ! ! ⋯𝑩 ℓ # #
, and the compact basis representation, 𝝍 -" ⋯-# .

Furthermore, these steps are performed on a loop which can be easily parallelized, because all operations are completely independent, except the first and the last one with the use of a mapping between the Smolyak scheme and the compact basis representation. In the current implementation, the loop on the Smolyak terms parallelized with OpenMP. However, due to the different sizes of the Smolyak terms, it is not easy to distribute the working load in a balanced way. Therefore, several aspects have to be considered: (i) The order of the Smolyak terms is important. In particular, it is not efficient to sort the Smolyak terms in an ascending order in |𝓵| (|𝓵| = ∑ ℓ , ), i.e. the graded lexicographic order, because all small terms would be at the beginning of the list. In our implementation, we use a simple lexicographic order with constraints where the small or large terms are spread everywhere in the list of Smolyak terms. (ii) From the previous list, the Smolyak terms are split among core equally. (iii) The Smolyak terms are split in c parts, such that the number of grid points are almost identical between each part.

In Figure 2, we present the speed-up (ratio between the real execution time with one core and c cores) obtained with OpenMP parallelization for Malonaldehyde (d=n=21) for the calculation of the Hamiltonian action on two quantum states with LB=5 and LG=6 (for the computational details, see next section). For 32 cores, the speed-up is 14 when the list of Smolyak terms is split among cores equally (point ii) and it increases to 20 when the list is split with respect to the grid point number (point iii). Of course, it is not perfect, nevertheless, it enables to get the results much faster than with a single core calculation.

VI. APPLICATION TO MALONALDEHYDE:

The aim of this study is to show the feasibility of the full dimensionality quantum calculation with the potential of Mizukami, Habershon and Tew 4 using the 21D-basis set and grid with the help of the second Smolyak scheme. To be able to perform the calculation of some vibrational levels, it is essential to select a good set of coordinates adapted to the isomerization process, which will enable us to use a compact basis set.

VI.a) Coordinates for malonaldehyde:

To simulate low excitation vibrational states (implied with fundamental transition), the curvilinear normal coordinates are well adapted since the quadratic coupling between the normal coordinates (curvilinear or rectilinear) are removed and the residual anharmonic coupling between the coordinates are reduced, with respect to rectilinear normal modes. However, for this proton transfer process, two different sets of curvilinear normal coordinates are required, one for each minimum. Therefore, one needs a set of coordinates which allows to switch from one minimum to the other one.

The coordinates associated to the reaction path Hamiltonian (RPH), [START_REF] Miller | Reaction path Hamiltonian for polyatomic molecules[END_REF] the reaction surface Hamiltonian 58 or similar approaches lead to using the local normal coordinates along a path [START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF] (or surface) and so that this kind of coordinates is well adapted for studying systems with several minima. In these methods, one needs to calculate a path (minimum energy path, steepest decent path ...) as well as the hessian along this path. However, in a way, this path is arbitrary and not necessarily the optimal one. Therefore, we want to use a simplified version of the RPH, which will not prevent performing an exact calculation in which the path and the hessian are obtained from an interpolation between the two minima. The curvilinear coordinates (Figure 3) adapted to this path are obtained after several coordinate transformations:

(i) Transformation to/from polyspherical coordinates, 𝑸 <WXY . This transformation is associated to the relation between Cartesian coordinates and 𝑸 <WXY . The eight vectors associated to 𝑸 <WXY are defined in the following order: the first vector, 𝑅 Š⃗ ) ()) , is the vector between the two oxygen atoms.

This vector is along the z-axis of the body fixed (BF) frame. 𝑅 Š⃗ # ()) vector between the center of mass of O-O and the central carbon atom. This vector is in the xy-plan of the BF frame. 𝑅 Š⃗ H ()) , vector between the central carbon atom and its bounded hydrogen (in the BF frame). 𝑅 Š⃗ I ()) , vector between the center of mass of O-O and the proton. The vector coordinates are cartesian in the BF frame and its z coordinate (z4) can be associated to a proton transfer coordinate. To avoid singularity during the dynamics, two sub-systems need to be defined for the two OCH fragments.

In the first subsystem (S1,1) on the left on the Figure 3 (ii) 1D-tranformation,𝑸 )[ : [START_REF] Sarka | Rovibrational energy levels of the F -(H 2 O) and F -(D 2 O) complexes[END_REF] All polyspherical coordinates will be used in linear combinations (the next two transformations) with range ]-∞, ∞[ since Harmonic oscillator basis set will be used. Therefore, for large grids, some polyspherical coordinates can be out of their mathematical range. In particular, it is the case for the coordinates 𝑢 , (Z) or 𝑢𝛽 , (Z) where the mathematical range is ]-1,1[. To avoid numerical trouble, those are transformed, so that the transformed coordinates are 𝑡𝑢 , (Z) = Arctanh (𝑢 , (Z) ) and with a mathematical range ]-∞, ∞[. For the other coordinates, the norms or the azimuthal angles, should be transformed but, numerically, it was not necessary. However, at the two reference geometries, the value of some Euler angles (𝛼 (),)) , 𝛾 (),)) and 𝛼 (#,)) ) is p. Therefore, to simplify the next transformations, these three angles are shifted by p.

(iii) Linear combinations, 𝑸 %YE : with this transformation, the previous coordinates (including the Euler angles) associated to the sub-systems, S1,1 and S2,1 are combined to form symmetrized coordinates. For instance, for the norm of the vectors 𝑅 Š⃗ ) (),)) and 𝑅 Š⃗ ) (#,)) , the two corresponding symmetrized coordinates are:

) # ˜𝑅) (),)) ± 𝑅 ) (#,)) š.
(iv) Approximate RPH:

Instead of an exact minimum energy path along the coordinate z4 (z-component of 𝑅 Š⃗

I ()) )
associated to the proton transfer, the approximate path is defined from the two equivalent minima, m1 and m2, as follows:

The Cartesian coordinates, 𝑸 \F$] , at the two minima (𝑸 \F$],E " and 𝑸 \F$],E $ ) are transformed with the help of TNUM into their 𝑸 %YE values (𝑸 %YE,E " and 𝑸 %YE,E $ ). Next, we perform a linear interpolation between (𝑸 %YE,E " and 𝑸 %YE,E $ ) in term of tz4, where tz4 is a function of z4 (𝑡𝑧 I = 𝑧 > • tanh (𝑧 I 𝑧 > ⁄ ), where z0=1.2 bohr). Except for z4, the linearly interpolated 𝑸 %YE,<F]^( 𝑧 I ) is given by:

𝑸 %YE,<F]^( 𝑧 I ) = #𝑡𝑧 I -𝑡𝑧 I,E " ' #𝑡𝑧 I,E $ -𝑡𝑧 I,E " ' • 𝑸 %YE,E $ + #𝑡𝑧 I -𝑡𝑧 I,E $ ' #𝑡𝑧 I,E " -𝑡𝑧 I,E $ ' • 𝑸 %YE,E " (27)
where, 𝑡𝑧 I,E " and 𝑡𝑧 I,E $ are, respectively, the 𝑡𝑧 I values at the two equivalent minima, m1 and m2.

It is important to note, that, in the linear interpolation, the use of 𝑡𝑧 I instead of 𝑧 I is essential to avoid unphysical geometries. Indeed, during the quantum calculation, the range of 𝑧 I can be so large that unphysical 𝑸 %YE,<F]^( 𝑧 I ) values (for instance negative distances ...) can be obtained from a linear interpolation directly performed in 𝑧 I . Now, with 𝑡𝑧 I , this trouble can be avoided, since the 𝑡𝑧 I range is ]-𝑧 > , 𝑧 > [ (and z0=1.2 bohr). The value of 𝑧 > has been chosen so as to obtain a smooth 𝑡𝑧 I function and 𝑸 %YE coordinates without unphysical ranges.

This procedure is used for the elements (except the ones implying 𝑧 I ) of the hessian matrix with respect to 𝑸 %YE . We obtain a hessian matrix, 𝒉 %YE,<F]^( 𝑧 I ), along the 𝑧 I . We also use 𝑡𝑧 I for the linear interpolation with the same parameter 𝑧 > , but in this case, it is to avoid unwanted imaginary frequencies.

At a given value of 𝑧 I , a quadratic 20D-Hamiltonian, parametrized by 𝑧 I , is obtained using, 𝑸 %YE,<F]^( 𝑧 I ), 𝒉 %YE,<F]^( 𝑧 I ) and 𝑮 %YE,<F]^( 𝑧 I ), where 𝑮 %YE,<F]^( 𝑧 I ) is the metric tensor at

𝑸 %YE,<F]^( 𝑧 I ). 𝑯 K𝑸 /12 #>[ (𝑧 I ) = -- ℏ # 2 𝐺 %YE,<F],
,D

(𝑧 I ) 𝜕 # 𝜕Δ𝑄 %YE , 𝜕Δ𝑄 %YE D #) ,,D(# + 1 2 ℎ %YE,<F],
,D

(𝑧 I )Δ𝑄 %YE , Δ𝑄 %YE D (28)
where,

Δ𝑄 %YE , is ˜𝑄%YE , -𝑄 %YE,<F], (𝑧 I )š.
Next, this coupled quadratic Hamiltonian along the path is transformed using a standard procedure to get the curvilinear normal coordinates to yield 20 uncoupled quadratic Hamiltonians in 𝑄 _`a , coordinate. In other words, each 𝑄 _`a , is expressed as a linear combination of the Δ𝑄 %YE D

.

𝑯 𝑸 345 #>[ (𝑧 I ) = - ℏ # • 𝜔 , (𝑧 I ) 2 - 𝜕 # 𝜕𝑄 _`a , # + 𝑄 _`a , # ¡ #) ,(# (29)
The 𝜔 , (𝑧 I ) are harmonic angular frequencies along the proton transfer path (z4). For some systems, like malonaldehyde, angular frequencies may present a sharp avoided crossing. Formally, it is not a problem, but numerically, it is difficult to deal with this kind of situation. In fact, it is similar to the sharp avoided crossings encountered in the Born-Oppenheimer potential energy curves and in our case, we also use some kind of quasi-diabatic curves to overcome this problem.

For malonaldehyde, two CH stretching modes present this sharp avoided crossing for z4=0 (w19(0)=3121.6 cm -1 and w20(0)=3122.7 cm -1 , see Figure 4) and at this point the energy difference between these two angular frequencies is 1.1 cm -1 (see Figure 4b). Like for Born-Oppenheimer potential energy curves, to avoid the numerical problems due to large non-adiabatic couplings, we use quasi-diabatic transformation. However, in our case it is very easy to bypass the sharp avoided crossing between w19 and w20. Indeed, starting from negative z4 value, we just need to swap the w19 and w20 values and also the normal coordinates associated to these two frequencies when z4 is positive. However, an important aspect must be considered here: the grid points along z4 should not be present in a small domain around zero, otherwise, it would be impossible to get smooth quasidiabatic frequencies. In particular, the point where z4=0 must be avoided. This aspect will have consequences for the basis set and grid associated to the z4 coordinate.

Remark: The quadratic Hamiltonians (the coupled or the uncoupled ones) are not used to perform the full quantum calculation, but they are used only to get the coefficients along the path (with smooth variation) of the linear combinations.

All the transformations (schematized in Figure 3) can be used to numerically compute the Cartesian coordinates from the active ones (here the RPH ones) and the full metric tensor (in 21D) used in the kinetic energy operator. The ordering of the coordinates, Q, used in the quantum calculation is the following: Q 1 =z4 and 𝑄 , = 𝑄 _`a , (i=2...21).

VI.b) Basis set and grid:

With this set of coordinates, we are able to choose the primitive basis sets (𝑩 , ), their corresponding grids (𝑮 , ) and their parameters (𝑛𝑏 ℓ ! , and 𝑛𝑞 ℓ ! , ) associated to the Smolyak scheme. For this system, all the primitive basis sets are in 1D (with one coordinate) and therefore n=d. Furthermore, two kinds of coordinates have to be considered: (i) The one associated to the proton transfer, Q 1 (or z4). Then, the corresponding basis set, 𝑩 ) , needs to be selected with care since it must describe a double-minimum. Here, we have chosen a scaled Harmonic Oscillator basis set, so that the unnormalized i th basis function is given by:

𝑏 , ) (𝑄 ) ) = 𝐻 ,B) (𝑥)𝐸𝑥𝑝 ¥- 𝑥 # 2 ¦ (30)
with 𝑥 = 𝑠. 𝑄 ) and where Hi-1 is the Hermite polynomial of degree i-1.

The value of the scaling parameter, s, is 3.8 for the normal malonaldehyde (with hydrogen atoms) and 4.5 for deuterated malonaldehyde. These values have been obtained by minimizing the zeropoint energy (the first energy level) in 21D using a small basis set (B0, see below).

The grid points associated to this basis set are the usual scaled gauss-Hermite quadrature grid points.

For this basis set, associated to z4, the integer sequences (𝑛𝑏 ℓ "

) and 𝑛𝑞 ℓ "

) ) need to be defined with extreme care. Indeed, the point, z4=0, must not be present in 𝑛𝑞 ℓ " ) grid points for all value of ℓ ) .

Therefore, 𝑛𝑞 ℓ " ) must be even numbers. Furthermore, for a sufficient large value of ℓ ) (about 7 or 8), the number a basis functions must be around 30 in order to be able to describe the eigenfunctions of the double minima potential. Therefore, for both 𝑛𝑏 ℓ " ) and 𝑛𝑞 ℓ " ) , we choose the following even sequence: 𝑛𝑏 ℓ " ) = 𝑛𝑞 ℓ " ) = 2 + 4ℓ ) . which gives, respectively, 30 and 34 basis functions and grid points for ℓ ) equal to 7 and 8.

(ii) the other coordinates, Q 2 to Q [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF] , where the corresponding basis sets are the unscaled Harmonic Oscillator basis set. Those basis sets do not need to be scaled, because the corresponding coordinates are the curvilinear normal coordinates along the path. Furthermore, for all those primitive basis sets, we used the usual Gauss Hermite quadrature points. For a given i, two integer sequences for both 𝑛𝑏 ℓ ! , and 𝑛𝑞 ℓ ! , are used:

𝑛𝑏 ℓ ! , = 𝑛𝑞 ℓ ! , = 1 + ℓ , or 𝑛𝑏 ℓ ! , = 𝑛𝑞 ℓ ! , = 1 + 2ℓ , .
With the parameters described above, three 21D-basis sets (or grid) are defined (i) an extremely small basis set (named B0) used to optimize the scaling parameter of 𝑩 ) . This 21D-basis set is expressed as a direct-product of a 1D-basis-set, 𝑩 ) (with 𝑛𝑏 ) basis functions and 𝑛𝑞 ) grid points) and a 20D-Smolyak basis set where all sequences are 1 + ℓ , for the basis sets or grids. For the optimization of the scaling parameter, the parameters associated to the basis set size are: 𝑛𝑏 ) = 8 and 𝑛𝑞 ) = 12, LB=0 and LG=2 for the normal malonaldehyde and 𝑛𝑏 ) = 12 and 𝑛𝑞 ) = 16, LB=1 and LG=3 for the mono-deuterated malonaldehyde. (ii) a small basis set (named B1), with a sequence equal to 2 + 4ℓ ) for 𝑩 ℓ " ) and 𝑮 ℓ "

) and where all other sequences are 1 + ℓ , for the basis sets or grids associated to the coordinates, Q 2 to Q 21 . (iii) a large basis set (named B2), with a sequence equal to 2 + 4ℓ ) for 𝑩 ℓ " ) and 𝑮 ℓ "

) and where some sequences are 1 + ℓ , and others are 1 + 2ℓ , for the basis sets or grids associated to the coordinates, Q 2 to Q 21 . The selection of the basis sets with the 1 + 2ℓ , sequence (4 for the normal malonaldehyde and 3 for deuterated one) is performed with the analysis of the eigenfunctions obtained with the small basis set, B1, with LB=5 (see below).

VI.b) Ground state and tunneling splitting

For this system with 21 degrees of freedom, the size of the basis set is too large (several millions of basis functions), so that the construction of the full Hamiltonian matrix is not possible. Therefore, one has to use an iterative diagonalization procedure to get the eigenstates. For the present study, the block Davidson technics is used, in which the Hamiltonian action is performed as described in section IV.c (with an on-the-fly computation of the KEO).

The first calculations were done for the normal malonaldehyde with the basis set B1 from LB=2 to LB=8 and the zero-point energy (ZPE) shows a nice convergence as function of LB (see Table V and Figure 5). However, the ZPE difference between the calculations at LB=7 and LB=8 is 4.2 cm -1 . Normally, it will require calculations at LB=9 and probably at LB=10, but with such parameters, the calculations will too large and too time-consuming on our computer. One difficulty with this basis set is that, the numbers of 1D-basis functions along all coordinates Q 2 to Q 21 are too small (𝑛𝑏 b , = 9

for LB=8). To improve that, one can use other sequences 1 + 2ℓ , . However, doing that for basis sets, 𝑩 , (i=2 to 21) will give extremely large grids. Thus, we must use some kind of selection.

Here, we use the following strategy: (i) A calculation with the small basis set B1 at LB=5. To get the first two levels, the calculation with a Davidson procedure takes about 3 hours and a half on a machine with 12 cores. (ii) Analysis of the average population on the two states of the diagonal elements of the 1Dreduced density matrix projected on the basis functions of 𝑩 , (i=2 to 21). Then, for all basis sets, 𝑩 , , with the population of the last basis functions (𝑏 c , ) larger than 0.3 10 -4 , the sequence 1 + 2ℓ , will be used for the basis set B2. Here, the basis sets (𝑩 d , 𝑩 )I , 𝑩 )b and 𝑩 #) ) are in this case.

As expected, with this new basis set, B2, the convergence in term of LB is improved (Figure 5), now the ZPE reaches 14719.2 cm -1 which is close to the DMC value, 4 14716.6 cm -1 (without extrapolation). Furthermore, the ZPE difference between the calculations at LB=7 and LB=8 is 2.1 cm -1 showing that the basis B2 gives a better description than with the B1 basis set.

From the calculation with both basis sets (B1 and B2), the tunneling splitting converges around 21.7 cm -1 . The agreement with the experimental value, 21.6 cm -1 , is remarkable. However, the difference between the DMC tunneling splitting (21.0 cm -1 ) is larger than expected.

(a) (b) For the mono-deuterated malonaldehyde, we have used the same procedure. The selection of the 𝑩 , to be used with the sequence 1 + 2ℓ , gives 3 1D-basis sets: 𝑩 d , 𝑩 )H and 𝑩 )e . The ZPE convergence (Table VI and Figure 6) with the basis set B2 as function of LB is not as efficient as the one with the normal malonaldehyde. For LB=8, we get 14004.8 cm -1 while the DMC value 4 is 13998.9 cm -1 (without extrapolation). However, the tunneling splitting for LB=8 is 2.9 cm -1 , which is very close to the DMC one (3.1 cm -1 ).

(a) (b)

VII. CONCLUSION AND PERSPECTIVES:

With the combination of several numerical approaches used in the present study, we were able to calculate the tunneling splitting of the normal isotoplogue of malonaldehyde, 21.7±0.3 cm -1 and of the mono-deuterated one, 2.9±0.1 cm -1 , in full dimensionality i.e. 21 degrees of freedom. The 0.1 incertitude is estimated from the difference between the small (B1) and the large basis set (B2) at LB=8. Furthermore, those values are in very good agreement (around 0.1 cm -1 difference) with the experimental ones, 21.6 cm -1 and 2.9 cm -1 for the normal isotopologue [START_REF] Turner | Microwave spectroscopic study of malonaldehyde. 4. Vibration-rotation interaction in parent species[END_REF][START_REF] Varella | Real-time observation of intramolecular proton transfer in the electronic ground state of chloromalonaldehyde: an ab initio study of time-resolved photoelectron spectra[END_REF] and the mono-deuterated one, [START_REF] Baughcum | Microwave spectroscopic study of malonaldehyde. 3. Vibration-rotation interaction and one-dimensional model for proton tunneling[END_REF] respectively. This shows the excellent quality of the potential of Mizukami et al. [START_REF] Mizukami | A compact and accurate semi-global potential energy surface for malonaldehyde from constrained least squares regression[END_REF] However, the calculations are not fully converged (see the right panels of Figure 5 and Figure 6) and at least LB=9 is required. Unfortunately, the memory requirement is too large for our computer. For instance, with LB=9 with the basis B2 for the normal isotopologue, the size of the basis set is not too large, (NB≃102.10 6 ) but the number of elements for the second Smolyak scheme 𝑁𝐸 -for the basis set and 𝑁𝐸 S for the grid are about 75.10 9 , hence, the need of 278 GB just for the mapping, see below.

Anyway, the present calculations in full dimensionality (21D) were possible, only because we have managed to combine several numerical procedures: (i) a well-adapted set of coordinates based on an approximated reaction path coordinate. The path is simpler than an intrinsic reaction path or minimum reaction path, but nevertheless, describes well the path between the two minima. However, we believe that an instanton path [START_REF] Benderskii | Tunneling splittings in vibrational spectra of non-rigid molecules. X. Reaction path Hamiltonian as zero-order approximation[END_REF][START_REF] Mil | Practical implementation of the instanton theory for the ground-state tunneling splitting[END_REF] may improve this description. Whatever the path, this set of coordinates implies the use of numerical but exact kinetic energy operators such as the one implemented in TNUM. [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF] (ii) on-the-fly: the different terms of the Hamiltonian (or other operators) are not necessarily stored on the full grid. Instead, they are recalculated on-the-fly, i.e. each time they are required. In the present study, we used this procedure for the metric tensor [START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF] only, the potential being stored in memory, although this can be avoided too.

(iii) the second Smolyak scheme. To perform such a large calculation with 21 degrees of freedom, we need an efficient parallelization (realized with OpenMP) and we need to avoid, as much as possible, large vectors. More precisely, the vectors of the full grid size, 𝑁𝐸 S , can be avoided completely. For instance, the operators are recalculated (see point ii) and furthermore, the wave functions are never stored on the full grid but rather, on the basis set in the compact representation (size NB), which cannot be avoided. However, NB is not that large and therefore the wave functions can be stored easily. However, we need to store a large vector of integers (length: 𝑁𝐸 -) for the mapping between the 𝑁𝐸 -elements (the basis functions) of the Smolyak scheme and the NB basis functions in the compact representation. This large storage cannot be avoided in the current implementation. All these features, including the coordinates and KEO aspects, are freely available in the fortran code, ELVIBROT-TNUMTANA. [START_REF] Lauvergnat | ElVibRot-Tnum[END_REF] The present implementation leaves room for improvements: (i) parallelization: presently, the loop on the Smolyak terms is parallelized with OpenMP on a single computer node with shared memory, but with the structure of our Smolyak algorithm (see section IV.b), a parallelization on a large number of computer nodes with MPI should be straightforward. Indeed, the number of Smolyak terms can be very large (14307150 for LG=9 with 21 degrees of freedom) and there is almost no communication between threads associated to the calculations on the Smolyak terms. (ii) memory requirement: to push the limit to larger systems or to have a better convergence, the large vector of integers required for the mapping between the basis functions on the Smolyak scheme and the compact representation has to be removed. It is probably the main bottle-neck, but we think it should possible using the relation between the first and the second Smolyak scheme. (iii) Presently, the Smolyak directproducts are selected in terms of excitations which works well when the couplings between the modes are of the same order of magnitude. However, it will be interesting to have several procedures for the selection of the Smolyak terms either automatically to add the important one only or using the physics of the system. For instance, we can imagine a selection with several L parameters, the main one (LG) used for the main part of the system and another one, L1, for a weakly coupled environment. The value of L1 will remain small (2 or 3) and in any case smaller than LG (around 9 or 10). With this approach, we should be able to study systems with a large number of degrees of freedom (much larger than 21) providing the system can be split into a main part and a weakly coupled environment.

where k and @𝑏 𝒌 )B/ B represents, respectively, a multidimensional index (𝑘 ) , … , 𝑘 , , … 𝑘 / ), the multidimensional basis function associated to k. The 𝜒 -")# (𝒌) and 𝜓 -")# (𝒌) are the coefficients of the basis expansion of |𝜒⟩ and |𝜓⟩, respectively.

Let us show that, in the case of the second Smolyak scheme, the scalar product as defined by Eq. ( 18) coincides with the usual scalar product.

In the second Smolyak scheme, the scalar product given by Eq. ( 18) can be compactly rewritten as follows ⟨𝜒|𝜓⟩ = -𝐷 𝓵 • J𝜒 𝓵 @𝜓 𝓵 B 𝓵 (32) where 𝓵 is the multidimensional index (ℓ ) , ℓ # , … ℓ / ), 𝜒 𝓵 and 𝜓 𝓵 are, respectively, the contributions on the Smolyak term associated to 𝓵 of 𝜒 and 𝜓.

Then, for a given 𝓵, the Smolyak contribution of 𝜓 and 𝜒 are given by: @𝜒 𝓵 B = -𝜒 -")# (𝒌′) • 𝜀 𝒌g 𝓵 • @𝑏 𝒌g )B/ B 𝒌g @𝜓 𝓵 B = ∑ 𝜓 -")# (𝒌) • 𝜀 𝒌 𝓵 • @𝑏 𝒌 )B/ B 𝒌 (33) where 𝜀 𝒌 𝓵 is the short notation for 𝜀 + " ⋯+ ! ⋯+ # ℓ " ⋯ℓ ! ⋯ℓ # . In Eq (33), the sums on k and k' run over all NB basis functions. Therefore, the coefficients, 𝜓 -")# (𝒌) are identical for all Smolyak terms (the same holds true for the 𝜒 -")# (𝒌′)).

Then by inserting Eq (33) into Eq (32), one gets: ⟨𝜒|𝜓⟩ = -𝐷 𝓵 • -𝜒 -")# (𝒌 g ) * • 𝜓 -")# (𝒌) • 𝜀 𝒌 𝓵 • 𝜀 𝒌 . 𝓵 J𝑏 𝒌 . )B/ @𝑏 𝒌 )B/ B 𝒌 . ,𝒌 𝓵 = -𝐷 𝓵 • -𝜒 -")# (𝒌′) * • 𝜓 -")# (𝒌) • 𝜀 𝒌 𝓵 • 𝜀 𝒌g 𝓵 • 𝛿 𝒌g,𝒌 𝒌g,𝒌 𝓵 = -𝐷 𝓵 • -𝜒 -")# (𝒌) * • 𝜓 -")# (𝒌) • 𝜀 𝒌 𝓵 • 𝜀 𝒌 𝓵 𝒌 𝓵 (34) Then, the two sums (in 𝓵 and 𝒌) can be swapped and since 𝜀 𝒌 𝓵 • 𝜀 𝒌 𝓵 = 𝜀 𝒌 𝓵 (because 𝜀 𝒌 𝓵 is equal to 0 or 1), so that Eq (34) gives: ⟨𝜒|𝜓⟩ = --𝐷 𝓵 • 𝜒 -")# (𝒌) * • 𝜓 -")# (𝒌) • 𝜀 𝒌 𝓵 𝓵 𝒌 = -𝜒 -")# (𝒌) * • 𝜓 -")# (𝒌) -𝐷 𝓵 • 𝜀 𝒌 𝓵 𝓵 𝒌 (35) Then, using the fact that the last sum in 𝓵 is equal to one (using Eq (14)), Eq (35) comes down to the usual scalar product expression: ⟨𝜒|𝜓⟩ = -𝜒 -")# (𝒌) * • 𝜓 -")# (𝒌) 𝒌 (36)

Figure 1 :

 1 Figure 1 : scheme of proton transfer isomerization of malonaldehyde.

  iv) Compute of 𝐺 $% (𝑸), 𝜌(𝑸) and 𝐽(𝑸) on the ∏ 𝑛𝑞 ℓ ! , / ,() grid points for the (ℓ ) ⋯ ℓ , ⋯ ℓ / ) Smolyak term and only this one. Since the ∏ 𝑛𝑞 ℓ ! , / ,() is small compare to 𝑁𝐸 S all these functions can be stored in memory without difficulty. (v) KEO action on 𝜓 action of <𝜌(𝑸)/𝐽(𝑸) on 𝜓 𝑮 ℓ " " ⋯𝑮 ℓ ! ! ⋯𝑮 ℓ # # , the computation is done as for the potential. (b) On the previous result, action of all first order derivatives, U.

  U. US / is not local in Q and it represented as a small matrix. (c) On the d previous results (s=1...d), multiplication by 𝐺 $% (𝑸) and sum on s. (d) On the d previous results (r=1...d), multiplication by 𝐽(𝑸). (e) On the d previous results (r=1...d), action of U.

Figure 2 :

 2 Figure 2 : Speed up with respect to the number of cores, c for the computation the Hamiltonian action on two states (with n=d=21, LB=5 and LG=6). The black line represents the ideal speed-up. For the dashed blue line, the list of Smolyak terms is split among thread equally. For the full red line, the list is split in c parts, such that the number of grid points are almost identical between each part.

Figure 3 :

 3 Figure 3 : scheme of the coordinate transformations used in the quantum calculation.

  )) are vectors along the C-O and C-H bounds, respectively. We use the same kind of parametrization for the second subsystem (S2,1) on the right on the Figure3.All eight vectors, 𝑅 Š⃗ , (Z) , (the upper index (s) stands for a particular sub-system or frame) have 3 coordinates: the norm of the vector, 𝑅 , (Z) , (the bond distance for all vectors except for 𝑅 Š⃗ ), the cosine of the polar angle (𝑢 , (Z) ) or the cosine of the second Euler angle (𝑢𝛽(Z) ) and the azimuthal angle (𝜑 , (Z) ) or the first (𝛼 (Z) ) or third Euler (𝛾 (Z) ) angles. As already mentioned, the 𝑅 Š⃗ I ()) vector is described by 3 cartesian coordinates (x4, y4, z4). For the first two vectors (𝑅 Š⃗ ) ()) and 𝑅 Š⃗ # ()) ), their norms and the polar angle of 𝑅 Š⃗ # ()) are internal coordinates, the three other angles are the Euler angles describing the overall rotation.

Figure 4 :

 4 Figure 4 : adiabatic angular frequencies associated to the CH stretching modes along z4. (a) z4 in the [-1.9 , 1.9] range. (b) zoom of the previous frequencies with z4 in the [-0.1 , 0.1] range

Figure 5 :

 5 Figure 5 :for two basis sets, B1 and B2, convergence as function of LB of the normal malonaldehyde: (a) ZPE (b) tunneling splitting, DE.

Figure 6 :

 6 Figure 6 :for two basis sets, B1 and B2, convergence as function of LB of the mono-deuterated malonaldehyde: (a) ZPE (b) tunneling splitting, DE.

  One of us (DL) gratefully acknowledges the hospitality of the School of Physics and the Institute of Condensed Matter and Nanosciences (NAPS) of the Université Catholique de Louvain, where part of the present work was carried out.APPENDIX, Scalar product in the second Smolyak scheme:Let consider two quantum states described by the kets, |𝜒⟩ and |𝜓⟩ which are expanded as follow:|𝜒⟩ = -𝝌 -")# (𝒌) • @𝑏 𝒌 )B/ B 𝒌 |𝜓⟩ = -𝜓 -")# (𝒌) • @𝑏 𝒌 )B/ B 𝒌

Table I .

 I Experimental and computed tunneling splitting of normal and mono-deuterated malonaldehyde.

		Exp.		PES2008 3		PES2014 4
			DMC 3	MCTDH 16 MCTDH 46	DMC 4
	Normal	21.583 a	21.6-22.6	23.8	23.4	21.0
	malonaldehyde
	mono-deutereted	2.915 b	3.0-3.1			3.2
	malonaldehyde

a) refs

  1 ! , ⋯ 𝑸 /2 ! , P and a set of weights H𝜔 1 " , , 𝜔 1 $ , ⋯ 𝜔 1 ! , ⋯ I where the 𝑸 1 !

,

and 𝜔 1 !

Table IV

 IV

					L or LB or LG
		2	3	4	5	6	7	8
	NB	253	2 024 12 650	65 780	296 010	1 184 040	4 292 145
	𝑁𝐸 -	946	14 190 163 185 1 533 939 12 271 512 85 900 584 536 878 650
	𝑁𝐸 S	946	14 190 163 185 1 533 939 12 271 512 85 900 584 536 878 650
	𝑛𝑞 ℓ " & = 1 + ℓ &
	𝑁𝐸 S	1.15 10 9 13. 10 9
	𝑛𝑞 ℓ " & = 2 + ℓ &

. As function of L (or LB or LG): (i) In line 2, the number of basis functions, NB, in a compact representation, (ii) the number of elements of the second Smolyak scheme, 𝑁𝐸 ) (basis functions) in line 3 or 𝑁𝐸 * (grid points) in lines 4 and 5. The Smolyak parameters are: 𝑛 = 21 and for all basis sets, the sequences in ℓ + are 𝑛𝑏 ℓ $ + = 1 + ℓ + . For the grids, the sequences in ℓ + are given in the first column of the table.

  1 + ℓ & for 𝑛𝑏 ℓ ! , and L=2, 𝑁𝐸 S increases from 946 with the sequence 1 + ℓ & for 𝑛𝑞 ℓ ! , to more than 10 9 with the 2 + ℓ & for 𝑛𝑞 ℓ !

, (see Table

Table V

 V

						LB/LG
			2/4	3/5	4/6	5/6	6/7	7/8	8/9
		NB	550	4 554	29 348	156 860	723 580	2 960 100	10 952 370
	Small	𝑁𝐸 *	354 750	3 10 6	28 10 6	28 10 6	196 10 6	1.246 10 9	7.158 10 9
	basis set,
	B1	ZPE	15063.7	14885.2	14793.6	14753.2	14734.9	14725.8	14721.6
		splitting	14.5	8.7	11.7	17.5	20.5	21.3	21.6
		NB	746	7 038	51 050	303 822	1 545 874	6 917 526	27 791 038
	basis set,	𝑁𝐸 *	488 060	5 10 6	43 10 6	43 10 6	329 10 6	2.210 10 9	13.4 10 9
	B2	ZPE	15050.9	14857.4	14774.1	14740.4	14726.9	14721.3	14719.2
		splitting	10.9	12.4	18.1	21.5	22.3	22.0	21.7

. For normal malonaldehyde, zero-point energy (ZPE), tunneling splitting and basis set or grid sizes as function of LB/LG for both basis set B1 and B2.

Table VI

 VI

						LB/LG
			2/4	3/5	4/6	5/6	6/7	7/8	8/9
		NB	550	4 554	29 348	156 860	723 580	2 960 100	10 952 370
	Small	𝑁𝐸 *	354 750	3 10 6	28 10 6	28 10 6	196 10 6	1.246 10 9	7.158 10 9
	basis set,
	B1	ZPE	14337.7	14173.4	14091.3	14041.7	14021.5	14012.0	14706.8
		splitting	51.6	14.7	4.8	0.5	2.0	2.6	2.8
		NB	694	6 344	44 706	259 116	1 286 758	5 630 768	22 160 270
	basis set,	𝑁𝐸 *	451 588	4.6 10 6	39 10 6	39 10 6	290 10 6	1.92 10 9	11.5 10 9
	B2	ZPE	14330.5	14161.3	14081.9	14034.8	14017.0	14008.9	14004.8
		splitting	51.8	13.7	2.6	1.8	2.7	2.9	2.9

. For mono-deutereted malonaldehyde, zero-point energy (ZPE), tunneling splitting and basis set or grid sizes as function of LB/LG for both basis set B1 and B2.

Author to whom correspondence should be addressed. Electronic mail: david.lauvergnat@universite-paris-saclay.fr