Benjamin Gonon

Benjamin Lasorne

Gabriel Karras

Loïc Joubert-Doriol

David Lauvergnat

Frank Billard

Bruno Lavorel

Olivier Faucher

St Éphane Gu Érin

Edouard Hertz

Fabien Gatti

A generalized vibronic-coupling Hamiltonian for molecules without symmetry: application to the photoisomerization of benzopyran

We present a model for the lowest two potential energy surfaces (PESs) that describe the photoinduced ring-opening reaction of benzopyran taken as a model compound to study the photochromic ring-opening reaction of indolinobenzospiropyran and its evolution toward its open-chain analogue. The PESs are expressed in terms of three effective rectilinear coordinates. One corresponds to the direction between the equilibrium geometry in the electronic ground state, referred to as the Franck-Condon (FC) geometry, and the minimum of conical intersection (CI); while the other two span the two-dimensional branching space at the CI. The model reproduces correctly the topography of the PESs. The ab initio calculations are performed with the extended multi-configuration quasi-degenerate perturbation theory at second order (XMCQDPT2) method. We demonstrate that accounting for electron dynamic correlation drastically changes the global energy landscape since some zwitterionic states become strongly stabilized. Quantum dynamics calculations using this PES model produce an absorption spectrum that matches the experimental one to a good accuracy.

I. INTRODUCTION

Fast and reproducible photochromism provides the potential for many optical switching processes from ophthalmic glasses to optical memories, and non-linear optical devices [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF][2][3] . In this paper, we focus on a particular class of photochromic molecules: spiropyrans, where a ring-opening reaction is induced by a π → π * excitation, as shown on Fig. 1 (a), from the electronic ground state S 0 to the first singlet excited state S 1 . Irradiating a spiropyran molecule with Ultra-Violet (UV) light leads to the breaking of the bond between the spiro carbon atom and the oxygen atom of the chromene moiety. After this first step, a cis-trans isomerization takes place around the three central C-C bonds leading to a planar structure called merocyanine. This planar form allows the π-system to be delocalized over the whole structure, resulting in an absorption band in the visible range for the merocyanine isomer 4,5 : see Fig. 1 (a). The two isomeric forms thus have quite different structures, and the nature of the mechanism that allows spiropyran/merocyanine systems to switch from one isomeric form to another is far from being fully understood.

Over the past decades spiropyrans have been the focus of many theoretical and experimental studies [START_REF] Celani | Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds[END_REF][7][START_REF] Broo | Theoretical characterization of a multifunctional electrooptical molecular device: Photochemical ring-opening mechanism of indolinospirobenzopyran[END_REF][9][10][11][12][13][START_REF] Rini | Ultrafast uv-midir investigation of the ring opening reaction of a photochromic spiropyran[END_REF] . In particular, time-resolved spectroscopic experiments 9,10 have provided insights into the dynamics of the photochromic process in benzospiropyran compounds. For instance, these have been shown to be sensitive to the initial chirp and polarisation state of the femtosecond pulses initiating the reaction 11,12 . We assume as in Ref. [START_REF] Celani | Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds[END_REF] that benzopyran (see Fig. 1 (c)) can be used a) E-mail: edouard.hertz@u-bourgogne.fr b) E-mail: fabien.gatti@u-psud.fr The π-system of the chromene moiety is orthogonal to the indoline moiety, thus explaining that the absorption spectrum corresponding to the π → π * excitation is, to a first approximation, the sum of the spectra of the two orthogonal fragments. (c) Ring opening of the benzopyran molecule. We assume as in Ref. [START_REF] Celani | Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds[END_REF] that benzopyran can be used as a model compound to study the photochromic ring opening of indolinobenzospiropyran.

as a model compound to study the photochromic ring-opening reaction of indolinobenzospiropyran and its evolution toward its open-chain analogue (see Fig. 1). This approximation is valid since only the benzopyran moiety initially absorbs before the isomerization 9,10 , which originates from the quasi orthogonality between the benzopyran and the indoline moiety (see Fig. 1 (b)) 4,15 . Note that calculations on the full system, exposed in Ref. [START_REF] Liu | Multiple pathways for the primary step of the spiropyran photochromic reaction: A caspt2//casscf study[END_REF][START_REF] Liui | Multireference ab initio density matrix renormalization group (dmrg)-casscf and drmg-caspt2 study on the photochromic ring opening of spiropyran[END_REF] , indicate the possibility of an alternative deactivation pathway involving C-N dissociation on the indoline moiety. Our model system is specifically designed to address deactivation via C-O dissociation on the benzopyran chromophore. The extent to which both pathways may compete will not be addressed in the present work.

Describing the corresponding photochemical process is a very difficult task. Indeed, simulating internal conversion between vibronically coupled electronic states requires a nonadiabatic quantum dynamics treatment (see, e.g. Ref. [START_REF] Lasorne | Excited-state dynamics[END_REF] and references therein). The extreme case occurs at conical intersections vicinities [START_REF]Conical Intersections: Electronic Structure, Dynamics & Spectroscopy[END_REF][START_REF]Conical Intersections, Theory, Computation and Experiment[END_REF] , namely molecular geometries where two potential energy surfaces cross (the corresponding states are degenerate). In the vicinity of a conical intersection, the Born-Oppenheimer approximation breaks down, allowing non-adiabatic processes to take place. Conical intersections were predicted by John von Neumann and Eugene Wigner in 1929 and seem to be ubiquitous nowadays in polyatomic molecules [START_REF]Conical Intersections, Theory, Computation and Experiment[END_REF][START_REF]Conical Intersections, Electronic Strucutre, Dynamics and Spectroscopy[END_REF] . To describe the corresponding dynamics, two typical strategies can be used: direct dynamics, where the potential energies and vibronic couplings are calculated on the fly along trajectories [START_REF] Worth | Solving the time-dependent Schrödinger equation in one step: Direct dynamics of non-adiabatic systems[END_REF] , and dynamics on grids, where an analytical model for the potential energies and vibronic couplings is a prerequisite [START_REF] Worth | Using the MCTDH wavepacket propagation method to describe multimode nonadiabatic dynamics[END_REF] . Trajectory-based methods -either semiclassical, e.g., trajectory surface hopping (TSH) [START_REF] Tully | Mixed quantum-classical dynamics[END_REF] , or quantum, e.g., ab initio multiple spawning (AIMS) [START_REF] Ben-Nun | Ab initio quantum molecular dynamics[END_REF] or direct dynamics variational multiconfiguration Gaussian (DD-vMCG) [START_REF] Worth | Solving the time-dependent Schrödinger equation in one step: Direct dynamics of non-adiabatic systems[END_REF][START_REF] Richings | Quantum dynamics simulations using Gaussian wavepackets: the vMCG method[END_REF] -have many advantages but are bound to lower accuracy than grid-based quantum methods such as multiconfiguration time-dependent Hartree (MCTDH) [START_REF] Meyer | Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF] , which will be adopted here.

For grid-based methods, one has to solve the timedependent Schrödinger equation for the nuclei after having solved the Schrödinger equation for the electronic Hamiltonian that provides the potential energy surfaces (PES)s and the couplings between the different electronic states. Obtaining such PESs and the couplings expressed as multidimensional functions of the nuclear coordinates is a significant bottleneck for the quantum-mechanical treatment of the vibration and/or reaction. In addition, quantum chemistry calculations produce adiabatic energies and non-adiabatic coupling vectors. However, quantum dynamics simulations are more easily run using a so-called quasi-diabatic representation, which induces less numerical instabilities in quantum dynamics simulations [START_REF] Mead | Conditions for the definition of a strictly diabatic electronic basis for molecular systems[END_REF][START_REF] Pacher | Adiabatic and quasidiabatic states in a gauge theoretical framework[END_REF][START_REF] Yarkony | On the adiabatic to diabatic states transformation near intersections of conical intersections[END_REF] . Instead of vector functions of the nuclear coordinates that are singular at conical intersections, quasi-diabatic representations lead to a Hamiltonian matrix where the off-diagonal elements are smooth-varying scalar functions. The quasi-diabatic potential energy surfaces (the diagonal elements of the Hamiltonian matrix) also are simple functions, as opposed to the two-dimensional cusp at the conical intersection between two adiabatic surfaces. Quasidiabatic functions are thus much easier to fit to regular mathematical expressions.

