Nataša Nadoveza

Ramón L Panadés-Barrueta

Lei Shi

Fabien Gatti

Daniel Peláez
email: daniel.pelaez-ruiz@universite-paris-saclay.fr

Analytical high-dimensional operators in Canonical Polyadic Finite Basis Representation (CP-FBR)

In the present work, we introduce a simple means of obtaining an analytical (i.e. grid-free) canonical polyadic (CP) representation of a multidimensional function which is expressed in terms of a set of discrete data. For this, we make use of an initial CP guess, even not fully converged, and a set of auxiliary basis functions (finite basis representation, FBR). The resulting CP-FBR expression constitutes the CP counterpart of our previous Tucker SOP-FBR approach. However, as is well-known, CP expressions are much more compact. This has obvious advantages in high-dimensional quantum dynamics. The power of CP-FBR lies in the fact that it requires a grid much coarser than the one needed for the dynamics. In a subsequent step, the basis functions can be interpolated to any desired density of grid points. This is useful, for instance, when different initial conditions (e.g. energy content) of a system are to be considered. We show the application of the method to bound systems of increased dimensionality: H 2 (3D), HONO (6D) and CH 4 (9D).

Introduction

Tensors are ordered collections of data which naturally arise in scientific and engineering problems when grid representations are employed. In computational science tensors are commonly referred to as multidimensional arrays. Tools to deal with these ubiquitous data structures are of high relevance in any real application in all branches of science and engineering [START_REF] Harris | Array programming with Numpy[END_REF]. Hereafter, we will focus on the grid-based representation of quantum operators for molecular quantum dynamics. However, the ideas we will put across are straightforwardly applicable to any grid-based problem.

Let us consider an f -dimensional (f D) discrete (grid, I) representation of the configura-tion space of a molecule (Q). One can represent this as the mapping: Q ≡ (q 1 , . . . , q f ) → (i 1 , . . . , i f ) ≡ I

(1)

where we have defined the composite indices (I and Q) as a tuple of f indices (i 1 , . . . , i f ) associated to f physical coordinates, (q 1 , . . . , q f ), e.g. distances, valence or dihedral angles.

In the case of a non-linear isolated molecular system f =3N-6 with N being the number of atoms. Hereafter, the number f will also be referred to as the number of degrees of freedom (DOF). We shall represent the dimensionality of the problem with the letter f for the sake of consistency with the MCTDH literature (f from German Freiheitsgrade, literally DOF).

Quantities expressed in configuration space read:

V (q 1 , . . . , q f ) ≡ V (Q) (2)
and when mapped onto a grid

V (Q) → V (I) (3)
In other words, they are represented by tensors:

V (I) ∈ R N 1 ×N 2 ...×N f (4)
here N κ represents the number of grid points employed in the discretisation of the κ-th DOF (κ = 1, 2, . . . f ). Furthermore, in our specific application, we shall concentrate in problems for which the function underlying the tensor elements is smooth and well-behaved. This is the most common situation when dealing with physical properties. However, it should be also noticed that this behaviour will strongly depend on the choice of the underlying set of coordinates. At this point, it is convenient to highlight that the minimum grid size or, alternatively, the minimum density of the grid is given a priori by the specific needs of the problem under consideration, e.g. the total energy content of the system as well as the availability of computer resources. Obviously, the size of the grid will imply an exponential growth of data-points and the concomitant number of operations to perform with them, the so-called curse of dimensionality. [START_REF] Bellman | Adaptive Control Processes: A Guided Tour[END_REF] To illustrate this, consider, for instance, the direct mapping of the Coulomb potential on a Cartesian grid for the explicit quantum dynamics of two electrons coupled to a continuum. This system would require of the order of ∼ 10 12

grid points [START_REF] Haller | Inter-Coulombic Decay in Laterally-Arranged Quantum Dots Controlled by Polarized Lasers[END_REF], or equivalently ∼ 10 3 GB of memory, which is beyond the capabilities of small workstations and most common-sized cluster nodes. This problem illustrates a major computational bottleneck in the field of grid-based (molecular) quantum dynamics.

In this context, the Multiconfiguration Time-Dependent Hartree (MCTDH) method [START_REF] Manthe | Wave-Packet Dynamics within the Multiconfiguration Hartree Framework: General Aspects and application to NOCl[END_REF][START_REF] Beck | An efficient and robust integration scheme for the equations of motion of the multiconfiguration time-dependent Hartree (MCTDH) method[END_REF][START_REF] Beck | The multi-configuration timedependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wave packets[END_REF][START_REF] Meyer | Quantum molecular dynamics: Propagating wavepackets and density operators using the multiconfiguration time-dependent Hartree (MCTDH) method[END_REF] for the resolution of the Time-Dependent Schrödinger Equation (TDSE) (in any of its associated software implementations, for instance [START_REF] Worth | Version 8.4.12[END_REF] or [START_REF] Worth | The QUANTICS Package, Version 1.1[END_REF]) has played a major role in this field. MCTDH is a variatonal method and thus tends to the exact result when the size of the MCTDH basis set tends to infinity. Hence, for decades now, it has served both as tool of choice and as method of reference for the development of other methods for multidimensional wavepacket propagations. In this way, all quantities involved in MCTDH calculations have to be formally represented on a grid, that is, every quantity has to have a value associated to each one of the system configurations (grid points). These raw grid representations suffer from aforementioned exponential scaling and consequently tensor decomposition methods are mandatory whenever high-dimensionality studies (effectively for f > 3) are aimed for.

Arguably the main bottleneck in this field is the computation and suitable representation of a potential energy surface (PES) or, more formally, the PES operator ( V ) . [START_REF] Jäckle | Product representation of potential energy surfaces II[END_REF][START_REF] Jäckle | Reaction cross sections for the H+D 2 (ν = 0, 1) system for collision up to 2.5 eV: A multiconfiguration time-dependent Hartree wavepacket propagation study[END_REF][START_REF] Gatti | Intramolecular Vibrational Energy Redistribution in Toluene: A nine dimensional Quantum mechanical study using the MCTDH algorithm[END_REF][START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF][START_REF] Otto | Multi-Layer Potfit: An accurate potential representation for efficient highdimensional quantum dynamics[END_REF][START_REF] Panadés-Barrueta | Specific Reaction Parameter Multigrid POTFIT (SRP-MGPF): Automatic Generation of Sum-of-Products Form Potential Energy Surfaces for Quantum Dynamical Calculations[END_REF] The latter, contrary to the case of elementary interactions, is specific for a molecular system (and, generally, this also holds for its isotopologues).

