Francesco Talotta

David Lauvergnat

Federica Agostini

    Describing the photo-isomerization of a retinal chromophore model with coupled and quantum trajectories

Describing the photo-isomerization of a retinal chromophore model with coupled and quantum trajectories

Francesco Talotta, 1 David Lauvergnat, 1 and Federica Agostini 1

1 Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405, Orsay, France

The exact factorization of the electron-nuclear wavefunction is applied to the study of the photoisomerization of a retinal chromophore model. We describe such an ultrafast nonadiabatic process by analyzing the time-dependent potentials of the theory and by mimicking nuclear dynamics with quantum and coupled trajectories. The time-dependent vector and scalar potentials are the signature of the exact factorization, as they guide nuclear dynamics by encoding the complete electronic dynamics and including excited-state effects. Analysis of the potentials is, thus, essential -when possible -to predict the time-dependent behavior of the system of interest. In this work, we employ the exact time-dependent potentials, available for the numerically-exactly solvable model used here, to propagate quantum nuclear trajectories representing the isomerization reaction of the retinal chromophore. The quantum trajectories are the best possible trajectory-based description of the reaction when using the exact-factorization formalism, and thus allow us to assess the performance of the coupled-trajectory, fully approximate, schemes derived from the exact-factorization equations.

I. INTRODUCTION

Trajectory-based schemes are widely employed to simulate excited-state nonadiabatic dynamics in molecular systems [START_REF] González | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Crespo-Otero | [END_REF][3][4][START_REF] Bonfanti | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Agostini | [END_REF][7]. The main reason for their success is twofold: first, when using trajectories it is possible to perform on-the-fly electronic structure calculations for the visited geometries along the nuclear dynamics [8][9][10], thus avoiding constructing the electronic potential energy surfaces and related quantities over the whole configuration space [START_REF] Bonfanti | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF]11]; second, evolving trajectories rather than the molecular wavefunction allows one to overcome the exponential scaling of the computational cost with the number of degrees of freedom inherent a fully quantum-mechanical calculation. While the outcome of the dynamics simulation is clearly affected by the accuracy of the underlying electronic-structure method [12], the overall evolution of the system, i.e., coupled nuclei and electrons, strongly depends on the forces used to propagate the trajectories and on the electron-nuclear coupling -intended here as the nuclear feedback on electronic dynamics.

Among the trajectory-based methods derived from the molecular time-dependent Schrödinger equation are, for instance, the direct-dynamics variational multiconfigurational Gaussian (DD-vMCG) [START_REF] Worth | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Richings | [END_REF], ab initio multiple spawning (AIMS) [3,15,[START_REF] Curchod | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF], and ab initio multiple cloning (AIMC) [START_REF] Makhov | [END_REF]18]. There, nuclear trajectories are evolved according to variationally-determined (DD-vMCG), or adiabatic (AIMS), or mean-field (AIMC with multiconfigurational Ehrenfest) forces, and are used to guide the basis functions appearing in the expansion of the molecular wavefunction. Electron-nuclear coupling is obtained by solving, alongside the trajectories, the quantum evolution equation for the expansion coefficients. In addition, spawning or cloning procedures aim to improve over time the sampling of the nuclear configuration space. An approximate solution of the molecular time-dependent Schrödinger equation can be obtained via quantum-classical approaches as well, like Ehrenfest dynamics [19,[START_REF] Kirrander | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF] and trajectory surface hopping (TSH) [START_REF] Tully | [END_REF]. The former is a mean-field method, where the classical nuclear force and the electron-nuclear coupling are given by averaging either over the electronic or over the nuclear time-dependent states of the system. The latter evolves classical trajectories according to adiabatic forces, and the electronic state of the system according to a Schrödinger equation where the nuclear trajectory appears as a time-dependent parameter. To account for nonadiabatic effects TSH is supplemented with a hopping procedure allowing the trajectories to change electronic state over time.

Trajectory-based schemes [22][23][24][25][START_REF] Agostini | Quantum chemistry and dynamics of excited states: Methods and applications[END_REF][START_REF] Ha | [END_REF] have also been derived from the exact factorization [28,29], which is essentially a rewriting of the molecular Schrödinger equation in terms of coupled electronic and nuclear equations. This formalism introduces the new concepts of time-dependent vector potential [30][31][32] and timedependent scalar potential, or time-dependent potential energy surface [33][34][35]. They govern the nuclear quantum-mechanical evolution and encode the full dynamical effect of the electrons, including nonadiabatic effects [36][37][38][39]. In turn, the electron-nuclear coupling is explicitly given in terms of the spatial variation of the nuclear wavefunction [40,41]. Evolving classicallike nuclei according to the force determined from the time-dependent potentials and reconstructing the nuclear distribution from an ensemble of trajectories yields the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm [22]. Simplification of this numerical procedure have been proposed, by combining the adiabatic evolution of the trajectories and the hopping idea of TSH with the reconstruction of the nuclear distribution, either by using coupled trajectories, i.e., CT-TSH [25], or by using auxiliary trajectories, i.e., SHXF [START_REF] Ha | [END_REF][42][43][44][45][46].

