Ari Pereira

Joachim Knapik

Ahai Chen

David Lauvergnat

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

Quantum molecular dynamics simulations of the effect of secondary modes on the photoisomerization of a retinal chromophore model

teaching and research institutions in France or abroad, or from public or private research centers.

Quantum molecular dynamics simulations of the effect of secondary modes on the photoisomerization of a retinal chromophore model

I. INTRODUCTION

Quantum molecular dynamics simulations provide extremely detailed information on the real-time evolution of molecular electron-nuclear systems, which are subject, for instance, to photo-excitation and nonadiabatic effects. There exist various flavors of simulation techniques used to propagate nuclear dynamics under the effect of electrons in the ground state as well as in the excited states [START_REF] González | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Crespo-Otero | [END_REF][3]. Those techniques range from wavepacket propagations with time-independent and time-dependent basis sets [START_REF] Bonfanti | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Worth | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Richings | [END_REF][7], to trajectory-guided, or trajectory-based, methods [8][START_REF] Kirrander | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Vacher | [END_REF][START_REF] Curchod | Quantum Chemistry and Dynamics of Excited States: Methods and Applications[END_REF][START_REF] Makhov | [END_REF][13][14][15][16] that employ (quasi-)classical, semi-classical or quantum trajectories. While the computational cost for solving the timedependent Schrödinger equation numerically exactly scales exponentially with the number of degrees of freedom, the problem can be mitigated with suitable choices of sparse grids and associated basis sets for purely quantum techniques or by introducing a classical-like description of the nuclei while maintaining a quantum perspective on the electronic degrees of freedom.

The purpose of the present work is to focus on some of those methods, discussing their performance, their limitations and their advantages, by applying them to the simulation of the photoisomerization process of a retinal chromophore model [17] including the effect of an environment of growing complexity [18]. The model Hamiltonian of the system represents the photo-induced 11-cis to all-trans isomerization of retinal in rhodopsin [17] (a simplified, and wellestablished, precursor of the model [19] studied in Ref. [20]) using two electronic states and two effective nuclear modes. The effect of the environment is accounted for by including secondary Raman-active modes [18,21] represented by a bath of harmonic oscillators coupled to the electronic states.

Our quantum molecular dynamics simulations are performed by propagating vibronic wavepackets and classical-like trajectories. For the former, we will focus on the use of different ideas to define the time-independent grid, i.e., direct product versus Smolyak [22,23] grids, in order to be able to include efficiently at a reasonable computational cost the secondary modes representing the environment. The calculations are performed with the ElVibRot code [START_REF] Lauvergnat | ElVibRot-TnumTana Quantum Dynamics Code[END_REF]. For the latter, we will test and compare Tully surface hopping [START_REF] Tully | [END_REF]26] (TSH), an independent-trajectory method, to its coupled-trajectory formulation [27][28][29] derived from the exact factorization of the electron-nuclear wavefunction [16,[30][31][START_REF] Agostini | Quantum chemistry and dynamics of excited states: Methods and applications[END_REF][START_REF] Agostini | Handbook of Materials Modeling[END_REF], namely the coupled-trajectory Tully surface hopping [START_REF] Pieroni | [END_REF]35] (CT-TSH) algorithm. In particular, we will present an in-depth analysis of the computational cost of independent-trajectory and coupled-trajectory schemes as the number of degrees of freedom increases. The trajectory-based methods employed in this work are implemented in the G-CTMQC code [36]. Both ElVibRot and G-CTMQC are interfaced with QuantumModelLib [START_REF] Lauvergnat | QuantumModelLib[END_REF], a library of analytical model potentials which includes the model system-bath Hamiltonian [18] used in this work.

The paper is organized as follows. In Section II, we introduce the theory by comparing the Born-Huang and the exact-factorization formulation of nonadiabatic dynamics, while in Section III, the numerical procedures used in our quantum dynamics simulations are presented. The model system used as test case in this work is described in Section IV. Section V provides information on the computational aspects of our simulations, whose results are presented in Section VI. Our conclusions are stated in Section VII.

II. THE TIME-DEPENDENT SCHRÖDINGER EQUATION

The quantum-mechanical time evolution of a system of interacting electrons and nuclei is governed by the timedependent Schrödinger equation (TDSE)

i ∂ t Ψ(r, R, t) = Ĥ(r, R)Ψ(r, R, t) (1)
with the non-relativistic electron-nuclear Hamiltonian Ĥ(r, R), and whose solution is the time-dependent molecular wavefunction Ψ(r, R, t). Electronic and nuclear coordinates are collectively indicated with the symbols r and R, respectively. The Hamiltonian is the sum of a nuclear kinetic energy operator Tn (R) = Nn ν=1 (-2 ∇ 2 )/(2M ν ) -in Cartesian coordinates -for the N n nuclei with masses M ν , and of the electronic Hamiltonian Ĥel (r, R), which is, in turn, the sum of the electronic kinetic energy and all interactions.

When looking at nonadiabatic processes, thus including the effect of electronic states and their couplings, it is convenient working with the Born-Huang representation or with exact factorization [3,30,31], by writing the molecular wavefunction as

Ψ(r, R, t) = l χ l (R, t)ϕ l (r; R) (2) = χ(R, t)Φ(r, t; R) (3)
Equation ( 2) is the Born-Huang expansion, and expresses the time-dependent molecular wavefunction in terms of electronic adiabatic states ϕ l (r; R), i.e., the eigenstates of the electronic Hamiltonian Ĥel (r, R), with parametric dependence on the nuclear positions R. The expansion coefficients are the time-dependent nuclear amplitudes χ l (R, t).

When inserted into the TDSE (1), the evolution equation for each nuclear amplitude reads

i ∂ t χ l (R, t) = m Tn (R) + E l (R) δ lm - ν 2 2M ν Ĉlm,ν (R) χ m (R, t) (4)
where E l (R) is the R-dependent eigenvalue of Ĥel (r, R) for state l, also known as adiabatic, or Born-Oppenheimer (BO), potential energy surface (PES), satisfying the electronic time-independent Schrödinger equation

Ĥel (r, R)ϕ l (r; R) = E l (R)ϕ l (r; R) (5)
and Ĉlm,ν (R) is a nonadiabatic operator that couples the dynamics of the l-th nuclear amplitude to that of the m-th nuclear amplitude. The nonadiabatic operator is Ĉlm,ν (R) = D lm,ν (R) + 2d lm,ν (R) • ∇ ν , and depends on the nonadiabatic coupling D lm,ν (R) = ∇ ν ϕ l (R)|∇ ν ϕ m (R) r as well as the nonadiabatic coupling vectors (NACVs)