Fortunately, dynamical processes in complex molecular systems often occur around a reactive core made of relatively few active degrees of freedom that are responsible of most of the dynamics. It is thus reasonable to invoke simplifying approximations for the inactive degrees of freedom when calculating the PES and simulating quantum dynamics [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF] . A widely used approach is the rigid-constraint one, which consists in freezing some nuclear degrees of freedom. Additional approximations may be considered. Indeed, since the quasidiabatic states depend smoothly on the nuclear coordinates, so do the quasi-diabatic Hamiltonian matrix elements that can thus be described by simple functions. For the treatment of processes with conical intersections, a very fruitful approach, termed vibronic-coupling Hamiltonian (VCH) model, has been formulated by Köppel and coworkers [START_REF] Cederbaum | An analogue to Koopman's theorem for the problem of vibrational structure in electron detachment spectra[END_REF][START_REF] Domcke | A many-body approach to the vibrational structure in molecular electronic spectra. ii application to nitrogen, carbon monoxide and formaldehyde[END_REF][START_REF] Cederbaum | Strong vibronic coupling effects in ionization spectra: The "mystery band" of butatriene[END_REF][START_REF] Köppel | Multimode molecular dynamics beyond the Born-Oppenheimer approximation[END_REF] . In this, the diabatisation is energy-based; it is done by ansatz: the PES (diagonal terms of the quasi-diabatic Hamiltonian matrix) and the potential-like interstate couplings (off-diagonal terms) are given as low-order Taylor expansions around the Franck-Condon geometry in terms of nuclear coordinates. The eigenvalues of the quasi-diabatic Hamiltonian matrix are then fitted to the adiabatic energies calculated with an electronic structure method over a range of geometries. This method is very efficient if the conical intersection (CI) is not far from the Franck-Condon (FC) geometry. In other cases, when motions of large amplitude are involved, Köppel and coworkers have devised models, where the diabatisation is also energy-based, exploiting the concept of regularized diabatic states [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF][START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] . In the latter models, symmetry plays a key role, in particular to define the "coupling" mode in the quasi-diabatic picture.

In our case, for the benzopyran molecule, the situation is particularly difficult since not only the CI lies far from the FC geometry but the molecule has no specific symmetry either. Our ultimate goal here is to explain and simulate experiments where the ring-opening process is controlled by laser pulses so that the system either goes selectively to the open form (photoreactivity) or back to the closed form (photostability). The present paper is rather methodology-oriented: we present a model inspired by the approach of Köppel and coworkers [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF][START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] but extended to a case where symmetry is not present. In addition, we highlight the dramatic impact of accounting for electron dynamic correlation on the topography of the PESs and on the physical process with respect to the original complete active space self-consistent field (CASSCF) calculations used in a recent work on the same system [START_REF] Joubert-Doriol | A generalized vibronic-coupling Hamiltonian model for benzopyran[END_REF] . The outline of this paper is as follows. Section II is devoted to the numerical details for the quantum chemistry calculations. Section III presents the ideas underlying the proposed approach. Section IV deals with the results, and the paper concludes presenting outlooks for the future.

II. COMPUTATIONAL DETAILS

Preliminary investigations have been performed with the complete active space self-consistent field (CASSCF) method. Following strategies used in previous works [START_REF] Celani | Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds[END_REF][START_REF] Joubert-Doriol | A generalized vibronic-coupling Hamiltonian model for benzopyran[END_REF] , our CAS, denoted (12,11) from now onward, was built with 12 electrons in 11 orbitals. We used state-averaged (SA) molecular orbitals with respect to the two electronic states S 0 and S 1 . We assume that the processes of absorption and excitation of the system mainly depend on the delocalized π system on the benzene and pyran parts of the molecule : the corresponding π/π * orbitals must thus be included. In order to describe the breaking of the CO bond in the pyran part, we also add the two σ/σ * orbitals of this bond. For the description of the ring opening, the non-bonding out-of-plane orbital on the oxygen also has to be included [START_REF] Joubert-Doriol | A generalized vibronic-coupling Hamiltonian model for benzopyran[END_REF] . The description of the active space is provided in section I of the supporting material.

The S 0 and S 1 stationary points were optimized at the SA-2-CASSCF/6-31G * level of theory using the GAUSSIAN 03 [START_REF] Frisch | Gaussian 03, Revision C.02[END_REF] and the MOLCAS [START_REF] Aquilante | Molcas 7: the next generation[END_REF] programs. The energy scans are very similar to what is shown on Fig. 9 of Ref. 40: the verticaltransition energy between the two electronic states is 4.65 eV (267 nm). On the excited electronic state, there is a transition state at 4.69 eV before a conical intersection (CI) point. The energy of the minimum of the conical intersection seam is 2.64 eV. The minimum of the cis open form of S 0 lies at 0.91 eV. The stationary points are thus M 0 (also denoted FC, for Franck-Condon), the minimum of S 0 for the closed form, M 1 , the minimum of S 1 , TS, the transition state on S 1 , CI, the minimum of the conical intersection seam, and cis, the minimum of S 0 for the open form. In the following, we will present several interpolated energy scans to compare the values obtained at the various levels of theory. The x-axis, corresponding to the scan coordinate, is in arbitrary unit: x = 0 corresponds to M 0 , x = 5 to M 1 , x = 10 to TS, x= 15 to CI, and x = 20 to cis. Intermediate points correspond to linear interpolations in internal coordinates between these reference geometries (note that from Fig. 2 to 5, all such five points are identical and were optimized at the SA-2-CASSCF/6-31G level of theory).

Several basis sets implemented in GAUSSIAN 03 and MOL-CAS have been compared: see Fig. 2. By visual inspection of the figure, we can conclude that the scans obtained with the 6-31G * and the valence triple-zeta polarization (TZVP) basis sets are indistinguishable. As expected, the small atomic natural orbital with valence triple zeta accuracy and with polarization functions, ANO-S-VTZP, basis set of MOLCAS gives different energies, but the positions and the nature of the stationary geometries remain very close. This is the reason why we have used the geometries optimized with the 6-31G * basis set for all the scans.

The ring-opening mechanism is described in detail in Ref. 40; the reactivity can be split into two steps: first, a π → π * excitation, mainly localized on the benzene moiety, exciting the system to the first electonic state and, second, an energy transfer to a π → σ * state, leading to the bond breaking. The change of dominant configuration within S 1 explains the existence of a transition state due to the presence of an avoided crossing between two diabatic states. The CASSCF method essentially accounts for static correlation. In Ref. 40, it helped us to rationalise the electronic states in terms of configuration interactions when building a quasi-diabatic PES model, since we assumed that such a picture would provide the chemical meaning of the quasi-diabatic states, and that further including dynamic correlation would be a mere refinement. On the contrary, we now show that the introduction of electron dynamic correlation leads to dramatic effects, changing the mechanism completely. In Ref. 40 When switching to complete active space with secondorder perturbation theory (CASPT2), the stabilization of zwitterionic states drastically changes the global picture: Fig. 3 (a) depicts the effect of CASPT2 calculations (obtained with the MOLCAS program). First, the degeneracy at the conical intersection (using the CASSCF geometries) is lifted significantly (about 1 eV) indicating that the geometries corresponding to the seam of the CI are affected by the perturbative treatment. In addition, using the CASSCF notation, we clearly see that, due to the perturbation, S 3 now lies between S 1 and S 2 . The difference between S 1 and S 3 is only around 0.2 eV. This is due to the zwitterionic character of S 3 40 : S 3 is lowered and interacts strongly with the ππ * -type electronic state. Not only are the energies strongly affected, but this also makes necessary to adopt a Multi-State (MS) procedure [START_REF] Finley | The multistate caspt2 method[END_REF] involving at least four electronic states. This means that electron dynamic correlation cannot be treated as a mere perturbation.

Fig. 3 (b) presents the results of our calculations with a multi-state (MS) approach over the first four electronic states: the scans with SA-2-CASSCF (12,11)/ANO-S-VTZP and SA-2-CASPT2 (12,11)/ANO-S-VTZP are depicted in black and red, respectively. They serve as a reference for the other cal- (12,11)/ANO-S-VTZP (blue). Note that SA-2-CASPT2 and MS-2-SA-2-CASPT2 are indistinguishable from M0 to CI. Note also that the new angular feature around x = 11 for MS-2-CASPT2 energies is due to a discontinuity of the second reference wavefunction from one to the other side of the avoided crossing marked by the S1 TS for CASSCF energies. This effect disappears when increasing the size of the reference space (MS-4-CASPT2 energies).

culations. The scan with MS-2-SA-2-CASPT2(12,11)/ANO-S-VTZP is shown in green: the difference between this curve and the curve with SA-2-CASPT2 in red is hardly visible. An important change with respect to CASSCF calculations is the fact that, with MS-4-SA-4-CASPT2, the S 1 and S 3 states interact and mix, leading to a strong lowering of the first excited state. Two qualitative features are interesting: first, the degeneracy at the conical intersection is still strongly lifted (we remind here that, for the CI, we use the geometries optimized with CASSCF), confirming that the position of the CI is affected by the introduction of dynamic correlation. Second, the TS disappears and is replaced by a mere shoulder. This should have a impact on the quantum dynamics of the system. A comprehensive study of the influence of the higher states has also been achieved: up to MS-8 calculations have been performed, indicating that the energies of the first two states are sensitive to S 6 for the cis geometry. The results are not shown here since we will return to this problem with extended multiconfiguration quasi-degenerate perturbation theory at second order (XMCQDPT2) calculations.