It is fair to say that nowadays the exact factorization is a reasonably well-known approach which provides an alternative perspective to the Born-Huang representa-tion [31,32] on the chatacterization of both static [47][48][49][50][51] and dynamic [24,[START_REF] Ha | [END_REF]34,52,53] properties of molecular -but not only molecular [54] -systems. Many authors have contributed in different ways to various aspects of the theory [39,41,[55][56][57][58][59][60][61][62][63]. Regarding, in particular, the trajectory-based picture introduced starting with the exact factorization of the time-dependent electron-nuclear wavefunction, abundant literature already reported studies on the time-dependent potentials [30,31,[33][34][35]64], on classical trajectories evolved using the (numericallyexact) time-dependent potentials [37,65], and on the CT-MQC [22][23][24] and CT-TSH [25] algorithms presented above. Despite those studies, some questions still remain unaddressed. For instance, can the classical trajectories, evolved on the time-dependent potentials, be considered as the best possible approximation to the exact-factorization dynamics? Is it plausible to expect that introducing the concept of quantum trajectories [32,66,67] within the exact factorization will significantly improve the performance of a trajectory-based description of nonadiabatic dynamics? How well CT-MQC and CT-TSH compare to quantum trajectories? Is there an actual gain in combining the coupled-trajectory idea of the exact factorization with surface hopping, which introduces a severe approximation on the force used to evolve the trajectories? All these questions touch upon very fundamental aspects of the theory. However, the success of the exact factorization, and related algorithms, is partially due to the fact that it is a rigorous framework, where approximations are systematically derived and tested. Therefore, we believe that answering those questions to evaluate the capabilities of the method will provide useful information on its current status.

It is worth mentioning that Ref. [32] already presented a comparison between classical and quantum trajectories in the framework of the exact factorization. However, the purpose of that study was different from the present work, as it aimed to use trajectories to clarify the effect of the Condon approximation in the treatment of a photo-excitation by an ultrashort laser pulse. Furthermore, in Ref. [25] the comparison between CT-MQC and CT-TSH has been reported, but the work mainly focused on the detailed derivation of the CT-TSH algorithm with a proof-of-principle illustration on a onedimensional model.

In this work we aim to answer the above questions by presenting a study the photo-isomerization reaction in a retinal chromophore model [68] completely based on the exact factorization: (i) we perform quantum vibronic wavepacket dynamics with ElVibRot [START_REF] Lauvergnat | ElVibRot-TnumTana Quantum Dynamics Code[END_REF], which has been adapted to output the necessary information to compute the time-dependent potentials, (ii) we extract the time-dependent potentials from this dynamics using the Exact Factorization Analysis Code (EFAC) [START_REF] Bessone | [END_REF], (iii) we generate quantum trajectories by evolving them using the Bohmian representation of the nuclear velocity field [START_REF] Holland | The Quantum Theory of Motion -An Account of the de Broglie-Bohm Causal Interpretation of Quantum Mechanics[END_REF] in combination with the time-dependent potentials, (iv) we compare the performance of CT-MQC and CT-TSH using the G-CTMQC code [START_REF] Agostini | G-CTMQC[END_REF] with the quantum trajectories by calculating observable properties usually of interest in a nonadiabatic simulation. Note that, the potentials constructed from the solution of the time-dependent Schrödinger equation (point (ii) above) can be viewed as incorporating all quantum dynamical effects of the electrons on the nuclei in nonadiabatic conditions. Our quantum trajectories are, in this sense, the best possible trajectory-based description of the photoisomerization process, thus allowing us to assess how well the coupled-trajectory approximate methods (could) perform. In addition, we will analyze the time-dependent potentials themselves in order to show that despite the used retinal model is very simple (two-dimensional in nuclear space with two electronic states) their behavior at long times is complex and, consequently, difficult to capture in a trajectory-based description of the dynamics, as was pointed out in previous work [START_REF] Pieroni | [END_REF].

The paper is organized as follows. In Section II we briefly recall the exact factorization and, in Section II A, the coupled-trajectory algorithms derived from it. Our numerical results are presented in Section III, focusing first on the time-dependent potentials, in Section III A, and then on the analysis of the trajectories, in Section III B. Our conclusions are stated in Section IV.

II. THEORETICAL FRAMEWORK

The non-relativistic molecular Hamiltonian

Ĥ(r, R) = Nn ν=1 -2 ∇ 2 ν 2M ν + Ĥel (r, R) (1)
describes a system of interacting electrons and nuclei, whose positions are collectively indicated as r and R, respectively. The nuclear kinetic energy operator is expressed in Cartesian coordinates and contains a sum over the N n nuclei, each labeled with the index ν, with spatial derivatives ∇ ν with respect to nuclear positions; M ν are the nuclear masses. The electronic Hamiltonian Ĥel (r, R) is the sum of the electronic kinetic energy and of all interactions. The time evolution of the electron-nuclear system is dictated by the timedependent Schrödinger equation (TDSE)

i ∂ t Ψ(r, R, t) = Ĥ(r, R)Ψ(r, R, t) (2)
whose solution yields the time-dependent molecular wavefunction Ψ(r, R, t).

When introducing the exact factorization [4,28], the molecular wavefunction is written as the product

Ψ(r, R, t) = χ(R, t)Φ(r, t; R) (3)
where the nuclear wavefunction χ(R, t) evolves according to the nuclear TDSE

i ∂ t χ(R, t) = (4) Nn ν=1 [-i ∇ ν + A ν (R, t)] 2 2M ν + (R, t) χ(R, t)
and the electronic conditional factor Φ(r, t; R), parametrically depending on R, evolves according to the electronic equation

i ∂ t Φ(r, t; R) = (5) ĤBO (r, R) + Ûen [Φ, χ] -(R, t) Φ R (r, t)
Both evolution equations ( 4) and ( 5) contain the time-dependent vector potential (TDVP) A ν (R, t) and the time-dependent potential energy surface (TDPES) (R, t) defined as

A ν (R, t) = Φ(t; R)| -i ∇ ν Φ(t; R) r (6) (R, t) = Φ(t; R)| ĤBO (R) + Ûen [Φ, χ] -i ∂ t |Φ(t; R) r (7)
The symbol • r stands for an integration over the electronic degrees of freedom, and we removed all dependencies on r within • r to imply that this variable is integrated out. The TDVP and the TDPES are electronic quantities that evolve in time as effect of electronic dynamics, which is coupled to the nuclear dynamics in nonadiabatic conditions, and completely determine the nuclear evolution via the nuclear TDSE (4). It is interesting to note that the TDVP is related to the nuclear momentum field, as

A ν (R, t) = Im[ Ψ(R, t)| ∇ ν |Ψ(R, t) r |χ(R, t)| 2 -∇ ν S(R, t) (8)
a property that will be used in Section II A to define our "quantum trajectories". The symbol S(R, t) is used to indicate the phase of the nuclear wavefunction χ(R, t) and it has the dimensions of an action. The electronic equation ( 5) contains the electron-nuclear coupling operator

Ûen [Φ, χ] = Nn ν=1 1 M ν [-i ∇ ν -A ν (R, t)] 2 2 (9)
+ -i ∇ ν χ(R, t) χ(R, t) + A ν (R, t) (-i ∇ ν -A ν (R, t))
which depends explicitly on the nuclear wavefunction, on the electronic wavefunction via the presence of the TDVP, and acts on the parametric dependence of Φ(r, t; R) as a spatial derivative. The ambiguity of the product form of the molecular wavefunction in Eq. ( 3) is partially eliminated by imposing the normalization condition |Φ(r, t; R)| 2 dr = 1 ∀R, t, which allows to identify |Φ(r, t; R)| 2 as a conditional probability density and |χ(R, t)| 2 = |Ψ(r, R, t)| 2 dr as a marginal probability density that yields the nuclear density from the full wavefunction. Equation ( 3) is thus unique up to a gauge encoded in a phase factor e (i/ )θ(R,t) (with θ(R, t) a real function): multiplying the nuclear wavefunction by this factor and the electronic term by its complex conjugate, Eq. (3) remains unaffected and the time-dependent potentials transform as standard gauge potentials. Therefore, this ambiguity has to be eliminated by imposing a choice of gauge as discussed in Section II A.

A. Trajectory-based description

In a trajectory-based description of the coupled electron-nuclear dynamics in the presence of nonadiabatic effects and described within the exact-factorization framework, nuclear dynamics from Eq. ( 4) is approximated in terms of an ensemble of coupled trajectories. Their positions and momenta evolve according to Hamilton equations under the effect of the TDVP and TD-PES [66]. In turn, the TDVP and TDPES are computed from the electronic wavefunction, which is expanded in the adiabatic basis, i.e., the basis of eigenstates of Ĥel (r, R). Equation ( 5) yields evolution equations for the expansion coefficients, which are solved coupled to the nuclear Hamilton equations. This is the essence of the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm [22]. It has been abundantly described in the literature, therefore, we refer the interested reader to Refs. [4,[74][START_REF] Agostini | Handbook of Materials Modeling[END_REF][START_REF] Gossel | [END_REF] for in-depth discussions, and we recall here only the main ideas and equations.

A nuclear trajectory, indicated with the symbol R α (t), is intended as a collection of 3N n positions that evolve in time. When the electronic wavefunction Φ(r, t; R α (t)) is expressed as a linear combination of adiabatic states, ϕ (m) (r; R α (t)), their parametric dependence on the nuclear position/trajectory induces an implicit time dependence. The expansion coefficients in Φ(r, t;

R α (t)) = m C m (R α (t), t) ϕ (m) (r; R α (t))
depend on the nuclear trajectory as well. In the following, the dependence on R α (t) will be indicated only via the index α for simplicity.

For the trajectory α, the classical force on the nucleus ν at time t is F α ν (t) and is written as the sum of three terms

F α ν (t) = F α ν,mf (t) + F α ν,na (t) + F α ν,ct (t) (10)
The mean-field contribution F α ν,mf (t) is the average of adiabatic forces -∇ ν E α m over the electronic states weighted by the population of the states |C α m (t)| 2 , namely

F α ν,mf (t) = m |C α m (t)| 2 (-∇ ν E α m ) (11)
The adiabatic energies are indicated as E m . The nonadiabatic contribution F α ν,na (t) depends explicitly on the nonadiabatic coupling vectors d α ν,ml =

ϕ (m) (R α (t))|∇ ν ϕ (l) (R α (t)) r and is F α ν,na (t) = m,l Cα m (t)C α l (t) (E α m -E α l ) d α ν,ml (12)
with Cα m (t) the complex conjugate of C α m (t). The last contribution F α ν,ct (t) arises as effect of the coupling among the trajectories

F α ν,ct (t) = 2 m |C α m (t)| 2 Nn µ=1 P α µ (t) • f α µ,m f α ν,m -A α ν (t) (13)
and depends on the TDVP A α ν (t), on the, so-called, accumulated adiabatic force

f α ν,m = t 0 (-∇ ν E α m
) dτ , and on the quantum momentum P α µ (t). Calculation of the quantum momentum requires knowledge of the positions of all trajectories at a given time t, thus the trajectories cannot be evolved independently [77]. The quantum momentum follows from the term in the expression of the electron-nuclear coupling operator (9) that depends explicitly on the nuclear wavefunction; its expression is related to the spatial variation of the nuclear density as

P α µ (t) = -∇ µ |χ(R α (t), t)| 2 2 |χ(R α (t), t)| 2 (14)
In CT-MQC the nuclear density is reconstructed as a sum of frozen Gaussians centered at the positions of the trajectories, and the quantum momentum is determined from Eq. ( 14). The electronic equation ( 5) yields evolution equations for the expansion coefficients C α m (t) in the form

Ċα m (t) = Ċα m,TSH (t) + Ċα m,ct (t) (15)
where the first term is identical as in a standard trajectory surface hopping (TSH) procedure [START_REF] Tully | [END_REF]

Ċα m,TSH (t) = - i E α m C α m (t) - l Nn ν=1 Ṙα ν (t) • d α ν,ml C α l (t) (16)
with Ṙα ν (t) the nuclear velocity, and the additional term Ċα m,ct (t) depends on the quantum momentum

Ċα m,ct (t) = Nn ν=1 P α ν (t) M ν • f α ν,m -A α ν (t) C α m (t) (17)
similarly to the expression of the force.

In deriving CT-MQC equations, the gauge has been chosen such that α (t) + ν Ṙα ν (t)•A α ν (t) = 0, with α (t) the value of the TDPES at the position of the trajectory.

Recently, the coupled-trajectory idea of CT-MQC, encoded in the quantum momentum, was combined with a TSH procedure to derive the CT-TSH scheme. CT-TSH dynamics promises to be more stable than CT-MQC at long times, at least when, as it was shown for a onedimensional system [25], the TDPES develops an oscillatory behavior. This can happen if the trajectories cross a nonadiabatic region several times and interfere. In addition, the computational cost for the force calculation is reduced if compared to CT-MQC because the nonadiabatic force, depending on the nonadiabatic coupling vectors, does not appear. In CT-TSH, the force driving a trajectory is simply determined from the gradient of the adiabatic energy of the running state, i.e., the state associated to that trajectory at a given time. The running state changes over time as a trajectory can hop from one state to another according to the standard fewest switches procedure. The probability for a hop from state m to state l to occur at time t is determined as [9]

P α ml (t) = max 0, -2 Cα m (t)C α l (t) ν d α ν,lm • Ṙα ν (t) |C α m (t)| 2 dt (18)
with dt the integration time step. Being P α ml (t) a real quantity, one has to take the real part of the numerator in brackets. CT-TSH trajectories remain, however, coupled, because Eq. ( 15) is used to propagate the electronic coefficients, as in CT-MQC, thus the quantum momentum has to be computed. When the electronic evolution equation of standard TSH is complemented with the coupled-trajectory term of Eq. ( 17) the overcoherence problem of TSH is overcome, since the quantum momentum allows one to account for decoherence effects, as shown in previous work [22,24,25,77].

In CT-MQC and CT-TSH both the electronic and the nuclear equations of the exact factorization are approximated, as described above. Therefore, a question that is indeed interesting to address is: How well can one hope to reproduce nuclear dynamics using trajectories? To answer this question, we introduce the concept of quantum trajectories, as trajectories that are propagated under the effect of the exact TDVP and TDPES [32]. For a model system, the TDVP and the TDPES can be computed numerically exactly from the output of vibronic wavepacket dynamics. Furthermore, if a gauge is chosen such that the TDVP is the nuclear momentum field [30][31][32], then the simple equation

Ṙα ν (t) = A α ν (t) M ν (19)
can be integrated to determine the positions of the trajectories at all times that closely follow the evolution of the (exact) nuclear density. Here, A α ν (t) stands for the value of the exact TDVP at the position R α ν (t). In the gauge where the nuclear wavefunction is always real and non-negative, the second term on the right-hand side of Eq. ( 8) identically vanishes and the TDVP equals the nuclear momentum field calculated from the molecular wavefunction. Therefore, to determine the (numerically exact) TDVP, the solution of the TDSE (2) is integrated over electronic positions to get the nuclear density, whose positive square root gives the nuclear wavefunction in this gauge. Dividing the molecular wavefunction by the nuclear wavefunction yields the electronic wavefunction, from Eq. ( 3), and Eq. ( 6) can be directly used to determine the TDVP.

Note that we refer to the trajectories determined according to Eq. ( 19) as quantum trajectories in a (Bohmian) hydrodynamic sense. If one focuses only on the nuclear TDSE (4), then the TDVP is the Bohmian definition of the momentum field for quantummechanical particles described by a real wavefunction in the presence of an external vector potential [START_REF] Holland | The Quantum Theory of Motion -An Account of the de Broglie-Bohm Causal Interpretation of Quantum Mechanics[END_REF]. It was recently proven [32] that this procedure is equivalent to propagating trajectories according to forces that include the combined effect of the TDPES and of the quantum potential -usually arising from the Bohmian formulation of quantum dynamics. In addition, the analysis reported here differs from the setup of Ref. [37] where the dynamics of the trajectories was driven only by the TDPES (without the quantum potential) and, thus, the trajectories were purely "classical".

III. NUMERICAL RESULTS

The model Hamiltonian [68] used in this work allows us to describe the 11-cis to all-trans isomerization of the retinal chromophore in rhodopsin [30,[78][79][80][81][82][83][84][85][86][87][88][89]. Photoabsorption takes the cis conformer from the electronic ground state S 0 to the first excited state S 1 , thus initiating an isomerization reaction towards the trans conformer. The nuclear coordinates defining the model are an angular reaction coordinate ϕ and a collective vibration q characterizing "a delocalized stretching motion of the polyene chain, whereby single and double bonds interchange" [68].

The Hamiltonian is given in the diabatic basis

Ĥ(ϕ, q) = P 2 ϕ 2m + P 2 q 2ω -1 + V 00 (ϕ, q) V 01 (q) V 10 (q) V 11 (ϕ, q) (20)
supposing here that R = ϕ, q (dimensionless coordinates). The nuclear kinetic energy operator is expressed via the momentum operators along ϕ and q, where P 2 ϕ and P 2 q are second-order derivatives with respect to ϕ and q, respectively. The elements of the 2 × 2 electronic matrix are

V 00 (ϕ, q) = 1 2 W 0 (1 -cos ϕ) + 1 2 ωq 2 (21) V 11 (ϕ, q) = E 1 - 1 2 W 1 (1 -cos ϕ) + 1 2 ωq 2 + κq (22) V 01 (q) = V 10 (q) = λq (23)
The stretching mode q is also called coupling mode, as the potential coupling V 10 only depends on q. The effective mass of the reaction coordinate m -1 = 4.84•10 -4 eV is chosen so as to achieve isomerization within 200 fs; the frequency of the coupling mode is ω = 0.19 eV and the interstate coupling is λ = 0.19 eV, with κ = 0.1 eV indicating a gradient in the excited state; the parameters W 0 = 3.6 eV, W 1 = 1.09 eV and E 1 = 2.48 eV are chosen in order to match the S 0 /S 1 gap to the center frequency of the absorption bands of the cis and trans isomers. The adiabatic potential energy surfaces (PESs) S 0 and S 1 are obtained by diagonalization of the electronic Hamiltonian in Eq. [START_REF] Kirrander | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF]. The PESs present conical intersections (CIs) at ϕ CI ± π 2 , q CI = 0. The cis conformer is characterized by the value of the reactive coordinate ϕ cis = 0, and at this geometry for q = 0 the excitation energy from S 0 to S 1 lies in the visible domain. The adiabatic PESs are shown in Fig. 1 as colormaps, where the positions of the CIs are indicated with crosses. At the cis geometry and for slightly negative values of the coupling mode, S 0 is a potential well while S 1 has a saddle shape; at the trans geometry (ϕ = ±π), only S 0 presents stable configurations.

The instantaneous excitation from S 0 to S 1 promotes the ground-state wavepacket to the excited state. The ground-state wavepacket is chosen as a two-dimensional Gaussian

χ S0 (ϕ, q, 0) = 4 1 πσ 2 ϕ e -ϕ 2 2σ 2 ϕ 4 1 πσ 2 q e -q 2 2σ 2 q (24)
with σ ϕ = 0.128 and σ q = 1.0, which corresponds to the vibrational ground state of S 0 at the cis geometry. The nuclear (probability) density at all times is given as the sum of the S 0 and S 1 (probability) densities

|χ(ϕ, q, t)| 2 = |χ S0 (ϕ, q, t)| 2 + |χ S1 (ϕ, q, t)| 2 (25)
Electronic populations of the ground (S 0 ) or excited (S 1 ) states are determined as the integral of the corresponding nuclear density,

ρ S k (t) = |χ S k (ϕ, q, t)| 2 dϕdq with k = 0, 1 (26)
Vibronic wavepacket dynamics after photo-excitation is simulated in the diabatic basis using ElVibRot [START_REF] Lauvergnat | ElVibRot-TnumTana Quantum Dynamics Code[END_REF]. The diabatic densities are transformed to the adiabatic representation for the comparison with the trajectorybased approaches. The nuclear wavepacket has been expanded in a two-dimensional basis set formed by the direct product of two uni-dimensional basis sets for each diabatic state. The harmonic-oscillator basis set has been used for q (40 basis functions), and a Fourier-series, periodic basis set has been used for ϕ (512 basis functions); the number of grid points along each coordinate is 80 and 552, respectively for q and ϕ. The propagation was performed with the Chebychev scheme for which the evolution operator is expanded on Chebychev polynomials [90] with time step 0.1 fs. The Hamiltonian needs to be renormalized so that its spectral range lies in the interval [-1, 1].

Figure 1 reports snapshots of the nuclear dynamics at times t = 10, 50, 100 fs, where the contour lines represent the adiabatic components of the nuclear density (in S 1 in the upper panels and in S 0 in the lower panels) and the colormaps represent the adiabatic PESs (S 1 in the upper panels and S 0 in the lower panels). The initially excited density on S 1 (t = 10 fs) extends over time towards the CIs (t = 50 fs) and finally funnels via the CIs to S 0 (t = 100 fs). While the dynamics shows coherent oscillations along q, the isomerization reaction from the cis to the trans conformer is characterized by the wavepacket reaching the values |ϕ| > π 2 . The isomerization process just described will be analyzed in detail employing the exact factorization, in Section III A by looking at TDPES and TDVP, and in Section III B by using quantum and coupled trajectories.

A. Time-dependent potentials

The time-dependent potentials emerging from the exact factorization can be constructed from the output of a vibronic wavepacket quantum simulation. If the electron-nuclear state of the system |Ψ(ϕ, q, t) is expanded in the diabatic basis with coefficients χ (d) k (ϕ, q, t) (k = 0, 1), the nuclear wavefunction in the chosen gauge is χ(ϕ, q, t) = |χ

(d) 0 (ϕ, q, t)| 2 + |χ (d)
1 (ϕ, q, t)| 2 . Therefore, when the electronic conditional state |Φ(t; ϕ, q) is expanded in the same basis, the expansion coefficients are χ (d) k (ϕ, q, t)/χ(ϕ, q, t) by virtue of the exact factorization (3). Using these coefficients and replacing the integration over r of Eqs. ( 6) and ( 7) with a sum over the diabatic states yields the (ϕ, q, t)-dependent TDPES and TDVP. Studying the evolution of their behavior in time allows one to qualitatively interpret and justify the dynamics of the nuclear wavefunction. To compute the TDPES and the TDVP we used the Exact Factorization Analysis Code (EFAC) [START_REF] Bessone | [END_REF].

Time-dependent potential energy surface

In Fig. 2 we report the TDPES (colormap) as function of the nuclear coordinates ϕ and q at the same time steps as in Fig. 1, superimposed with the nuclear density |χ(ϕ, q, t)| 2 (contour lines) and indicating the CIs as crosses. Note that the TDPES is only calculated, and shown, in regions where the nuclear density is larger than 10 -7 for numerical convenience.

At the beginning of the simulated dynamics, t = 10 fs in Fig. 2 (left panel), the TDPES along q is basically a potential well with minimum between q = -2 and q = 0 where the nuclear density is localized; along ϕ and in the region q ∈ [-2, 0] the TDPES is slightly more flat than along q, so as to allow the density to spread towards the CIs. This shape resembles the adiabatic S 1 PES, since at this early time, nuclear dynamics fully takes place in the excited state.

At later times, t = 50 fs in Fig. 2 (central panel), the nuclear density starts spreading symmetrically towards the CIs as consequence of the appearance of a saddle shape along ϕ with maximum in ϕ = 0 in the region q ∈ [-3, 0]. In the regions q < -3 and q > 0, the increase of the TDPES localizes the nuclear density along the q direction. In particular, the shape of the TDPES around ϕ = 0, q = 2 indicates the presence of a small portion of the nuclear density in the ground state. We recall that quantum dynamics is performed in the diabatic basis and, thus, the nuclear state is initialized with full population in the diabatic state k = 1. This means that a very small amount of density is found in the (adiabatic) ground state even at time t = 0, whose evolution is also dictated by the TDPES. In the portion of configuration space q > 0 at t = 50 fs the (adiabatic) excited-state density is much smaller than the ground-state density, thus the TDPES follows the shape of the S 0 PES.

At the final time shown in Fig. 2, t = 100 fs (right panel), the central portion of the TDPES where the nu-clear density is mainly located has, once again, a saddle shape along ϕ with maximum in ϕ = 0, allowing in this way the nuclear density to reach the CIs and to "funnel" to the ground state, thus accessing the trans configuration. However, note that in a description based on the TDPES, we do not need to invoke any funnelling process to justify the transition from the cis to the trans region, which is, in fact, simply accompanied by the shape of the TDPES. The highly oscillatory features that develop in the colored areas of the TDPES arise from the interferences between the S 1 density transferred in S 0 and the S 0 density. After t = 100 fs and up to t = 200 fs, the shape of the TDPES in the trans region is mainly characterized by interferences between the S 1 portion of the density that continuously funnels through the CIs and the S 0 portion that accumulates there.

It is worth noting that periodic boundary conditions are used along the reactive coordinate ϕ. Therefore, interferences are produced as well by S 0 portions of the nuclear wavepacket that reach the trans configuration at +π (-π) and continue their evolution "reentering" the simulation box from -π (+π).

Time-dependent vector potential

In Fig. 3 the TDVP offers a complementary perspective on the nuclear dynamics to the TDPES, and, in the used gauge, it is identified with the nuclear momentum field. As for the TDPES, the TDVP is only shown in regions where the nuclear density is larger than 10 -7 . At t = 10 fs (left panel) the TDVP is small in magnitude because the dynamics starts with zero average momentum, and appears symmetric in the ϕ direction with respect to ϕ = 0. This feature is maintained all along the simulated dynamics, as the nuclear density moves symmetrically towards the CIs. The TDVP reported in all panels in Fig. 3 changes sign in the ϕ direction at ϕ = 0.

As the nuclear density approaches the CIs, at t = 50 fs (central panel), the magnitude of the TDVP increases. This behavior can be easily explained having in mind the double-cone shape of the adiabatic PESs, with the slope of the S 1 PES inducing the nuclear density to efficiently relax towards S 0 via the CIs.

At the final time shown in Fig. 3, t = 100 fs (right panel), the magnitude of the TDVP reflects the interference patterns already observed in the TDPES, with a clear separation between the smooth behavior in the cis region |ϕ| < π 2 and a highly oscillatory behavior mainly identified in the trans region |ϕ| > π 2 . In all panels in Fig. 3, the vector field representing the direction of the TDVP is nearly parallel to the horizontal axis, indicating that the q component of the momentum field is much smaller than the ϕ component. This results in dynamics mainly characterized by motion from the cis to the trans configuration.

B. Quantum and coupled trajectories

Quantum trajectories are determined by intergrating Eq. ( 19) where the TDVP is calculated as described in Section III A from the output of our quantum-dynamics simulation. We refer to A ϕ (ϕ, q, t) and to A q (ϕ, q, t) as the two components of the TDVP. Initial conditions for the quantum trajectories are determined by sampling 1000 positions R α (0) = ϕ α (0), q α (0) from the Gaussian probability distribution |χ S0 (ϕ, q, 0)| 2 associated to the chosen initial condition for the quantum simulation. Once the initial positions are given, the initial velocities are simply determined as φα (0), qα (0) = A ϕ (ϕ α (0), q α (0), 0)/m, A q (ϕ α (0), q α (0), 0)/ω -1 .

Coupled trajectories are propagated according to the CT-MQC and CT-TSH algorithms. We used 1000 initial conditions, whose positions and momenta are independently sampled from the Wigner distribution determined from |χ S0 (ϕ, q, 0)| 2 (the initial positions are identical to the initial quantum trajectories). In the case of CT-MQC, we extract nuclear and electronic observables by averaging over 1000 trajectories. In the case of CT-TSH, we perform 4 independent runs in order to account for different histories of hops, thus yielding a total of 4000 trajectories whose results are averaged to extract the observables that will be discussed below.

CT-MQC and CT-TSH dynamics are performed in the adiabatic basis, with adiabatic PESs determined by diagonalizing the electronic Hamiltonian in Eq. [START_REF] Kirrander | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF]. The nonadiabatic couplings are computed analytically with the usual Hellmann-Feynman theorem [91] from the eigenvectors of the electronic Hamiltonian in Eq. ( 20). This procedure is implemented in the ModelLib library [START_REF] Lauvergnat | QuantumModelLib[END_REF]. The velocity-Verlet algorithm is used to integrate classical nuclear equations, and the fourth-order Runge-Kutta method is used for the electronic evolution; nuclear and electronic dynamics are integrated with at time step of 0.0024 fs (0.1 au). We checked the stability of CT-TSH results by reducing/increasing by half the integration time step, since the hopping probability of Eq. ( 18) depends on it. The electronic initial conditions

Distributions of trajectories

In Fig. 4 we compare the distributions of nuclear trajectories with the quantum nuclear density (green areas) at times t = 10, 50, 100 fs along the dynamics. Quantum trajectories are indicated as lilac dots in the upper panels, CT-MQC trajectories as blue dots in the middle panels and CT-TSH trajectories as orange dots in the lower panels. In general, we observe that CT-MQC and CT-TSH trajectories resemble closely the distribution of quantum trajectories, in agreement with the initialization of all trajectories ensembles from the same distribution in configuration space.

In Fig. 5 the distributions of momenta are shown at times t = 10, 50, 100 fs along the dynamics. In the upper panels, CT-MQC distributions (blue dots) are compared to the quantum trajectories (lilac dots), whereas, in the lower panels, CT-TSH distributions (orange dots) are compared to the quantum trajectories. In both cases a clear qualitative difference is observed, which follows from the different initial momenta used for the quantum and for the coupled trajectories. In CT-MQC and CT-TSH positions and momenta are independent variables, and thus they are randomly sampled from uncorrelated (Wigner-transformed) positions-momenta distributions. Note that the distributions are uncorrelated since the (quantum) initial density is a two-dimensional Gaussian in the variables ϕ and q. Figure 5 shows that the momenta distribution of quantum trajectories is much more localized along P q than the distributions of coupled trajectories, whereas they have comparable widths along P ϕ at time t = 10 fs (left panels). The momenta of the quantum trajectories rapidly increase only along P ϕ , and the distribution in this direction at time t = 50 fs (central panels) is as broad as CT-MQC and CT-TSH. At later times, t = 100 fs in the figure (right panel), the distributions of all ensembles of trajectories have comparable broadening along P ϕ . However, quantum trajectories remain extremely localized in momentum space along P q if compared to CT-MQC and CT-TSH.

The results shown in Fig. 5 are consistent with the previously observed property of the TDVP of being much larger along the ϕ direction than along the q direction. In Fig. 3 the vector field representing the direction of the TDVP is basically parallel to the horizontal axis. Fig. 5 shows, in fact, that P ϕ attains values up to approximately 50 times the values of P q .

In general, Figs. 4 and5 show that the initial dynamics simulated with CT-MQC and CT-TSH is basically identical, and this is because the trajectories evolve adiabatically on S 1 . When they reach the CIs, however, the treatment of the dynamics differs regarding both the nuclear forces and the evolution of the electronic coefficients. When and after the nonadiabatic event occurs, CT-MQC and CT-TSH yield (slightly) different results. The analysis of Figs. 4 and 5 also suggests that we can directly compare average configurational properties extracted from quantum trajectories and from CT-MQC/CT-TSH trajectories. On the other hand, the behavior of quantum and coupled (or classical-like) trajectories in momentum space is qualitatively different as consequence of the different relations between initial positions and initial mo-menta in the two schemes. This point will be developed further below.

Averaged observables

Figure 6 shows the population of the electronic excited state S 1 as function of time (left panel), where reference results are provided by Eq. ( 26) and are indicated in red. For CT-MQC (blue) and CT-TSH (orange) this same observable is estimated as the average over the trajectories of |C α S1 (t)| 2 , namely as

ρ q. S1 (t) = 1 N traj Ntraj α=1 |C α S1 (t)| 2 (27)
We refer to this average as "quantum" and we indicate it as (q.) in Fig. 6. In addition, in CT-TSH, one can also count how many trajectories are "on" the S 1 PES at time t, i.e., N S1 (t), and divide this number by the total number of trajectories (green line in Fig. 6), namely

ρ cl. S1 (t) = N S1 (t) N traj (28)
We refer to this estimate of the electronic population as "classical" and we indicate it as (cl.) in Fig. 6. In Fig. 6 we report as well TSH results (black) using the energy decoherence correction (with parameter 0.1 Ha), indicated as TSH-ED; only the classical estimate of the electronic population obtained from TSH-ED is shown in the figure. Figure 6 reports the probability of the trans configuration as function of time as well (right panel). Exact results (red) are determined by integrating at all times the total nuclear density between ϕ = -π 2 and ϕ = + π 2 (and over all q), thus obtaining P cis (t), and then by defining P trans (t) = 1 -P cis (t). Trajectory-based results are determined by counting at each time how many trajectories have a value of |ϕ α (t)| > π 2 , independently of q α (t), and by dividing this number by the total number of trajectories. CT-MQC results are shown in blue, CT-TSH in orange, TSH-ED in black, quantum trajectories (evolved with the TDVP) in lilac in Fig. 6 (right panel). Note that, when the trajectories are evolved with the TDVP, only their distribution in configuration and momentum space is available. For this reason, we can only compute the probability of the trans configuration using quantum trajectories and we do not have access to the electronic population.

Quantum (exact) results from Fig. 6 show that after around 80 fs of dynamics in the excited state, the photoexcited wavepacket reaches the CIs and starts transferring population to the ground state. At the same time, the trans configuration is being formed as the wavepacket tends towards values of the reactive coordinates |ϕ| > π 2 . While the S 1 population keeps decreasing until 200 fs, i.e., the end of the simulated dynamics, the trans configuration reaches a maximum at around 180 fs and decreases FIG. 6. Left: Population of the electronic excited state S1 as function of time from quantum dynamics (Exact, in red), CT-MQC (in blue), CT-TSH (in orange and in green), and TSH-ED (in black). For CT-MQC only the quantum (q.) estimate of the population is accessible, whereas for CT-TSH we compare the quantum (q.) and the classical (cl.) expressions, as defined in the text. For TSH-ED, only the classical (cl.) estimate is shown. Right: Probability of finding the system in the trans configuration as function of time from quantum dynamics (Exact, in red), CT-MQC (in blue), CT-TSH (in orange), TSH-ED (in black), and quantum trajectories (in lilac).

afterwards. This behavior suggests transfer of density from the trans basins to the cis basin on the S 0 PES without population transfer, since the population of S 1 does not increase after 180 fs. An additional feature that we can observe in the plot of the population of S 1 is that the curve hints towards a stabilization after 200 fs.

The formation of the trans conformer is reproduced by the quantum trajectories (lilac) in extremely good agreement with exact results, indicating that a trajectorybased representation of nuclear dynamics for this retinal model is adequate to study the trans photo-product formation after excitation by light. However, TDVP results slightly deviate from the reference towards the end of the simulated dynamics, due to numerical errors in the time integration of Eq. ( 19) when the TDVP develops a highly oscillatory behavior (see Fig. 3 at time 100 fs).

CT-MQC results (blue) reproduce qualitatively and quantitatively well the expected behavior of the S 1 population, eve though it seems that the time trace of the population reaches a plateau at around 180-190 fs. The trans probability as function of time is also reproduced extremely well, with the maximum only slightly overestimated if compared to the reference.

Interestingly, the quantum (orange) and classical (green) estimates of the electronic population in CT-TSH are very close all along the dynamics, even though this has not been imposed in the derivation of the algorithm. The quantum estimate is, however, closer to CT-MQC, as the population reaches a plateau at around 180-190 fs, which is not observed in the classical estimate. The classical estimate of the S 1 population resembles closely the reference curve as well as TSH-ED results. In addition, the prediction of the trans probability by CT-TSH is very close to the exact curve, and shows a slight (quantitative) improvement with respect to both CT-MQC and TSH-ED. FIG. 7. Average kinetic energy along the coupling mode q (left panel) and along the reactive coordinate ϕ (right panel) as function of time. The color code is the same as in Fig. 6.

In Fig. 7 we analyze the nuclear kinetic energies along the coupling mode q (left panel) and along the reactive coordinate ϕ (right panel) as function of time for all simulation methods as in Fig. 6.

The comparison clearly highlights the fact that the kinetic energy of the quantum trajectories along q is shifted at the initial time with respect to classical-like coupled trajectories, due to the different initial conditions. More specifically, quantum-mechanically the expectation value of the kinetic energy operator is the sum of a pseudoclassical term related to the momentum and of an additional contribution arising from the delocalization of the wavefunction. The latter remains unaccessible to our quantum-trajectory approach and, consequently, the shift observed at the initial time persists all along the simulated dynamics. Another aspect that appears clearly in the TDVP results shown in Fig. 7 is the appearance of numerical instabilities after 130 fs, producing (unphysical) spikes in the kinetic energy of quantum trajectories along both coordinates. As observed for Figs. 2 and3, the time-dependent potentials, and here in particular the TDVP, develop a highly oscillatory behavior due to interfering portions of the nuclear wavepackets (and not due to numerical issues). Those interferences are indeed difficult to capture within a quantum-trajectory approach, which might largely limit the possibility of combining quantum trajectories with the exact factorization beyond the proof-of-principle studies [66,67] reported so far. However, studies in this direction are currently ongoing.

The instabilities developed at long time by the timedependent potentials seem to affect as well CT-MQC trajectories. Although the instabilities appear less severe, they manifest themselves in the kinetic energy along the coupling coordinate as an increase of the energy with respect to reference (exact) results and as coherent oscilla-tions persisting until the end of the simulated dynamics. Some oscillations are also visible in the kinetic energy behavior along the reactive coordinate.

In general, Fig. 7 attests to an actual improvement of CT-TSH over CT-MQC in capturing the kinetic energy along q and along ϕ, with qualitative and quantitative good agreement if compared to quantum results. The energy along q is slightly underestimated while the energy along ϕ is slightly overestimated, but the qualitative behavior is extremely good, especially in capturing the loss of coherence over time along the coupling mode. In CT-TSH the force used to propagate the trajectory is purely adiabatic, as in standard surface hopping, which helps in completely overlooking the interferences appearing in the time-dependent potentials, thus producing a more stable dynamics. It is worth noting that such a detailed analysis of the energy for CT-TSH was not reported before in Ref. [25]. As in Fig. 6, we observe that CT-TSH improves as well, even though mainly quantitatively, on TSH-ED.

IV. CONCLUSIONS

In this paper we reported the analysis of the relaxation dynamics through a conical intersection in a twodimensional retinal model and fully based on the exact factorization. The aim of the work is to validate approximate trajectory-based results obtained by applying CT-MQC and CT-TSH against quantum wavepacket dynamics and against quantum trajectories, which represent for us the best trajectory-based results accessible within the exact factorization. To propagate quantum trajectories, we calculated numerically the time-dependent vector potential and time-dependent potential energy surface from the solution of the coupled electron-nuclear problem.

The agreement between quantum trajectoires and quantum wavepacket dynamics is excellent all along the simulated dynamics, even though the time-dependent potentials show complex interferences at long times due to the reduced dimensionality of the model. Despite this complex dynamical features, the approximate trajectorybased schemes CT-MQC and CT-TSH are in close agreement to quantum results, i.e., trajectories and wavepackets. Furthermore, CT-TSH shows a slightly improved quantitative agreement with the reference over CT-MQC.

FIG. 1 .

 1 FIG. 1. Quantum nuclear dynamics in terms of adiabatic wavepackets contributions (contour lines) superimposed to the adiabatic PESs (colormaps). In each panel the value Maxiso is indicated, which is the largest iso-value of the density represented in the figure. The upper panels refer to S1 and the lower panels to S0. The color bars report the energies in Ha. The nuclear densities in S0 and in S1 are shown at times t = 10, 50, 100 fs as indicated in the figure. The CIs are indicated by the crosses.

FIG. 2 .

 2 FIG. 2. TDPES (colormaps) at times t = 10, 50, 100 fs as indicated in the figure. The color bar reports the energy in Ha. The nuclear density (contour lines) at the same times is superimposed for reference. The CIs are indicated by the crosses.

FIG. 3 .

 3 FIG. 3. Magnitude of the TDVP (colormaps) and its direction (unit vector) at times t = 10, 50, 100 fs as indicated in the figure. The nuclear density (contour lines) at the same times is superimposed for reference. The CIs are indicated by the crosses.

FIG. 4 .

 4 FIG. 4. Exact nuclear density (colormap) at times t = 10, 50, 100 fs as indicated in the figure, compared with the positions of the quantum trajectories (lilac dots in the upper panel), CT-MQC trajectories (blue dots in the middle panels) and CT-SH trajectories (orange dots in the lower panels) at the same times. The CIs are indicated by the crosses.

FIG. 5 .

 5 FIG. 5. Distribution of the trajectories in momentum space at times t = 10, 50, 100 fs as indicated in the figure. Upper panels: Comparison between quantum trajectories (lilac dots) and CT-MQC trajectories (blue dots). Lower panels: Comparison between quantum trajectories (lilac dots) and CT-TSH trajectories (orange dots).

ACKNOWLEDGEMENTS

This work was supported by the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014 and by the ANR HYTRAJ project, Grant No. ANR-19-CE30-0039 of the French Agence Nationale de la Recherche.


André Nauts

David Lauvergnat
email: david.lauvergnat@lcp.u-psud.fr

Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol

Keywords: Torsional level, Methanol, Kinetic energy operator, On-the-fly

In quantum dynamics, physically well-adapted curvilinear coordinates are coordinates that lead to a Hamiltonian operator as separable as possible, in order to simplify the resolution of the corresponding time-independent or time-dependent Schrödinger equations. Various equivalent curvilinear expressions of the kinetic energy operator (KEO) are well known. They can be used in either an analytical or a numerical approach. The latter has the feature of allowing to straightforwardly compute the KEO in terms of sophisticated (yet easy to define) physically welladapted curvilinear coordinates. Nevertheless, the number of terms to be computed on a full grid, scales as n 2 /2 (n being the number of degrees of freedom), so that, for systems with n>10, the memory storage of the KEO's becomes extremely demanding and therefore often unrealistic. We show here that it is possible, starting from the basic quantum expression of the KEO as a curvilinear Laplacian operator, to reduce the memory storage bottleneck by numerically computing the KEO on-the-fly i.e. each time it is required without computing the extrapotential term. This new approach opens the way to rigorous quantum studies of systems with many degrees of freedom.

The comparison of torsional levels of methanol obtained by the present on-the-fly method with our previous results shows excellent agreement.

I Introduction

In molecular physics it is often essential for solving problems efficiently to use curvilinear coordinates (such as bond lengths, bond angles, torsional angles…) physically well adapted to the molecular systems and the processes under study. Indeed, an optimal set of curvilinear coordinates has to lead to a Hamiltonian operator that is as separable as possible, i.e. that reduces as much as possible the coupling terms between the coordinates, in order to simplify the resolution of the corresponding time-independent or time-dependent Schrödinger equations.

The quantum kinetic energy operator, KEO, has various well-known equivalent expressions in terms of curvilinear coordinates [START_REF] Podolsky | Quantum-mechanically correct form of Hamiltonian function for conservative systems[END_REF][2] [START_REF] Wilson | Molecular vibrations : the theory of infrared and Raman vibrational spectra[END_REF][4] [START_REF] Gatti | Chpt. 6, Applications of Quantum Dynamics in Chemistry[END_REF]. These expressions lead to two main approaches: (i) an analytical one, adapted to particular sets of coordinates, such as Jacobi or polyspherical coordinates [START_REF] Chapuisat | Vector parametrization of the N-body problem in quantum mechanics: Polyspherical coordinates[END_REF][7] [START_REF] Ndong | Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach[END_REF]. (ii) a numerical one, [START_REF] Meyer | General Internal Motion of Molecules, Classical and Quantum-Mechanical Hamiltonian*[END_REF][10] [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][13] [START_REF] Strobusch | A general nuclear motion Hamiltonian and non-internal curvilinear coordinates[END_REF] particularly well-suited when sophisticated physically well-adapted coordinates are needed. For instance, with TNUM, [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][15] [START_REF] Sarka | Rovibrational energy levels of the F -(H 2 O) and F -(D 2 O) complexes[END_REF] it is possible to successively perform several coordinate transformations (symmetrized coordinates, flexible coordinates, curvilinear normal modes, coordinates associated with a reaction-path Hamiltonian [START_REF] Miller | Reaction path Hamiltonian for polyatomic molecules[END_REF] [START_REF] Carrington | Reaction surface Hamiltoian for the dynamics of reactions in polyatomic systems[END_REF][16]…).

Let the configuration of an N-atom molecular system be described by means of n curvilinear coordinates 𝑸 = [𝑄 ! , ⋯ 𝑄 " , ⋯ 𝑄 # ]. The standard quantum expression of the KEO in terms of n curvilinear coordinates is given by the following Laplacian operator [START_REF] Podolsky | Quantum-mechanically correct form of Hamiltonian function for conservative systems[END_REF] Eq.1b where the 𝐺 "' (𝑸) (i,j=1...n) are the contravariant components of the metric tensor, 𝑮(𝑸), and where 𝐽(𝑸) = det (𝑮(𝑸)) )!/& . 𝐽(𝑸) is also used to define the volume element, 𝑑𝜏 $ , corresponding to the given curvilinear coordinates (see Eq. 1b). It is worth noting that (i) when n<3N, the system is subjected to constraints as, for instance, in the case of reduceddimensionality models. (ii) when n=3N, the system is, of course, not constrained and, in addition, 𝐽(𝑸) is the Jacobian determinant for changing from 3N Laboratory-Fixed Cartesian coordinates to 3N curvilinear coordinates (for more details see, for instance, ref [START_REF] Gatti | Chpt. 6, Applications of Quantum Dynamics in Chemistry[END_REF]), in which case the volume element is Euclidean.

In most cases, the volume element, 𝑑𝜏 $ , is not adapted to the basis set into which the wave function is to be expanded. Indeed, the choice of a set of orthonormal basis functions very often requires to define a new volume element, 𝑑𝜏 + = 𝜌(𝑸)𝑑𝑄 ! ⋯ 𝑑𝑄 # (Eq. 2e), where 𝜌(𝑸) is a (real positive) weight function. In that case, which may also be regarded as a change in the normalization convention of the wave functions, the corresponding KEO has to be modified [START_REF] Podolsky | Quantum-mechanically correct form of Hamiltonian function for conservative systems[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] as follows:

[2][3][4][5]
𝑇 * + +𝑸, 𝜕 𝑸 - = ? 𝐽(𝑸) 𝜌(𝑸) 𝑇 * $ +𝑸, 𝜕 𝑸 -? 𝜌(𝑸) 𝐽(𝑸) = - ℏ & 2 ? 𝐽(𝑸) 𝜌(𝑸) 1 𝐽(𝑸) 1 𝜕. 𝜕𝑄 " 𝐽(𝑸) 1 𝐺 "' (𝑸) 𝜕. 𝜕𝑄 ' ? 𝜌(𝑸) 𝐽(𝑸) # '(! # "(! = - ℏ & 2 1 𝜌(𝑸) 1 𝜕. 𝜕𝑄 " 𝜌(𝑸) 1 𝐺 "' (𝑸) 𝜕. 𝜕𝑄 ' + 𝑉 ,-(𝑸) # '(! # "(! = 1 𝑓 & "' (𝑸) 𝜕 & . 𝜕𝑄 " 𝜕𝑄 ' # ".' + 1 𝑓 ! " (𝑸) 𝜕. 𝜕𝑄 " # " + 𝑉 ,-(𝑸)
and 𝑑𝜏 + = 𝜌(𝑸)𝑑𝑄 ! ⋯ 𝑑𝑄 # Eq. 2a

Eq. 2b

Eq. 2c Eq. 2d

Eq. 2e

The derivation of Eq. 2d has been known for many years. [START_REF] Nauts | Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians[END_REF]

[5][9][10][11][12][13][14]
However, to be complete and to ease the analysis of the complexity in terms of numerical operations (multiplications), an overview of the derivation of Eq. 2d is given in the appendix A.

As mentioned above, curvilinear coordinates are essential to compute ro-vibrational levels and to propagate wave packets of a molecular system efficiently. The numerical procedures for the timedependent approach (Chebychev [START_REF] Tal-Ezer | An accurate and efficient scheme for propagating the time dependent Schrödinger equation[END_REF], Runge-Kutta , Taylor expansion [START_REF] Lauvergnat | A simple and efficient evolution operator for time-dependent Hamiltonians: the Taylor expansion[END_REF]) or the timeindependent approach (iterative diagonalization techniques, such as Davidson [START_REF] Davidson | The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices[END_REF] [START_REF] Ribeiro | A Jacobi-Wilson description coupled to a block-Davidson algorithm: an efficient scheme to calculate highly excited vibrational levels[END_REF] or Lanczos [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF]), require repeated computation of the action of the potential and the KEO on a wave function or a wave packet.

When the number of degrees of freedom, n, is small (n<10) either for small systems (tetra or penta-atomic molecules) or for reduced dimensionality models, the numerical approach is efficient and easy to implement for treating the action of the KEO (as given by Eq. 2d) on a wave function. However, when the number of degrees of freedom is large, this numerical approach has to face two limitations:

(i) The number of terms becomes very large (i.e. the number of 𝑓 ! "# (𝑸), 𝑓 $ " (𝑸) and 𝑉 %& (𝑸)) in the KEO (Eq. 2d) is about n 2 /2 (more precisely, (n+1)(n+2)/2) ), and each term has to be computed and stored on the full grid. So, in our previous study on methanol in 12D [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF], the KEO is composed of 91 terms and for the largest calculation, the number of grid points was about 0.9 billion. On our computer, the KEO on the full grid could not be stored in memory (about 605 GB) and was therefore stored on a disk. (ii) The numerical calculation of the KEO can be time-consuming. Indeed, the computations of 𝑓 ! "# (𝑸), 𝑓 $ " (𝑸) and 𝑉 %& (𝑸) scale as n 3 , n 4 and n 5 , respectively (see appendix A for the expressions of these terms (Eq. A1) and a discussion of their numerical cost).

So, for instance, with a 12-core computer, about 9 days were needed to calculate the Hamiltonian (potential and KEO) on the full grid for the largest calculation of our previous 12D study on methanol.

Difficulties of the same kind are also present in the standard DVR approach where usually Eq. 2c is used [START_REF] Light | Discrete-variable representations and their utilization[END_REF][13], instead of Eq. 2d, where the number of terms to be stored on the full grid is also about n 2 /2 and where the 𝑉 %& (𝑸) has also to be calculated.

The unfavorable scaling (n 5 ) of 𝑉 %& (𝑸) and the large number of terms in the KEO's prevent an efficient use of the numerical kinetic energy operator for systems with a large number of degrees of freedom.

The main objective of the present work is to show that the use of a numerical approach is possible without storing the full on-the-grid KEO in memory or on disk, provided that the KEO is repeatedly computed on-the-fly i.e. each time the KEO action on a wave function is required. In addition, we will show that the unfavorable (n 5 ) scaling of 𝑉 %& (𝑸) can be avoided by using Eq. 2b directly [START_REF] Császár | The fourth age of quantum chemistry: molecules in motion[END_REF].

II KEO action on a wave function

The action of the KEO, 𝑇 * + +𝑸, 𝜕 𝑸 -, on a wave function, Ψ(𝑸), is the most time-consuming part of the quantum dynamics approaches (time-dependent or time-independent). Several steps are required to perform this operation: (i) Calculation of the derivatives (first and second order) of a wave function known on the grid (ii) taking into account the kinetic energy terms (the 𝑓 ! "# (𝑸), 𝑓 $ " (𝑸) and 𝑉 %& (𝑸) for Eq. 2d or G ij (Q) and 𝑉 %& (𝑸) for Eq. 2c).

IIa) derivatives of a wave function:

Here, we assume that the wave function, Ψ(𝑸), is known on the full grid, i.e. this wave function is expressed as a column vector, C ⋮ Ψ(𝑸 𝑼 ) ⋮ F, where 𝑸 𝑼 is the value of the n coordinates for the U th grid point and U runs from 1 to the total number, NQ. Then, to calculate its derivatives, ∂ " Ψ = '( ') ! , with respect to coordinate Q i on the full grid, we have to compute 01 02 ! for all 𝑸 𝑼 points and it will be noted as a column vector, C

⋮ ∂ 𝑖 Ψ(𝑸 𝑼 ) ⋮ F.
To simplify the discussion, let assume that the wave function has two coordinates only, x and y, and is expressed as a direct product with nqx and nqy grid points along the x and y coordinates, respectively (see also the discussion in the review by Light and Carrington on the DVR approach [START_REF] Light | Discrete-variable representations and their utilization[END_REF]). The wave function is thus represented on the 2D-grid as a matrix,

Ψ 33 (𝑢 4 , 𝑢 5 )
, where the double upper index, gg, means that Ψ is expanded on the grid for both coordinates, x and y. The indices, ux and uy, run from one to nqx and nqy grid points, respectively.

To compute the derivative with respect to x, 𝜕 4 Ψ = 𝜕Ψ 𝜕𝑥 , the wave function has to be partially expanded on the orthonormal basis set, B x , associated to the coordinate, x. This 1D-basis-set is composed of nbx basis functions: H𝐵 ! 4 (𝑥) ⋯ 𝐵 6 " 4 (𝑥) ⋯ 𝐵 #7 " 4 (𝑥) K. This wave function, Ψ 73 (𝑘 4 , 𝑢 5 ), where the upper indices are b (for the basis expansion on x) and g (for the grid expansion on y), is given by:

Ψ 73 (𝑘 4 , 𝑢 5 ) = 1 Ψ 33 (𝑢 4 , 𝑢 5 ) • 𝐵 6 " 4 (𝑢 4 ) • 𝜔 4 (𝑢 4 ) #8 " 9 " (! ∂ 4 Ψ 33 (𝑢 4 , 𝑢 5 ) = 1 Ψ 73 (𝑘 4 , 𝑢 5 ) • ∂ 4 𝐵 6 " 4 (𝑢 4 ) #7 "
6 " (! Eq. 3a

Eq. 3b

The first equation (Eq. 3a) is simply the partial integration over x of the product P𝐵 6 " 4 (𝑥) • Ψ(𝑥, 𝑦)R computed by means of a Gaussian quadrature. In Eq. 3a, the 𝜔 4 (𝑢 4 ) are the weights associated to the grid point 𝑢 4 . Next, the Ψ 73 (𝑘 4 , 𝑢 5 ), is transformed back on the grid (Eq. 3b), but taking into account the derivative with respect to x, hence, the use of ∂ 4 𝐵 6 " 4 (𝑢 4 ).

These two transformations can be merged into a single one using a square matrix, 𝑫 4 (!) (size nqx.nqx), which represents the first derivative operator in x. A similar procedure is also used in the well-known FBR-DVR approaches, where, T𝜔 4 (𝑢 4 ) is taken into account in 𝐵 6 " 4 (𝑢 4 )

and ∂ 4 𝐵 6 " 4 (𝑢 4 ) (see ref [START_REF] Light | Discrete-variable representations and their utilization[END_REF]). Furthermore, our approach can be easily modified when second order derivatives are required (one has then to use

∂ 44 𝐵 6 " 4 (𝑢 4 ) instead of ∂ 4 𝐵 6 " 4 (𝑢 4 ))
and the matrix representation is noted 𝑫 4 (& . Of course, it can be extended to systems with more than two degrees of freedom. In that case, the coordinate y and its corresponding index, uy, have to be viewed as a collective coordinate and a collective index, respectively.

Furthermore, this matrix, 𝑫 4 (!) , can be obtained using a pseudoinverse procedure, [START_REF] Penrose | A generalized inverse for matrices[END_REF] which is numerically more stable when nbx is much smaller than nqx and in particular when a contracted basis set is used (see the appendix B, for an explanation of the procedure).

IIb) KEO action:

In the previous subsection, we have briefly recalled how to compute the first derivatives,

∂ " Ψ = '( ') ! , and second derivatives, ∂ "# Ψ = ' " ( ') ! ') #
, of a wave function on the full grid. In other words, we know the following column vectors:

𝚿 = C ⋮ Ψ(𝑸 𝑼 ) ⋮ F, ∂ " 𝚿 = * ⋮ ∂ " Ψ(𝑸 𝑼 ) ⋮ +, ∂ "# 𝚿 = * ⋮ ∂ "# Ψ(𝑸 𝑼 ) ⋮ + Eq. 4
where, 𝑸 𝑼 is the value of the n coordinates for the U th grid point. Thus, on the grid, the KEO action on Y gives a new vector, c:

𝛘 = C ⋮ χ(𝑸 𝑼 ) ⋮ F = 1 C𝑓 & "' (𝑸 𝑼 ) ⋮ ∂ 𝑖𝑗 Ψ(𝑸 𝑼 ) ⋮ F # ".' + 1 C ⋮ 𝑓 ! " (𝑸 𝑼 ) ∂ 𝑖 Ψ(𝑸 𝑼 ) ⋮ F # " + C ⋮ 𝑉 ,-(𝑸 𝑼 )Ψ(𝑸 𝑼 ) ⋮ F Eq. 5
This is close to the procedure used to calculate the action of a potential function on a wave function in the standard FBR-DVR approach [START_REF] Light | Discrete-variable representations and their utilization[END_REF][13], yet more terms being present here.