d lm,ν (R) = ϕ l (R)|∇ ν ϕ m (R) r .
The symbol • | • r indicates an integration over the electronic positions, and for this reason we removed all dependencies on r within the brackets. Equation ( 3) is the exact factorization [30,31], which expresses the time-dependent molecular wavefunction as a product of the time-dependent nuclear wavefunction χ(R, t) and of the time-dependent conditional electronic amplitude Φ(r, t; R) -with parametric dependence on nuclear positions R. Note that, conversely to the Born-Huang expansion, the electronic sub-system is explicitly time dependent in the exact factorization. One can easily recover Eq. (2) from Eq. (3), by representing the electronic amplitude using a Born-Huang-like expansion, i.e., Φ(r, t; R) = l C l (R, t)ϕ l (r; R), and identifying the nuclear amplitudes as χ l (R, t) = χ(R, t)C l (R, t). When inserted into the TDSE (1), and using the partial normalization condition Φ(t; R)|Φ(t; R) r = 1 ∀ R, t, the coupled nuclear and electronic evolution equations are derived

i ∂ t χ(R, t) = Nn ν [-i ∇ ν + A ν (R, t)] 2 2M ν + (R, t) χ(R, t) (6)
i ∂ t Φ(r, t; R) = Ĥel (r, R) + Ûen [Φ, χ] -(R, t) Φ(r, t; R) (7)
We introduced here the time-dependent vector potential,

A ν (R, t) = Φ(t; R)| -i ∇ ν Φ(t; R) r , and the time- dependent potential energy surface, (R, t) = Φ(t; R)| Ĥel (R) + Ûen [Φ, χ] -i ∂ t |Φ(t; R) r
, driving the evolution of the nuclear wavefunction via a nuclear TDSE [START_REF] Richings | [END_REF]. The electron-nuclear coupling operator,