Finally, we have achieved calculations with the XMC-QDPT2 method [START_REF] Granovsky | Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory[END_REF] using the FIREFLY8 program [START_REF] Granovsky | Firefly version 8[END_REF][START_REF] Schmidt | [END_REF] . This approach has been proven to provide a more efficient and robust description than MS-CASPT2 approaches for getting correct energies near a conical intersection [START_REF] Ioffe | Photoisomerization of Stilbene: The Detailed XMCQDPT2 Treatment[END_REF][START_REF] Quick | Photoisomerization Dynamics of Stiff-Stilbene in Solution[END_REF][START_REF] Gozem | Dynamic Electron Correlation Effects on the Ground State Potential Energy Surface of a Retinal Chromophore Model[END_REF][START_REF] Gozem | Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection[END_REF] . The convergence of XMCQDPT2 calculations with respect to the number of electronic states in the multi-state treatment is presented in Fig. 4 (a): MS-4-SA-2-XMCQDPT2 (12,11)/TVZP (in green) and MS-8-SA-2-XMCQDPT2 (12,11)/TVZP (in blue) curves for the first two adiabatic states are almost indistinguishable, proving that MS-4 is sufficient. Interestingly enough, the convergence is faster than with MS-CASPT2 methods: as aforementioned, MS-CASPT2 energies were sensitive to S 6 for the cis geometry, which is no longer the case here. In Fig. 4 (b), we present scans where the molecular orbitals have been state-averaged over more than two states: the MS-4-SA-2-XMCQDPT2(12,11)/TVZP curve from before (in red) along with the MS-4-SA-4-XMCQDPT2 (12,11)/TVZP (in green), and the MS-5-SA-5-XMCQDPT2/TVZP (in blue) curves. Some points are missing at the MS-5-SA-5-XMCQDPT2/TVZP level of theory due to some problems of convergence (the problem was even more severe with MS-5-SA-5-CASPT2). We observe a perfect agreement for the first two adiabatic states, indicating that the calculations have converged. We notice that the transition state has disappeared at all levels of theory.

For completeness, we compare XMCQDPT2 and CASPT2 calculations in Fig. 5. The geometries of reference are slightly different due to the different atomic orbitals used in the two programs (MOLCAS and FIREFLY8). Before the CI, the curves are almost identical. Around the CI, the situation is very different. This difference can be attributed to the better convergence and behavior of the XMCQDPT2 method to describe regions of strong non-adiabatic couplings.

We thus have chosen the MS-4-SA-4-XMCQDPT2 (12,11)/TVZP level of theory for our model, to be used in quantum dynamics simulations. The stationary points were optimized at this level of theory using FIREFLY8. We never found any minimum in S 1 before rapidly reaching the CI seam. The stationary points are thus M 0 , the minimum of S 0 for the closed form, CI, the minimum of the seam of the conical intersection, and cis, the minimum of S 0 for the open form. The corresponding energies are given in Tab. I.

The experimental absorption spectrum 15 is the one of 2,2diethylchromene (where the CH 2 group is replaced by C-(CH 2 -CH 3 ) 2 ) at 77 K in 3-methylpentane. Our calculation at the MS-4-SA-4-XMCQDPT2 (12,11)/TVZP level of theory gives a vertical transition energy at 4.02 eV (308 nm), very close to the absorption maximum of the experimental spectrum at 310 nm.