It is worth noting, that, when the full KEO-vectors (𝑓 ! "# ,𝑸 𝑼 -, 𝑓 $ " ,𝑸 𝑼 -and 𝑉 %& (𝑸 / )), can be stored in memory, these operations are very fast because everything is known on the grid.

However, when it is not possible to store them in memory, we can calculate them on-the-fly, but this too will be time-consuming for large systems due to the n 5 scaling of 𝑉 %& (𝑸 / ), and therefore unrealistic.

To circumvent those difficulties, we now show that Eq. 2b can be used, directly and easily, instead of Eqs. 2c or 2d (see ref [START_REF] Császár | The fourth age of quantum chemistry: molecules in motion[END_REF]) and so, in particular, that the computation of the extrapotential term can be avoided.

The following steps are required, starting from the right-hand side of Eq. 2b:

(i) Multiplication by !𝜌(𝑸) 𝐽(𝑸) ⁄ :

𝛘 𝟏 = C ⋮ χ ! (𝑸 𝑼 ) ⋮ F = C ⋮ T𝜌(𝑸 𝑼 ) 𝐽(𝑸 𝑼 ) ⁄ Ψ(𝑸 𝑼 ) ⋮ F
Eq. 6

(ii) Differentiation with respect to Q j : we use the procedure of the previous subsection on 𝛘 ! , to get ∂ # 𝛘 1 .