Ûen [Φ, χ] = ν ( [-i ∇ν -Aν (R,t)] 2 2Mν + 1 Mν [ -i ∇ν χ(R,t) χ(R,t) + A ν (R, t)][-i ∇ ν -A ν (R, t)],
induces the coupling to the nuclear dynamics in the electronic evolution equation (7), as it is explicitly dependent on the nuclear wavefunction χ(R, t), and via the time-dependent vector potential, A ν (R, t), also contains implicitly the electronic wavefunction Φ(r, t; R). Quantum molecular dynamics simulations performed in this work are based either on the Born-Huang representation or on the exact factorization. More specifically, vibronic wavepacket dynamics is carried out using the electronic diabatic representation, rather than the adiabatic representation, where essentially in Eq. ( 4), the BO PESs are replaced by the diabatic PESs and the nonadiabatic operator is replaced by the diabatic potential coupling. The transformation from diabatic-to-adiabatic basis allows us to reconstruct information along the dynamics in the adiabatic representation, for a direct comparison with trajectory-based results. Trajectory-based calculations are performed in the adiabatic basis, and even when working with the exact factorization, the electronic conditional amplitude is represented as a Born-Huang-like expansion as discussed above. More details on the trajectory-based algorithms will be given in Section III.

III. QUANTUM AND QUANTUM-CLASSICAL ALGORITHMS

A. Vibronic wavepacket dynamics For a system described by n degrees of freedom and in a standard quantum approach (either time-dependent or time-independent), the wavefunction or wavepacket are expanded in a single direct-product:

S DP λ = S 1 λ1 ⊗ S 2 λ2 ⊗ • • • ⊗ S n λn (8)
In Eq. ( 8), the

S i λi (with i ∈ [1 • • • n]
) can be a finite basis set (B i λi with nb i λi basis functions) or a finite quadrature grid (G i λi with nq i λi grid points and weights). The index, λ i , enables to control the sizes of the i th basis set and grid (see below for their relations). Tipically, nb i λi and nq i λi values are larger than 10, so that the total number of basis functions, N b or grid point N q are larger than 10 n .

In 1963 [START_REF] Smolyak | [END_REF], Smolyak proposed two equivalent schemes to overcome the exponential scaling of the direct-product approach when n increases. With respect to contraction schemes [39,40] or others basis set selections or pruning approaches [41], the Smolyak schemes enable to use selected basis sets and grids in a consistent way. [22,[42][43][44][45][46][47][48][49][50][51][52][53][54][55] Furthermore, in the second scheme, Smolyak has shown that a large direct-product can be substituted by a sum of small direct-products:

S Srep L = R(λ) D λ • S DP λ (9)
where R(λ) are constraints on the λ i which enable to control the number of the small direct-products in the Smolyak expansion. The usual constraints are defined as L-n+1 ≤ |λ|≤ L where L is the Smolyak parameter and |λ|= n i=1 λ i . Furthermore, the D λ coefficients are (-1) L-|λ| C L-|λ| n-1 and C j i are binomial coefficients. However, in the present study, we have used constraints adapted to a system-bath separation [55] with three Smolyak parameters: (i) L b , a parameter associated to the bath modes (i

∈ [3 • • • n]), such i=3,n λ i ≤ L b . (ii) L s ,
a parameter associated to the system modes (the first and the second modes), such i=1,2 λ i ≤ L s . (iii) L m , a parameter which controls the coupling between the system and bath modes, such i=1,n λ i ≤ L m . Furthermore, the Smolyak parameters are subject to other constraints: L m ≥ L s ≥ L b and L m ≤ L s + L b . This approach enables us to tune the basis set expansion of the system part (the first two modes) and the bath part (the remaining modes). With this system-bath separation, D λ coefficients are obtained from the transformation between the first and the second Smolyak scheme. Furthermore, to improve the integration, the value of Smolyak parameter associated to the grid, L G , can be larger than L m and typically,

L G = L m + 1 or L G = L m + 2.
Finnally, nb i λi and nq i λi can be defined by any increasing sequences in λ i . In the present study, their expressions are given by:

nb i λi = nq i λi = A i + B i • λ i (10)
Then, a state vector |Ψ can be represented as a "weighted" sum of Smolyak contributions Ψ λ and each contribution is expanded on a small direct-product,

S DP λ |Ψ = R(λ) D λ • Ψ λ (11)
This representation as a sum of Smolyak contributions ( Ψ λ ) of a state vector is fundamental so that the action of an operator, Ô, on a quantum state, |Ψ , can be split on the different Smolyak contributions.

To summarize, in the Smolyak scheme, several direct-products, S DP λ , are used and due to the constraint on λ, those direct-products are small. Indeed, when one λ i is large (for instance equal to L s ), all the other ones are equal to zero. In the direct-product approach, one large direct-product is used, S DP Ls with all λ i being equal to L s . Furthermore, the Smolyak scheme overcomes the usual exponential scaling with system size (it becomes polynomial of degree L), which is the main advantage of the approach with respect to the standard direct-product one. As the direct-product approach, the Smolyak schemes are able to use the full grid, therefore, the diabatic-to-adiabatic transformation of the density is possible. In the present study, this feature will be used to analyse the wavepacket density along the propagation. In particular, the adiabatic populations and other properties will be calculated.

B. Trajectory-based approaches

Trajectory-based -quantum-classical -simulations are performed using the independent-trajectory method Tully surface hopping (TSH) [START_REF] Tully | [END_REF], originally derived starting from the Born-Huang representation of the molecular wavefunction, and its coupled-trajectory formulation, namely the coupled-trajectory Tully surface hopping (CT-TSH) [START_REF] Pieroni | [END_REF]35], derived from the exact factorization.

These methods are characterized by the fact that classical trajectories are used to mimic the quantum nuclear dynamics. A trajectory α is identified by a collection of 3N n coordinates, indicated henceforth with the symbol R α (t). Along each trajectory, and coupled to the nuclear evolution, the electronic time-dependent state of the system evolves nonadiabatically as a superposition of adiabatic states. In such a superposition, the expansion coefficients are indicated as C α l (t), for each electronic state l included and along each trajectory α. In TSH, the electronic coefficients evolve according to

Ċα l,TSH (t) = - i E α l C α l (t) - m Nn ν=1 Ṙα ν (t) • d α lm,ν C α m (t) (12)
where E α l is the value of the l-th BO PES at the position R α (t), d α lm,ν is the value of the NACV at the position R α (t), and Ṙα ν (t) is the velocity of the trajectory α at time t. In CT-TSH, the electronic evolution equation ( 12) is corrected by a term accounting for the coupling among the trajectories, namely

Ċα l,CT-TSH (t) = Ċα l,TSH (t) + Ċα l,qm (t) (13)
The additional term depends on the so-called quantum momentum (qm)

P α ν (t) = (-∇ ν |χ(R α (t), t)| 2 )/(2∇ ν |χ(R α (t), t)| 2 ) as Ċα l,qm (t) = Nn ν=1 P α ν (t) M ν • f α ν,l -A α ν (t) C α m (t) (14)
with f α l,ν = t 0 (-∇ ν E α l ) dτ the accumulated adiabatic force and the time-dependent vector potential approximated as

A α ν (t) = m |C α m (t)| 2 f α ν,m .
Given the expression of the quantum momentum in terms of the nuclear density, it appears clear that the nuclear density has to be reconstructed along the dynamics at each time from the distribution of the classical trajectories. Therefore, the trajectories need to be propagated altogether and are, thus, coupled.

In TSH, the trajectories are propagated adiabatically, namely

F α ν,TSH (t) = -∇ ν E α a ( 15
)
where the active state a is selected stochastically at each time step of dynamics according to the fewest switches probability [14]

P α m→l (t) = max   0, - 2Re Cα m (t)C α l (t) ν d α lm,ν • Ṙα ν (t) |C α m (t)| 2 dt   (16)
which gives the probability for the trajectory α to hop from the BO PES m to the BO PES l from time t to time t+dt. TSH results presented in Section VI are obtained including the so-called energy decoherence (ED) correction [26], we will thus referred to TSH-ED henceforth. In this case, the electronic coefficients along the trajectory α associated to the non-active states, say m, are exponentially damped over a typical decoherence time

τ α m = ( /|E α m -E α a |)(1+C/T α n ), where C = 0.1 E h and T α
n is the instantaneous nuclear kinetic energy of the trajectory α. In CT-TSH, the trajectories are propagated adiabatically as in TSH-ED, thus the force is the same as in Eq. ( 15), with hopping probability given by Eq. ( 16). In CT-TSH, decoherence effects are automatically taken care of by the quantum-momentum term in Eq. ( 13), therefore, as discussed above, the trajectories have to propagated simultaneously in order to evaluate the nuclear density from the distribution of the classical trajectories.

IV. PHOTOISOMERIZATION MODEL INCLUDING SECONDARY MODES

Quantum molecular dynamics simulations employing the algorithms described in Section III are performed on a model describing the photoisomerization process of the retinal in rhodopsin, which accounts for the effects of the environment via the inclusion of secondary modes in a system-bath formulation, where the bath is treated harmonically [17,18,21].

The original model developed by Hahn and Stock [17] is a a two-electronic-state two-dimensional model for the photo-induced 11-cis to all-trans isomerization of retinal in rhodopsin. The effective nuclear modes defining the model are an angular reaction coordinate ϕ and a collective vibration q characterizing "a delocalized stretching motion of the polyene chain, whereby single and double bonds interchange" [17]. Electronic ground (S 0 ) and first-excited (S 1 ) singlet states present a conical intersection (CI) at ϕ CI ± π 2 , q CI = 0, and are represented in Fig. 1 (upper panel). These BO PESs are obtained by diagonalization of the electronic Hamiltonian Ĥel (ϕ, q) given in the diabatic basis

Ĥel (ϕ, q) = V 00 (ϕ, q) V 01 (q) V 10 (q) V 11 (ϕ, q) (17)
The full Hamiltonian of the system reads

ĤS = P 2 ϕ 2m + P 2 q 2ω -1 + Ĥel (ϕ, q) (18)
where Pϕ and Pq are the momentum operators along the two nuclear coordinates, and m and ω -1 the corresponding masses. Note that ϕ and q are dimensionless variables. The elements of the electronic Hamiltonian in the diabatic basis are

V 00 (ϕ, q) = 1 2 W 0 (1 -cos ϕ) + 1 2 ωq 2 (19) V 11 (ϕ, q) = E 1 - 1 2 W 1 (1 -cos ϕ) + 1 2 ωq 2 + κq (20)
V 01 (q) = V 10 (q) = λq

The potential coupling V 10 only depends on q, thus the stretching mode is also called coupling mode. The parameters of the model are: the effective mass of the reaction coordinate m -1 = 4.84 • 10 -4 eV is chosen so as to achieve isomerization within 200 fs; the frequency of the coupling mode is ω = 0.19 eV and the interstate coupling is λ = 0.19 eV, with κ = 0.1 eV indicating a gradient in the excited state; the parameters W 0 = 3.6 eV, W 1 = 1.09 eV and E 1 = 2.48 eV are chosen in order to match the S 0 /S 1 gap to the center frequency of the absorption bands of the cis and trans isomers. At the cis geometry (ϕ cis = 0), the excitation energy ∆E ex to induce a S 0 -to-S 1 transition at q = 0 lies in the visible domain, being ∆E ex = 2.48 eV (λ ex = 500 nm). The model is parametrized to reproduce features observed in Raman and time-resolved experiments [17], and was for long time the preferred choice to study the ultrafast relaxation process of the photo-excited retinal based on quantum-dynamical simulations [18,[56][57][58][59][60][61][62][63][64][65][66][67].

The secondary modes representing the effect of the environment are included via a bath Hamiltonian [18,21]

ĤB = Nm j=1 1 2 ω j p 2 j + x 2 j 0 0 Nm j=1 1 2 ω j p 2 j + x 2 j + c j x j (22)
such that the full system-bath Hamiltonian is

Ĥ = ĤS + ĤB (23)
Note that, even though there is no direct coupling between the modes of the system ϕ, q and the modes of the bath x j , the bath affects the energies of the electronic states in different ways. Therefore, via the nonadiabatic coupling

-π -π/2 0 π/2 π -2 -1 0 1 2 0 1 2 3 4 E S 0 E S 1 ϕ q Energy (eV) -π -π/2 0 π/2 π -2 -1 0 1 2 0 1 2 3 4 ϕ q Energy (eV) 0 1 2 3 4 -π -π/2 0 π/2 π CI position (system+bath) CI position (system) Energy (eV) ϕ E S 0 (system) E S 1 (system) E S 0 (system+bath) E S 1 (system+bath)
FIG. 1: Upper panel: adiabatic BO PESs (S0 in blue and S1 in red) as functions of ϕ, q when the system is not coupled to the bath. The contour plot represent the initial Gaussian nuclear density in ϕ, q. Lower panel: adiabatic BO potential energy curves as functions of ϕ obtained for q = 0 when the system is not coupled to the bath (S0 in blue and S1 in red) and when the system is coupled to the modes 1, 8 and 23 fixing the bath positions at x1 = 0.5, x8 = 1.0, x23 = 1.5. Vertical dashed lines highlight the positions of the CIs in the two cases.

between electronic and nuclear motion, a certain degree of correlation is established over time within the full system.

In our simulations, we use up to N m = 23 modes with frequencies ω j (in cm -1 ) and dimensionless excited-state shifts c j /ω j (in parenthesis) from Ref. [18] 1 (lower panel) compares the adiabatic BO potential energies as functions of the variable ϕ for q = 0 considering only the system (in blue and red) and considering the system coupled to only 3 of the modes listed here (in orange and dark green), i.e., mode 1, 8 and 23, at an arbitrary choice of the bath positions

x 1 = 0.5, x 8 = 1.0, x 23 = 1.5. The figure highlights with the black and gray dashed lines the displacement of the position of the CI along the ϕ coordinate. The evolution of the Gaussian nuclear wavepacket -prepared in S 0 before photo-excitation -is initiated by an instantaneous S 0 -to-S 1 transition. Its expression is

χ S0 (ϕ, q, {x j }, t = 0) = 1 πσ ϕ σ q e -ϕ 2 2σ 2 ϕ + q 2 2σ 2 q Nm j=1 4 1 πσ 2 j e - x 2 j 2σ 2 j (24)
with σ ϕ = 0.128, σ q = 1.0, σ j = 1.0 ∀ j. The two-dimensional Gaussian density in ϕ, q space is represented as contour lines in Fig. 1 (upper panel).

V. COMPUTATIONAL DETAILS

A. Quantum dynamics propagation

The propagation is performed in the diabatic representation of the wavefunction with the Short Iterative Lanczos (SIL) scheme. [68,69] With this scheme, the evolution operator, Û (∆t), is decomposed into spectral projectors

Û (∆t) = k Exp - i E k ∆t |Φ k Φ k | ( 25
)
where the |Φ k and E k are, respectively, the eigenvectors and the eigenvalues of the Hamiltonian in the Krylov subspace. The first vector of this subspace is the wavepacket at time t and other vectors are added until convergence (the contribution of the last |Φ k is smaller than 10 -12 ). Furthermore, to avoid numerical instabilities, a vector orthogonalization is enforced via a Gram-Schmidt procedure at each iteration k.

For our application, the number vectors in the Krylov subspace remains more or less constant for all basis set sizes and from 2 to 25 degrees of freedom, but it depends on the time step. For 1 or 5 fs time step, the number of vectors are about 20 or 50, respectively. This is clearly an advantage over Chebyshev scheme [69,70] for which, the number of vectors depends on the full Hamiltonian spectral range which increases with the basis set size. Furthermore, for the SIL scheme and by construction, the wavepacket norm is conserved and the energy is virtually conserved. [69] The basis sets are defined as follows. (i) A real Fourier series for the first mode, associated to the coordinate ϕ in a [-π, π] range. The quadrature grid is defined by equidistant points and the weights are constant. (ii) A scaled Harmonic oscillator basis set, HO, for the other modes. The not-normalized expression of the j th basis function of the i th mode, q i , is b i j (q i ) = H j-1 (x(q i )) • Exp(-1/2 • x(q i ) 2 ) where x(q i ) = s HO i • (q i -q HO 0,i ). Each parameter s HO i is adapted to the curvature and the mass of the 1D harmonic oscillator model associated to the mode q i . For the second mode of the system, s HO 2 = 0.9592722 while for all the bath modes s HO i = 1.0. The parameters q HO 0,i are defined, usually, at the corresponding equilibrium positions of the potential. However, since they are different for the two diabatic potentials, q HO 0,i are defined in the middle of the two minima. From the expressions of the diabatic potentials (Eqs. ( 17) and ( 22)), q HO 0,2 = -κ/(2ω) for second mode of the system and q HO 0,i = -c i /(2ω i ) for the bath modes. The size of the basis sets, nb i λi , has been checked carefully with respect to the convergence of the diabatic populations. Therefore, for the direct-product calculations, nb 1 , nb 2 and the nb bath (for all bath modes) are, respectively, 256, 30 and 8. With those values, the diabatic and adiabatic populations up to 250 fs, the total energy and the other properties of interest are well converged. This means that, when the nb i values are increased, respectively, to 300, 40 and 12, the analysed properties remain unchanged (see Table I).

For the calculation with 4, 8 and 25 degrees of freedom, the Smolyak scheme with the system-bath separation is used. The relation between, nb i λi and λ i are: (i) nb 1 λ1 = 1 + 36λ 1 ; (ii) nb 2 λ2 = 1 + 4λ 2 ; (iii) nb i λi = 1 + 2λ i for the bath modes. The same relations are used for nq i λi . The convergence of the properties will be checked against L s , the system Smolyak parameter. L b , the bath parameter, is fixed to 2, the coupling parameter, L m , is defined as L s + 1 and the Smolyak grib parameter, L G is L m + 1. In Table II, we can see that the energy is already converged for L s = 6. However, L s = 7 is needed to converge the populations up to 200 fs with an accuracy of 0.1 %. Furthermore, it is worth noting that all properties remain almost unchanged when L b = 3 is used. For the 4D-model (n = 4), the comparison between the direct-product scheme and the Smolyak scheme shows an excellent agreement for the energy and the populations up to 200 fs. However, the computation with Smolyak scheme is about 30 time faster than the one with the direct-product approach. For a larger number of degrees of freedom, only calculations with the Smolayk scheme are possible in a reasonable time.

B. Quantum-classical propagation

Trajectory-based calculations are performed with the G-CTMQC code [36], where the algorithms presented in Section III B are implemented. In all cases, an ensemble of N traj = 1000 classical trajectories have been propagated starting from Wigner-sampled initial conditions. The Gaussian Wigner density in positions and momenta is determined as the Wigner-transform of the quantum nuclear density |χ S0 (ϕ, q, {x j }, t = 0)| 2 from Eq. ( 24). Electronic and nuclear evolution equations are integrated with the Runge-Kutta-Gill algorithm and the velocity-Verlet algorithm, respectively, with a time step of 0.0024 fs (0.1 a.t.u.). In Section VI, the results are presented for simulations of duration 300 fs (for the study of the numerical costs of independent-and coupled-trajectory procedures when increasing the system size) and 1 ps (for the study of the effect of the environment on the system dynamics).

C. Observables of interest

In Section VI, the dynamics of the system w/o the effect of the environment is studied by reporting on some key nuclear and electronic observables. We provide here a list of these observables, together with information on how to determine them based on quantum wavepacket calculations and on trajectory-based calculations.

The ϕ, q nuclear (probability) density of the system when N m = 0 is given as the sum of the S 0 and S 1 (probability) densities

|χ(ϕ, q, t)| 2 = |χ S0 (ϕ, q, t)| 2 + |χ S1 (ϕ, q, t)| 2 (26)
and becomes the reduced ϕ, q nuclear (probability) density in the case of the system-bath dynamics by integrating over all bath {x j } configurations

|χ(ϕ, q, t)| 2 = |χ(ϕ, q, {x j }, t)| 2 Nm j=1 dx j (27) = |χ S0 (ϕ, q, {x j }, t)| 2 + |χ S1 (ϕ, q, {x j }, t)| 2 Nm j=1 dx j (28)
The probability for the system to be in the electronic ground (S 0 ) or excited (S 1 ) state is determined as the integral of the corresponding nuclear (reduced) density,

ρ S l (t) = |χ S l (ϕ, q, t)| 2 dϕdq with l = 0, 1 (29)
which we refer to as "electronic population". This electronic population can be decomposed into two contributions as ρ S l (t) = P cis S l (t) + P trans S l (t), from the cis and trans conformers,

P cis S l (t) = |ϕ|< π 2 |χ S l (ϕ, q, t)| 2 dϕdq (30)
P trans S l (t) = |ϕ|> π 2 |χ S l (ϕ, q, t)| 2 dϕdq, (31)
yielding the probability of each conformer in each state.

Given the time-dependent vibronic state of the system, which can be expressed in vector notation as

|Ψ(ϕ, q, {x j }, t) = χ S0 (ϕ, q, {x j }, t) χ S1 (ϕ, q, {x j }, t) (32)
average values of any operator, say Ô, can be determined as

O(t) = Ô (t) = Ψ(ϕ, q, {x j }, t)| Ô |Ψ(ϕ, q, {x j }, t) dϕdq Nm j=1 dx j (33)
In Section VI, we report results for various choices of the operators, namely

Ô = P 2 ϕ /(2m), P 2 q /(2ω -1 ), q, φ, x1 , x2 , qx 1 , qx 2 , φx 1 , φx 2 .
The corresponding quantum-classical observables, to be compared to the quantum calculations, can be usually determined as averages over the trajectories distribution. Quantum-classical electronic populations, corresponding to Eq. ( 29), are computed using the TSH-based methods as

ρ qc S l (t) = N S l (t) N traj (34)
where N S l (t) is the number of trajectories on state S l at time t (see Ref. [71] for a detailed discussion on how to estimate electronic populations with CT-TSH). The cis and trans probabilities of Eqs. ( 30) and ( 31) as functions of time can be simply estimated by counting the number of trajectories with |ϕ α (t)| < π 2 (for cis) or |ϕ α (t)| > π 2 (for trans), irrespective of the values of q α (t) and of {x α j (t)}. Finally, since the observables of interest can be expressed as functions of classical phase space variables for each trajectory, i.e., O α (t) = O(P α ϕ (t), P α q (t), ϕ α (t), q α (t), {p α j (t)}, {x α j (t)}), the expression corresponding to Eq. ( 33) in the quantum-classical case is a simple average over the trajectories

O qc (t) = 1 N traj Ntraj α=1 O α (t) (35)

VI. NUMERICAL RESULTS

The scope of this section is to compare the performance of CT-TSH and TSH-ED in simulating the photoisomerization dynamics of the 2D retinal model, i.e., the system, embedded in an environment of harmonic bath modes (BMs) and to benchmark some results against quantum dynamics calculations. We report on the population dynamics in Section VI A, on the correlation between system and bath degrees of freedom in Section VI B and on the nuclear kinetic energy in Section VI C. We recall here that both trajectory-based methods include the effect of quantum decoherence within an independent-trajectory (TSH-ED) or a coupled-trajectory (CT-TSH) framework. Therefore, we find instructive to analyze their computational cost as well, when using the G-CTMQC code in its serial implementation (in Section VI D). Similar analysis has been the focus of previous work [START_REF] Pieroni | [END_REF]35], where CT-TSH and TSH-ED have been compared and benchmarked against some reference (quantum) results in one-dimensional and two-dimensional models. However, in this work we increase the complexity of the model up to 25 nuclear coordinates (2D + 23BM) and we simulate short-time dynamics as well as long-time dynamics in Section VI E (not accessible to quantum vibronic wavepackets techniques).

Quantum calculations are performed over 200 fs, since we already observed a sizeable effect on the electronic populations due to the environment on this time scale. Since this is not the case for the quantum-classical dynamics, we report CT-TSH and TSH-ED results up to 300 fs. Long-time dynamics is simulated up to 1 ps.

A. Population analysis

In Fig. 2 (left) we show the population of the first diabatic state, V 00 , as function of the propagation time with 0, 1, 2, 6 and 23 bath modes. For the 2D, 2D + 1BM and 2D + 2BM models, the standard direct-product approach has been used with a 1 fs time step. With this short time step, we can observe that the population presents small oscillations up to 150 fs with a period of about 1-2 fs. For 2D + 6BM and 2D + 23BM, the Smolyak scheme has been used with a 5 fs time step and the small oscillations observed with the 1 fs time step are not visible anymore although they are present in the dynamics. It is worth noting, that this large time step does not prevent us from performing an exact propagation. However, it enables us to reduce the number of Hamiltonian actions required to propagate by a factor two for full propagation time. As it appears clearly from the results just discussed, the time-dependent quantum treatment with a sparse grid scheme is feasible for the studied vibronic model with a large number of degrees of freedom (up to 25). The computational time necessary to obtain the results in Fig. 2 employing the present Smolyak FIG. 2: Left: Population of the first diabatic state, V00, using a quantum dynamics propagation approach for different model sizes, from the 2D system without the environment to the 2D + 23BM system+bath model which includes all modes listed in Section IV. Right: Population of the S1 adiabatic state using a quantum dynamics propagation approach for different model sizes, from the 2D to the 2D + 23BM (Ls = 6 has been used for 2D + 23BM model).

scheme is longer than using the multi-configuration time-dependent Hartree method (MCTDH) [START_REF] Meyer | Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF] or its multilayer generalization [START_REF] Wang | [END_REF][74][75][76], although we demonstrate that calculations are indeed doable. Specifically, for the 2D + 23BM model, more than 9 days are required to obtain the full propagation with a MPI parallelization scheme on a 32 core computer (see Table II). Note that with this Smolyak scheme and since the wavepacket is known on the full sparse grid, the adiabatic populations and other adiabatic properties can be computed easily.