We notice that the geometries of the stationary points are a little bit more "twisted" (i.e. less planar) for XMCQDPT2 cal- (12,11)/TVZP level of theory. M0 denotes the minimum of S0 for the closed form, CI the minimum of the seam of the conical intersection, and cis the the minimum of S0 for the open form. We also give relative energies (in eV) at the SA-2-CASSCF(12,11/6-31G * level of theory (the corresponding absolute value for M0 (S0) is -420.43712 hartree). The difference of energy between S0 and S1 at CI is 7. 10 -5 hartree = 0.002 eV.

(a) coworkers [START_REF] Liu | Multiple pathways for the primary step of the spiropyran photochromic reaction: A caspt2//casscf study[END_REF][START_REF] Liui | Multireference ab initio density matrix renormalization group (dmrg)-casscf and drmg-caspt2 study on the photochromic ring opening of spiropyran[END_REF] have performed quantum chemistry simulations at the CASSCF//CASPT2 levels on the full spiropyran system and not on the benzopyran only. These static calculations have been followed by semiclassical simulations at a similar level of theory for the electronic part [START_REF] Zhai | Detailed molecular dynamics of the photochromic reaction of spiropyran: A semiclassical dynamics study[END_REF] . They claim that their CASPT2//CASSCF study does not support the widely accepted conical-intersection mechanism (adopted in the present work) since the CASSCF-optimized CI on the C-O bond-cleavage path corresponds to an avoided crossing with an energy gap (> 20 kcal/mol) at the CASPT2 level. They propose another internal conversion funnel. These calculations do not include the zwitterionic states in the active space and the conical intersection is optimized at the CASSCF level only. In view of the present study, it is clear that such a descritpion cannot be reliable for our system. Accounting for the zwitterionic states is mandatory for a correct description of the mechanism. In addition, our CASPT2 calculations at the CASSCF-optimized geometries also give a significant energy gap (about 1 eV and thus > 20 kcal/mol): see the red curves on Fig. 3 (a) and (b). This does not prove that there is no conical intersection but simply that is is displaced somewhere in the vicinity of the latter and that its geometry must be obtained by a re-optimization at a post-CASSCF level, such as performed in the present work. To conclude, we cannot completely rule out the competitive pathway proposed in Refs. [START_REF] Liu | Multiple pathways for the primary step of the spiropyran photochromic reaction: A caspt2//casscf study[END_REF][START_REF] Liui | Multireference ab initio density matrix renormalization group (dmrg)-casscf and drmg-caspt2 study on the photochromic ring opening of spiropyran[END_REF] , but our calculations prove that the reaction leading to the ground state can be understood only at a post-CASSCF level including the lowest zwitterionic states in the reference space.

III. CONCEPTUAL DEVELOPMENTS

A. Choice of the active coordinates

In Ref. 40, a first model for the photodynamics of benzopyran based on CASSCF calculations was proposed. To this end, we chose very carefully a set of twelve active polyspherical coordinates [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF][START_REF] Gatti | Vector parametrization of the n-atom problem in quantum mechanics. I. Jacobi vectors[END_REF] and described the two adiabatic states of interest in terms of three underlying quasi-diabatic states that vary smoothly along this set of "active" curvilinear coordinates of the polyspherical kind. The polyspherical coordinates are curvilinear, i.e. coordinates involving angles, which, in principle, are better adapted to describe motions of large amplitude. The method of diabatisation was energy-based [START_REF] Köppel | Diabatic representation: methods for the construction of diabatic states[END_REF] . In other words, it was done by ansatz : we optimized the parameters defining the quasi-diabatic Hamiltonian matrix elements expressed as functions of the polyspherical coordinates by fitting the eigenvalues of this matrix to energies calculated on a grid of points.

This approach was devised as a generalization of the VCH model proposed by Köppel and coworkers [START_REF] Köppel | Multimode molecular dynamics beyond the Born-Oppenheimer approximation[END_REF][START_REF] Köppel | Diabatic representation: methods for the construction of diabatic states[END_REF][START_REF] Cederbaum | Born-oppenheimer approximation and beyond[END_REF] . In contrast to the original formulation, the number of diabatic states in our case was larger than the number of adiabatic ones (three and two, respectively).

Since the quasi-diabatic states depend smoothly on the nuclear coordinates, so do the quasi-diabatic Hamiltonian matrix elements that can thus be described by simple functions (typically linear or quadratic). The ability to fit the adiabatic data by diagonalisation of the quasi-diabatic Hamiltonian matrix and adjustment of its parameters ensures the quasi-diabatic character of the states involved in the model. In this approach, each quasi-diabatic state is associated with the dominant configuration around its own equilibrium geometry. The key point was the choice of additional quasi-diabatic states to be introduced in the model, which is determined from an analysis of static electron correlation. In other words, this was achieved at the CASSCF stage of the study.

Much care must be taken to understand the chemical meaning of the quasi-diabatic states in terms of dominant electronic configurations that can become quasi-degenerate at some point along the reaction path. For our model, we found out that the level of theory had to be further refined so as to include electron dynamic correlation. It is only the ability to fit at this higher level of theory that proves that the model is ultimately correct. This idea had already been tested successfully on a two-dimensional model of ethylene [START_REF] Lasorne | The electronic excited states of ethylene with large-amplitude deformations: A dynamical symmetry group investigation[END_REF][START_REF] Somoza | A generalised 17-state vibronic-coupling Hamiltonian model for ethylene[END_REF] : we used a model including 17 quasi-adiabatic states to fit five adiabatic states calculated with multireference configuration interaction (MRCI). Note that, in the model developed for ethylene [START_REF] Lasorne | The electronic excited states of ethylene with large-amplitude deformations: A dynamical symmetry group investigation[END_REF][START_REF] Somoza | A generalised 17-state vibronic-coupling Hamiltonian model for ethylene[END_REF] , one could assign to each quasi-diabatic state a single electronic configuration. For benzopyrane, each quasi-diabatic state corresponds to a family of many electronic configurations with a similar chemical nature. The difference comes from the fact that for ethylene we could consider two active electrons only in the reference space.

As already explained, three quasi-diabatic were selected : (ref ) bp corresponding to the closed-form molecule denoted with index bp (for benzopyran), (ππ * ), and (ref ) mr corresponding to the open-form molecule denoted with index mr (for merocyanine). The third state, (ππ * ), was added because the absorption is induced by a π → π * transition mainly localized on the benzene ring of the molecule so that the TS point can be interpreted as an avoided crossing between the two quasi-diabatic states (ππ * ) and (ref ) mr .

From the detailed discussion of Sec. II, we can already deduce that the addition of electron dynamic correlation is not a mere refinement of the picture at the CASSCF level as we originally thought in Ref. 40. On the contrary, the stabilization of zwitterionic states drastically modifies the rationalization of the process: not only the transition state (TS) disappears but also a larger anharmonicity appears along the path leading from M 0 to CI on S 0 . This is well summarized by the comparison of the two scans with either CASSCF or XMCQDPT2 depicted on Fig. 7. The energies at the SA-4-CASSCF (12,11)/TZVP and MS-4-SA-4-XMCQDPT2 (12,11)/TZVP levels of theory are given along two different energy cuts, where the S 0 minimum correspond to x = 12 and the minimum of the CI seam to x = 34. The arrows highlight the main changes for the post-CASSCF treatment: the disappearance of TS and the high anharmonicity along the path leading from M 0 to CI on S 0 . The analysis of the dominant configurations confirms that these changes are mainly due to the stabilization of the zwitterionic state, S 3 . The latter has another physical effect. The absorption of light corresponds to a π → π * excitation, but whereas, for CASSCF calculations, the excited orbitals at the FC geometry on S 1 are mainly the π orbitals on the benzene moiety, for XMCQDPT2 calculations, these are the π orbitals on the pyran moiety that are mainly involved. In other words, we observe a π pyr → π pyr * excitation, which partly explains the disappearance of the transition barrier, since no energy trans-fer is required from the benzene to the pyran ring. To illustrate this more precisely, Fig. 8 shows the potential energy gradient at the FC geometry on S 1 . This figure can be compared with the Fig. 6 (top) of Ref. 40 : at the CASSCF level, the gradient is very different and looks similar to the breathing mode in benzene. Of course, such changes will have an impact on the model originally developed in Ref. 40 : after several attempts, it became obvious that we could no longer fit the two adiabatic PESs with a Hamiltonian matrix based on three quasi-diabatic states as was achieved at the CASSCF level even with our 12 polyspherical coordinates. In Ref. 40, only harmonic functions had been used for the diagonal elements of the potential matrix. But even with higher-order expressions or Morse functions, we could not fit the adiabatic PESs. The most severe problem we had to face was the difficulty to reproduce the topography around the CI owing to the strong anharmonicity of S 0 (although the CI could be stabilized at the right position). In a energy-based diabatisation, the ability to reproduce the topography around the CI is an important test to ensure that the non-adiabatic coupling (NAC) is described correctly in the model [START_REF] Köppel | Diabatic representation: methods for the construction of diabatic states[END_REF] . The NAC is the basic quantity that governs the internal conversion between same-spin electronic states around conical intersections [START_REF] Lengsfield | Non-adiabatic interactions between potential energy surfaces: Theory and applications[END_REF] and the possibility to reproduce the splitting of the two adiabatic states within the so-called branching plane is a signature, at least at the first order, that the NAC is correctly taken into account in the model. This plane is formed by the two directions, the two branching space (BS) vectors, that lift degeneracy to first order. These vectors are known as the gradient-difference (GD) and derivative-coupling (DC) vectors [START_REF] Davidson | Global topology of triatomic potential surfaces[END_REF][START_REF] Atchity | Potential energy surfaces near intersections[END_REF][START_REF] Yarkony | Diabolical conical intersections[END_REF] . The first conclusion is thus the following. Determining the quasi-diabatic states to be included in the model for the potential energies at the CASSCF stage and assuming that post-CASSCF calculations will merely refine the model cannot work for a complex system such as benzopyran. This probably has a general scope, except for systems with a small number of active electrons such as ethylene [START_REF] Lasorne | The electronic excited states of ethylene with large-amplitude deformations: A dynamical symmetry group investigation[END_REF][START_REF] Somoza | A generalised 17-state vibronic-coupling Hamiltonian model for ethylene[END_REF] . If the number of active electrons increases, it becomes more likely that spurious zwitterionic states may strongly interact with the reference states for post-CASSCF treatments. One solution could be to add even more quasi-diabatic states (associated with the zwitterionic configurations) to reproduce the anharmonicity observed on the adiabatic PESs. However, the procedure to choose and to locate theses new quasi-diabatic states in the configuration space is not straightforward since the electronic configurations are intrinsically mixed almost everywhere. This is the reason why we have adopted a new strategy here.

Since it is very difficult to reproduce the topography around the CI, one solution, following the philosophy of Yarkony and coworkers [START_REF] Schuurman | On the vibronic coupling approximation: A generally applicable approach for determining fully quadratic quasidiabatic coupled electronic state hamiltonians[END_REF] , is to build an effective model that enforces this topography around the CI to be correct from the onset. This guarantees, at least to some order and locally, that the NAC is correctly described. The next step is thus to generate the BS vectors at the conical intersection (here, we mean the minimum of the seam of the CI). The definition of the first-order non-adiabatic couplings among electronic states and thus of the BS vectors relies on the availability of the electronic wavefunctions and the knowledge of how they vary with respect to nuclear coordinates. In particular, the DC vector requires in principle to produce the first derivatives of the electronic wavefunctions with respect to the nuclear coordinates. Even though much progress has been achieved in the analyticgradient techniques to the production of DC vectors [START_REF] Shiozaki | Communication: Extended multi-state complete active space second-order perturbation theory: Energy and nuclear gradients[END_REF][START_REF] Chen | Analytical derivative techniques for TDDFT excited-state properties: Theory and application[END_REF] , this approach is tractable only for relatively small molecules and is not available yet for the XMCQDPT2 method.

Fortunately, the branching space can be obtained in the vicinity of a conical intersection from energies only. The energy difference is not smooth along both BS vectors around a CI but its square is regular [START_REF] Bearpark | A direct method for the location of the lowest energy point on a potential surface crossing[END_REF] . Based on an eigenvector analysis of the square energy difference of the CI, this has led to the design of a procedure for generating BS vectors, from energies and their derivatives only, with no direct involvement of wavefunctions [START_REF] Maeda | Updated Branching Plane for Finding Conical Intersections without Coupling Derivative Vectors[END_REF] . This strategy was first proposed by Köppel and Schubert for a two-dimensional model [START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF] and further exploited in more general situations in Ref. [START_REF] Maeda | Updated Branching Plane for Finding Conical Intersections without Coupling Derivative Vectors[END_REF] . More recently, we have reinvestigated various formulations of this approach and x (1) (GDV)

x (2) (DCV) FIG. 9.

The two vectors of the BS at the MS-4-SA-4-XMCQDPT2 (12,11)/TZVP level of theory: the gradient difference vector (GDV) and the derivative coupling vector (DVC). The vectors were obtained with a numerical approach 66 via finite differences with a step of 0.005 Å (note that this procedure produces orthogonal BS vectors by construction).

have shown that if one relies on the calculation of a large number of single points, the method is very robust with respect to numerical differentiation [START_REF] Gonon | On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections[END_REF] . We have then used the numerical approach of Ref. 66 to obtain the BS of benzopyran : the two vectors, x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) , shown on Fig. 9.

Let us stress again here (see also Appendix A) that, in a general situation, both BS vectors (x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) , also often denoted g and h, see Ref. [START_REF] Yarkony | On the adiabatic to diabatic states transformation near intersections of conical intersections[END_REF] ) have no reason to be orthogonal nor to have the same length (except for a Jahn-Teller problem). They are defined only up to some gauge angle related to the arbitrary rotation within a degenerate pair among all possible choices of so-called crude adiabatic states. There will always exist a specific rotation providing orthogonal BS vectors, hence made unique up to arbitrary permutation and signs, as shown by Yarkony in Ref. [START_REF] Yarkony | On the adiabatic to diabatic states transformation near intersections of conical intersections[END_REF] (they will coincide with symmetry-specific ones in the case of asymmetry-induced conical intersection). Reversely, rotating two vectors of different lengths through an angle that is not a multiple of 90 • produces nonorthogonal vectors. As first pointed out by Köppel and Schubert in Ref. [START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF] , which we further addressed recently in Ref. [START_REF] Gonon | On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections[END_REF] , a possible pair of BS vectors can be generated via the diagonalization of the Hessian of the square energy difference. Being eigenvectors of a symmetric matrix, they are orthogonal and thus coincide with the unique (up to permutation and signs) pair of orthogonal BS vectors. From a geometric point of view, they correspond to the semi-major and semi-minor axes of the ellipse corresponding to isovalues of the energy difference in the vicinity of the conical intersection.

Interestingly, the GD vector corresponds mainly to an "outof-plane" motion and the DC vector corresponds mainly to an "in-plane" motion, although the molecule has no real planar geometry. This "quasi-symmetry" is all the more surprising given that the XMCQDPT2 geometry of CI is rather twisted (i.e. even less planar than the CASSCF one: see Fig. 6).

