Emanuele Marsili

Federica Agostini
email: federica.agostini@universite-paris-saclay.fr

André Nauts

David Lauvergnat
email: david.lauvergnat@universite-paris-saclay.fr

Quantum Dynamics with curvilinear coordinates: Models and Kinetic Energy Operator

Keywords: Quantum dynamics, curvilinear coordinates, kinetic energy operator

In order to simplify the numerical solution of the time-dependent or time-independent Schrödinger equations associated with the atomic and molecular motions, the use of well-adapted coordinates is essential. Usually, this set of curvilinear coordinates leads to a Hamiltonian operator as separable as possible. Although their corresponding kinetic energy operator (KEO) expressions can be derived analytically for small systems or special kinds of coordinates, a numerical and exact approach allows to compute them in terms of sophisticated curvilinear coordinates. Furthermore, the numerical approach enables to easily define reduced dimensionality or constrained models. We present here a recent implementation of this numerical approach that allows nested coordinates transformations and therefore leading to a large flexibility in the definition of the curvilinear coordinates. Furthermore, this implementation has no-limitations in terms of number of atoms or coordinate transformations. The quantum dynamics of the cis-trans photoisomerization of part of the retinal chromophore illustrates the building of the coordinates and KEO part of a 3D-model.

Introduction

From the general point of view of quantum mechanics any set of coordinates can be used to construct the Hamiltonian operator and, in particular, Cartesian coordinates, in which case the kinetic energy operator (KEO) has a simple form. Although the expression of the potential energy surfaces can be involved, its complexity can be overcome easily. On the one hand, for fitted force-fields, the Cartesian coordinates are transformed into geometrical parameters (distances, angles, etc.) required to calculate the potential. [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF] This transformation disappears completely with direct-dynamics or ab initio dynamics, with or without quantum effects, for which the potential is easily calculated on-the-fly in terms of Cartesian coordinates. [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF][START_REF] Iftimie | Ab initio molecular dynamics: Concepts, recent developments, and future trends[END_REF] On the other hand, standard quantum simulations of atomic motion in molecules (with time dependent or independent approaches) are usually performed by means of well-chosen coordinates, such as normal modes or curvilinear coordinates. There are several reasons for this choice:

-To analyze experimental data of a spectrum or a chemical (physical) process: The analysis of particular motions along some coordinates (distances, angles, etc.) is simpler if the dynamics is performed directly using these coordinates. Of course, one can always transform the coordinates used for the dynamics into those required for the analysis as done in classical molecular dynamics. However, since the wavepacket representing the quantum mechanical state of the system is delocalized, the coordinate transformation has to be applied for all grid points and at each time step. Therefore, this operation could be more time consuming than performing a coordinate transformation only once for the Hamiltonian or, more precisely, for the KEO.

-To use an effective model: In this case, the effective model is often a reduced dimensionality model with a well-chosen set of coordinates. For instance, models with a single torsional angle can be used to describe a hindered rotation of two fragments. Most of the times, the parameters of the model are adjusted based on experimental data. This approach has been intensively used in the quantum treatment of atomic motion in ro-vibrational spectroscopy of isolated molecules. [START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF][START_REF] Perrin | Recent progress in the analysis of HNO3 spectra[END_REF][START_REF] Lauvergnat | New assignments in the torsional spectrum of[END_REF] -To reduce the basis set or grid sizes: Most of the times, the multidimensional basis set is expressed as a direct-product of 1D-basis sets or grids (or low dimensional basis sets). Therefore, to reduce the size of the full basis, the correlation or the coupling between coordinates must be as small as possible. In a way, it is equivalent to finding coordinates that give an almost separable Hamiltonian. [START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF] Indeed, a well-chosen set of curvilinear coordinates reduces the coupling between modes, even though one must find a balance between the kinetic and the potential couplings. It is important to note, that when the multidimensional basis set is not expressed as a direct-product but as a sum of localized multidimensional functions (such as "moving" gaussians), one can work efficiently in Cartesian coordinates as done with DD-vMCG [START_REF] Worth | A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets[END_REF] (direct dynamics variational multi-configuration Gaussian wavepacket), multiple spawning approaches [START_REF] Martinez | Multi-electronic-state molecular dynamics: A wave function approach with applications[END_REF][START_REF] Curchod | Ab Initio Nonadiabatic Quantum Molecular Dynamics[END_REF] and other similar approaches.

Unfortunately, the analytical expression of the KEO can be derived relatively easily only for some kinds of coordinates such as Jacobi or polyspherical coordinates. [START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF][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] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[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][START_REF] Wang | A simple method for deriving kinetic energy operators[END_REF] For more sophisticated coordinates, the analytical derivation of the KEO can be cumbersome, in particular for molecular systems with constraints. [START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF][START_REF] Justum | One-dimensional quantum description of the bending vibrations of HCN/CNH[END_REF] Therefore, some authors have proposed to derive the KEO numerically. [START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][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] Luckhaus | 6D vibrational quantum dynamics: Generalized coordinate discrete variable representation and (a)diabatic contraction[END_REF][START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Yurchenko | Theoretical ROVibrational Energies (TROVE): A robust numerical approach to the calculation of rovibrational energies for polyatomic molecules[END_REF][START_REF] Matyus | Toward black-box-type full-and reduced-dimensional variational (ro)vibrational computations[END_REF] In the next section, we will present the general strategy to obtain KEO in terms of any set of curvilinear coordinates with or without constraints. Then, we will show how to adapt this procedure to get a numerical but exact KEO with several nested coordinate transformations. This work can be considered as an overview of some of our previous studies on numerical KEO, [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Lauvergnat | Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates[END_REF][START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF] where a more detailed presentation of the coordinate transformations is given. As an illustration, the cis-trans photo-isomerization process of the 2cis-penta-2,4-dieniminium cation (PSB3), which is a minimal model for the retinal chromophore of Rhodopsin, will be simulated with a quantum dynamics propagation scheme using a 3D-reduced dimensionality model. [START_REF] Marsili | Two-State, Three-Mode Parametrization of the Force Field of a Retinal Chromophore Model[END_REF][START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF] 2. Coordinates and kinetic energy operator 2a) General expressions Let a molecular system be composed of N atoms, each denoted Al, of mass, ml, and in the laboratory (L) frame, its Cartesian coordinates being 𝑿

""⃗ ! " = [𝑥 ! " 𝑦 ! " 𝑧 ! "
]. These Cartesian atomic coordinates can be projected onto the body-fixed (BF) or molecular frame. The relation between the BF Cartesian coordinates, 𝑿 ""⃗ #$ " , and the L coordinates is given below (see also fig 1):

𝑿 ""⃗ ! " = 𝑿 ""⃗ ! %&' + 𝓡(𝛉)𝑿 ""⃗ #$ "
Eq. 1

where 𝑿 ""⃗ ! %&' is the molecular center-of-mass (COM) and 𝓡(𝛉) is the Euler rotation matrix associated to the three Euler angles, 𝛉 = [𝛼, 𝛽, 𝛾], or two angles in the case of a diatomic molecule. For a given molecular system, this expression enables us to split the motion of a molecule into a deformation part (𝑿 ""⃗ #$ " ), a rotation part (𝛉) and the COM motion (𝑿 ""⃗ ! %&' ). [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] Chapuisat | Exact quantum molecular Hamiltonians[END_REF] Furthermore, the total angular momentum can be easily expressed in terms of the three Euler angles [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] Chapuisat | Exact quantum molecular Hamiltonians[END_REF] and, when the molecule is isolated, the COM can be separated out.