Figure 2 (right) reports the S 1 adiadatic population for the 2D, 2D + 1BM and 2D + 2BM, 2D + 6BM and 2D + 23BM models and, as for the diabatic populations, the bath modes have a small effect on the populations. The fast oscillations observed up to 150 fs in the diabatic populations are absent in the adiabatic populations. Such small oscillations arise from the non-zero nonadiabatic coupling contributions around the initial diabatic wavepacket; the corresponding initial adiabatic wavepacket is slightly delocalized on the S 0 adiabatic state and its initial S 0 population is about 0.3 %. This delocalization creates population oscillations with a corresponding period, τ 1 (τ 1 = 1.67 f s), associated to the energy difference between the two diabatic states or to the parameter E 1 of the model (see Section IV).

Figure 3 (left) reports the time evolution of the population of the excited state S 1 after photo-excitation up to 300 fs of dynamics simulated with CT-TSH (upper panel) and with TSH-ED (lower panel) and computed according to Eq. [START_REF] Pieroni | [END_REF]. The colored lines represent different simulation setups, from the 2D model without inclusion of the environment (dark red thick lines) to the 2D + 23BM model that includes all modes listed in Section IV (dark magenta thick lines). The bath modes are included one by one following the order given in Section IV. In the same figure (right), we show the probability of forming the trans conformer, by counting the number of trajectories such that the absolute value of the angle ϕ is larger than π/2, irrespective of the occupied electronic state.