Traditionally, in the context of MCTDH, an already existing PES needs to be pre-processed (refitted ) in order to achieve the so-called sum-of-products (SOP) form:

V = r c r κ v (κ) r (5)
where the v (κ) are mono-(or low-dimensional) basis functions (aka factors) expressed on a grid. In this work, we will solely consider the Tucker (Eq.6) and the canonic polyadic (CP) forms (Eq.8) which are commonly used in MCTDH. The CP decomposition (CPD) is also known as CANDECOMP or PARAFAC in the literature. In both ansätze, the PES value for a configuration V i 1 ,...,i f is approximated by a sum of products of so-called factor matrices (v (κ) ) or single-particle functions (SPP) auxiliary in the MCTDH literature. The set of weights associated to each of the configurations (Hartree product) is known as core tensor (C j 1 ...j f ) in Tucker form and simply weights (c r ) in CP form. For comparable accuracy, CP requires a much smaller number of terms than Tucker and this renders CP as the form of choice to represent a high-dimensional operator on a grid. However, as it will be argued later, the optimisation of the CP form might lead to well-known (at least in the mathematical community) bottlenecks owing to the non-orthogonal nature of its bases (vide infra). On the other hand, the corresponding algebraic process in Tucker is smooth and painless.

In any case, as a general statement, it should be evident that any reduction on the grid necessary to compute a given decomposition will be highly beneficial. In the case of grid-based methods, this was the underlying idea behind all different POTFIT-related flavours. [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF][START_REF] Otto | Multi-Layer Potfit: An accurate potential representation for efficient highdimensional quantum dynamics[END_REF][START_REF] Schröder | Transforming high-dimensional potential energy surfaces into sum-of-products form using Monte Carlo methods[END_REF] Moreover, during the last decade, several groups have also proposed powerful fitting methods in this very same spirit, such as neural-network based ones, [START_REF] Manzhos | Neural Network Potential Energy Surfaces for Small Molecules and Reactions[END_REF].

The interested reader is referred to our previous publication in which an extensive discussion on this topic and the extended list of references is provided. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] At this point, it is necessary to mention that it has been recently shown the possibility of obtaining an analytical CP expansion through a Gaussian Regression process. [START_REF] Song | Direct Canonical-Polyadic-Decomposition of the Potential Energy Surface from Discrete Data by Decoupled Gaussian Process Regression[END_REF] In the context of MCTDH, a major breakthrough has been the introduction of a Monte Carlo based Canonical Polyadic Decom-position (MCCPD) of the potential by Schröder. [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into a canonical polyadic decomposition using Monte Carlo methods[END_REF] MCCPD has been shown to be able to refit potential energy surfaces with a surprisingly good accuracy with a much larger number of degrees of freedom than with the different variants of POTFIT.

In this same line of thought, some of us have recently shown that it is also possible to achieve this by turning these tensor-decomposed expressions into an analytical form through the use of a set of basis functions, hence the name of sum-of-products finite-basisrepresentation (SOP-FBR). [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] In this respect, we have only recently find out that a related idea was suggested by Long and Long [START_REF] Long | Surface Approximation and Interpolation via Matrix SVD[END_REF] for the specific case of 3D SVD problems. Unfortunately no follow-up was proposed and, hence, our contribution is independent and general.

Back to SOP-FBR, in order to achieve this form one can use a set of ab initio points (geometry and associated energy) directly. SOP-FBR yields an analytical Tucker-form expression for multidimensional PES with an accuracy that outperforms a direct (algebraic) POTFIT algorithm. Furthermore, in SOP-FBR the use of grids much coarser than the primitive ones, required by POTFIT, leads to a much lower computational cost, both in terms of memory requirements and CPU time. In view of these advantageous features, in this work we propose the extension of the SOP-FBR concept to the CP form, hence the name of CP-FBR ansatz.

The basic idea is presented in Figure 1. We start from an initial CPD guess which in the case of large systems will be obtained from MCCPD on a coarse grid (which is not dense enough for quantum dynamics). The CPD guess will be turned into an analytical form by the use of a set of basis sets (adapted to the problem), as we did in SOP-FBR. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] Having an analytical expression for the PES has the great advantage that there is no need to refit (e.g. PF or MCCPD) many times the potential whenever the density of the grid has to be changed (vide infra). Consequently, CP-FBR further extends the power of MCCPD (and hence MCTDH) to deal with high-dimensional systems by: (i) providing its analytical CP expression; and thus (ii) by reducing the computational cost of achieving a high-dimensional CP representation of the, say, PES. Note in passing that our approach can be used with any potential-like operator, not only the PES. This includes numerically computed KEO and extrapotentials. [START_REF] Peláez | Full-dimensional MCTDH/MGPF study of the ground and lowest lying vibrational states of the bihydroxide H 3 O -2 complex[END_REF] This article is structured as follows. First, we introduce the notation together with the basic theory around the aforementioned grid-based tensor decomposition methods (Tucker and CP). Then we discuss the computation of the Finite Basis Representation (FBR) in the context of the the Canonical Polyadic (CP) form. Finally, we shall illustrate its power and discuss its behaviour through its application to fully general 3D, 6D and 9D problems in spectroscopy. This article concludes with a summary of our results and future perspectives.

r = 1 r = 1 r = 1 r = R r = R r = R λ 1 λ R λ 1 λ 1 λ R λ R . .

Theory

In the following subsections, we will present and succinctly discuss the most common tensor decomposition algorithms in the context of the MCTDH method for quantum dynamics and the necessary tensor notation used in this paper.

Grid-based Tucker and CP representations