Besides the fact that the molecule has no specific symmetry, the CI point is located far from the FC point. To illustrate this, let us define an effective "reaction coordinate" along the vector x (0) joining FC to CI,

x (0) = x CI -x FC , (1)
where x CI and x FC are the Cartesian coordinates of the CI and FC (identical to M 0 ) points in the body-fixed (BF) frame, respectively. Tab. II provides the normalized overlaps between the four vectors, x (0) , g FC (the S 1 potential energy gradient at FC), x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] (the gradient difference vector at CI), and x (2) (the derivative coupling vector at CI). Each vector here is canonically identified to the collection of its Cartesian components in the BF frame. The small overlap between x (0) and the two vectors x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) defining the two-dimensional branching space indicates that the path leading from FC to CI must be somewhat curved. The vector g FC is related to the velocity of the wavepacket when leaving the Franck-Condon region, and the width of the absorption spectrum is determined by the slope of the PES along this direction [START_REF] Schinke | Photodissociation Dynamics[END_REF] . The strong overlap between x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and g FC allows us to come to the conclusion that the subspace built by x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] , x (2) and x (0) is sufficient to produce a minimal model, relevant for the first step of the dynamics, which is the main goal of the present work.

Normalized overlap x (0) gFC x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] x (2) x (0)

1.0000 -0.0439 0.1240 -0.2130 gFC -0.0439 1.0000 -0.7448 -0.1767 x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] 0.1240 -0.7448 1.0000 0.0000 x (2) -0.2130 -0.1767 0.0000 1.0000 TABLE II. Normalized overlaps (scalar products between normalized vectors) between four important vectors: x (0) , the vector joining FC to CI, gFC, the S1 potential energy gradient at FC, x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] , the gradient difference vector at CI, and x (2) the derivative coupling vector at CI.

As already pointed out, the dynamics of ultra-fast photochemical processes is governed by non-adiabatic transitions that are intense around conical intersections where adiabatic PESs cross [START_REF] Lasorne | Excited-state dynamics[END_REF][START_REF]Conical Intersections, Theory, Computation and Experiment[END_REF][START_REF]Conical Intersections, Electronic Strucutre, Dynamics and Spectroscopy[END_REF] . Quantum-chemistry calculations produce adiabatic energies and non-adiabatic coupling vectors, but a quasi-diabatic representation is often more suitable for the numerical treatment in quantum-dynamics simulations [START_REF] Mead | Conditions for the definition of a strictly diabatic electronic basis for molecular systems[END_REF] . In diabatic representations, the large non-adiabatic or derivative couplings that occur in the adiabatic representation, with possible singularities, are replaced by smooth potential coupling terms. In Ref. 40, we had developed a model based on CASSCF calculations that could be viewed as a generalization of the vibronic-coupling model (VCH) model formulated by Köppel and coworkers [START_REF] Köppel | Multimode molecular dynamics beyond the Born-Oppenheimer approximation[END_REF][START_REF] Köppel | Diabatic representation: methods for the construction of diabatic states[END_REF][START_REF] Cederbaum | Born-oppenheimer approximation and beyond[END_REF] . In this approach, we had included more quasi-diabatic than adiabatic states to introduce more anharmonicity and simply describe the transition state between FC and CI. The presence of additional anharmonicity due to the stabilization of the zwitterionic states at the post-CASSCF level prompted us to switch to another approach. The latter is closer, to some extent, to the approach developed by Köppel and coworkers for instance for the photoinduced ring-opening of furan [START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] : here, the diagonal parts of the diabatic potential matrix are expressed as Taylor expansions along some active coordinates. This allows one to describe the strong anharmonicity of the potential and the complexity of the path leading from FC to CI. The off-diagonal element of the diabatic potential matrix is simply a linear term with respect to a unique "coupling" coordinate. This approach can be viewed as an approximate treatment of the concept of regularized diabatic states [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF] . However, there is in our problem a noticeable difference : the molecule has no specific symmetry, while this usually helps to determine the "coupling" coordinate.

Our goal is to build a new type of model that describes the first part of the photochemical process from FC to CI, preserving the subspace built by the three vectors {x (0) , x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] , x (2) } and simplifying as much as possible the expression of the Hamiltonian operator. Since symmetry cannot help us to define the coupling coordinate, we have to make a choice that will also specify the two quasi-diabatic states (see Appendix A). In the aim of designing the simplest possible form for the diabatic Hamiltonian, we apply two transformations: 1) rotation of the diabatic states to impose zero coupling at FC, and 2) rotation of the rectilinear coordinates to confine the coupling only along one coordinate in the spirit of the effective mode transformation [START_REF] Gindensperger | Short-time dynamics through conical intersections in macrosystems. I. Theory: Effective-mode formulation[END_REF] (note that the third vector, namely the gradient average, denoted s by Yarkony [START_REF] Yarkony | On the adiabatic to diabatic states transformation near intersections of conical intersections[END_REF] , complementing the gh branching plane plane with the gradient of the intersection space so as to obtain the three directions required for a linear expansion of the potential energy surfaces around the conical intersection, is absent from our description because we consider the optimized conical intersection, i.e., the minimum of the intersection seam, which requires s be parallel to g-h). The procedure is described in appendix A. It results in the construction of a new set of rectilinear coordinates, which form a basis for the 3-dimensional subspace B = {X 0 , X 1 , X 2 } (the three basis vectors are now normalized and orthogonal; see Appendix A for their detailed definitions). By construction, the diabatic coupling only depends on X 2 . Also by construction, B contains the BS at CI and the effective vector leading from FC to CI. The three figures are shown on Fig. 10. From now onwards, we will express the nuclear geometries in the reduced basis with FC as the origin of the frame:

X = (X 0 , X 1 , X 2 ) , X 0 = X 0 • x , X 1 = X 1 • x , X 2 = X 2 • x , Q FC = (Q FC, 0 , Q FC, 1 , Q FC, 2 ) = (0, 0, 0) Å , Q CI = (Q CI, 0 , Q CI, 1 , Q CI, 2 ) = (2.54, 0, 0) Å , (2)
where x represents any point in the full-dimensional space, and X its projection onto the three-dimensional subspace. The notation Q is used instead of X when it corresponds to the coordinates of a specific point such as FC or CI. Note that the origin of the frame attached to B = {X 0 , X 1 , X 2 } is the FC point.

B. Hamiltonian operator

The choice of our coordinates allows us to split our active subspace into one "diabatic subspace" built by X 0 and X 1 and one coupling coordinate, X 2 . It leads to write the off-diagonal element of the diabatic potential matrix as a single linear term in X 2 : (2) identifies to x (2) , i.e. the length of

H 1 2 (X) = x (2) • (X -Q CI ) = x (2) 2 (X 2 -Q CI, 2 ) , (3) where x (2) 2 = X 2 •x
x (2) (since x (2) a = X a • x (2) = 0 for a = 2 by construction; see Appendix A).
The expression is identical to the term in the off-diagonal element in Eq. ( 5) of Ref. 39 except that in Ref. 39, the choice of the coupling coordinate has been made following symmetry considerations, which is not the case here. Our specific choice of X 0 , which is perpendicular to the coupling coordinate, X 2 , allows us to neglect the impact of the coupling H 1 2 (X) along the coordinate X 0 . We assume the following expression for the two diagonal matrix elements, i = 1, 2,

H i i (X) = P 4 i (X 0 , X 1 ) + p i,5 (X 0 -Q CI, 0 ) 5 + p i,6 (X 0 -Q CI, 0 ) 6 + C i (X 1 , X 2 ) . ( 4
)
P 4 i (X 0 , X 1
) is a fourth-order polynomial in both X 0 and X 1 . We added higher-order terms in X 0 , p i,5 (X 0 -Q CI,0 ) 5 and p i,6 (X 0 -Q CI,0 ) 6 in order to account for the strong anharmonicity we observed with post-CASSCF calculations. The first two lines correspond to the two PESs for the subspace {X 0 , X 1 } without any coupling between the two electronic states (X 2 = 0). The last term, C i (X 1 , X 2 ), contains terms in X 2 as well as the correlation between X 1 and X 2 . This results in 17 parameters for each quasi-diabatic PES. The parameters defining the polynomial P 4 i , as well as the parameters p i,5 and p i,6 , are obtained by a fit of the ab initio adiabatic energies on a two-dimensional grid (containing 315 points) in X 0 and X 1 , where X 2 = 0. On this grid, the terms C i (X 1 , X 2 ) are equal to zero, therefore the adiabatic and the diabatic potentials coincide and a direct linear fit of the ab initio adiabatic energies on the 2D-grid can be used (a simple method of least squares on a two-dimensional grid can be used). The parameters of the grid are given in Tab. III. Note that, at this level, we do not use the grid along X 2 . We give higher weights to the points in the most relevant regions, especially, around FC and CI. In addition, p i,6 and the parameters for the fourth-order terms are enforced to remain positive, in order to avoid holes in the fitted PESs. The functions C 1 (X 1 , X 2 ) and C 2 (X 1 , X 2 ) are extracted from local derivatives around FC by a semi-analytical procedure. The basis vector X 2 is totally included in the branching space: it is parallel to x (2) but not orthogonal to x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] because x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) are not orthogonal (while x (1) ∝ X 1 and x (2) ∝ X 2 are; see Appendix). Thus, the first-order terms along X 2 are directly obtained from the vectors of the branching space. More precisely, keeping first-and secondorder terms only and assuming that there is no coupling term between X 0 and X 2 (as in Ref. 39), the functions