The analysis of the S 1 population up to 300 fs shows some effect of the environment, which tends to slightly decrease the probability that the system remains in, or actually re-populates, the excited state. Such an effect is stronger in TSH-ED than in CT-TSH, with an abrupt jump observed for TSH-ED between the 2D + 11BM and the 2D + 12BM simulations. However, this effect appears to be in disagreement with quantum dynamics results, that predict an increase of the final S 1 population as an effect of the environment. The stabilization effect due to the environment

2D 2D + 1BM 2D + 2BM 2D + 3BM 2D + 4BM 2D + 5BM 2D + 6BM 2D + 7BM 2D + 8BM 2D + 9BM 2D + 10BM 2D + 11BM 2D + 12BM 2D + 13BM 2D + 14BM 2D + 15BM 2D + 16BM 2D + 17BM 2D + 18BM 2D + 19BM 2D + 20BM 2D + 21BM 2D + 22BM 2D + 23BM

CT-TSH

Trans conformer

2D 2D + 1BM 2D + 2BM 2D + 3BM 2D + 4BM 2D + 5BM 2D + 6BM 2D + 7BM 2D + 8BM 2D + 9BM 2D + 10BM 2D + 11BM 2D + 12BM 2D + 13BM 2D + 14BM 2D + 15BM 2D + 16BM 2D + 17BM 2D + 18BM 2D + 19BM 2D + 20BM 2D + 21BM 2D + 22BM 2D + 23BM
FIG. 3: Left: population of the S1 state over 300 fs of dynamics using CT-TSH (upper panel) and TSH-ED (lower panel). Right: probability of forming the trans conformer over 300 fs of dynamics using CT-TSH (upper panel) and TSH-ED (lower panel). The colored lines represent different models, from the 2D system without the environment to the 2D + 23BM system+bath model which includes all modes listed in Section IV.