As discussed in the Introduction section, MCTDH grid-based quantities, particularly the PES, are expressed in SOP form. This form has two main advantages, it enables the expression of (ubiquitous) multidimensional integrals into sums of product of low-dimensional integrals and it reduces the exponential scaling. [START_REF] Beck | The multi-configuration timedependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wave packets[END_REF] We will not discuss here the powerful multilayer potfit [START_REF] Otto | Multi-Layer Potfit: An accurate potential representation for efficient highdimensional quantum dynamics[END_REF] or the tree-tensor network states approaches [START_REF] Larsson | Computing vibrational eigenstates with tree tensor network states (TTNS)[END_REF] which are out of the scope of the present work. The most common of the non-recursive approaches make use of the so-called Tucker form which reads:

V i 1 ,...,i f ≈ V Tucker i 1 ,...,i f = V Tucker I = m 1 j 1 =1 • • • m f j f =1 C j 1 ...j f f κ=1 v (κ) iκ,jκ = J C J Ω IJ (6)
In the expression above, we have utilised the composite indices I ≡ (i 1 , . . . , i f ) and J ≡ (j 1 , . . . , j f ) and have defined the configurations Ω IJ (Hartree products). Then the expression Ω (ν) I ν J ν refers to a configuration in which the factor corresponding to the ν-th DOF has been integrated out. In that case, one can write the equivalent, but computationally more advantageous, [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF][START_REF] Jäckle | Time-dependent calculation of reactive flux employing complex absorbing potentials: General aspects and application within MCTDH[END_REF] expression:

V Tucker I = J ν D (ν) J ν Ω (ν) I ν J ν (7)
From all of the above, it should be clear that it is always possible to re-express any PES in Tucker form on a grid. Let us discuss some of its features. As observed, the original tensor (V) is re-expressed into another tensor (core tensor, C) of the same dimensionality which contains the weights of the configurations or Hartree products (products of factor matrices, SPPs, v (κ) ). Several algorithms have been proposed for this POTFIT (PF) being the pioneer. [START_REF] Jäckle | Time-dependent calculation of reactive flux employing complex absorbing potentials: General aspects and application within MCTDH[END_REF] The interested reader is deferred to the abundant literature, in particular, to the variational PF variants which enable the use of larger grids than PF. [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF][START_REF] Schröder | Transforming high-dimensional potential energy surfaces into sum-of-products form using Monte Carlo methods[END_REF] The Tucker factor matrices are orthonormal (both on grid and basis spaces). [START_REF] Jäckle | Time-dependent calculation of reactive flux employing complex absorbing potentials: General aspects and application within MCTDH[END_REF] Provided that a sufficiently dense grid is considered their shape is smooth, at least for the set of j κ lowest functions (this value depends on the case). This fact was observed by the original authors [START_REF] Jäckle | Time-dependent calculation of reactive flux employing complex absorbing potentials: General aspects and application within MCTDH[END_REF] as well as others. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF][START_REF] Manzhos | A random-sampling high dimensional model representation neural network for building potential energy surfaces[END_REF] With respect to the core tensor elements, their value typically spans a very range of values making if feasible its pruning. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] Despite these nice properties, the core tensor suffers dramatically from the curse of dimensionality. As mentioned above, it is convenient to recall that several POTFIT flavors exist that partially circumvent this limitation. [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF][START_REF] Schröder | Transforming high-dimensional potential energy surfaces into sum-of-products form using Monte Carlo methods[END_REF]. This is also true for some types of neural networks as shown by Manzhos

and Carrington [START_REF] Manzhos | A random-sampling high dimensional model representation neural network for building potential energy surfaces[END_REF], exploited by Brown and coworkers [START_REF] Pradhan | A ground state potential energy surface for HONO based on a neural network with exponential fitting functions[END_REF] as well as Koch and Zhang, the latter in the context of Gaussian expansion of wave packets. [START_REF] Koch | Communication: Separable potential energy surfaces from multiplicative artificial neural networks[END_REF] To alleviate the dimensionality issues associated to the exponential growth of the core, Schröder [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into a canonical polyadic decomposition using Monte Carlo methods[END_REF] has recently introduced a (Monte Carlo based) CP decomposition scheme (MC-CPD) which is also of SOP form. In its normalised version a CPD reads:

V CP i 1 ,...,i f = R r=1 λ r f κ=1 v (κ) iκ,r (8)
where {λ} R r=1 is a vector of weights and v (κ) are the factor matrices. As it can be readily observed, CP is much more computationally convenient than Tucker in terms of memory.

Indeed, λ is a vector and C is a tensor and typically R ≪ f κ=1 m κ . Furthermore, CP expansions are surprisingly accurate even for low ranks (R). This is due to the non-orthogonal character of their factor matrices. This very last feature is actually the source of some difficulties (multilinearity, swamping) when considering the optimisation of the CP ansatz. [START_REF] Battaglino | A Practical Randomized CP Tensor Decomposition[END_REF] Sum-of-products finite basis set representation (X-FBR, X=Tucker/SOP,

CP)

To finalize this section, we will briefly discuss our recent grid-free approach leading to analytical SOP PES. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] A SOP-FBR expression reads:

V (q 1 , . . . , q f ) = m 1 j 1 • • • m f j f C j 1 ,...,j f f κ ( tκ µ c (κ) µ,jκ T µ (q κ )) (9)
Several observations common to all Tucker-related methods have led to some of us to suggest this ansatz. First, the PES SOP representation is typically of a much lower rank than the actual rank of the full primitive tensor, that is f κ m κ ≪ f κ N κ (see Eqs.4 and 6).

In other words, one needs much less factors (SPPs) than grid points to expand the potential.

Second, the (discrete) values of the lowest-j factor matrices (SPPs) can be easily fitted to a suitable family of, say, orthogonal polynomials. We call the latter Schmidt basis sets in analogy to electronic structure methods. The choice of the latter depends on the type of PES to fit: topologically similar PES (e.g. without singularity, such as in the Coulomb potential or in the presence of conical intersections) require similar basis sets. Finally, the core can be pruned as it is effectively sparse. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] This means that every single potential configuration (Hartree product, Ω IJ ) is not necessary for most common applications (reactivity or vibrational spectroscopy). This, in turn, is simply the same idea underlying multiconfigurational methods such as MCSCF. At this point, it is convenient to point out that our experience indicates that the pruning cannot be done through a simple index summation rule.