C 1 (X 1 , X 2 ) and C 2 (X 1 , X 2 ) read C 1 (X 1 , X 2 ) = f1, 12 X 1 X 2 + f1, 22 X 2 2 , (5)
and

C 2 (X 1 , X 2 ) = 2x (1) 2 X 2 + f2, 12 X 1 X 2 + f2, 22 X 2 2 , (6) where x (1) 2
= X 2 • x (1) ∝ x (2) • x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] . In other words, the linear term along X 2 in C 2 (X 1 , X 2 ) reflects the nonorthogonality of x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) . It corresponds to an expansion along the gradient difference at both CI and FC

(Q CI, 1 = Q FC, 1 = 0 and Q CI, 2 = Q FC, 2 = 0). Since Q FC corresponds to a minimum of H 1 1 (X)
, in particular along X 1 , we have enforced the linear term to be carried by H 2 2 (X), which is the reason why there is a 2x

(1) 2 X 2 term in C 2 (X 1 , X 2 ) rather than ±x (1) 2 X 2 terms in both C 1 (X 1 , X 2 ) and C 2 (X 1 , X 2 ).
Since we made the choice of Q ref = Q FC to determine the angle θ and thus the quasi-diabatic states (see Appendix A), the coupling between the two electronic states is zero by construction at this geometry. This allows us to obtain the parameters fi, ab analytically, the corresponding terms being understood as a second-order Jahn-Teller correction (see for instance Eq. (4.5) in Ref. 69):

f1, ab = 1 2 ∂ 2 H 11 (Q FC ) ∂X a ∂X b = 1 2 ∂ 2 1 (Q FC ) ∂X a ∂X b -2 x (2) a x (2) b 2 (Q FC ) -1 (Q FC ) , (7)
f2, ab = 1 2

∂ 2 H 22 (Q FC ) ∂X a ∂X b = 1 2 ∂ 2 2 (Q FC ) ∂X a ∂X b + 2 x (2) a x (2) b 2 (Q FC ) -1 (Q FC ) , (8)
where 1 and 2 denote the adiabatic potential energies. The adiabatic second derivatives at FC appearing in the above equations were obtained in a first step by finite difference using the ab initio adiabatic energies around this point. Since we assumed that a linear model is valid for the coupling (see Eq.

(

)), x 3
is directly provided by the numerical calculation of the branching space at CI 66 (let us remind that x

(2) a = 0 for a = 2 by construction; see Appendix A).

At this level, we stress that, for a molecule with a specific symmetry (like the planar molecule in Ref. 39), the linear and cross terms in Eqs. ( 5) and ( 6) would have disappeared. The second-order Jahn-Teller "corrections" in Eqs. ( 5) and ( 6) can be viewed to some extent as a sort of local regularization. On the other hand, in the approach based on regularized diabatic states [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF] by Köppel and coworkers, the correction due to the regularization scheme has a more global character. The main difference between the present approach and the model developed in Ref. 39 is, besides the problem of symmetry, the fact that in Ref. 39 the coupling term x (symbolized by λ in Ref. 39) is obtained by a non-linear fit (see Eq. (7) in Ref. 39). Here, the non-linear fit of Ref. 39 is avoided. This nonlinear fit would have been very complicated due to the lack of symmetry. We do not enforce a global correction (or regularization) due to the presence of the coupling, rather we enforce a correct topography at FC and around the minimum of the conical intersection.

Our model is based on the idea that the quasi-diabatic potential energy surfaces can be fitted by simple functions in X 2 and that there is a relative decoupling between X 2 and the subspace {X 0 , X 1 } along the path from FC to CI.

IV. RESULTS

A. Potential energy surfaces