appears to be less severe when analyzing the probability of formation of the trans conformer in S 0 and in S 1 , being nearly absent at the end of the simulated dynamics for TSH-ED. Nonetheless, overall we conclude that even at short times the photo-product formation, i.e., the trans conformer in S 0 , is stabilized by the presence of the environment.

B. Correlation analysis

In Fig. 4, we analyze the effect of the coupling between the system and the bath in the 2D + 2BM simulation. From the expression of the Hamiltonian given in Eq. ( 23), it is clear that there is no direct coupling between the modes of the system ϕ, q and the modes of the bath {x j }. However, the bath affects the energies of the electronic states in different ways (see Eq. ( 22)), and thus, due to the nonadiabatic electron-nuclear coupling within the system, correlation is established along the dynamics between the nuclear degrees of freedom of the system and the bath degrees of freedom.

We demonstrate this observation in Fig. 4 by reporting measures of correlation between ϕ, q and x 1 , x 2 calculated in the quantum dynamics simulations (left) and in the quantum-classical simulations (right), distinguishing CT-TSH (upper panel) and TSH-ED (lower panel) results. For instance, for the degrees of freedom ϕ and x 1 , such a measure is estimated as Corr(ϕ, x 1 ) = ϕ • x 1 -ϕ • x 1 in the quantum case and as Corr(ϕ, x 1 ) = ϕ qc (t)x qc 1 (t) -(ϕx 1 ) qc (t) in the quantum-classical case. Each term with superscript qc is evaluated as an average over the trajectories as given in Eq. (35).