A SOP-FBR expression can be obtained through a combination of global and local optimisations of an objective function. This is typically defined as the rmse of the energy in the SOP-FBR parameter space ({C j 1 ,...,j f }, {c µ,jκ }). Considering the analysis above, the formally large number of SOP-FBR parameters can be dramatically reduced and, moreover, we have observed that it is possible to freeze the factors (SPPs). Consequently, the whole SOP-FBR process would consist on fine tuning (local optimisation) the core coefficients which have survived the pruning. A good guess for both core (and factors, if needed) can be readily obtained through the use of an approximated level of electronic structure theory such as a semiempirical Hamiltonian (or a reparametrisation thereof). [START_REF] Panadés-Barrueta | Specific Reaction Parameter Multigrid POTFIT (SRP-MGPF): Automatic Generation of Sum-of-Products Form Potential Energy Surfaces for Quantum Dynamical Calculations[END_REF] The CP-FBR approach

In the following, we will consider a high-dimensional problem that might be beyond the reach of both grid-based (PF and related methods) and grid-free (SOP-FBR) Tucker format. In other words, a system for which the Tucker rank ( f κ=1 m κ ) is such that the core tensor exceeds or requires a large fraction of the available memory. Note that for this estimate one should also take into account memory required for the storage of several other quantities, for instance, the wave function. Also the computation of the reference data is not a negligible aspect either. In those cases, the CP format as introduced in MCCPD is of clear advantage over Tucker since typically R ≪ f κ=1 m κ .

Hereafter, we will consider the standard Alternating Least Squares CP algorithm (CP-ALS) for small systems, for which the full grid can be taken into account. More specifically, we have made use of the Python [START_REF] Van Rossum | Python 3 Reference Manual[END_REF] and Numpy [START_REF] Harris | Array programming with Numpy[END_REF] based Tensorly implementation for this. [START_REF] Kossaifi | Tensor Learning in Python[END_REF] For larger systems beyond 6D, we use the recent and powerful Monte Carlo flavour of CP-ALS, the so-called MCCPD. [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into a canonical polyadic decomposition using Monte Carlo methods[END_REF] By construction, MCCPD lives on the primitive grid and the required Monte Carlo sampling points are randomly chosen from it according to a distribution law (e.g. Boltzmann). This becomes an expensive task for large dimensional systems where the number of required sampling points, that is, the set of configurations which is taken as representative of the volume of configuration space accessible for a given set of conditions, becomes large. In addition to this, an integration over all DOF following 1D-cuts [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into a canonical polyadic decomposition using Monte Carlo methods[END_REF] must be performed for each and every sampling point. As a rule of thumb one could consider that the cuts increase the number of points by a factor ∼ 10 × f so that for a 10D problem, on would need to multiply the number of PES evaluations on the sampling points by a factor 100. These aspects are very important whenever modifications to the density of the grid are necessary since this would imply new (full) MCCPD calculation. This is typically the case in two common situations: (i) non converged results for propagation for one (or several) DOF or (ii) a series of simulations requiring new physical, say, initial conditions (e.g. higher energy content of the system). The former aspect is universal to all QD calculations and the latter is more common in scattering studies, e.g. bimolecular collisions or molecule-surface scattering. To circumvent these issues and to achieve an analytical representation of the PES, we have devised a simple two-grid scheme which leads to an analytical PES in Canonical Polyadic form. For this, we use a set of auxiliary basis functions (FBR), hence the name CP-FBR. In what follows, we shall adopt most of the ideas and conventions presented in the previous Section on SOP-FBR. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] In a nutshell, first, it should be emphasized that the CP-FBR is performed on a grid much less dense than the primitive one (aka Discrete Variable Representation, DVR [START_REF] Bacić | Highly excited vibrational levels of "floppy" triatomic molecules: A discrete variable representation -Distributed Gaussian approach[END_REF][32][START_REF] Bacić | Accurate localized and delocalized vibrational states of HCN/HNC[END_REF][34][START_REF] Light | Time-Dependent Quantum Molecular Dynamics[END_REF][START_REF] Light | Discrete variable representations and their utilization[END_REF]). Second, CP-FBR delivers an analytical CP expression from an initial CPD guess computed on this coarse grid. And third, the CP-FBR can be interpolated to primitive grids of any density, thus circumventing the aforementioned problems.

Let us consider a primitive grid (DVR), denoted by Ĩ, as the one needed for numerically converged QD calculations. This grid will consist of f κ=1 N κ grid points. Let us also define a coarser one (I) as in MGPF [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF]. It should be highlighted that this coarse grid is not dense enough for the convergence of the sought QD calculation. Moreover, it will consist of f κ=1 n κ grid points with typically N/n > 2. The underlying hypothesis in the CP-FBR approach is that it is possible to obtain an analytical expression in CP format valid for both fine and coarse grids (see Figure 1):

V CP (Q) = R r=1 λ r f κ=1 v (κ) r (q (κ) ) Q = (q (1)
i 1 , . . . , q

(f ) i f ) (10)
where q

(κ)
iκ is the value of the physical coordinate corresponding to the κ-th DOF at the i κ grid point and Q is a collective representation thereof. More specifically, the CP-FBR factor matrices are expressed by expansion in terms of t κ basis sets (T µ (q (κ) )):

v (κ) r (q (κ) ) = tκ µ=1 c µ,r T µ (q (κ) ) (11)
as it was done for the Tucker SOP-FBR form. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] In other words, it should be possible to express the target PES within a given set of boundary conditions (DVR limits, thus problem dependent) as a weighted sum of R Hartree products of analytical factor matrices (SPPs) expressed in a given finite basis. As previously discussed, the latter are universal for topologically equivalent potentials so that one can apply the same type of basis functions (T µ (q (κ) ))

to topologically equivalent PES (or more generally, tensors). [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] In this sense, a Coulomb potential (r -1 ), owing to its singularity, would not be topologically equivalent to the PES of, say, an isolated molecular system. In what follows, we will describe our current algorithm for the computation of the CP-FBR expression. Note in passing that this analytical CP-FBR form is differentiable ad infinitum.