(iii) Multiplication by 𝐺 "# (𝑸), summation over j and multiplication by 𝐽(𝑸 𝑼 ):

𝛘 &" = C ⋮ χ &" (𝑸 𝑼 ) ⋮ F = ⎣ ⎢ ⎢ ⎢ ⎡ 𝐽(𝑸 𝑼 ) ⋮ 1 𝐺 "' (𝑸 𝑼 ) • ∂ 𝑗 𝛘 ! (𝑸 𝑼 ) # '(! ⋮ ⎦ ⎥ ⎥ ⎥ ⎤ Eq. 7
For large systems, 𝐺 "# (𝑸 𝑼 ) is calculated on-the-fly, but only once, for all vectors 𝛘 &" . So, for instance, during the block-Davidson procedure, the G matrix is re-computed for each iteration.

(iv) Differentiation of 𝛘 &" with respect to Q i , by means of the procedure given in the previous subsection, to get ∂ " 𝛘 2𝑖 .

(v) Summation over i and multiplication by -ℏ $ /$4𝐽(𝑸)𝜌(𝑸) to get the full action of the KEO on the wave function:

𝛘 = C ⋮ χ(𝑸 𝑼 ) ⋮ F = ⎣ ⎢ ⎢ ⎢ ⎡ -ℏ & /T4𝐽(𝑸 𝑼 )𝜌(𝑸 𝑼 ) ⋮ 1 ∂ " χ &" (𝑸 𝑼 ) # "(! ⋮ ⎦ ⎥ ⎥ ⎥ ⎤ Eq. 8
In this approach, we need to compute only the metric tensor, G(Q), and J(Q). There is no need to compute the extrapotential term, Vep, anymore, which is one of the main differences with the standard-DVR approaches previously implemented, where Eq. 2c is used instead of Eq. 2b, so that the action of the extrapotential term, 𝑉 ,-, on the wave function has to be calculated.