To test the validity of the model, we will compare the two potentials obtained after diagonalizing the Hamiltonian matrix with the two ab initio adiabatic PESs. The criterion that the model is correct is the ability to reproduce the ab initio data for the adiabatic PESs. One of the main problems we had to face was the presence of a very strong anharmonicity in particular along the path leading from FC to CI. Fig. 11, which gives the cuts through the ab initio and model S 0 and S 1 PESs along the effective X 0 coordinate, proves that this anharmonicity is now correctly described. In addition, Fig. 12 provides similar figures but along X 1 and X 2 and around FC (Fig. 12 (a) and (b)) and around CI ((Fig. 12 (c) and(d)). The most important result is the fact that we can reproduce correctly the topography around CI since it is linked to the NAC. A more global picture is provided by the twodimensional cuts depicted in Figs. 2, 3 and 4 of section II of the supporting material. The figures show that the very complex and anharmonic topographies of the to ab initio PESs are correctly described. The very sinuous valley from FC (at X 0 =0 and X 1 = 0) to CI (at X 0 =2.54 a.u. and X 1 = 0 ) on S 0 along X 0 and X 1 = 0 is very well reproduced (see Figs. 2 (a) and (c) in section II of the supporting material). We also see that the local minimum of the open form of the molecule after CI corresponds to the second quasi-diabatic states (the quasi-diabatic states are depicted on Figs. 2, 3 and 4 (e) and (f)). In view of Figs. 3 (a), (b), (c), and (d) in section II of the supporting material, we observe that the curvature of the paths from FC to CI on both electronic states are also correctly described. The curvature is generated by the coupling between the two quasi-diabatic states only since there is no curvature along X 2 = 0 for the quasi-diabatic states (Figs. 3 (e),(f)). This justifies retrospectively our assumptions in the model, i.e. a relative decoupling between {X 0 ; X 1 } and X 2 , as well as a simple treatment for X 2 . Finally, Figs. 3 (a), (b), (c), and (d) in section II of the supporting material are complementary to (Fig. 12 (c) and (d) : we see that the topography around the CI point is correct even in two-dimensions, which is the most difficult region to reproduce. The root-mean square deviation between the ab initio and model PESs is equal to 0.8580 and 0.6876 eV for S 0 and S 1 , respectively. Now, if we define a "relevant" region, where the process takes place and defined by {[-0.4; 2.8], [-0.15; 0.15], [-0.15; 0.15]}, the root-mean square deviation is equal to 0.0880 and 0.1230 eV for S 0 and S 1 , respectively: a difference that is perfectly acceptable for what we need.

Figs. 5 in the supporting material provide the cuts through the model PESs on a large range along the different coordinates: here we depict values of the potential that are very high in energy and that cannot be reached by the system. On the other hand, we see no hole in the PESs. To calculate the vibrational ground in the electronic ground state, we have applied the improved relaxation method [START_REF] Meyer | Calculation and selective population of vibrational levels with the Multiconfiguration Time-Dependent Hartree (MCTDH) algorithm[END_REF][START_REF] Doriol | Computation of vibrational energy levels and eigenstates of fluoroform using the multiconfiguration time-dependent Hartree method[END_REF] of the Heidelberg Multi-Configuration Time-Dependent Hartree (MCTDH) package [START_REF] Worth | The MCTDH Package[END_REF] . This method is very sensitive to the presence of any hole in the PESs. Since we could converge correctly the vibrational ground state, we are confident that no hole is present at least on the coordinate range used for the dynamics.

B. Absorption spectrum

The 3D theoretical absorption spectrum of benzopyran has been calculated using the PES model developed in the present paper, the Heidelberg MCTDH package [START_REF] Meyer | Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Meyer | The multi-configurational time-dependent Hartree approach[END_REF][START_REF] Beck | The multiconfiguration time-dependent Hartree (MCTDH) method: A highly effi-cient algorithm for propagating wave packets[END_REF] , and the kinetic energy operator (KEO) obtained after a simple transformation from the Cartesian coordinates to the three active coordinates of the model (the masses that appear in the KEO are provided automatically by the TANA program [START_REF] Ndong | Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach[END_REF][START_REF] Ndong | Automatic computer procedure for generating exact and analytical kinetic energy operators based on the polyspherical approach: general formulation and removal of singularities[END_REF] : see Appendix B for their actual values). The Multi-Configuration Time Dependent Hartree method is a time-dependent method in which each degree of freedom is associated with a small number of single particle functions (SPFs) which, through their time de-pendence, allow efficient description of the molecular dynamical processes. The total MCTDH wave function is expanded in products of single-particle functions

Ψ (Q 1 , . . . , Q f , t)= n1 j1=1 • • • n f j f =1 A j1•••j f (t) f κ=1 φ (κ) jκ (Q κ , t) = J A J Φ J , (9)
where f is the number of degree of freedom of the system, Q 1 , . . . , Q f are the nuclear coordinates, A j1•••j f denotes the MCTDH expansion coefficients, and φ

ϕ (κ) jκ (Q κ , t) = Nκ iκ=1 c (κ) iκjκ (t)χ (κ) iκ (Q κ ) (10)
where in general the n κ of Equation ( 9) is such that n κ < N κ . Thus, the MCTDH method propagates a wave function on a small time-dependent, variationally optimized basis set of single-particle functions, which in turn are defined on a timeindependent primitive basis set. The primitive basis, its range and the number of the SPFs used for the calculations of the rovibrational states are summarized in The theoretical spectrum is depicted in Fig. 13 (in blue). The experimental absorption spectrum (in purple) is the one of 2,2-diethylchromene (where the CH 2 group is replaced by C-(CH 2 -CH 3 ) 2 ) at 77 K in 3-methylpenthane. It presents two bands. The band between 280 and 340 nm corresponds to the absorption studied in the present work. The second one between 240 and 280 nm corresponds to the excitation to another electronic state, higher in energy, which is not included in our model. Our spectrum appears in the correct domain of the electronic state of interest: we did not have to shift artificially the spectrum for the comparison with the experiment. In yellow, the spectrum obtained with the 12D model of Ref. 40 is depicted after having been shifted to lower transition energies to account for the lack of electron dynamic correlation of CASSCF calculations. Not only the spectrum obtained with the present model is now at the correct place, it presents several similitudes with the experimental spectrum : there are four peaks, two being shoulders on the left and right ; the widths between the peaks are rather similar. Only the global width of the theoretical spectrum is a little bit too large. This is not surprising since the subspace built by our three active coordinates does not include completely the potential energy gradient at the FC geometry on S 1 (see Tab. II), which is the most important quantity related to the width of the absorption spectrum. To obtain the correct width, it would be necessary to add a fourth coordinate to make sure that this gradient is fully included in the active subspace.

V. CONCLUSIONS

We have developed a new model for the PESs that describe the ring-opening of benzopyran taken as a model compound to study the photochromic ring-opening reaction of indolinobenzospiropyran and its evolution toward its open-chain analogue. This model has been improved with respect to some aspects compared to what we had proposed in Ref. [START_REF] Joubert-Doriol | A generalized vibronic-coupling Hamiltonian model for benzopyran[END_REF] .

The process we are looking at is very complex due to the size of the system, the presence of motions of large ampli-tude, the fact that we have to deal with its reactivity and the lack of any symmetry. Whereas in Ref. 40 we worked at the CASSCF level of theory, we develop here a new model at the XMCQDPT2 level.

Several important points are addressed. Accounting for electron dynamic correlation drastically changes the global picture originally sketched at the CASSCF level: zwitterionic states are stabilized and one of them strongly interacts with the first excited state S 1 of the CASSCF picture. The transition state, that was present between FC and CI, disappears due to this interaction. A multi-state treatment including the zwitterionic states becomes unavoidable. In other words, the usual CASSCF//CASPT2 approach cannot give even the correct qualitative picture of the process. Many other differences are observed: the stationary geometries are even more "twisted" (i.e. less planar) for XMCQDPT2 calculations than for CASSCF ones. The absorption of light corresponds to a π → π * excitation as expected, but whereas, for CASSCF calculations, the excited orbitals at the FC geometry on S 1 are mainly the π orbitals on the benzene moiety, for XMCQDPT2 calculations, it is the π pyr , i.e. the π orbitals on the pyran moiety, that is mainly involved. In other words, we observe a π pyr → π pyr * excitation, which partly explains the disappearance of the TS, since no energy transfer is required from the benzene to the pyran ring. It is expected that the fact that the dynamic correlation is not a perturbative correction for the process is not specific to this system. The same situation may occur in many other molecules.

After having adopted the MS-4-SA-4-XMCQDPT2 (12,11)/TZVP level of theory, we developed an effective model for the PESs that can describe the absorption, the path leading to the CI, found at the post-CASSCF level, and the non-adiabatic couplings at CI. The approach of the present work has some similarities with the treatment based on regularized diabatic states [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF][START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] , but we had to face the major problem of the lack of symmetry, which makes the choice of the diabatic states and the coupling mode less straightforward. To overcome this difficulty, we have proposed to choose the coupling mode so that it is as much decoupled as possible from the sinuous path from FC to CI. We have also enforced the model to reproduce correctly the topography at FC and around CI, which can be seen as a kind of "local" regularization. The results are very encouraging, the main weakness of the present model being probably the low number of degrees of freedom (three here), which cannot describe the open form precisely and the dissipation of energy after reaction. However, the general strategy proposed here may be of general interest and could be used for other systems, several improvements of the model being easily possible in particular by adding more degrees of freedom. For the ring-opening process of spiropyran, we can already predict, from our model, that the first step of the mechanism, corresponding to the nonadiabatic decay via the conical intersection, takes about 50 fs, in agreement with the shortest time constant observed in previous experiments [START_REF] Poisson | Gas-phase dynamics of spiropyran and spirooxazine molecules[END_REF] and in recent experiments including the coherent control of the process that will give rise to a future publication 78 .

VI. APPENDIX A

The branching-space (BS) vectors at a conical intersection are defined up to a gauge rotation angle, which corresponds to the arbitrary choice of a certain pair of degenerate adiabatic states, Φ 1 and Φ 2 , at the geometry of the conical intersection, Q CI . They are defined in principle as

x (1) = 1 2 Φ 2 ; Q CI |∂ Q H el (Q CI )|Φ 2 ; Q CI - 1 2 Φ 1 ; Q CI |∂ Q H el (Q CI )|Φ 1 ; Q CI , (11)
and

x (2) = Φ 1 ; Q CI |∂ Q H el (Q CI )|Φ 2 ; Q CI . (12)
In general situations, such vectors have no reason to be orthogonal in practice, except sometimes for symmetry reasons.

In our case, however, the two BS vectors of the branching space calculated numerically using the procedure exposed in Ref. [START_REF] Gonon | On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections[END_REF] and shown on Fig. 9 have been produced in such a way that they are orthogonal. Now, at the the conical intersection, Q CI , equally acceptable adiabatic states, Φ i , can be defined up to a rotation with an arbitrary angle, θ, independent from the coordinates,

|Φ 1 ; Q CI = cos θ |Φ 1 ; Q CI + sin θ |Φ 2 ; Q CI |Φ 2 ; Q CI = -sin θ |Φ 1 ; Q CI + cos θ |Φ 2 ; Q CI . (13)
The new pair of BS vectors are obtained from the following similarity transformation, -x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] x (2) x (2) x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] = cos θ sin θ sin θ cos θ • -x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] x (2) x (2) x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • cos θsin θ sin θ cos θ , (14) leading to

x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] = cos 2θ x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] sin 2θ x (2) x (2) = sin 2θ x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] + cos 2θ x (2) .

Any value of θ yields a possible pair of BS vectors that lift degeneracy to first order. θ = 0 corresponds only to the result given by the numerical calculation providing arbitrary in principle, but determined in practice, degenerate states. We now consider diabatic states that coincide with the rotated adiabatic states, Φ i , at the conical intersection but not necessarily elsewhere (crude adiabatic states). At the first order, the off-diagonal term of the Hamiltonian matrix in this diabatic basis set (by definition of the derivative coupling vector) reads

H 1 2 (X) = x (2) • (X -Q CI ) , (16)
where X denotes the set of coordinates of a running point. In the applications of the vibronic-coupling (VCH) model [START_REF] Köppel | Multimode molecular dynamics beyond the Born-Oppenheimer approximation[END_REF][START_REF] Köppel | Diabatic representation: methods for the construction of diabatic states[END_REF][START_REF] Cederbaum | Born-oppenheimer approximation and beyond[END_REF] or in the approach based on regularized diabatic states [START_REF] Köppel | Construction scheme for regularized diabatic states[END_REF][START_REF] Köppel | The concept of regularized diabatic states for a general conical intersection[END_REF][START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] , due to using the symmetry of the systems, the term in Eq. ( 16), often involves a single (normal or curvilinear) coordinate only (for instance a coordinate associated with an out-ofplane mode of vibration for a planar molecule), or at least a single irreducible representation. This simplifies significantly the definition of the quasi-diabatic states upon fixing the aforementioned angle θ.

Here, as already explained, we cannot exploit symmetry for our molecule. In practice, choosing the rotating angle, θ, is equivalent to setting a geometry, Q ref , for which we wish the coupling to vanish,

H 1 2 (Q ref ) = 0 . ( 17
)
This is equivalent to defining one direction, x (0) defined such as

x (0) = Q CI -Q ref , (18)
that is orthogonal to the derivative coupling vector for our choice of quasi-diabatic states. In other words,

x (2) • x (0) = 0 .

Thus, in view of Eq. ( 15), this yields sin 2θ x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • x (0) + cos 2θ x (2) • x (0) = 0 ,

sin 2θ cos 2θ = -x (2) • x (0)

x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • x (0) ,

leading to 2θ = arctanx (2) • x (0)

x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • x (0) .