The actual CP-FBR workflow is presented in Figure 2. It should be noted that when using MCCPD the reference points are usually given as a (DVR) grid-based Metropolis MC trajectory. As observed in the flowchart, one needs first to define a coarse grid (I). For the sake of efficiency, the choice of the coarse grid points should be such that it is significantly smaller than the one needed for the dynamics but sufficiently large to capture a good fraction of the correlation between the different DOFs (PES topography). This condition was achieved in MGPF via an algebraic process (see optimisation of MGPF grids). [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF] In CP-FBR, we obtain a Monte Carlo estimate of the physically meaningful subset of points within our coarse grid using MCCPD. This is analogous to the guess in Tucker SOP-FBR.

Using this information, an initial CP guess on the coarse grid (using either CP-ALS [START_REF] Kossaifi | Tensor Learning in Python[END_REF] or MCCPD [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into a canonical polyadic decomposition using Monte Carlo methods[END_REF]) can be obtained. It should be noted that the use of CP-ALS would imply sampling the coarse grid in full, hence it is limited to small-medium systems. In both cases, after fixing a rank (R), one obtains a series of weights and factors, which will be considered as normalised hereafter. The latter are matrices whose columns correspond to the mapping of a basis function on the coarse grid (see the upper panel of Figure 3). The quality of a CPD decomposition depends on number of configurations in the expansion (rank, R) used to represent the tensor. However a large rank may lead to overfactoring (multilinearity) and, unfortunately, determining the exact rank of a given tensor is not a trivial task and it still constitutes an open question. [START_REF] Kolda | Tensor Decompositions and Applications[END_REF] Furthermore, in the specific case of CP-FBR, an important aspect to consider is the smoothness of the factors on the coarse grid. In this regard, the choice of the rank is also relevant. Indeed, the behaviour of factors may vary dramatically when going beyond the actual rank (see Fig. 3). This is illustrated in Figure 3 where we display the variation of the CP factors, corresponding to the valence angle in the water molecule, upon an increase in the rank (from R=15 to R=30). As observed, they become noisy (highly peaked) thus hindering the fitting step. In fact, none of them is smooth anymore and this would even prevent a simple pruning. Concerning this possibility, the increase in the rank leads to weights presenting a distribution qualitatively similar to that of those from a lower rank expansion (see Figure 4). However, the number of weights above a certain pruning threshold increases with the rank, thus also hindering the aforementioned direct pruning. Fortunately for us, as it will be illustrated, it is possible to obtain analytical CP-FBR PES for relatively small values of the rank (see Results section). The task at hand now is to find a suitable number of fitting functions (e.g. Chebyshev polynomials) that will lead to a satisfactory representation of the factors on the coarse grid while, at the same time, providing a good quality interpolation towards the fine grid (DVR). This was achieved by performing a series of 1D fits for which we used available numpy implementation of least square fits to Chebyshev series of a first kind, numpy.polynomial.chebyshev.chebfit.

It should be clear that by simply increasing the number of fitting polynomials, we consistently improve the quality of the representation on the coarse grid. However, from a certain point, one reaches saturation (excellent fit on the coarse grid) and a further increase will lead to a poorer interpolation (aka overfitting). This is illustrated in Figure 5. The quality of interpolation is assessed through a validation set generated from the fine grid (hence different from the coarse one). For the sake of benchmarking, we have used validation sets which are slightly denser than the ones used for the initial CP decomposition. The objective function to be minimized is rmse on the validation set. For this, we increase sequentially for each DOF (or combinations thereof) the number of fitting functions (the rest been kept fixed) until improvement in the rmse is below, say, 0.01 cm -1 or when it starts to increase thus signaling overfitting (see Fig. 5) The actual process will be illustrated in each of the study systems presented in the Results section. Note that the use of a grid (coarse/fine) is not necessary, one can directly use a list of reference points. As a final remark, it will be shown that the CP-FBR approach preserves to a good degree the rmse on the coarse grid when transferring to the fine grid.

Results and discussion

In the following section we will present the application of our method to several systems of increasing size. As a proof of concept we have considered two benchmark systems: H 2 O (3D) and HONO (6D) as well as a higher-dimensionality application, methane (9D). The Figure 5: Variation of rmse (cm -1 ) upon increase of order of fitting polynomials. Overfitting can be clearly observed.

assessment of the quality of our representation is done through measurements of the quality of the fit (e.g. correlation plots, cumulative rmse plots, etc.) as well as calculations of vibrational eigenstates using the Improved Relaxation algorithm. [START_REF] Meyer | Calculation and selective population of vibrational levels with the Multiconfiguration Time-Dependent Hartree (MCTDH) algorithm[END_REF] Proof of concept, single-well: H 2 O (3D)

As first application, we have considered the water molecule (H 2 O) through the use of the Polyanski-Jensen-Tennyson H 16 2 O spectroscopically refined PES known as PJT2. [START_REF] Polyansky | The potential energy surface of H 2 16 O[END_REF] The definition of coarse and fine grids, together with the the DVR functions is presented in Table 1. A PES tensor mapped on the coarse grid was initially decomposed using the CP routine from the Tensorly library [START_REF] Kossaifi | Tensor Learning in Python[END_REF] using a rank of R = 15. The so-obtained factors were fitted using coarse grid as reference geometries and fine as validation. As observed, the fine grid is about 13 times the size of the coarse one.

The two O-H bonds were taken as equivalent and, consequently, the order of the fitting Chebyshev polynomials for both DOF was kept identical. The angular basis was optimized independently in a subsequent step. The variation of the rmse at every iteration of the optimisation is displayed in Figure 6. It can be observed the differential behaviour of the rmse upon increase of the of order of the fitting polynomials for different DOF. This has In the upper part of the graph is specified which DOF is being optimized within a given iteration. Final values are shown in Table 1 been observed in other grid based methods such as MGPF (see Table IV in Ref. [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF]. The global rmse of the CP-FBR PES (0.686 cm -1 ) compares well to the value from our initial CPD decomposition (0.281 cm -1 ). A more thorough test of the quality of the CP-FBR PES is obtained through the correlation and cumulative rmse plots, presented in Fig. 7 and Fig. 8, respectively. To further assess the quality of the PES, we computed the geometry corresponding to the minimum of the CP-FBR PES using steepest descend method. The obtained geometry compares very well with the one obtained using original PES, with an agreement up to the fourth decimal (Table 2). It should be highlighted that in our SOP-FBR approach, CP-FBR in this case, no a priori explicit knowledge of the physics of the problem (e.g. energy barriers, minima) is utilised in the optimisation process. Further assessment of the quality of our CP-FBR PES has been achieved through the calculation of vibrational eigenstates and their comparison with the exact ones (exact expansion) calculated under the same conditions (SPFs, integrator, etc.). It should be mentioned that our grid definition is not meant to yield spectroscopically accurate eigenstates. Our aim here is simply to compare the results obtained using our analytical CP-FBR approach with an equivalent grid-based tensor-decomposed (refitted) PES. First, the ZPE value from CP-FBR of 4644.238 cm -1 compares very well with the reference one of 4644.234 cm -1 as well as the vibrational eigenstates which are in excellent agreement with the reference (exact) calculations (see Table 3). Benchmark two-well system: HONO (6D)

H O O N r OH r ON r N=O θ 2 θ 1 φ
We focus now on a more general problem, the 6D two-well system HONO molecule in the cis-trans izomerization region. In recent years this has become a standard benchmark. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF][START_REF] Pradhan | A ground state potential energy surface for HONO based on a neural network with exponential fitting functions[END_REF] As source of reference data, we use the CCSD(T)/cc-pVQZ-quality analytical PES of Richter et al. [START_REF] Richter | A study of mode-selective trans-cis isomerisation in HONO using ab initio methodology[END_REF] In Figure 9 we present the definition of the polyspherical coordinates [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF] as those used in Refs. [START_REF] Richter | Time-dependent wavepacket study on trans-cis isomerisation of HONO[END_REF][START_REF] Richter | Time-dependent wave packet study on trans-cis isomerisation of HONO driven by an external field[END_REF]. Concerning the angles, Θ 1 and Θ 2 they were used as their cosines, u 1 and u 2 , respectively. The definition of the coarse grid is given in Table 4. In this case, the guess decomposition (on the whole coarse grid) was done using the CP-ALS routine from the Tensorly library. This is comparable to our previous SOP-FBR calculation on the same system [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF]. The CP rank was set to R=200. The validation set consisted of 10 6

geometries uniformly sampled from a slightly bigger grid defined in such a way to avoid the coincidence with the coarse grid. The initial order of Chebyshev fitting basis was 4 for all DOFs. In Figure 10, we present the variation (decrease) of the rmse upon optimization (fitting to higher order) for each DOF.

As initial remark, we observe that the rmse is already of chemical accuracy (∼ 350 cm -1 )

for the lowest number of basis sets. The final rmse is 2.213 and 2.929 cm -1 on the reference and the validation sets, respectively. A proof of the stability of our optimisation scheme is that the final expansion order (see Table 4) nicely agrees (within a numerical error) to the one obtained with the grid based method MGPF. [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF] To further prove our hypothesis of interpolation to any fine grid, we defined the (primitive) grid shown in Table 4 and uniformly sampled 10 6 geometries from it. From the correlation plot given in Figure 11 can be seen that CP-FBR values are in good agreement with reference ones in entire range of energies (up to 10 5 cm -1 ). The cumulative rmse plot given in Figure 12 confirms that value of 5.345 cm -1 , even though slightly bigger, compares well with values on reference and validation sets and is in good agreement with other methods. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF][START_REF] Pradhan | A ground state potential energy surface for HONO based on a neural network with exponential fitting functions[END_REF] E CP_FBR / cm -1 To further assess the topography of our CP-FBR PES, we located the stationary points.