All quantities shown in Fig. 4 (right) are non-zero, even though quantum calculations show that the correlations between ϕ and all bath modes are zero. This is due to the fact that the initial wave packet is symmetric with respect to ϕ = 0 as well as the Hamiltonian, thus producing a wavepacket that remains symmetric during the propagation. We expect that in the limit of infinite number of classical trajectories to reproduce the distribution along ϕ = 0, the quantum-classical correlations involving ϕ would tend to zero as well.

We find interesting to analyze the correlation between q and the x i (i = 1, 2) in quantum results. Figure 4 (left)

shows the indirect and weak coupling between the system coordinate, q, and the bath coordinates, which presents oscillation with a period between 19 and 24 fs. This period is close to the intrinsic period of 21.8 fs associated to the vibrations along q of the system. Corr(q,x 1 ) Corr(q,x 2 )

FIG. 4: Left: Measure of correlation between the nuclear degrees of freedom of the system ϕ,q and the modes of the bath x1,x2 along the dynamics simulated for the 2D + 2BM model with the quantum dynamics propagation approach. The curves are: Corr(q, x1) (blue), Corr(q, x2) (purple). Right: Measure of correlation between the nuclear degrees of freedom of the system ϕ, q and the modes of the bath x1, x2 along the dynamics simulated for the 2D + 2BM model with CT-TSH (upper panel) and the TSH-ED (lower panel). In each plot, the curves are: Corr(ϕ, x1) = ϕ qc (t)x qc 1 (t) -(ϕx1) qc (t) (red), Corr(ϕ, x2) = ϕ qc (t)x qc 2 (t) -(ϕx2) qc (t) (green), Corr(q, x1) = q qc (t)x qc 1 (t) -(qx1) qc (t) (blue), Corr(q, x2) = q qc (t)x qc 2 (t) -(qx2) qc (t) (purple).

C. Kinetic energy analysis

Figure 5 shows the behavior of the nuclear kinetic energy along the reaction coordinate ϕ (left) and along the coupling mode q (right) as functions of time computed according to Eq. ( 35) and based on CT-TSH (upper panels) and on TSH-ED (lower panels).

The effect of the environment on the torsion of the system is to slightly reduce its kinetic energy from the 2D to the 2D + 23BM simulation, which is expected since the additional degrees of freedom, modelled here as a bath of harmonic oscillators, is that of stabilizing the photo-product in the trans configuration once it is formed. As it is observed in general in all results reported in Figs. 3 and5, including a large number of bath modes reduces the deviation between CT-TSH and TSH-ED (at least, up to 300 fs): when comparing the kinetic energy along the coupling mode for a small number of bath modes, CT-TSH and TSH-ED show a slightly different behavior, which is, instead, "washed out" by the presence of the environment. The coupled-trajectory approach, by contrast to the independent-trajectory approach, can be thought as preserving some degree of quantumness in its description of nuclear dynamics, since the trajectories, being coupled, behave as belonging to a single entity, i.e., the quantum nuclear density. This is the feature that, when accounted for in the electronic evolution equation, correctly accounts for quantum decoherence. When the number of degrees of freedom of the bath increases, the system is coupled to a large number of modes and encounters a large number of energy dissipation and decoherence channels, which make the system somehow more "classical" than when it is isolated.

Proceeding on a similar analysis, Fig. 6 reports the kinetic energy of the harmonic bath modes as function of time for the 2D + 23BM simulation with TSH-ED (red) and CT-TSH (blue). The two methods yield basically identical results for each mode, thus, we deduce that the differences observed in Fig. 5 for the 2D + 23BM simulation between

(E h ) 2D 2D + 1BM 2D + 2BM 2D + 3BM 2D + 4BM 2D + 5BM 2D + 6BM 2D + 7BM 2D + 8BM 2D + 9BM 2D + 10BM 2D + 11BM 2D + 12BM 2D + 13BM 2D + 14BM 2D + 15BM 2D + 16BM 2D + 17BM 2D + 18BM 2D + 19BM 2D + 20BM 2D + 21BM 2D + 22BM 2D + 23BM 0
(E h ) 2D 2D + 1BM 2D + 2BM 2D + 3BM 2D + 4BM 2D + 5BM 2D + 6BM 2D + 7BM 2D + 8BM 2D + 9BM 2D + 10BM 2D + 11BM 2D + 12BM 2D + 13BM 2D + 14BM 2D + 15BM 2D + 16BM 2D + 17BM 2D + 18BM 2D + 19BM 2D + 20BM 2D + 21BM 2D + 22BM 2D + 23BM
FIG. 5: Nuclear kinetic energy along the torsion coordinate ϕ (left) and along the coupling mode q (right) expressed in Hartree (E h ) as function of time using CT-TSH (upper panels) and TSH-ED (lower panels). The color code is the same as in Fig. 3.

TSH-ED and CT-TSH can be solely ascribed to a different treatment of the nonadiabatic coupling within the 2D (sub-)system.

D. Computational time analysis

In general, the analysis reported so far on short-time system+bath dynamics up to 300 fs has shown small differences between CT-TSH and TSH-ED, that even tend to disappear/reduce when increasing the number of degrees of freedom of the environment. Therefore, we find informative to investigate the dependence of the computational cost of the coupled-trajectory (CT-TSH) and of the independent-trajectory (TSH-ED) methods. We recall that all calculations are performed with the G-CTMQC code [36] interfaced with QuantumModelLib [START_REF] Lauvergnat | QuantumModelLib[END_REF], on a node comprising of two 16-core 2.30GHz Intel(R) Xeon(R) E5-2698 v3 CPUs with 192 GB of DDR4 2133 MHz ECC DRAM.

Figure 7 shows the dependence of the calculation time on the dimensionality of the model, from the 2D model (with 0BM) to the 2D + 23BM model. Aside from a slight deviation for the 2D + 12BM simulation, which is probably due to input/output at a maximum total shared bandwidth of 119 MB/s, it appears clear that the computational cost of CT-TSH grows linearly with the number of (nuclear/classical) degrees of freedom, while TSH-ED is basically independent. Therefore, the limiting step in CT-TSH appears to be the calculation of the quantum momentum, which is not performed in TSH-ED.