This determines the value of θ and provides the two vectors x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) using Eq. (15). Since x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) are not normalized, x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) have no reason to be orthogonal, even if the former were. A new vector, x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] (orthogonal to x (2) ) can be generated by with the Gram-Schmidt procedure,

x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] = x (1) -x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • x (2) x (2) • x (2) x (2) -x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] • x (2) x (0) • x (0) x (0) . ( 23)

Now, a choice for the vector x (0) in Eq. ( 18) must be made. Following the approach of Köppel and coworkers [START_REF] Gromov | Theoretical study of photoinduced ring-opening in furan[END_REF][START_REF] Gromov | Ab initio quantum dynamical study of photoinduced ring-opening in furan[END_REF] , our goal is to decouple as much as possible the direction of the coupling in the Hamiltonian operator between the two quasidiabatic states, defined by x (2) , from the rest of the "active" subspace in which we will perform the dynamics. At the first order, there is no coupling betweeen x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] and x (2) at the CI point by construction. A natural choice is thus to choose x (0) parallel to the vector X 0 from FC to CI. In other words, we choose the coupling vector, x (2) , and thus the quasi-diabatic states, so that by construction there is as little coupling as possible between X 0 and x (2) . Our goal is to build a twodimensional subspace, {X 0 , X 1 } that is as decoupled as possible from the direction X 2 ∝ x (2) , so that H 1 2 (X) can be written as a linear term along X 2 only. Mathematically, this means that, in Eq. ( 18), Q ref = Q FC . Physically, it also means that the adiabatic and quasi-diabatic states are identical at FC. We thus build the orthonormal basis set B = {X 0 , X 1 , X 2 } such that

X 0 = x (0)
x (0) , X 1 =

x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] , X 2 = x (2) x (2) . ( 24)

VII. APPENDIX B

The Body-Fixed (BF) KEO in terms of the three active (BF) coordinates {X 0 , X 1 , X 2 } reads TBF = -

M i j ∂ 2 ∂X i ∂X j , (25)
in uma -1 . As explained in Appendix A, the three vectors {X 0 , X 1 , X 2 } were obtained from orthonormalizing the three original vectors denoted {x (0) , x [START_REF] Liu | Photoelectrochemical information storage using an azobenzene derivative[END_REF] , x (2) }. They correspond to linear combinations of Cartesian basis vectors in the BF frame, thus defining a rectangular transformation matrix with constant coefficients. The standard full-dimensional KEO has a similar expression to the one above, but with derivatives with respect to all Cartesian coordinates and a reciprocal mass matrix containing the inverse of the masses of all atoms (repeated three times, once for each direction). The three-dimensional KEO was obtained via a mere projection of the full-dimensional one onto the three-dimensional space, using the rectangular transformation matrix that connects both types of derivatives.

VIII. SUPPLEMENTARY MATERIAL

See supplementary material for the description of the active orbitals and a thorough comparison between the ab initio and model potential energy surfaces.

FIG. 1 .

 1 FIG. 1. (a)Photo-induced ring opening of indolinobenzospiropyran into its merocyanine analogue. (b) The π-system of the chromene moiety is orthogonal to the indoline moiety, thus explaining that the absorption spectrum corresponding to the π → π * excitation is, to a first approximation, the sum of the spectra of the two orthogonal fragments. (c) Ring opening of the benzopyran molecule. We assume as in Ref.[START_REF] Celani | Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds[END_REF] that benzopyran can be used as a model compound to study the photochromic ring opening of indolinobenzospiropyran.

  FIG. 3. (a) Effect of dynamic correlation in CASPT2 calculations on the first four electronic states: scans with SA-2-CASSCF(12,11)/ANO-S-VTZP (black) and with SA-2-CASPT2(12,11)/ANO-S-VTZP (red). (b) MS-CASPT2 applied to the first four electronic states: SA-2-CASSCF(12,11)/ANO-S-VTZP (black) as a reference, along with SA-2-CASPT2(12,11)/ANO-S-VTZP (red), MS-2-SA-2-CASPT2(12,11)/ANO-S-VTZP (green), and MS-4-SA-4-CASPT2(12,11)/ANO-S-VTZP (blue). Note that SA-2-CASPT2 and MS-2-SA-2-CASPT2 are indistinguishable from M0 to CI. Note also that the new angular feature around x = 11 for MS-2-CASPT2 energies is due to a discontinuity of the second reference wavefunction from one to the other side of the avoided crossing marked by the S1 TS for CASSCF energies. This effect disappears when increasing the size of the reference space (MS-4-CASPT2 energies).

FIG. 5 .FIG. 6 .

 56 FIG. 4. (a) MS-XMCQDPT2 energies: the curve SA-2-CASSCF(12,11)/TVZP in black serves as a reference. The other curves are with MS-XMCQDPT2: MS-2-SA-2-XMCQDPT2(12,11)/TVZP in red, MS-4-SA-2-XMCQDPT2(12,11)/TVZP in green, MS-8-SA-2-XMCQDPT2(12,11)/TVZP in blue. (b) Influence of the number of states in the state-averaged procedure: SA-2-CASSCF(12,11)/TVZP in black serves again as a reference. The MS-4-SA-2-XMCQDPT2(12,11)/TVZP curve is in red, the MS-4-SA-4-XMCQDPT2(12,11)/TVZP curve in green, and the MS-5(12347)-SA-5-XMCQDPT2(12,11)/TVZP curve in blue.

FIG. 7 .

 7 FIG. 7. Comparison between scans with CASSCF and XMC-QDPT2 calculations. In red are given the energies at the SA-4-CASSCF(12,11)/TZVP level of theory, in yellow those at the MS-4-SA-4-XMCQDPT2(12,11)/TZVP one. The global S0 minimum corresponds to x = 12, the CI minimum to x = 34 (the corresponding geometries are optimized consistently with the levels of theory used for the energies).

FIG. 10 .

 10 FIG. 10. Two different view of the vectors X0, X1, and X2 at CI (from above and side).

FIG. 11 .FIG. 12 .

 1112 FIG.11. Cuts through the ab initio and model S0 and S1 PESs along the effective X0 coordinate (X1 = 0 ; X2 = 0) : in green and blue the ab initio adiabatic energies at the XMCQDPT2 level of theory ; in red and light blue the energies given by the model.

  jκ (Q κ , t) are the n κ SPFs associated with each degree of freedom κ. The subsequent equation of motion for the coefficients and single particle functions are derived after replacing the wave function ansatz into the time dependent Schrödinger equation. To solve the equations of motion, the κ SPF are represented on a primitive basis or discrete variable representation (DVR) grid of N κ points,

FIG. 13 .

 13 FIG. 13. Theoretical absorption spectra of benzopyran (blue line). Purple line: experimental absorption spectrum of 2,2diethylchromene (the CH2 group is replaced by C-(CH2-CH3)2) at 77 K in 3-methylpenthane. In yellow the spectrum with the 12D model of Ref. 40 at the CASSCF level of theory shifted to smaller to lower transition energies.

TABLE I .

 I Energy values (absolute values in hartree and relative values in eV) at the four relevant points optimized at the MS-4-SA-4-XMCQDPT2

	Stationary point (state) XMCQDPT2 (hartree) XMCQDPT2 (eV) CASSCF (eV)
	M0 (S0)	-421.82563	0.00	0.00
	FC (S1)	-421.67799	4.02	4.65
	FC (S2)	-421.65001	4.78	6.53
	CI (S0/S1)	-421.73801	2.38	2.64
	cis (S0)	-421.78991	0.97	0.95

TABLE III .

 III Parameters of the (X0, X1, X2)-grid for the potential.

		X0	X1	X2
	Number of grid points	21	15	15
	Grid spacing ( Å)	0.2	0.05	0.05
	Grid length ( Å)	[-0.6;3.4] [-0.35;0.35] [-0.35;0.35]

Table IV

 IV

					. The spec-
	trum is obtained by first calculating the vibrational ground
	state in the electronic ground state and by a Fourier transform
	of the autocorrelation function of the vibrational wavepacket
	propagated on the excited state (Franck-Condon transition)
	during T = 100 fs (nonadiabatic decay from S 1 to S 0 via
	the conical intersection starts around 50 fs and is complete
	in about 100 fs). We have also used a damping function, more
	precisely a Gaussian function, cos(π/2 * T ) exp (-t/τ ) with
	τ = 10 fs, multiplying to the autocorrelation function in the
	Fourier transform, to convoluate the spectrum in order to com-
	pare the theoretical spectrum with the experimental data.
	Coords. Primitive Number of	Range	SPF
		Basis	Points		basis
	X0	HO	150	[-2 ; 90, 90
	X1	HO	70	[ -0.8 ; 0.8 ] 65 , 65
	X2	HO	60	[-0.7 ; 0.7 ] 50 , 50

TABLE IV .

 IV Parameters of the primitive basis used for the calculation of absorption spectrum. HO stands for the Harmonic Oscillator DVR. The units for the coordinates are angströms.

ACKNOWLEDGMENTS

Numerous discussions with Prof. Hans-Dieter Meyer and Prof. Horst Köppel are gratefully acknowledged. We acknowledge support from the CoConicS Project (ANR-13-BS08-0013-03).