For the sake of efficiency, we imposed the torsion to its optimal value and optimized other coordinates with steepest descent method. Final gradients were all below 10 -7 hartree/(bohr/radian).

As it can be seen in Table 5, the obtained values agree very well with the reference ones, up to the fourth decimal for all DOFs, even though none of these geometries was explicitly included in the reference data sets. Finally, we used the primitive DVR grid (analogous to that of Refs. [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF][START_REF] Pradhan | A ground state potential energy surface for HONO based on a neural network with exponential fitting functions[END_REF]) to calculate the ZPE and several vibrational eigenstates. The obtained ZPE of 4365.39 cm -1 compares well the reference (exact) value of 4367.73 cm -1 . As for the vibrational states, values given in Table 6 also exhibit a good agreement.

Vibrational eigenstates of CH 4 using an analytical full-dimensional (9D) CP-FBR PES As a final example, we consider a 9D system, the methane molecule. We used the PES from Zhang et al. [START_REF] Zhang | Full-dimensional quantum mechanical calculations of the reaction probability of the H + CH 4 reaction based on a mixed Jacobi and Radau description[END_REF] as reference. The coordinates system used is based on a mixed Radau (for non reacting CH 3 group) and Jacobi (for the remaining C-H bond) polyspherical coordinates [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF]. This follows reference [START_REF] Zhang | Full-dimensional quantum mechanical calculations of the reaction probability of the H + CH 4 reaction based on a mixed Jacobi and Radau description[END_REF] (see Fig. 13). Again, for the angles β s , θ 1 and θ 2 we use the corresponding cosines, labeled as u b , u 1 and u 2 , respectively.

Due to the size of a system, decomposition of a tensor taking into account the full grid, even for a coarse grid, had memory requirements which rendered it impractical, at least for the current CP-ALS implementation of Tensorly. For this reason, here we made use of the powerful MCCPD algorithm to perform the decomposition on the coarse grid.