Furthermore, and since the calculation of G(Q) scales as n 3 , it is possible to compute it with an on-the-fly procedure i.e. without storing it in memory on the full grid.

In the present procedure and for a single wave function, we need the following vectors on the full grid (Y or c 1 or c, ∂jc1, ∂ic2i and J, r) and one more vector for the potential. Altogether, we need 2n+4 vectors on the full grid, whereas, in our previous implementation, about n 2 /2 vectors were required.

III Application to methanol

In the present paper, we have used the same potential [START_REF] Bowman | Vibrational Levels of Methanol Calculated by the Reaction Path Version of MULTIMODE , Using an ab initio, Full-Dimensional Potential[END_REF] as in our previous study. [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF] The coordinates are defined as a series of transformations between coordinate sets [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF] (see figure 1) and are given as follows: (i) z-matrix coordinates, Qzmat: the first set of curvilinear coordinates is defined through a z-matrix (fig. 1). The corresponding transformation relates the Cartesian coordinates to z-matrix ones, Qzmat. (ii) linear transformation, giving Qsym: the three dihedral angles associated to the hydrogen of the methyl group of the z-matrix, j1, j2 and j3 are combined to give three angles, fa, fe1 and fe2 (see fig. 1). The torsion of the methyl group with respect to the OH bond is given by the first angle fa. It is important to note that, at the end of all these transformations, the coordinates used in the dynamics are closely related to coordinates associated with a reaction path Hamiltonian, RPH, or similar approaches in which the curvilinear normal modes evolve along the path.