The way the BF frame is attached to the molecule is crucial. (i) One can use two vectors, 𝑹 ""⃗ ( and 𝑹 ""⃗ ) , so that the vector 𝑹 ""⃗ ( is along the z-axis in the BF frame and (𝛽, 𝛼) are associated to the two spherical angles of this vector in the L frame. The vector 𝑹 ""⃗ ) lies in the xz-BF plane, and 𝛾, the third Euler angle, completes the 3Drotation of the BF frame with respect to the laboratory frame. These two vectors can be defined as atomic vectors (𝐴 ( 𝐴 ) """"""""""⃗ , and 𝐴 ( 𝐴 * """"""""""⃗ or 𝐴 ) 𝐴 * """"""""""⃗ ) as in a z-matrix or as vectors between atoms or center-of-mass as in Jacobi or polyspherical 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] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF] Fig. 1 shows an illustration of the orientation a triatomic molecule in the different frames. (ii) One can use an Eckart frame [START_REF] Eckart | Some Studies Concerning Rotating Axes and Polyatomic Molecules[END_REF][START_REF] Mátyus | Vibrational energy levels with arbitrary potentials using the Eckart-Watson Hamiltonians and the discrete variable representation[END_REF][START_REF] Pesonen | Eckart frame vibration-rotation Hamiltonians: Contravariant metric tensor[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][START_REF] Szalay | Understanding nuclear motions in molecules: Derivation of Eckart frame rovibrational Hamiltonian operators via a gateway Hamiltonian operator[END_REF][START_REF] Lauvergnat | Numerical and exact kinetic energy operator using Eckart conditions with one or several reference geometries: Application to HONO[END_REF], and in this case all atomic vectors are involved in the definition of the Euler rotation matrix. This enables one to separate more efficiently the vibrational from the rotational motions, although in general the analytical expression of 𝓡(𝛉) is unknown. (iii) One can use a more unusual embedding in which the z-axis of the BF frame bisects or trisects several vectors. [START_REF] Szidarovszky | The role of axis embedding on rigid rotor decomposition analysis of variational rovibrational wave functions[END_REF]

Fig 1: Illustration of the body-fixed frame orientation with respect to the laboratory frame for a triatomic molecule (A1-A2-A3) defined with the help of a z-matrix. The space-fixed (SF) frame is an intermediate frame parallel to the LF frame but centered on the center-of-mass.

Generally, the BF Cartesian coordinates are expressed as functions of n-internal curvilinear coordinates, 𝑸, with 𝑛 = 3𝑁 -6 (or 𝑛 = 1 for a diatomic molecule). Usually, the 𝑸s are atomic distances or vector norms (for Jacobi like coordinates), valence angles (range ]0, p[) and dihedral angles (2p periodic). However, more sophisticated coordinates can be used (see below). In this framework, we can define 3𝑁-curvilinear coordinates, q, which include the n-internal curvilinear coordinates, the Euler angles and the three coordinates associated with the molecular center-of-mass:

𝒒 = >𝑸, 𝛉, 𝑿 ""⃗ ! %&' ? Eq. 2
One of the major difficulties in the use of curvilinear coordinates is the potentially cumbersome analytical expression of the kinetic energy operator,𝑻 A . Fortunately, for some particular kinds of curvilinear coordinates, such as Jacobi or polyspherical coordinates, one can find analytical expressions for molecular systems of any size. In those types of coordinates, it is assumed that the BF frame is attached to two vectors only (𝑹 ""⃗ ( and 𝑹 ""⃗ ) ). For polyspherical coordinates (which include Jacobi or Radau coordinates ...), the KEO analytical expression is always expressed as a sum of products of one-dimensional functions or operators. This feature is numerically very interesting in quantum dynamics because the KEO action on a wave function or a wavepacket does not require a multidimensional grid explicitly. Therefore, those coordinates are heavily used in many applications or approaches such as the multi-configurational time-dependent Hartree (MCTDH) method

A 2 z BF x BF z SF x SF Θ z L x L X L COM X L 1 X BF 1 A 3
A 1 [START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF][START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF][START_REF] Meyer | Multidimensional Quantum Dynamics[END_REF], the Vibrational Configuration Interaction (VCI) or Vibrational Coupled Cluster (VCC) methods. [START_REF] Klinting | Vibrational Coupled Cluster Computations in Polyspherical Coordinates with the Exact Analytical Kinetic Energy Operator[END_REF] Furthermore, some of the authors were able to develop a numerical code, TANA, [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][START_REF] Lauvergnat | ElVibRot-TnumTana fortran quantum dynamics and kinetic energy operator codes[END_REF] to get analytical expressions for molecular systems of any size with a lot of freedom in the selection of the vector parametrization and the sub-systems (local body fixed frames).

For other embeddings or other sets of coordinates, the analytical expressions are more involved. It is therefore that, in the literature, they have been obtained for small molecular systems only. The KEO, 𝑻 A + , is related to the massweighted Laplacian associated to the 3N coordinates 𝒒 and with a Euclidean (E) volume element (𝑑𝜏 + ), 𝑻 A + reads:

𝑻 A + = - ℏ ! )
-. / (𝒒) "# 3 34 $ . / (𝒒)5 6 $% (𝒒)

3 34 % &' $,%)# 𝑑𝜏 + = 𝐽 E (𝒒)𝑑𝑞 ( ⋯ 𝑑𝑞 *7
Eq.

where, 𝐽 E (𝒒) and 𝐺 I 89 are, respectively the Jacobian of the coordinate transformation from the curvilinear to the Cartesian coordinates and the contravariant components of the metric tensor associated to the coordinates 𝑞 8 and 𝑞 9 .

When the separation in terms of internal coordinates, Euler angles and position of the center-of-mass is considered (Eq. 2), the metric tensor, 𝑮 K , can be factorized into a product of three matrices: [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] Chapuisat | Exact quantum molecular Hamiltonians[END_REF][START_REF] Meyer | General Internal Motion of Molecules, Classical and Quantum-Mechanical Hamiltonian*[END_REF]]

𝑮 K (𝒒) = L 𝟏 :×: 𝛀(𝛉) <,>( 𝟏 *×* O • 𝑮(𝑸) • L 𝟏 :×: 𝛀(𝛉) >( 𝟏 *×* O Eq.
The matrix-vector product between (-𝑖ℏ)𝛀(𝛉) >( and L 𝜕 ? 𝜕 @ 𝜕 A

O gives the three components of the total angular momentum operator

L 𝐽 T B 𝐽 T C 𝐽 T D O.
From the factorization of the metric tensor (Eq. 4), the Jacobian can be factorized as a product of deformation part, 𝐽(𝑸), a rotational part, 𝑠𝑖𝑛(𝛽), and a constant for the translation part. Furthermore, the KEO can be expressed in terms of 𝑮(𝑸), the components of the total angular, the 3 3E $ , and the kinetic energy operator of the center-of-mass. For instance, the deformation part of the KEO or more precisely its expression at J=0 is given by:

𝑻 A FGH + = - ℏ ! ) -.(𝑸) "# 3 3E $ .(𝑸)5 $% (𝑸) 3 3E % * $,%)# 𝑑𝜏 FGH + = 𝐽(𝑸)𝑑𝑄 ( ⋯ 𝑑𝑄 :
Eq.

The KEO relations (Eq. 3 and Eq. 5) look similar, but in the latter expression the double sum runs over the n-internal coordinates, 𝑸, and 𝑮(𝑸) depends only on the internal coordinates.

There are several ways to get the expression of 𝑮. One can start with 𝑸(𝑿 ""⃗ #$ " ), which gives directly the deformation part, however, in this case, 𝑮 is expressed as a function of 𝑿 ""⃗ #$ " instead of 𝑸. Otherwise, using the expression 𝑿 ""⃗ #$ " (𝐐), one calculates the covariant components, 𝑔 89 , of the metric tensor, rather than 𝑮, but directly as function of 𝑸:

𝒈 = L 𝑺 J×J 𝑪 J×* 𝒕 𝟎 𝑪 *×J 𝑰 *×* 𝟎 𝟎 𝟎 𝑴 *×* O Eq.
where, 𝑺 J×J , 𝑪 *×J , 𝑰 *×* and 𝑴 *×* are, respectively, the deformation, the Coriolis, the inertia, and the translational parts of the metric tensor. The expressions of those matrices are functions of 𝑿 ""⃗ #$ " (𝐐) and their derivatives (

3𝑿 MM⃗ +, - (𝐐)
3E $ ) with respect to the internal coordinates, 𝑄 8 (see below in the next section for their expressions). Then, the expression of 𝑮 is simply the inverse of 𝒈. Furthermore, the deformation part of the Jacobian, 𝐽(𝑸), can be obtained from 𝒈 as well: 𝑮 = 𝒈 >𝟏 𝐽(𝑸) = _det(𝒈) Eq. 7

The expression of the deformation KEO (Eq. 5) can be used directly, however, most of the times, the Euclidean volume element is not adapted to the basis set used to solve the Schrödinger equation. Therefore, one uses a volume element, 𝑑𝜏 FGH Q = 𝜌(𝑸)𝑑𝑄 ( ⋯ 𝑑𝑄 : , so that the multidimensional basis set is orthonormal. The expression of the KEO must be accordingly adapted to this new volume element: [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] Chapuisat | Exact quantum molecular Hamiltonians[END_REF]

𝑻 A FGH Q = d 𝐽(𝒒) 𝜌(𝒒) 𝑻 A FGH + d 𝜌(𝒒) 𝐽(𝒒) = - ℏ ! ) -Q(𝑸) "# 3 3E $ Q(𝑸)5 $% (𝑸) 3 3E % * $,%)# + 𝑉 GR (𝒒) = - ℏ ! ) -5 $% 3 ! 3E $ 3E % * $,%)# - ℏ ! ) -S 3T:(Q) 3E $ 5 $% U 35 $% 3E $ V 3 3E % * $,%)# + 𝑉 GR (𝒒) 𝑑𝜏 Q = 𝜌(𝑸)𝑑𝑄 ( ⋯ 𝑑𝑄 :
Eq. 8

In the previous expression, 𝑉 GR (𝒒) is a scalar term, part of the KEO. It is often called extra-potential [START_REF] Chapuisat | Exact quantum molecular Hamiltonians[END_REF] or pseudo-potential term. [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] It is worth noting, that with this new volume element, the scalar operators (such as the potential or the dipole moments) remain unchanged.

One of the major interests in the use of curvilinear coordinates is the possibility to a define reduced dimensionality model for which only a subset of coordinates (the active ones) are present in the dynamics, while the other coordinates (the inactive ones) are neglected. Clearly, the number of active coordinates (n), is smaller than 3N-6. Those kinds of models have been used intensively in ro-vibrational spectroscopy of floppy molecular systems to calculate or analyze spectra. Usually, the definition of the potential is not a problem. However, one has to derive the KEO with extreme care to avoid non-Hermitian issues. For instance, in the third line of Eq. 8, one cannot neglect terms associated to the inactive coordinates. Instead, one has to define a reduced dimensionality metric tensor, 𝑮 𝑹𝑫 , and plug it into Eq. 8 in which n is the number of active coordinates. There are several ways to define the 𝑮 𝑹𝑫 from: (i) a mathematically constrained model. [START_REF] Justum | One-dimensional quantum description of the bending vibrations of HCN/CNH[END_REF][START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][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] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF][START_REF] Nauts | HAMILTONIANS FOR CONSTRAINED N-PARTICLE SYSTEMS[END_REF][START_REF] Chapuisat | A general property of the quantum mechanical Hamiltonians for constrained systems[END_REF] The constraints have to be implemented rigorously in the Cartesian coordinates. Then, the covariant components of the reduced dimensionality metric tensor, 𝒈 𝑹𝑫 , can be obtained from the usual expressions, and 𝑮 𝑹𝑫 is just the inverse of 𝒈 𝑹𝑫 . Usually, the overall rotation is included in the active part of the metric tensor. Several types of constraints can be defined: (a) rigid constraints, where inactive coordinates are held fixed, or (b) flexible constraints, where inactive coordinates are defined as functions of the active coordinates.

(ii) a submatrix of 𝑮. [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF][START_REF] Avila | Nonproduct quadrature grids for solving the vibrational Schrödinger equation[END_REF] In this case, one can calculate only the exact contravariant metric components associated to the active coordinates. For instance, in the well-known book by Wilson-Decius and Cross, [START_REF] Wilson | Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra[END_REF] one can find tabulated contravariant metric components associated to usual coordinates (distances, valence angles, dihedral angles). However, when the coordinates are more involved, the full metric tensor has to be computed and 𝑮 𝑹𝑫 has to be extracted from it.

(iii) ad hoc expressions. For instance, the diagonal contravariant component associated to a torsional degree of freedom can be defined as the inverse of a moment of inertia.

Several comments regarding the reduced dimensionality metric tensor are in order: (i) Other approximations could be considered for 𝑮 𝑹𝑫 , such as using a constant expression or a diagonal expression or both of them. [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF][START_REF] Scribano | Fast vibrational configuration interaction using generalized curvilinear coordinates and self-consistent basis[END_REF] (ii) For a molecular system, the usual KEO deformation part is obtained at J=0 (zero total angular momentum). So, formally, the overall rotation is included in the model. (iii) In most of the studies, the center-of-mass (COM) motion is not included in the model. Indeed, for isolated molecule, the wave function can be factorized rigorously into a COM contribution and a rovibronic one. However, when the molecular system is coupled to an external field or an environment (a cavity, a surface …), the COM motion can be added easily to the model. [START_REF] Bahel | Six-dimensional quantum treatment of the vibrations of diatomic adsorbates on solid surfaces: CO on Cu(100)[END_REF][START_REF] Bačić | Perspective: Accurate treatment of the quantum dynamics of light molecules inside fullerene cages: Translation-rotation states, spectroscopy, and symmetry breaking[END_REF][START_REF] Lauvergnat | H 2 , HD, and D 2 in the small cage of structure II clathrate hydrate: Vibrational frequency shifts from fully coupled quantum sixdimensional calculations of the vibration-translation-rotation eigenstates[END_REF] (iv) Compared to the constrained approach, the one using the submatrix of 𝑮 gives the smoothest 𝑮 𝑹𝑫 [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF] and the simplest analytical expression. For instance, when polyspherical coordinates are used, the metric tensor components are expressed as a sum of products of one-dimensional functions. This feature is lost with the constrained approach.

(iv) The constrained approach seems to be the most rigorous one, since it gives the correct mathematical constrained metric tensor. However, from a physical or chemical point of view, a molecular system is never strictly constrained (for instance, a distance in a molecule cannot be held fixed) and the "best" reduced dimensional model is not necessary the mathematically constrained model.

(vi) The constrained approach and the approach with a submatrix of 𝑮 are related. [START_REF] Justum | One-dimensional quantum description of the bending vibrations of HCN/CNH[END_REF][START_REF] Lauvergnat | Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates[END_REF][START_REF] Nauts | HAMILTONIANS FOR CONSTRAINED N-PARTICLE SYSTEMS[END_REF] Indeed, the constrained approach is equivalent to extracting the active covariant components, 𝒈 𝒂𝒂 , from the exact one, 𝒈 𝒆𝒙𝒂𝒄𝒕 . Then 𝑮 𝑹𝑫 (𝒊) is just the inverse of 𝒈 𝒂𝒂 (upper part of fig. 2). For the second approach, with a submatrix of 𝑮, 𝑮 𝑹𝑫 (𝒊𝒊) is extracted from 𝑮 𝒆𝒙𝒂𝒄𝒕 (lower part of fig. 2). Using block matrix inversion, one can add a correction to 𝑮 𝑹𝑫 (𝒊𝒊) to recover 𝑮 𝑹𝑫 (𝒊) . [START_REF] Brill | Photoinduced nonadiabatic dynamics of ethene: Six-dimensional wave packet propagations using two different approximations of the kinetic energy operator[END_REF] This equivalence works only with rigid constraints (i.e., without flexible constraints). However, the flexible constraints can be transformed into rigid constraints by the use of flexible coordinates. [START_REF] Gatti | Quantum-Mechanical Description of Rigidly or Adiabatically Constrained Molecular Systems[END_REF] This approach, with or without constraints, is completely general and it has been used in the past to get analytical expressions of the KEO for some small molecular systems; it has been used intensively within a numerical scheme, where the numerical values of the metric tensor (or similar expression) and the Jacobian are calculated at given geometries (i.e. on a grid).

2b) Numerical approach

The approach to numerically calculate the KEO, or the elements of the metric tensor, has been known for many years. [START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][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] Luckhaus | 6D vibrational quantum dynamics: Generalized coordinate discrete variable representation and (a)diabatic contraction[END_REF][START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Yurchenko | Theoretical ROVibrational Energies (TROVE): A robust numerical approach to the calculation of rovibrational energies for polyatomic molecules[END_REF][START_REF] Matyus | Toward black-box-type full-and reduced-dimensional variational (ro)vibrational computations[END_REF] It relies on the numerical calculations of the BF Cartesian coordinates, 𝑿 ""⃗ #$ " (𝑸), and their derivatives,

3𝑿 MM⃗ +, - (𝑸)
3E $ . Most of the times, those derivatives are obtained by means of a finite difference scheme. However, in our approach, we are using an Automatic Differentiation (AD) scheme, [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Wengert | A simple automatic derivative evaluation program[END_REF][START_REF] Bartholomew-Biggs | Automatic differentiation of algorithms[END_REF][START_REF] Yachmenev | Automatic differentiation method for numerical construction of the rotational-vibrational Hamiltonian as a power series in the curvilinear internal coordinates using the Eckart frame[END_REF] so that the derivatives, up to the third order, are calculated exactly. Our approach has been developed twenty years ago [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF] and resorts to the forward mode and to the use of an extension of the dual number approach to get high order derivatives. [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF][START_REF] Bartholomew-Biggs | Automatic differentiation of algorithms[END_REF] Note that, in Ref. [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF], even though the scheme Automatic Differentiation was employed, this exact wording was not used to identify the method. Therefore, one can use sophisticated coordinate transformations without losing accuracy. For instance, in our program TNUM, [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Lauvergnat | ElVibRot-TnumTana fortran quantum dynamics and kinetic energy operator codes[END_REF] it is particularly easy to use several coordinates transformations between the active coordinates, 𝑸 or 𝑸 𝒂𝒄𝒕 , and the BF Cartesian coordinates, 𝑿 ""⃗ #$ " (𝑸), (see the fig. 3 as a general illustration). In our implementation, for a given transformation, 𝐓𝐫 < , the coordinates 𝑸 𝒕 are transformed into the coordinates 𝑸 𝒕>𝟏 , and the required derivatives with respect to 𝑸 𝒂𝒄𝒕 are obtained through the chain rule. High order derivatives are easily obtained as well.

j 𝑄 <>( f = 𝑇𝑟 < f (𝑸 𝒕 ) 𝜕𝑄 <>( f 𝜕𝑄 gh< 8 = m 𝜕𝑄 < 9 𝜕𝑄 gh< 8 • 𝜕𝑇𝑟 < f 𝜕𝑄 < 9 : 9
Eq. 9

In Eq. 9, the derivatives $ are known from the previous transformation, 𝐓𝐫 <U( . Those expressions are completely general, but some special transformations need to be discussed in more detail:

• The active transformation: this transformation is essential when constraints (rigid or flexible) are used (𝑛 < 3𝑁 -6). In this case, one needs a relation between the 3𝑁 -6 coordinates, 𝑸 𝒅𝒚𝒏 , which describe the system, and the n-active coordinates, 𝑸 𝒂𝒄𝒕 , used in the dynamics. Three main cases have to be considered to simplify Eq. 9:

(i) When the coordinate 𝑄 FC:

f is an active coordinate and corresponds to 𝑄 gh< 9 , then the relations in Eq. 9 simplify as: Active and inactive coordinates

!" '&-# < = < >=
In terms of the AD approach, this first step is the seed of the AD procedure and, more precisely, the active coordinates are the n-independent variables among the 3N-6 ones and the remaining variables, the inactive coordinates, are the dependent variables.

• The linear transformation: this transformation is particularly simple and the coordinates 𝑸 𝒕>𝟏 are expressed as linear combinations of the coordinates 𝑸 𝒕 . Then, the relations in Eq. 9 simplify as:

⎩ ⎪ ⎨ ⎪ ⎧ 𝑄 <>( f = m 𝑀(𝑘, 𝑗) • 𝑄 < 9 9 𝜕𝑄 <>( f 𝜕𝑄 gh< 8 = m 𝑀(𝑘, 𝑗) • 𝜕𝑄 < 9 𝜕𝑄 gh< 8 9
Eq. 12 Formally, the full matrix M must be defined. However, it is possible to define a few linear combinations and, then, M is automatically computed. Moreover, it is often easier to define the inverse relation (with 𝑴 >( ), with 𝑸 𝒕 as functions of 𝑸 𝒕>𝟏 , or the symmetrized coordinates as functions of the unsymmetrized ones. With this transformation, we can easily define symmetrized coordinates, curvilinear normal modes and coordinates associated to concerted motions.

• The primitive transformation: This transformation allows one to get the Cartesian coordinates, 𝑿 ""⃗ $ " , in an intermediate frame from the so-called primitive coordinates, 𝑸 𝒑𝒓𝒊𝒎 . For instance, when the primitive coordinates are defined with the help of a z-matrix (other options are possible) and for a triatomic system (see fig, 1), 𝑿 ""⃗ $ " are defined as follows: For larger molecular system, the expressions are more involved but they can be obtained without difficulty.

⎩ ⎪ ⎨ ⎪ ⎧ 𝑿 ""⃗ $ ( =
• The Cartesian transformation: the main purpose of this transformation is to translate the previous Cartesian coordinates with respect to the center-of-mass to get the Cartesian coordinates in the BF frame. Furthermore, it enables one to take into account the rotational Eckart conditions. [START_REF] Lauvergnat | Numerical and exact kinetic energy operator using Eckart conditions with one or several reference geometries: Application to HONO[END_REF] The advantage of the approach presented above is its flexibility. Indeed, there is no built-in limitation in terms of number of atoms or coordinates, and number of coordinate transformations. [START_REF] Lauvergnat | ElVibRot-TnumTana fortran quantum dynamics and kinetic energy operator codes[END_REF] So far, only 𝑿 ""⃗ #$ " (𝑸) and its derivatives have been computed numerically and exactly for a given value of 𝑸. As for the Cartesian coordinates, the metric tensor, its derivatives, the Jacobian and the extra-potential term, they can also be computed numerically and exactly. The main part of the algorithm is presented in fig 4.

Application to reduced dimensionality photo-isomerization.

Recently, some of the authors have studied the cis-trans photoisomerization of PSB3 (see fig. 5), which is a minimal model of a retinal chromophore of Rhodopsin. [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF] Fig. 5 provides a schematic representation of PSB3, where the colours indicate the active coordinates that define our model (see description below). One of the main aspects of this study was the comparison of several quantum-classical methods with an exact quantum calculation for the cis-trans isomerization dynamics upon a photo-excitation. Therefore, it was important to use a model compatible for all approaches. The potential has been obtained by Marsili et al. [START_REF] Marsili | Two-State, Three-Mode Parametrization of the Force Field of a Retinal Chromophore Model[END_REF] and it consists of a two-states potential developed along three coordinates. Those coordinates, the torsion along the central double bond (Tors, depicted in blue in fig. 5), the bond-length-alternation stretching (BLA, represented in red in fig. 5), and the hydrogen-out-of-plane wagging of the central bond (HOOP, in green in fig. 5) capture the main features of the cis-trans photoisomerization. The quasi diabatic potential was fitted using points at XMCQDPT2//CASSCF/6-31G* ab initio level [START_REF] Marsili | Two-State, Three-Mode Parametrization of the Force Field of a Retinal Chromophore Model[END_REF] and was adjusted for Jacobian and its derivatives

S %& = ) #'( ) * # %! !" # %& $ + %! !" # %& & %S %& %& * = ) #'( ) * # % + ! !" # %& * %& $ + %! !" # %& , + %! !" # %& $ + % + ! !" # %& * %& & 8 ./ = ) #'( ) * # ! !" # + ! !" # 9 ./ -: !" #. + : !" #/ ; ., = ) #'( ) * # ! !" # × %! !" # %& $ . = ./ = * 010 + 9 ./ . = > 2×2 ? 2×4 5 0 ? 4×2 A 4×4 0 0 0 B 4×4 C 1 C 2 C 3 C 4 C 5 NH 2 H 8 H 7
the study reported in Ref. [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF]. The three coordinates were defined with the help of four transformations (see fig. 6 as an illustration): (i) A transformation of the primitive coordinates, 𝑸 𝒑𝒓𝒊𝒎 , from or to Cartesian coordinates, 𝑿 𝑭 , in the z-matrix frame.

(ii) Linear combinations of the primitive coordinates to form the 𝑸 𝒍𝒊𝒏 . This transformation helps to define the three main coordinates (BLA, Tors and HOOP).

(iii) A 1D-transformation to obtain 𝑸 𝟏𝑫 , namely a shift of the 𝑄 T8: * coordinate, so that the HOOP value is zero for the cis geometry (the reference geometry in our study). (iv) Definition of the active coordinates, 𝑸 𝒂𝒄𝒕 , from 𝑸 𝟏𝑫 or 𝑸 𝒅𝒚𝒏 . Rigid constraints are imposed at the reference cis geometry of the ground electronic state. Rigid constraints are imposed for this study, although other type of constraints could have been used. For instance, in a model with flexible constraints, the inactive coordinates are functions of the active coordinates, 𝑄 HTGB f (𝑸 𝒂𝒄𝒕 ) and they can be adapted to symmetry modifications along the active coordinates. In the majority of the studies that employed this flexible model, [START_REF] Justum | One-dimensional quantum description of the bending vibrations of HCN/CNH[END_REF][START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][START_REF] Senent | Determination of the kinetic energy parameters of non-rigid molecules[END_REF][START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF]] it has been used in combination with the adiabatic approximation, i.e., with a single electronic potential energy surface, and the 𝑄 HTGB f (𝑸 𝒂𝒄𝒕 ) functions have been obtained from a geometry optimization along 𝑸 𝒂𝒄𝒕 . The drawback is that, even when several electronic states are involved, as is the case for the PSB3 system, the geometry optimization has to be performed with one electronic state. The choice of such state is clearly arbitrary. Therefore, rigid constraints are preferred. [START_REF] Brill | Photoinduced nonadiabatic dynamics of ethene: Six-dimensional wave packet propagations using two different approximations of the kinetic energy operator[END_REF] The present implementation of the quantum-classical methods in G-CTMQC [START_REF] Agostini | G-CTMQC[END_REF] requires that the contravariant metric tensor is diagonal and constant, so that the KEO and the volume element are given by: Cartesian coordinates Primitive coordinates, (z-matrix):

C ! C " C ! " #$%& ! C ' C " " #$%& " C ! " #$%& ' C ( C ' " #$%& ( C " " #$%& ) C ! " #$%& * C ) C ( " #$%& + C ' " #$%& , C " " #$%& - # * C ) " #$%& !. C ( " #$%& !! C ' " #$%& !" H + C ( " #$%& !' C ' " #$%& !( C ) " #$%& !) H , C ' " #$%& !* C ( " #$%& !+ C " " #$%& !, ⋮ 3 active coordinates: ! ,-. , ! /012 , ! 3445 CD 6%"'#* CD "! ≡CD #$!%&' CD ( ≡CD )*%+ ! FGH ! TPK ! NOPQ
Active and inactive (rigid) coordinates defined as linear combinations of primitive ones:

" /01 = " 234 5 = " #$%& " + " #$%& 6 2 - " #$%& 5 + " #$%& 7 + & >?$@ (A 3 & B1?C = & D%E + = & >?$@ F & D%E G = -" #$%& 58 -" #$%& 59 " R
shift of ! 789 :

:

& HIIJ = & D%E G + 6V
CD (;

! MW ≡ ! IJK 𝑻 A FGH = - ℏ ! ) S5 +23,+23 3 ! 3E +23 ! U5 4567,4567 3 ! 3E 4567 ! U5 899:,899: 3 ! 3E 899: ! V 𝑑𝜏 = 𝑑𝑄 #!n • 𝑑𝑄 iojp • 𝑑𝑄 q&&r
Eq. 14

However, this constraint can be lifted for the exact quantum study, as, in fact, the metric tensor is neither constant nor diagonal. The constrained metric tensor along the Tors coordinate (see fig. 7) shows that, as function of the Tors coordinate, the variation amplitude of each component is relatively small. Furthermore, all components are symmetric with respect to the trans geometry (𝑇𝑜𝑟𝑠 = 180°), except 𝐺 #!n,iojp and 𝐺 #!n,q&&r which are anti-symmetric (see the right panel of fig. 7). In the present study, two metric tensors are compared: the constrained metric tensor (𝑮 𝑹𝑫 (𝒊) , see fig. 7) and the metric tensor obtained directly from the exact contravariant components of the metric tensor, 𝑮 𝑹𝑫 (𝒊𝒊) (see fig. 8). The comparison shows, as expected, that 𝑮 𝑹𝑫 (𝒊𝒊) presents less oscillations than 𝑮 𝑹𝑫 (𝒊) (see also reference [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF]), although the amplitude of 𝐺 q&&r,q&&r and 𝐺 iojp,iojp are larger for 𝑮 𝑹𝑫 (𝒊𝒊) than for 𝑮 𝑹𝑫 (𝒊) . From a computational point of view, the numerical calculation of the KEO (including the extra-potential term) obtained from exact metric tensor (𝑮 𝑹𝑫 (𝒊𝒊) ) is more time consuming than the constrained one, because the exact metric tensor in full dimensionality has to be computed. In fact, for the PSB3 molecule, the KEO calculation on the grid of 𝑮 𝑹𝑫 (𝒊) is about 1800 time faster than the calculation of 𝑮 𝑹𝑫 (𝒊𝒊) .

The analysis of the variations of the metric tensor components allows us to propose: [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF] (i) a constant metric tensor, 𝑮 𝑪𝒐𝒏𝒔𝒕 , which is obtained from an average of the constrained metric tensor, 𝑮 𝑹𝑫 (𝒊) , over 45 geometries, and it includes only the symmetric components, whose values are given in Table 1; (ii) a diagonal constant metric tensor, 𝑮 𝑫𝒊𝒂𝒈>𝑪𝒐𝒏𝒔𝒕 , in which 𝐺 iojp,q&&r is set to zero.

Contravariant components

Value in a.u.

𝐺 #!n,#!n 7.981 × 10 -5 𝐺 iojp,iojp
2.599 × 10 -5 𝐺 q&&r,q&&r 40.375 × 10 -5 𝐺 iojp,q&&r 4.025 × 10 -5 Table 1: Non-zero values of the averaged constant metric tensor components.

To further assess the influence of those approximations, exact quantum dynamics simulations have been performed with the three metric tensors, namely 𝑮 𝑫𝒊𝒂𝒈>𝑪𝒐𝒏𝒔𝒕 , the constrained metric tensor 𝑮 𝑹𝑫 (𝒊) , and the one extracted from the exact metric tensor, 𝑮 𝑹𝑫 (𝒊𝒊) .

The quantum dynamics calculations have been performed with the Chebyshev propagation scheme [START_REF] Tal-Ezer | An accurate and efficient scheme for propagating the time dependent Schrödinger equation[END_REF][START_REF] Leforestier | A comparison of different propagation schemes for the time dependent Schrödinger equation[END_REF] implemented in ELVIBROT fortran code [START_REF] Lauvergnat | ElVibRot-TnumTana fortran quantum dynamics and kinetic energy operator codes[END_REF] with time step of 1 fs and for 500 fs propagation duration. It is important to note, that the split operator scheme [START_REF] Feit | Solution of the Schrödinger Equation by a Spectral Method[END_REF] cannot be used for the present study. When the metric tensor is a function of all the active coordinates, one has to use a general scheme such as Chebyshev scheme, [START_REF] Tal-Ezer | An accurate and efficient scheme for propagating the time dependent Schrödinger equation[END_REF] Short Iterative Lanczos [START_REF] Leforestier | A comparison of different propagation schemes for the time dependent Schrödinger equation[END_REF] and even a Taylor expansion. [START_REF] Lauvergnat | A simple and efficient evolution operator for time-dependent Hamiltonians: the Taylor expansion[END_REF] The propagation is performed in the diabatic representation, and in each electronic state the wavepacket is expanded in a 3D-direct product grid or basis set. The primitive basis sets associated with the coordinates, 𝑄 #!n , 𝑄 iojp and 𝑄 q&&r , are, respectively, a Harmonic oscillator (HO) basis set, a Fourier series and another HO. The number of basis functions, 𝑛𝑏, and the number of grid points, 𝑛𝑞, and the HO parameters are given in Table 2 (see also Ref. [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF]). The (not-normalized) expression of the HO basis functions is 𝐻 l (𝑆𝑐(𝑄 -𝑄 k ))𝐸𝑥𝑝>-# ! (wh(E>E ; )) ! ?, where 𝐻 l are Hermite polynomials. The initial wavepacket starts in the second diabatic state at the cis configuration and it is expressed as product of three 1D-gaussians (𝑔(𝑄) = 𝐸𝑥𝑝(-(𝑄 -𝑄 h ) ) /∆𝑄 ) ) (notnormalized expression). The parameters, 𝑄 h and ∆𝑄, have been derived from the quadratic expansion of the cis ground state potential and the cis-metric tensor (see table 2 and also Ref. [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF]).

Coordinates

Primitive basis

Wavepacket Fig. 9 shows a comparison of the populations (full and restricted to the trans domain, i.e., when |𝑇𝑜𝑟𝑠| ∈ [𝜋/2, 𝜋]) of the 𝑆 k adiabatic state for different metric tensor models. As already shown in reference [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF] the populations with the 𝑮 𝑫𝒊𝒂𝒈>𝑪𝒐𝒏𝒔𝒕 and 𝑮 𝑹𝑫 (𝒊) give similar results. This is an expected result, since the diagonal constant metric tensor is obtained from an average operation of the constrained metric tensor. All models present the same trend for the 𝑆 ( -to-𝑆 k population transfer. After about a 30 fs delay from the initial photoexcitation, a fast transfer occurs, and then the 𝑆 k population increases slowly. However, the population transfer to the 𝑆 k trans isomer is two times less efficient with the model 𝑮 𝑹𝑫 (𝒊𝒊) than with the others. This feature is due to the trapping of the wavepacket around the cis configuration on the 𝑆 ( state which is stronger for the 𝑮 𝑹𝑫 (𝒊𝒊) model than for the 𝑮 𝑹𝑫 (𝒊) one. Indeed, at the end of the dynamics (500 fs), the reduced density in 𝑆 ( along the torsional angle presents a peak at around 0° (cis configuration) which is two times higher for the 𝑮 𝑹𝑫 (𝒊𝒊) model than for the 𝑮 𝑹𝑫 (𝒊) model (see fig. 10). Furthermore, the reduced density in 𝑆 k is almost completely delocalized (see fig. 10). Therefore, the population transfer from 𝑆 ( -to-𝑆 k around the trans configuration (when |𝑇𝑜𝑟𝑠| ∈ [𝜋/2, 𝜋]) is about half of the total population transfer (see fig. 9). Although the trend of the population transfer after photoexcitation is similar for the 𝑮 𝑹𝑫 (𝒊) and 𝑮 𝑹𝑫 (𝒊𝒊) models, the quantitative details are different. Probably, this effect is consequence of the strong kinetic couplings between the active and inactive coordinates, since the correction from the 𝑮 𝑹𝑫 (𝒊𝒊) model to 𝑮 𝑹𝑫 (𝒊) involves the metric tensor components between active and inactive coordinates, 𝒈 𝒂𝒊 , (see fig 2). [START_REF] Brill | Photoinduced nonadiabatic dynamics of ethene: Six-dimensional wave packet propagations using two different approximations of the kinetic energy operator[END_REF] From the photoisomerization point of view, it means that some coordinates should be explicitly added in the dynamics to reduce this coupling.

type 𝑛𝑏 𝑛𝑞 𝑄 k Sc 𝑄 h ∆𝑄 𝑄 #!n HO 20

Discussion

As mentioned earlier, the numerical implementation of the KEO calculation with generalized coordinates is not new. [START_REF] Meyer | Flexible Models for Intramolecular Motion, a Versatile Treatment and Its Application to Glyoxal[END_REF][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] Luckhaus | 6D vibrational quantum dynamics: Generalized coordinate discrete variable representation and (a)diabatic contraction[END_REF] However, our implementation in TNUM [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF][START_REF] Lauvergnat | ElVibRot-TnumTana fortran quantum dynamics and kinetic energy operator codes[END_REF] has several advantages: (i) It is exact. The derivatives are not calculated with a finite difference scheme, as is often the case in other codes, instead they are computed with an Automatic Differentiation scheme. [START_REF] Lauvergnat | Exact numerical computation of a kinetic energy operator in curvilinear coordinates[END_REF] (ii) No limitations. Our implementation has no built-in limitations in terms of number of atoms, degrees of freedom or coordinate transformations.

(iii) Its flexibility. Along with the usual z-matrix or polyspherical coordinates, several nested coordinates transformations can be added, and it becomes very easy to test new coordinates. Furthermore, if a coordinates transformation is not available, its implementation in TNUM is straightforward. (iv) Easy to interface. As a library, TNUM can be interfaced easily with other quantum dynamics codes for different kinds of applications, such as collisions, [START_REF] Baloïtcha | Cumulative reaction probability by constrained dynamics: H transfer in HCN, H2CO, and H3CO[END_REF][START_REF] Von Horsten | An efficient route to thermal rate constants in reduced dimensional quantum scattering simulations: Applications to the abstraction of hydrogen from alkanes[END_REF] spectroscopy (ro-vibrational, [START_REF] Lauvergnat | Numerical and exact kinetic energy operator using Eckart conditions with one or several reference geometries: Application to HONO[END_REF][START_REF] Sarka | Rovibrational energy levels of the F -(H 2 O) and F -(D 2 O) complexes[END_REF] non-linear optics [START_REF] Luis | Treatment of nonlinear optical properties due to large amplitude anharmonic vibrational motions: Umbrella motion in NH3[END_REF][START_REF] Garcia-Borràs | A Full Dimensionality Approach to Evaluate the Nonlinear Optical Properties of Molecules with Large Amplitude Anharmonic Tunneling Motions[END_REF]), quantum dynamics of photo-induced processes [START_REF] Marsili | Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model[END_REF][START_REF] Piechowska-Strumik | Quantum dynamics around a non planar conical intersection in vinoxy radical relaxation[END_REF] or quantum dynamics with laser pulses and optimal control. [START_REF] Lauvergnat | A simple and efficient evolution operator for time-dependent Hamiltonians: the Taylor expansion[END_REF][START_REF] Bomble | Controlled full adder or subtractor by vibrational quantum computing[END_REF][START_REF] Bomble | Vibrational computing: simulation of a full adder by optimal control[END_REF][START_REF] Bomble | Controlled full adder-subtractor by vibrational computing[END_REF][START_REF] Ndong | Optimal control simulation of the Deutsch-Jozsa algorithm in a two-dimensional double well coupled to an environment[END_REF] In the past, it has been interfaced with some "in-house" codes [START_REF] Lauvergnat | Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates[END_REF][START_REF] Scribano | Fast vibrational configuration interaction using generalized curvilinear coordinates and self-consistent basis[END_REF][START_REF] Baloïtcha | Cumulative reaction probability by constrained dynamics: H transfer in HCN, H2CO, and H3CO[END_REF][START_REF] Piechowska-Strumik | Quantum dynamics around a non planar conical intersection in vinoxy radical relaxation[END_REF][START_REF] Lauvergnat | Wave packet propagation for constrained molecular systems: spectroscopic applications to triatomic molecules[END_REF] and with codes developed by various groups. [START_REF] Von Horsten | An efficient route to thermal rate constants in reduced dimensional quantum scattering simulations: Applications to the abstraction of hydrogen from alkanes[END_REF][START_REF] Luis | Treatment of nonlinear optical properties due to large amplitude anharmonic vibrational motions: Umbrella motion in NH3[END_REF][START_REF] Von Horsten | Fingerprints of delocalized transition states in quantum dynamics[END_REF][START_REF] Brogaard | Wave Packet Simulation of Nonadiabatic Dynamics in Highly Excited 1, 3-Dibromopropane[END_REF] Curvilinear coordinates along with the corresponding kinetic contributions have been used with semiclassical dynamics (a single thawed gaussian propagation [START_REF] Lauvergnat | Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates[END_REF]) or with constrained classical dynamics. [START_REF] Lauvergnat | Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates[END_REF] For the latter, the SHAKE algorithm [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF] has not been used, and instead Hamilton's equations have been solved directly in curvilinear coordinates. This possibility to perform classical or quantum-classical dynamics with generalized coordinates with TNUM has been seldom used. However, interface with quantum-classical dynamics codes using Ehrenfest dynamics [START_REF] Tully | Mixed quantum-classical dynamics[END_REF] or Surface Hopping, [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] and TNUM should be straightforward even when constraints are required. An interesting development can be envisaged aiming at interfacing TNUM with the Coupled-Trajectory Mixed Quantum-Classical [START_REF] Agostini | Quantum-Classical Nonadiabatic Dynamics: Coupled-vs Independent-Trajectory Methods[END_REF] (CT-MQC) algorithm and G-CTMQC code [START_REF] Agostini | G-CTMQC[END_REF], even though it is likely that CT-MQC equations need to be (re)derived from the exact-factorisation equations by introducing generalised coordinates in the original formulation. This feature could be tested easily on the cis-trans photo-isomerization of PSB3, so that the exact reduced dimensionality metric tensor could be used without the diagonal and constant approximation.

Tables

Contravariant components

Value in a.u. * , so that the HOOP value is zero at the cis-configuration. 𝑸 𝟏𝑫 are also the dynamical coordinates 𝑸 𝒅𝒚𝒏 . (v) 𝑸 𝒂𝒄𝒕 : the three active coordinates are selected from the 3N-6 𝑸 𝒅𝒚𝒏 coordinates. The other coordinates, the inactive ones, are constrained to their values at the cis-configuration.

𝐺 #!n,#

Figure and table captions

Fig 2 :

 2 Fig 2: Scheme to obtain the reduced dimensionality metric tensors.

  Matrix inversion correction . 67895 = . 88 . 8% . %8 . %% -67895 = -88 -8% -%8 -%% . :; = . 88 -:; (%) = . 88 -( -:; (%%) = -88 --8% + -%% -( + -%8

Fig 3 :

 3 Fig 3: Illustration of the coordinate transformations in TNUM, from the active coordinates, 𝑸 ≡ 𝑸 𝒂𝒄𝒕 , to the body fixed Cartesian coordinates. Each transformation, 𝑻𝒓 < , is represented, by a pink arrow and a green label. The different coordinates, 𝑸 𝒂𝒄𝒕 , 𝑸 𝒅𝒚𝒏 , 𝑸 𝒕 , 𝑸 𝒑𝒓𝒊𝒎 , 𝑿 𝑭 , 𝑿 𝑩𝑭 are represented by blue boxes.

%

  of the transformation, 𝐓𝐫 < , are known from its definition and the derivatives 3E . % 3E 01.

  When the coordinate 𝑄 FC: f is inactive and equal to a constant value, 𝑄 FC:k f , all derivatives are zero. (iii) When 𝑄 FC: f is inactive and equal to a function of the active coordinates (flexible constraints), then the relations in Eq. 9 simplify as: j 𝑄 FC: f = 𝑄 HTGB f (𝑸 𝒂𝒄𝒕 )

Fig 4 :

 4 Fig 4: Illustration of TNUM implementation.

Fig 5 :

 5 Fig 5: PSB3 scheme and the corresponding z-matrix. The primitive coordinates highlighted in color are the three coordinates used in the model.

Fig 6 :

 6 Fig 6: PSB3 TNUM coordinates transformations (pink arrows and green boxes) with the five sets of coordinates (blue and yellow boxes): (i) 𝑿 𝑭 : The Cartesian coordinates in the z-matrix frame where the first three atoms (C1, C2 and C3) are in the xz plan and C1-C2 are along the z-axis. (ii) 𝑸 𝒑𝒓𝒊𝒎 : the primitive or z-matrix coordinates. The ones involve in the definitions of active coordinates are highlighted in color. (iii) 𝑸 𝒍𝒊𝒏 : the three linear combinations of 𝑸 𝒑𝒓𝒊𝒎 helping to define the three main coordinates (BLA, Tors and HOOP). (iv) 𝑸 𝟏𝑫 : The coordinates obtained after a 1D-transformation from 𝑸 𝒍𝒊𝒏 . Here, it is just a shift of the 𝑄 T8: * , so that the HOOP value is zero at the cisconfiguration. 𝑸 𝟏𝑫 are also the dynamical coordinates 𝑸 𝒅𝒚𝒏 . (v) 𝑸 𝒂𝒄𝒕 : the three active coordinates are selected from the 3N-6 𝑸 𝒅𝒚𝒏 coordinates. The other coordinates, the inactive ones, are constrained to their values at the cisconfiguration.

Fig 7 :

 7 Fig 7: Metric tensor components along the Tors for the constrained metric tensor, 𝑮 𝑹𝑫 (𝒊) .

Fig 8 :

 8 Fig 8: Metric tensor components, 𝑮 𝑹𝑫 (𝒊𝒊) , along the Tors from the contravariant components of the exact metric tensor.

Fig 9 :

 9 Fig 9: Population in the 𝑆 0 adiabatic state for several metric tensors. Left pane full population, right panel population around the trans isomer.

Fig 10 :

 10 Fig 10: Reduced density along the torsional angle on the 𝑆 k and the 𝑆 ( states at the end of the dynamics (500 fs). Left panel 𝑮 𝑹𝑫 (𝒊) model. Right panel 𝑮 𝑹𝑫 (𝒊𝒊) model.

Fig 1 :

 1 Fig 1: Illustration of the body-fixed frame orientation with respect to the laboratory frame for a triatomic molecule (A1-A2-A3) defined with the help of a z-matrix. The space-fixed (SF) frame is an intermediate frame parallel to the LF frame but centered on the center-of-mass.

Fig 2 :

 2 Fig 2: Scheme to obtain the reduced dimensionality metric tensors.

Fig 3 :

 3 Fig 3: Illustration of the coordinate transformations in TNUM, from the active coordinates, 𝑸 ≡ 𝑸 𝒂𝒄𝒕 , to the body fixed Cartesian coordinates. Each transformation, 𝑻𝒓 < , is represented, by a pink arrow and a green label. The different coordinates, 𝑸 𝒂𝒄𝒕 , 𝑸 𝒅𝒚𝒏 , 𝑸 𝒕 , 𝑸 𝒑𝒓𝒊𝒎 , 𝑿 𝑭 , 𝑿 𝑩𝑭 are represented by blue boxes.

Fig 4 :Fig 5 :

 45 Fig 4: Illustration of TNUM implementation.

Fig 6 :

 6 Fig 6: PSB3 TNUM coordinates transformations (pink arrows and green boxes) with the five sets of coordinates (blue and yellow boxes): (i) 𝑿 𝑭 : The Cartesian coordinates in the z-matrix frame where the first three atoms (C1, C2 and C3) are in the xz plan and C1-C2 are along the z-axis. (ii) 𝑸 𝒑𝒓𝒊𝒎 : the primitive or z-matrix coordinates. The ones involve in the definitions of active coordinates are highlighted in color. (iii) 𝑸 𝒍𝒊𝒏 : the three linear combinations of 𝑸 𝒑𝒓𝒊𝒎 helping to define the three main coordinates (BLA, Tors and HOOP). (iv) 𝑸 𝟏𝑫 : The coordinates obtained after a 1Dtransformation from 𝑸 𝒍𝒊𝒏 . Here, it is just a shift of the 𝑄 T8:* , so that the HOOP value is zero at the cis-configuration. 𝑸 𝟏𝑫 are also the dynamical coordinates 𝑸 𝒅𝒚𝒏 . (v) 𝑸 𝒂𝒄𝒕 : the three active coordinates are selected from the 3N-6 𝑸 𝒅𝒚𝒏 coordinates. The other coordinates, the inactive ones, are constrained to their values at the cis-configuration.

Fig 7 :Fig 8 :

 78 Fig 7: Metric tensor components along the Tors for the constrained metric tensor, 𝑮 𝑹𝑫 (𝒊) .

Fig 9 :

 9 Fig 9: Population in the 𝑆 0 adiabatic state for several metric tensors. Left pane full population, right panel population around the trans isomer.

Fig 10 :FiguresFig 1 :

 101 Fig 10: Reduced density along the torsional angle on the 𝑆 k and the 𝑆 ( states at the end of the dynamics (500 fs). Left panel 𝑮 𝑹𝑫 (𝒊) model. Right panel 𝑮 𝑹𝑫 (𝒊𝒊) model.

Fig 2 :Fig 3 :

 23 Fig 2: Scheme to obtain the reduced dimensionality metric tensors.

Fig 4 :Fig 5 :

 45 Fig 4: Illustration of TNUM implementation.

Fig 6 :

 6 Fig 6: PSB3 TNUM coordinates transformations (pink arrows and green boxes) with the five sets of coordinates (blue and yellow boxes): (i) 𝑿 𝑭 : The Cartesian coordinates in the z-matrix frame where the first three atoms (C1, C2 and C3) are in the xz plan and C1-C2 are along the z-axis. (ii) 𝑸 𝒑𝒓𝒊𝒎 : the primitive or z-matrix coordinates. The ones involve in the definitions of active coordinates are highlighted in color. (iii) 𝑸 𝒍𝒊𝒏 : the three linear combinations of 𝑸 𝒑𝒓𝒊𝒎 helping to define the three main coordinates (BLA, Tors and HOOP). (iv) 𝑸 𝟏𝑫 : The coordinates obtained after a 1D-transformation from 𝑸 𝒍𝒊𝒏 . Here, it is just a shift of the 𝑄 T8:* , so that the HOOP value is zero at the cisconfiguration. 𝑸 𝟏𝑫 are also the dynamical coordinates 𝑸 𝒅𝒚𝒏 . (v) 𝑸 𝒂𝒄𝒕 : the three active coordinates are selected from the 3N-6 𝑸 𝒅𝒚𝒏 coordinates. The other coordinates, the inactive ones, are constrained to their values at the cisconfiguration.

Fig 7 :

 7 Fig 7: Metric tensor components along the Tors for the constrained metric tensor, 𝑮 𝑹𝑫 (𝒊) .

Fig 8 :Fig 9 :

 89 Fig 8: Metric tensor components, 𝑮 𝑹𝑫 (𝒊𝒊) , along the Tors from the contravariant components of the exact metric tensor.

Fig 10 :

 10 Fig 10: Reduced density along the torsional angle on the 𝑆 k and the 𝑆 ( states at the end of the dynamics (500 fs). Left panel 𝑮 𝑹𝑫 (𝒊) model. Right panel 𝑮 𝑹𝑫 (𝒊𝒊) model.

Table 2 :

 2 Primitive basis set definitions and 1D-gaussian parameters of the initial wavepacket.

				32 0.1725 9.157 0.1725 0.1544
	𝑄 iojp	Fourier 128 136			0.	0.1833
	𝑄 q&&r	HO	24	60	0.	3.482	0.	0.4061

Table 1 :

 1 Non-zero values of the averaged constant metric tensor components.

			!n		7.981 × 10 -5
			𝐺 iojp,iojp		2.599 × 10 -5
			𝐺 q&&r,q&&r	40.375 × 10 -5
			𝐺 iojp,q&&r	4.025 × 10 -5
			Primitive basis		Wavepacket
	Coordinates	type	𝑛𝑏 𝑛𝑞	𝑄 k	Sc	𝑄 h	∆𝑄
	𝑄 #!n	HO	20	32 0.1725 9.157 0.1725 0.1544
	𝑄 iojp	Fourier 128 136			0.	0.1833
	𝑄 q&&r	HO	24	60	0.	3.482	0.	0.4061

Table 2 :

 2 Primitive basis set definitions and 1D-gaussian parameters of the initial wavepacket.

Table 1 :

 1 Non-zero values of the averaged constant metric tensor components.

Table 2 :

 2 Primitive basis set definitions and 1D-gaussian parameters of the initial wavepacket.

Acknowledgments

Please acknowledge anyone who contributed to the study but did not meet the authorship criteria.

Funding Statement E.M. and F.A. thank CNRS for financial support via the International Emerging Action.

Data Accessibility

The codes used for the present study are freely available from github repository.

-1-The quantum dynamics code, ELVIBROT, and the code for the KEO, TNUM: https://github.com/lauvergn/ElVibRot-TnumTana -2-Quantum Model Library (QML) including the PSB3 potential: https://github.com/lauvergn/QuantumModelLib

Additional Information

Information on the following should be included wherever relevant.

Competing Interests

We have no competing interests.

Authors' Contributions

EM carried out the calculations and implemented the PSB3 potential in the Quantum Model Library (QML). FA and DL supervised the work on PSB3. DL is the main developer of the QML and ELVIBROT-TNUM-TANA codes. AN and DL developed the TNUM algorithm. All authors contributed to writing the manuscript.