MCCPD used a set of Metropolis MC sampled grid-points of that grid. To focus on the low and 1500 cm -1 to generate set of 15000 points in total, from the coarse grid. According to our experience, an MCCPD performed on such a small number of points seems to lead to a not fully converged CP representation. By this, we mean that it predicts more accurately the reference energies than any other geometry out of this reference set, even within the same domain of energies. Interestingly, such a CPD can still be used as initial guess for our CP-FBR fitting/interpolation. For the sake of assessment of the quality of our approach, we have decided to define three different grids. Note that only two are necessary. These grids are: (i) coarse which refers to the sampling points required for the MCCPD as well as for the fitting of the factors; (ii) validation from which 11500 points was sampled and used to control overfitting; and (iii) primitive (DVR) grid which is the one needed for the quantum dynamical simulation to be converged. Relevantly, neither the coarse grid nor the validation are subsets of the primitive grid. It should be clear that the primitive grid is not used at all in the CP-FBR process and hence our CP-FBR ansatz will be simply interpolated onto the primitive grid, which is also the denser one. For the sake of comparison, we performed a MCCPD directly on the primitive grid as well, using 70000 sampling points and keeping the rank the same, at 250 Finally, as an extra (strictly speaking unnecessary) assessment, we sampled (same conditions) a 10-times bigger number of geometries just for the sake of further testing the obtained CP-FBR form. Concerning the factor fitting, we took into account the symmetry of the molecule by treating on the same footing the Radau-defined C-H bonds and, thus, optimizing them simultaneously. The same was applied to u 1 and u 2 . We estimated the quality of the resulting CP-FBR on the coarse, fine, and primitive grids. In Figure 14 we present the correlation plot for the primitive grid. It can be clearly observed that the agreement between CP-FBR and original surface is best for the lowest energies. The same behaviour is observed in Figure 15, where we present the comparison of the cumulative rmse plots for CP-FBR on all of the three grids together with MCCPD for the coarse and the primitive. Both methods show a good agreement in the targeted region up to 17000 cm -1 .

H 3 H 1 H 2 H 4 C R 1 R 2 R 3 R 4 θ 2 θ 1 φ β s γ s G s G CH3 G H3
As a final test, we calculated the ZPE and several vibrational eigenstates of methane and compared our CP-FBR results with the ones obtained with MCCPD (on the primitive grid), those from the reference values (see Ref. [START_REF] Zhang | Full-dimensional quantum mechanical calculations of the reaction probability of the H + CH 4 reaction based on a mixed Jacobi and Radau description[END_REF]), as well as experimental values. The values of energies obtained for vibrational states are presented in Table 8. Our ZPE value of 9681.495 cm -1 is slightly lower than the one from the reference (9689.43 cm -1 ) and compares well with the one from MCCPD (9681.07 cm -1 ). Concerning the vibrational eigenvalues, CP-FBR is in good agreement with all of our reference values, particularly with MCCPD. Finally, the CPU time for MCCPD on the primitive grid was 24m:06s whereas for the whole CP-FBR process (i.e. MCCPD decomposition on the coarse grid and subsequent FBR fitting) took 12m:25s. It should be noted that our current CP-FBR implement is sub-optimal and not fully automated. decomposition can be taken. An important advantage of CP-FBR is that it enables direct interpolation between a coarse grid and the target one (DVR) needed for dynamics. This is achieved through fit of the factors using a set of auxiliary basis functions. CP-FBR has been tested in systems of increasing size: water molecule (3D), HONO (6D) and methane (9D). In all three cases, it was observed that the rmse was well preserved when transferring between the coarse and primitive grids. The quality of the obtained expressions was carefully assessed through correlation and cumulative rmse plots, topographical analysis, as well as calculation of vibrational eigenstates with MCTDH. In sum, CP-FBR constitutes a further proof of the possibility of obtaining the high-dimensional analytical representation of potential-like operators, as epitomized by the PES. As shown previously in the case of a Tucker format, the CP SPPs (factor matrices) can also be interpolated from a coarse (cheap, small) grid to the primitive one (needed for QD). In a sense, this operation generalizes the seminal ideas underlying the chnpot method of MCTDH. It does so taking into account global information of coarse and fine grids, whereas chnpot is local and bound to a grid. Our (CP-)FBR approach leads to a fully general expression and multidimensional function with a topography that is in very good agreement with the reference level of electronic structure theory. In addition to this, our X-FBR approaches (X=Tucker, CP) provide a conceptual frame linking the generation of analytical (aka routine) high-dimensional PES in SOP form using as reference an amount of information even coarser than the MC sets generated from the primitive grid. We are currently using CP-FBR for the fit of high-dimensional unbound PES. A future development will concern the efficient and parallel computation of energies for single geometries using our CP-FBR form. This will enable the optimisation using lists of points generated from the full set of stationary points as some of us already suggested. [START_REF] Panadés-Barrueta | Specific Reaction Parameter Multigrid POTFIT (SRP-MGPF): Automatic Generation of Sum-of-Products Form Potential Energy Surfaces for Quantum Dynamical Calculations[END_REF] Our software will be freely available and interfaced to the MCTDH software package.

Figure 2 :

 2 Figure 2: CP-FBR flowchart for the generation of a PES representation for quantum dynamics

Figure 3 :

 3 Figure 3: Effect of rank in the smoothness of factors for the illustrative case of the valence angle in water. Upper panel, discrete values of the factors on the coarse grid points. Lower panel, interpolated lines on those points for the sake of better visualization. The phase of the factors is random and has not been corrected. Upon increase of the rank, the factors lose their smooth behaviour.

Figure 4 :

 4 Figure 4: Effect of rank in the distribution of weights. Upon increase of the rank, more weights lie above a given threshold thus preventing efficient pruning.

Figure 6 :

 6 Figure6: Variation of the rmse (cm -1 ) upon variation of the number of fitting functions. In the upper part of the graph is specified which DOF is being optimized within a given iteration. Final values are shown in Table1

Figure 7 : 2 OFigure 8 :

 728 Figure 7: Correlation plot on the fine grid for H 2 O

Figure 9 :

 9 Figure 9: Definition of the coordinates for HONO

Figure 10 :

 10 Figure10: Variation of the rmse (cm -1 ) upon variation of the number of fitting functions. In the upper part of the graph is specified which DOF is being optimized with a given iteration.

Figure 11 :Figure 12 :

 1112 Figure 11: Correlation plot on the primitive grid for HONO

Figure 13 :

 13 Figure 13: Coordinate system for CH 4 .

Figure 14 : 4 Figure 15 :

 14415 Figure 14: Correlation plot on the primitive grid for CH 4Figure 15: Cumulative rmse (cm -1 ). Note that the black trace (MCCPD on coarse grid) lies right below the red trace FBR on coarse grid).