E. Long-time dynamics

Figure 8 reports some results on the 1 ps dynamics simulated based on CT-TSH (upper panels) and on TSH-ED (lower panels) for the isolated system in 2D (green lines) and for the same system coupled to the bath of 23 harmonic oscillators (purple lines). In the figure, the dynamics is analyzed by showing the behavior of the excited-state S 1 population (left panels) and the formation probability of the trans conformer (right panels) as functions of time. In addition, we are able to compare the trans probability with the corresponding results obtained in Ref. [21] with the multi-configuration time-dependent Hartree method (MCTDH) [START_REF] Meyer | Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF] and its multilayer generalization [START_REF] Wang | [END_REF][74][75][76].

When including the environment, the system tends to approach the ground state S 0 slightly faster than when it is isolated, as the decrease in the S 1 population attests for both CT-TSH and TSH-ED. Furthermore, as observed previously, the inclusion of the environment reduces the deviation between the two methods even at long times. Conversely, both trajectory-based methods fail to correctly capture the decrease in the probability of forming the trans conformer at long times, that is expected when comparing reference MCTDH results from Ref. [21] for the 2D and for the 2D + 23BM models. However, CT-TSH reproduces slightly better than TSH-ED the damping of the oscillations of the trans probability over time when the environment is included in the simulation.

Finally, in Fig. 9, we compare the long-time dynamics up to 1 ps in the isolated 2D model and in the system+bath model by analyzing the nuclear kinetic energy as function of time along the reactive coordinate ϕ (left panels) and along the coupling mode q (right panels) using CT-TSH (upper panels) and TSH-ED (lower panels). While the inclusion of the environment does not strongly affect the behavior of the system along ϕ, we observe a striking damping of the nuclear kinetic energy along the high-frequency coupling mode q due to the bath of harmonic oscillators. The effect is severe when the dynamics is simulated with CT-TSH, but it is not negligible even with TSH-ED, and shows once again that the presence of the environment tends to reduce the differences in dynamics simulated with the two trajectory-based methods.

VII. CONCLUSIONS

In this work, we presented an in-depth analysis of the performance of various algorithms in simulating the photoisomerization process in the retinal of rhodopsin using a system-bath model of up to 25 degrees of freedom and two electronic states.

We reported numerical results based on quantum vibronic wavepacket dynamics, which required the use of special sparse methods to handle the computational cost required for the full system (25 degrees of freedom) and for 200 fs of dynamics. Quantum-classical calculations are, on the other hand, more affordable from the computational point of view, even though the coupled-trajectory Tully surface hopping scheme derived from the exact factorization, i.e., CT-TSH, is more expensive than standard fewest-switches surface hopping, i.e., TSH-ED, which uses independent trajectories. While some differences have been observed between quantum and quantum-classical results, that require additional investigation -as for instance the effect of the environment on the electronic excited-state population, it seems that the agreement between CT-TSH and TSH-ED improves as the number of environment modes increases, suggesting that the environment washes out some quantum-mechanical features that might be present in the coupled-trajectory scheme. FIG. 9: Nuclear kinetic energy along the torsion coordinate ϕ (left) and along the coupling mode q (right) expressed in Hartree (E h ) as function of time using CT-TSH (upper panels) and TSH-ED (lower panels) over 1 ps of dynamics. The color code is the same as in Fig. 8.

FIG. 6 :

 6 FIG.6: Kinetic energy along the coordinates {xj} as function of time expressed in Hartree (E h ) for the simulation including 2D + 23BM. TSH-ED results are reported in red and CT-TSH results are reported in blue. For visualization purposes, the curves are translated along the y-axis with increasing value of the index j, such that the bottom red and blue curves represent the kinetic energy of the mode x1 while the top red and blue curves represent the kinetic energy of the mode x23.

FIG. 7 :

 7 FIG.7: Dependence of the calculation time (in hours) as function of the dimensionality of the model for the coupled-trajectory (CT-TSH, circles) and for the independent-trajectory (TSH-ED, squares) methods. The x-axis is labeled in terms of the number of bath modes in each simulated dynamics.

TABLE I :

 I Size of the basis set, nbi, energy in eV (E), population of the first diabatic state (PopV 00 ) at 200 fs and computational time as functions of the number of degrees of freedom, n, for the direct-product scheme.

	n nb1 nb2 nb bath E (eV) PopV 00 Computational time
	2 256 30	/	0.09483 0.3992	6min
		300 40	/	0.09483 0.3992	11min
	3 256 30	8	2.62963 0.4034	1h 18min
		300 40	12 2.62963 0.4034	2h 41min
	4 256 30	8	2.68187 0.4095	18h 38min
		300 40	12 2.68187 0.4096	54h 28min
	n Ls	N b		Nq	E (eV) PopV 00 Computational time
	4 6 25 133	321 048 2.68187 0.3385	39min
	7 35 733	560 040 2.68187 0.4099	1h 6min
	8 48 205	912 210 2.68187 0.4096	1h 44min
	8 6 149 365 2 272 032 2.91089 0.3591	5h 45min
	7 215 837 3 963 360 2.91089 0.3863	12h 15min
	8 294 549 6 455 640 2.91089 0.3879	21h 1min
	25 4 633 249 7 324 450 4.32296 0.2229	1d 21h
	5 1 161 577 15 906 982 4.32296 0.2814	4d 15h
	7 2 695 593 53 031 480 4.32296 0.4184	9d 8h a

TABLE II :

 II Size of the basis set (N b ) and grid (Nq), energy (E), population of the first diabatic state (PopV 00 ) at 200 fs and computational time as functions of the number of degrees of freedom, n and the Smolyak parameter, Ls.

a MPI paralelization scheme with 32 cores.

Acknowledgments

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

CT-TSH

Trans conformer 2D 2D + 23BM 2D from Ref. [44] 2D + 23BM from Ref. [44] FIG. 8: Left: population of the S1 state over 1 ps of dynamics using CT-TSH (upper panel) and TSH-ED (lower panel). Right: probability of forming the trans conformer over 1 ps of dynamics using CT-TSH (upper panel) and TSH-ED (lower panel), compared to some reference MCTDH results extracted from Ref. [21] for the 2D model (continuous thin black lines) and for the 2D + 23BM model (dashed thin black lines). The colored lines represent different models, i.e., the 2D system without the environment (in green) and the 2D + 23BM system+bath model (in purple).