However, in our previous [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF] and in the present studies, the curvilinear normal modes remain constant along the torsional mode. This simplified RPH version does not prevent from performing exact quantum calculation for a given potential.

Furthermore, the full 12D-basis set is expressed as a direct product of (i) a 1D-contracted basis set for the torsion, made of 24 eigenfunctions obtained from the diagonalization of a 1D-Hamiltonian. The uncontracted basis set is made of 36 sine and cosine functions and the numerical integrations were performed using 48 quadrature grid points. (ii) a 11D-Smolyak basis set for the other coordinates. More precisely, the K th 11D-basis functions, 𝐵 G !!H (𝑸), is defined as a product of 11 Harmonic Oscillator basis functions as follows:

𝐵 G !!H (𝑸) = 𝐵 6 # ! (𝑄 ! ) … 𝐵 6 ! " +𝑄 " -… 𝐵 6 ## !! (𝑄 !! )
Eq. 8

z-matrix definition

Symmetrization: Linear transfo.

Flexible transfo.

Normal modes: Linear transfo.

Normal coordinates with all coordinates except f a For instance, for a given 11D-basis function defined with this set of 𝑘 " ,

Q zmat Q cart Q sym Q flex Q NM ≡Q ! " = $ % & $ + & ( + & % ! *$ = $
[𝑘 ! … 𝑘 " … 𝑘 !! ] = [1, 1, 𝟔, 1, 𝟏𝟏, 1, 1, 1, 1, 1, 1], the ℓ " set is [ℓ ! … ℓ " … ℓ !! ] = [
0, 0, 𝟑, 0, 𝟒, 0, 0, 0, 0, 0, 0] and the sum of the ℓ " is 7. Therefore, if the parameter LB, which defines the size of the basis set is equal or larger than 7, this 11D-basis function is kept and otherwise it is discarded. It is important to note for example, that if k3=7 (instead of 6), the value of ℓ J remains the same (3). In a certain way and for each index i, the ℓ " define the number of basis functions, 𝑛𝑏 " , such as 𝑛𝑏 " (ℓ " ) = 𝑛𝑏 " (ℓ " -1) + ℓ " and 𝑛𝑏 " (0) = 1 (see table 1). This kind of selection is similar to the pruned basis set of Carrington's group [31][32]. With this non-direct product basis set, a Smolyak sparse grid can be used [START_REF] Smolyak | Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions[END_REF][34] [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF][36] [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF]; the grid size is determined by means of a Smolyak parameter, LG, defined as LG=LB+1 and the number of quadrature grid points for each HO basis set is 𝑛𝑞 " (ℓ " ) = 𝑛𝑞 " (ℓ " -1) + ℓ " . To converge the energy levels within 0.01 cm -1 , the values of LB and LG are 7 and 8, respectively.

In our previous study, [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF] the Smolyak basis set was defined in a slightly different way i.e.

𝑛𝑏 " (ℓ " ) = 𝑛𝑏 " (ℓ " -1) + 1 and 𝑛𝑞 " (ℓ " ) = 𝑛𝑞 " (ℓ " -1) + 1. Therefore, LB and LG had to be larger (LB=9 and LG=10) than those of the present study in order to converge the energy levels. However, in both studies, the basis set size (NB) and the number of grid points (NQ) are of the same order of magnitude (see table 2).

It is important to note that the diagonalization was performed using a block Davidson procedure [START_REF] Ribeiro | A Jacobi-Wilson description coupled to a block-Davidson algorithm: an efficient scheme to calculate highly excited vibrational levels[END_REF]. However, during the Davidson iterations and with the present approach, the intermediate Hamiltonian matrix is not Hermitian (due to the numerical approximation of the differentiation in the steps ii and iv of the KEO action). This is why we have used a Lapack subroutine to diagonalize this non-symmetric matrix.

In table 2, we present some torsional energy levels from the present study and from our previous ones. The comparison is almost perfect; the energy differences are not larger than 0.01 cm -1 for all levels. Furthermore, due to the symmetry, some levels are degenerate (those corresponding to the following experimental values: 9.1, 208.9, 510.3 and 751.0 cm -1 ) and the computed energy components are not strictly equal due to approximations in the numerical integration. However, the largest difference between the two components is only 0.005 cm -1

for the level at 741.78 cm -1 . It will be noted that in our previous study, this was not possible (because it required too much time), except for the zero-point energy (ZPE) to obtain the energy levels, for the largest calculation (LB=9, LG=10), which already took 20 days. Now with the present approach, in which the metric tensor is computed for each Hamiltonian operator action on a wave function, the calculation takes only 8 days for all the treated states.

Furthermore, and as expected, the memory requirement is much smaller (about five times), for grids of similar sizes.

VI Conclusion and perspectives

From a numerical point of view, the approach developed in the present paper opens the way to using sophisticated coordinates with a numerical on-the-fly KEO for large systems without introducing approximations [START_REF] Stare | Fourier Grid Hamiltonian Method for Solving the Vibrational Schrödinger Equation in Internal Coordinates: Theory and Test Applications[END_REF][38] (such as using constant or Taylor-expanded metric tensors or neglecting the extra-potential term). Indeed, we show that the numerical procedure implemented here enables to reproduce the previous results exactly [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF] with less computational resources (memory and cpu time).

However, that the memory scales as 2n (n being the number of degrees of freedom) is clearly a limitation, although not as strong as the n 2 /2 scaling in our previous implementation. To address this problem, one can use the data structure of the Smolyak representation (the sum of small direct product bases or grids) in order to avoid almost completely, except for the potential, to store vectors on the full grid. In other words, we need to merge the Smolyak scheme with the on-the-fly numerical action of the KEO's on the wave functions.

There is also another issue worth analyzing in more detail. Indeed, in order to deal with operators such as ' '1 , the use of 𝑫 4 (!) (obtained from Eq. 3) may give rise to a non-Hermitian

Hamiltonian matrix [START_REF] Wei | Discrete variable representations of complicated kinetic energy operators[END_REF][25] when %. %' ! (Eq. 2b or Eq. 2c), is applied. Unfortunately, when using (Eq. 2b or Eq. 2c), the non-hermiticity cannot be reduced by increasing the number of grid points as usual. For the calculation of eigenlevels, the use of a diagonalization procedure for non-symmetric matrices solves this problem, although some spurious eigenvectors may crop up, which can be easily discarded upon inspection. Furthermore, the consequences of the loss of hermiticity of the Hamiltonian operators during the time-dependent approach are not clear. One might probably try differentiation on the left %. %' " ⃖,,,,, [START_REF] Wei | Discrete variable representations of complicated kinetic energy operators[END_REF] instead of the usual differentiation on the right, %. %' ! .

Moreover, it is worth pointing out that this non-hermicity issue is almost absent when all derivative operators act directly on the wave function (Eq. 2d). More precisely, the Hamiltonian matrix may be slightly non-Hermitian, but this can be numerically suppressed, by increasing the number of grid points without increasing the number of basis functions.

Anyway, and whatever the remaining difficulties, numerical on-the-fly implementations as developed in the present paper, will certainly contribute to open the way to rigorous quantum computations of systems with many degrees of freedom. The key ingredients are the contravariant components of the metric tensor, 𝐺 "' (𝑸). Its matrix representation, G can be obtained as follows:

𝑮 = H𝐺 "' (𝑸)K = H𝑔 "' (𝑸)K )! (i,j = 1, ....n) Eq. A2
where, 𝑔 "' (𝑸) are the covariant components of the metric tensor, associated to the matrix, g, with

𝑔 "' (𝑸) = 1 𝑚 L 𝜕𝑥 L (𝑸) 𝜕𝑄 " 𝜕𝑥 L (𝑸) 𝜕𝑄 ' JM L(!
Eq. A3

In the previous equation, N is the number of atoms and ml is the atomic mass associated with the l th Cartesian coordinates. Furthermore, when all degrees of freedom are considered (no constraints), n is equal to 3N.

From the two previous equations, one can see that the number of operations (i.e. multiplications) to build the G matrix and therefore the 𝑓 & "' terms (Eq. A1a), scales as O(n 3 ): indeed, 6N

operations to build each of the n 2 component, 𝑔 "' (Eq. A3), and about n 3 operations for the matrix inversion (Eq. A2).

Next, to compute the 𝑓 ! " , one needs several ingredients: (i) the G matrix which is already known. (ii) 𝜌(𝑸), the weight function associated with the chosen volume element and the derivative of 𝜌(𝑸). In all applications, the expression of 𝜌(𝑸) is really simple (product of 1D-functions). Therefore, the numerical cost associated to the 𝜌(𝑸) contributions is small. (iii) the n first derivatives of the G matrix, which are the costly contributions (see Eqs. A4, A5).

The derivation of the G matrix reads:

𝜕𝑮 𝜕𝑄 6 = -𝑮 𝜕𝒈 𝜕𝑄 6 𝑮 (k = 1, ....n) Eq. A4
and the n derivatives of the components of g are obtained through the differentiation of Eq. A.3 and are given as follows:

𝜕𝑔 "' (𝑸) 𝜕𝑄 6 = 1 𝑚 L { 𝜕 & 𝑥 L (𝑸) 𝜕𝑄 6 𝜕𝑄 " 𝜕𝑥 L (𝑸) 𝜕𝑄 ' + 𝜕𝑥 L (𝑸) 𝜕𝑄 " 𝜕 & 𝑥 L (𝑸) 𝜕𝑄 6 𝜕𝑄 ' | JM L(! Eq. A5
From the previous equations, one can see that the number of operations scales (i) as 2n 4 : 2n matrix multiplications (n 3 operations) for the n derivatives for Eq. A4. (ii) as 3x3N for each term

03 !$ (𝑸)
02 % (n 3 terms). Therefore, the numerical cost to get all the 𝑓 ! " scales as O(n 4 ).

The computation of the extrapotential term, 𝑉 ,-(𝑸), is the costliest among the Eq. A1 terms. Indeed, all elements from Eq. A1c are already known, except those implying 𝐽(𝑸) and its first and second derivatives with respect to the curvilinear coordinates. 𝐽(𝑸) is expressed as a function of the determinant of g as follows:

𝐽(𝑸) = Tdet (𝒈) Eq. A6

The numerical cost scales as O(n 3 ), for det(g). Furthermore, the computation of the first and the second derivatives of 𝐽(𝑸) requires, respectively, the first (Eq. A5) and the second (obtained from the differentiations of Eq. 5 with respect to 𝑄 ℓ ))-derivatives of g. Thus, for all i and j values, the numerical cost of

0 & 3(𝑸) 02 ! 02 ! and 0 & $(𝑸)
02 ! 02 $ scale as O(n 5 ). Therefore, the numerical cost to get the 𝑉 ,-(𝑸) scales as O(n 5 ).

The method developed in the present paper is valid for the quantum treatment of any system with n degrees of freedom. It goes without saying that it can be used in the particular case of a N-atom molecular system (constrained n<3N or not constrained n=3N) by introducing the following n generalized coordinates 𝑸 = [𝒒, 𝚯, 𝑿 𝑪𝑴 ], where 𝒒 = [𝑞 ! ⋯ 𝑞 #)Q ] are internal (shape) coordinates, 𝚯 are the three Euler angles that orient a Body-fixed frame with respect to the Space-fixed frame and 𝑿 𝑪𝑴 are the three Cartesian coordinates that determine the position of the center of mass of the molecular system in the Laboratory-fixed frame.

B) Derivatives of a wave function, pseudoinverse procedure:

The procedure described in section IIa (Eq. 3) works well in most cases and, when the numbers of basis function, nb, and grid points, nq, are equal, this procedure is equivalent to the one used in the DVR approach. However, in some circumstances and in particular for contracted basis sets, the numerical precision of the matrix representation on the grid of the first, 𝑫 4 (!) , and second, 𝑫 4 (&) , derivative operators is not satisfactory.

To show this feature, we will compute the differences between the direct computation of the first derivative of the k th basis function 𝐵 6 (𝑥) at a given point xu,

0I % (4) 04 ƒ 4 '
= 𝜕 4 𝐵 6 (𝑥 9 ) using the usual relations of the orthogonal polynomials or the basis functions and the two numerical approaches, the standard one using Eq. 3 and another one using a pseudoinverse procedure. [START_REF] Penrose | A generalized inverse for matrices[END_REF] Both approaches proceed in two steps, the second one is the same in both cases and corresponds to Eq. 3b, while in the first step, a vector known on the grid, Ψ for the basis expansion and for the grid expansion) as follows:

Ψ 7 = 𝑴 73 • Ψ 3 Eq. A7
where 𝑴 73 is a matrix (size nb.nq), which is computed by using two different procedures: (i) In the standard approach, 𝑴 73 is given as in Eq. 3a (𝑴 73 (𝑘, 𝑢) = 𝐵 6 (𝑥 9 ) • 𝜔 4 (𝑢)).

(ii) In the pseudoinverse approach, 𝑴 73 is obtained from solving a system of linear equations, where the k th -basis function on the grid, C ⋮ B 6 (𝑥 9 ) ⋮ F, is Ψ 3 so that the corresponding, Ψ 7 , is a vector for which all components are zero except C 6 , which is equal to one. All vectors C ⋮ B 6 (𝑥 9 ) ⋮ F can be put together to form a matrix 𝑴 37 (size nq.nb) so that all the corresponding Ψ 7 vectors form the identity matrix, 𝑰 #7•𝒏𝒃 (size nb.nb).

𝑰 #7•𝒏𝒃 = 𝑴 73 • 𝑴 37 (𝑴 37 ) E = 𝑴 73 • (𝑴 37 • (𝑴 37 ) E )

Eq. A8a

Eq. A8b

When nb < nq, 𝑴 37 is not invertible, so that the pseudoinverse approach is required. [START_REF] Penrose | A generalized inverse for matrices[END_REF] In the present study, we multiply on the right by (𝑴 37 ) E (the transpose of 𝑴 37 ) to get Eq A8b. Then, the matrix (𝑴 37 • (𝑴 37 ) E ) on the right-hand-side of Eq. A8b is diagonalized and expressed as (𝑃 E • 𝐷 • 𝑃). The pseudoinverse of (𝑴 37 • (𝑴 37 ) E ) is +𝑃 Š • 𝐷 ‹ )! • 𝑃 Š Ewhere the tilde implies that only the non-zero eigenvalues, 𝐷 ‹ , and their corresponding eigenvectors, 𝑃 Š are taken into account.

In our case, the number of non-zero eigenvalues is the number of basis functions, nb. In that case, the matrix 𝑴 73 , which is also the pseudoinverse of 𝑴 37 , reads:

𝑴 73 = (𝑴 37 ) E • +𝑃 Š • 𝐷 ‹ )! • 𝑃 Š E - Eq. A9
For instance, in the present study, the basis set associated with the torsion is a Fourier series with 36 terms (cosine and sine functions) contracted with 24 basis functions. Furthermore, the number of grid points is 48 for both contracted and un-contracted basis sets. The largest difference (for k=1...nb and u=1...nq) between 𝜕 4 𝐵 6 (𝑥 9 ) obtained from the usual relations of the orthogonal polynomials and the corresponding values obtained from the standard and the pseudoinverse approaches are given in table 3. The analysis of the results shows that when the Fourier basis set is not contracted the largest error on 𝜕 4 𝐵 6 (𝑥 9 ) is of the same order of magnitude for the standard and the pseudoinverse procedure, while for the contracted basis set, the error obtained from the standard procedure is about 10 9 larger than the one obtained from the pseudoinverse approach.

  [START_REF] Kemble | The Fundamental Principles of Quantum Mechanics[END_REF][3][4][START_REF] Gatti | Chpt. 6, Applications of Quantum Dynamics in Chemistry[END_REF]:𝑇 * $ +𝑸, 𝜕 𝑸 -=and 𝑑𝜏 $ = 𝐽(𝑸)𝑑𝑄 ! ⋯ 𝑑𝑄 # Eq. 1a

4 "

 4 This transformation is essential to preserve the symmetry of the molecule, [29][30]. (iii) flexible coordinates, Qflex: the i th flexible coordinates, 𝑄 @A,is given by 𝑄 B5C " -𝑄 D-E " (𝜙 F ), where 𝑄 D-E " (𝜙 F ) is the optimal value of the 𝑄 B5C "coordinates along the torsional path (in fa). Furthermore, these optimal values are expanded in Fourier series, using the three equivalent minima and three equivalent transition states. (iv) curvilinear normal modes, QNM: these coordinates are defined as a constant linear transformation between 11-flexible ones (all coordinates except fa). The 11x11 hessian and kinetic matrices used to get the curvilinear normal modes are obtained from averaging hessian and kinetic matrices at the three equivalent minima.

Fig 1 :

 1 Fig 1: Coordinate definitions adapted to the symmetry and reducing the coupling between modes.

  with nq grid points) is projected onto the nb basis functions,Ψ 7 = C ⋮ C 6 ⋮ F, (the upper indices b and g stand, respectively,

table 1 ) : Table 1 :

 1:1 where 𝐵 6 ! " +𝑄 "is the ki th Harmonic Oscillator basis function associated with the i th coordinate, Q i . Then, the 𝐵 G !!H (𝑸) are selected by means of the parameters, ℓ " and LB, such as 0 ≤ ≤ 𝐿 I and a relation (for all i) between the ℓ " and the 𝑘 " given in the following table (relation between the 𝑘 " , the ℓ " and 𝑛𝑏 " for all values of index i.

	!! ∑ ℓ " "(!
	𝑘 "	1	2	3 4 5 6 7 8 ... 11 12 ... 16 17 ... 22 23 ... 29 30 ...37
	𝑛𝑏 "	1	2	4	7	11	16	22	29	37
	ℓ "	0	ϕ 1 1	2	3	4	5	6	7	8
	C O C 2 ?9"< $		ϕ 3	ϕ 2
	H O 2 ?9"< (	C 2 ?9"< %
	H $ C 2 ?9"< A		O 2 ?9"< B	H & $
	H $ C 2 ?9"< C		O 2 ?9"< D	H & (
	H $ C 2 ?9"< $E		O 2 ?9"< $$	H & %

Table 2 :

 2

		Exp.	ELVIBROT previous work	ELVIBROT present work
		Exp.	LB=6, LG=8 LB=9, LG=10 LB=6, LG=7 LB=7, LG=8
		/	11088.76	11088.50	11088.55	11088.51
		9.1	9.14, 9.15		9.15	9.15
		208.9	205.34		205.34	205.34
	Energy levels	294.7	290.69		290.69	290.69
		353.0	347.49		347.49	347.49
		510.3	503.54		503.53	503.53
		751.0	741.80		741.78	741.78
	Basis set size: NB	/	297 024	4 031 040	920 064	3 196 008
	Grid size: NQ	/	88 584 912 892 360 032 134 737 200 588 600 720
	Memory requirement	/	60 GB	605 GB on disk	28 GB	123 GB

First column, torsional energy levels (in cm -1 ) of methanol; second column, experimental levels; third and fourth columns, levels and parameters from our previous study

[START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF]

; fifth and sixth columns, levels and parameters from the present work using ELVIBROT. The levels at 9.1, 208.9, 510.3 and 751.0 are degenerate.

Table 3 :

 3 Largest difference between 𝜕 4 𝐵 6 (𝑥 9 ) obtained from the usual relations of the orthogonal polynomials and the corresponding values obtained from the standard (columns 2, 4) and the pseudoinverse (columns 3, 5) for the un-contracted (columns 2,3) and contracted (columns 4,5) basis sets.

	largest error	nb=36 nq=48	contracted, nb=24 nq=48
		Standard	Pseudoinverse	Standard	Pseudoinverse
	𝜕 4 𝐵 6 (𝑥 9 )	5. 10 -14	2. 10 -14	6. 10 -5	2. 10 -14

Acknowledgements

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

A) Derivation of the KEO and its numerical cost:

The derivation of the numerical kinetic energy operator has been known for many years and there are several ways to obtain the expression given in Eq. 2d. [START_REF] Nauts | Momentum, quasi-momentum and hamiltonian operators in terms of arbitrary curvilinear coordinates, with special emphasis on molecular hamiltonians[END_REF][5] [START_REF] Meyer | General Internal Motion of Molecules, Classical and Quantum-Mechanical Hamiltonian*[END_REF][10] [11][12][13] [START_REF] Strobusch | A general nuclear motion Hamiltonian and non-internal curvilinear coordinates[END_REF] In the present study, the 𝒙(𝑸) approach is used, where Q is the set of n curvilinear coordinates and x is the set of Cartesian coordinates and in this case the 𝑓 & "' (𝑸), 𝑓 ! " (𝑸) and 𝑉 ,-(𝑸) are given by the following expressions:

𝜌(𝑸) y

Eq. A1a

Eq. A1b

Eq. A1c