Table 1 :

 1 Definition of the DVR grid for water molecule. O-H bond distances (r 1 and r 2 ) are expressed in bohrs and the H-O-H angle (θ) in radians. The number of DVR points is given for both coarse (n) and primitive (N) grids. Fit indicates the Chebyshev order.

	Coordinate DVR n N	Range	Fit
	r1	sin 15 34 1.4500, 2.4500 13
	r2	sin 15 34 1.4500, 2.4500 13
	θ	Leg 20 50 1.1400, 3.1416 16
		r 1 and r 2	Ѳ

Table 3 :

 3 Comparison of the 20 lowest vibrational eigenvalues of H 2 O for the original PES (POTFIT) and the CP-FBR one.

	State Eigenenergies (cm -1 )
		CP-FBR POTFIT
	ZPE	4644.24	4644.23
	1	1594.23	1594.23
	2	3151.68	3151.68
	3	3681.12	3681.14
	4	3779.57	3779.51
	5	4670.76	4670.75
	6	5256.67	5256.67
	7	5353.37	5353.33
	8	6149.39	6149.38
	9	6796.91	6796.90
	10	6894.15	6894.13
	11	7278.66	7278.63
	12	7333.69	7333.59
	13	7497.74	7497.64
	14	7581.88	7581.87
	15	8299.68	8299.65
	16	8402.48	8402.46
	17	8831.88	8831.83
	18	8884.26	8884.18
	19	8944.09	8944.10
	20	9049.27	9049.21

Table 4 :

 4 Definition of the DVR grid for the HONO molecule. Bonds are given in bohrs and angles in radians. The number of DVR points is given for both coarse (n) and primitive (N) grids. Fit represents the Chebyshev order

	Coordinate DVR n N	Range	Fit
	r N =O	HO	7 13 1.90, 2.60	5
	r ON	HO	7 16 2.10, 3.25	6
	r OH	HO	8 18 1.30, 2.45	7
	u 1 (cos θ 1 )	HO	7 13 -0.65, -0.10 5
	u 2 (cos θ 2 )	HO	7 18 -0.65, 0.25	5
	ϕ	cos 15 32	0,π	10

Table 5 :

 5 Comparison of the geometrical parameters optimized with the steepest descent method on the CP-FBR and original PESs. All final gradients were below 10 -7 .

				Coordinate
	PES	r N O	r ON	r OH	θ 1	θ 2	ϕ
				Trans
	reference 2.2133 2.6967 1.8229 1.9315 1.7776 3.1416
	CP-FBR 2.2134 2.6972 1.8235 1.9315 1.7776 3.1416
				Cis
	reference 2.2374 2.6312 1.8416 1.9752 1.8223	0
	CP-FBR 2.2376 2.6304 1.8427 1.9753 1.8223	0
				TS
	reference 2.2009 2.8475 1.8180 1.9287 1.7577 1.5080
	CP-FBR 2.2008 2.8475 1.8185 1.9286 1.7571 1.5080

Table 6 :

 6 Comparison of the lowest vibrational eigenstates of HONO for the original PES and the CP-FBR one. ZPE from CP-FBR is 4365.39 and from reference is 4367.73 cm -1 . The states are denoted with an ordinal and the label t (trans) or c (cis).

	State Eigenenergies (cm -1 )
		CP-FBR Reference
	ZPE	4365.39	4367.73
	1c	92.79	94.05
	2t	600.75	600.81
	3c	709.15	710.65
	4t	795.51	795.86
	5c	942.67	944.12
	6t	1055.43	1055.38
	7t	1187.77	1188.02
	8t	1264.30	1264.89
	9c	1305.91	1306.59
	10c	1311.22	1312.78
	11t	1384.82	1385.35
	12c	1403.28	1404.80
	13c	1546.33	1547.81
	14t	1574.27	1574.92
	15t	1640.79	1640.85
	16t	1690.37	1689.88
	17c	1725.52	1725.95
	18t	1762.06	1762.68
	19c	1778.10	1779.61
	20t	1828.65	1828.94

Table 7 :

 7 Coordinates and DVR functions used to define the primitive grid. Distances are given in Bohrs and angles in Radians. The number of DVR points on the coarse grid (n) and primitive grid (N). The column Fit presents the Chebyshev order

	Coordinate DVR n N	Range	Fit
	ub	sin	7 14 -0.85 0.3	6
	γ s	sin	7 17 4.5	6	16
	R 1	HO 7 10 1.5	2.8	6
	R 2	HO 7 10 1.5	2.8	6
	R 3	HO 7 10 1.5	2.8	6
	R 4	sin	7 14 1.6	2.9	6
	u 1	sin	7 14 -0.9 0.25 7
	u 2	sin	7 14 -0.9 0.25 7
	ϕ	sin	7 19 1.3	3	7
	energy region, we used Metropolis sampling at 3 different temperatures (k B T): 500, 1000

Table 8 :

 8 Obtained energies (in cm -1 ) for vibrational states (J=0), comparison with values from the reference and experimental ones

	State	CP-FBR MCCPD Reference[44] Experimental[45]
	ZPE	9681.50	9680.97	9689.43
		1308.92	1309.78	1309.19
	(0001)F 2	1309.85	1310.27	1309.83	1310.76
		1310.35	1310.91	1310.33
	(0100)E	1529.56 1530.87	1531.10 1531.22	1530.40 1530.83	1533.33
	(0002)A 1 2588.06	2588.85	-	2587.04
	(0002)F 2 2611.81	2612.78	-	2614.26

Acknowledgements

The authors are very thankful to H.-D. Meyer and M. Schröder (Heidelberg) for the careful reading of the manuscript, their useful remarks on MCCPD, and fruitful discussions in general. NN is grateful for its Ph.D. funding from the University Paris-Saclay. The authors are very thankful to J.-Y. Bazzara, A. Borissov, and P. C ¸arçabal (ISMO) for their computer support. As a final remark, the authors would like to publicly acknowledge the support and efforts of the Fortran, Python, Numpy, Tensorly, and Jupyter communities which have much facilitated our work.

Availability of data

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Summary and conclusions

We have presented an approach for obtaining the canonical polyadic counterpart of Tucker SOP-FBR, [START_REF] Panadés-Barrueta | Low-rank sum-of-products finite-basis-representation (SOP-FBR) of potential energy surfaces[END_REF] hence the name CP-FBR. Our method yields an analytical CP form for multidimensional tensors. CP-FBR requires a guess for which a non-converged initial CP