Evgeny Posenitskiy
email: posenitskiy@irsamc.ups-tlse.fr

Fernand Spiegelman

Didier Lemoine

ON APPLICATION OF DEEP LEARNING TO SIMPLIFIED QUANTUM-CLASSICAL DYNAMICS IN ELECTRONICALLY EXCITED STATES PREPRINT

Deep Learning (DL) is applied to simulate non-adiabatic molecular dynamics of phenanthrene, using the Time-Dependent Density Functional based Tight Binding (TD-DFTB) approach for excited states combined with mixed quantum-classical propagation. Reference calculations rely on Tully's Fewest-Switches Surface Hopping (FSSH) algorithm coupled to TD-DFTB, which provides electronic relaxation dynamics in fair agreement with various available experimental results. Aiming at describing the coupled electron-nuclei dynamics in large molecular systems, we then examine the combination of DL for excited-state potential energy surfaces with a simplified trajectory surface hopping propagation based on the Belyaev-Lebedev scheme. We start to assess the accuracy of the TD-DFTB approach upon comparison of the optical spectrum with experimental and higher-level theoretical results. Using the recently developed SchNetPack[1] for DL applications, we train several models and evaluate their performance in predicting excited-state energies and forces. Then, the main focus is given to the analysis of the electronic population of low-lying excited states computed with the aforementioned methods. We determine the relaxation timescales and compare them with experimental data. Our results show that DL demonstrates its ability to describe the excited-state potential energy surfaces. When coupled to the simplified Belyaev-Lebedev scheme considered in this study, it provides reliable description of the electronic relaxation in phenanthrene as compared with either the experimental data or the higher-level FSSH/TD-DFTB theoretical results. Furthermore, the Deep Learning performance allows high-throughput analysis at a negligible cost.

Introduction

Modelling the properties and evolution of large molecular systems is a challenging task in chemical and biological sciences, which usually requires solving the stationary or the time-dependent Schrödinger equation. Recently, Machine Learning (ML) has appeared as a promising tool that can be used to fully or partially avoid electronic structure calculations in atomistic simulations. [START_REF] Schütt | Quantumchemical insights from deep tensor neural networks[END_REF][START_REF] Volker | Machine learning interatomic potentials as emerging tools for materials science[END_REF][START_REF] Westermayr | Machine learning and excited-state molecular dynamics[END_REF][START_REF] Pavlo | Quantum chemistry in the age of machine learning[END_REF][START_REF] Oliver | Physnet: A neural network for predicting energies, forces, dipole moments, and partial charges[END_REF][START_REF] Meuwly | Transformative applications of machine learning for chemical reactions[END_REF] ML methodologies may differ depending on the objective and can be categorized in three main branches: (i) supervised learning, (ii) unsupervised learning and (iii) reinforcement learning. The major focus in reinforcement learning is given to finding a balance between exploration (of uncharted areas) and exploitation (of current knowledge). [START_REF] Pack Kaelbling | Reinforcement learning: A survey[END_REF] The main objective of the unsupervised learning is to look for previously undetected patterns in a dataset. The supervised learning aims to find optimal parameters for a function that maps an input to an output based on a given set of input-output pairs (also called the training set). Ideally, the supervised ML model should infer some knowledge from the training set in order to perform reasonably well on new data entries. For example, the global structure optimization can be enhanced by each of the aforementioned techniques. [START_REF] Sørensen | Accelerating atomic structure search with cluster regularization[END_REF][START_REF] Jørgensen | Combining evolutionary algorithms with clustering toward rational global structure optimization at the atomic scale[END_REF][START_REF] Jacobsen | On-the-fly machine learning of atomic potential in density functional theory structure optimization[END_REF][START_REF] Malthe | Efficient global structure optimization with a machine-learned surrogate model[END_REF][START_REF] Jørgensen | Atomistic structure learning[END_REF] Artificial Neural Networks (NNs) have been actively studied and used in various fields of science, thus forming a Deep Learning domain of ML. [START_REF] Gasteiger | Neural networks in chemistry[END_REF] In particular, it has been proven that multilayer feed-forward NNs are "universal approximators", [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] i.e. they can approximate an unknown multidimensional function with an arbitrary accuracy based on a set of known function values. The more hidden layers and nodes are included, the more flexible functional form of the network is obtained. The supervised Machine Learning models can be used to construct complex interatomic potentials [START_REF] Behler | Generalized neural-network representation of high-dimensional potentialenergy surfaces[END_REF][START_REF] Volker | Machine learning interatomic potentials as emerging tools for materials science[END_REF] that can in turn be used to perform extensive molecular dynamics simulations. [START_REF] Gastegger | Machine learning molecular dynamics for the simulation of infrared spectra[END_REF][START_REF] Shakouri | Accurate neural network description of surface phonons in reactive gas-surface dynamics: N2 + ru(0001)[END_REF][START_REF] Zhang | Bridging the gap between direct dynamics and globally accurate reactive potential energy surfaces using neural networks[END_REF] Ideally, this would correspond to an ab initio accuracy at the computational cost of a force field. Some supervised ML models include forces in the training process, which drastically improves the quality of a machine-learned Potential Energy Surface (PES). [START_REF] Chiriki | Neural network potentials for dynamics and thermodynamics of gold nanoparticles[END_REF][START_REF] Li | Molecular dynamics with on-the-fly machine learning of quantum-mechanical forces[END_REF][START_REF] Zhang | Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics[END_REF][START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] Alternatively, NNs have also been used to correct DFT [START_REF] Ramakrishnan | Electronic spectra from tddft and machine learning in chemical space[END_REF] and DFTB [START_REF] Zhu | Artificial neural network correction for density-functional tight-binding molecular dynamics simulations[END_REF] results based on a ∆-Machine Learning [START_REF] Ramakrishnan | Big data meets quantum chemistry approximations: The δ-machine learning approach[END_REF] approach. Several software packages are available for applications of supervised ML to atomistic simulations, e.g. MLatom [START_REF] Pavlo | Mlatom: A program package for quantum chemical research assisted by machine learning[END_REF], DeePMD-kit [START_REF] Wang | Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics[END_REF] and SchNetPack [START_REF] Schütt | Schnetpack: A deep learning toolbox for atomistic systems[END_REF].

Addressing electronically excited states and non-adiabatic dynamics in excited states obviously proves both methodologically and computationally much more demanding than merely ground electronic state properties and dynamics, which involves among other things complex shapes of PES and the occurrence of conical intersections. Several theoretical approaches have been devised to incorporate non-adiabatic effects in the dynamics. The most accurate type is fully ab initio and deals with quantum dynamics of nuclei, such as the Multi-Configuration Time-Dependent Hartree (MCTDH) method [START_REF] Meyer | The multi-configurational time-dependent hartree approach[END_REF][START_REF] Beck | The multiconfiguration time-dependent hartree (mctdh) method: a highly efficient algorithm for propagating wavepackets[END_REF] propagating time-dependent wavepackets. Alternatively, mixed quantum-classical schemes have been developed based on a classical description of the nuclear motion, e.g. the mean-field propagation [START_REF] Ehrenfest | Bemerkung über die angenäherte gültigkeit der klassischen mechanik innerhalb der quantenmechanik[END_REF] or the Trajectory Surface Hopping (TSH) approach. [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF][START_REF] Fabiano | Implementation of surface hopping molecular dynamics using semiempirical methods[END_REF][START_REF] Crespo | Recent advances and perspectives on nonadiabatic mixed quantum-classical dynamics[END_REF] Within the generalized TSH picture, the nuclear wavepacket evolution is simulated by an ensemble of independent classical trajectories and non-adiabatic effects are taken into account via a probability to switch from the current electronic state to another one. Recently, a Deep Learning architecture SchNet [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] has been interfaced with SHARC [START_REF] Mai | Nonadiabatic dynamics: The sharc approach[END_REF] code for non-adiabatic molecular dynamics. [START_REF] Westermayr | Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics[END_REF] Other studies have also reported excited-state dynamics using supervised ML schemes. [START_REF] Chen | Deep learning for nonadiabatic excited-state dynamics[END_REF][START_REF] Hu | Inclusion of machine learning kernel ridge regression potential energy surfaces in on-the-fly nonadiabatic molecular dynamics simulation[END_REF][START_REF] Pavlo | Nonadiabatic excited-state dynamics with machine learning[END_REF][START_REF] Westermayr | Machine learning enables long time scale molecular photodynamics simulations[END_REF][START_REF] Westermayr | Neural networks and kernel ridge regression for excited states dynamics of CH 2 NH + 2 : From single-state to multi-state representations and multi-property machine learning models[END_REF] Some of them require an a priori knowledge about the location of conical intersections in order to include more training points in this region [START_REF] Chen | Deep learning for nonadiabatic excited-state dynamics[END_REF][START_REF] Hu | Inclusion of machine learning kernel ridge regression potential energy surfaces in on-the-fly nonadiabatic molecular dynamics simulation[END_REF], while others attempt to include the non-adiabatic couplings in the ML model, [START_REF] Pavlo | Nonadiabatic excited-state dynamics with machine learning[END_REF][START_REF] Westermayr | Machine learning enables long time scale molecular photodynamics simulations[END_REF] which is another challenging and computationally demanding task. Alternatively, direct quantum dynamics can be performed within the Bayesian framework [START_REF] Krems | Bayesian machine learning for quantum molecular dynamics[END_REF] or using the MCTDH propagation on a machine-learned PES [START_REF] Gareth | Direct quantum dynamics using grid-based wave function propagation and machine-learned potential energy surfaces[END_REF]. Recent progress regarding ML applications to electronically excited states has been reviewed by Westermayr et al. [START_REF] Westermayr | Machine learning for electronically excited states of molecules[END_REF] The main motivation of the present work is to investigate the ability of Deep Learning schemes at describing nonadiabatic dynamics in relatively large molecular systems with no a priori knowledge about the topology of the excited-state PESs, based on on-the-fly blind propagation involving all nuclear degrees of freedom and as many excited states as needed. Electronic structure calculations can become computationally extremely demanding for medium and large polyatomic molecules even within the commonly used linear response Time-Dependent Density Functional Theory (TD-DFT) approach. [START_REF] Casida | Time-Dependent Density Functional Response Theory for Molecules[END_REF] This has led us to choose the approximate Density Functional based Tight-Binding (DFTB) [START_REF] Porezag | Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon[END_REF][START_REF] Seifert | Calculations of molecules, clusters, and solids with a simplified lcao-dft-lda scheme[END_REF][START_REF] Elstner | Density functional tight binding[END_REF][START_REF] Elstner | Selfconsistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] scheme that can be applied to large systems, due to the use of parameterized (precomputed) integrals and of a minimal valence basis set. The linear response TD-DFTB approach has been developed in order to access excited states. [START_REF] Niehaus | Tight-binding approach to time-dependent density-functional response theory[END_REF] It has been coupled to nuclear dynamics via Tully's Fewest-Switches Surface Hopping (FSSH) scheme and has been applied to medium-sized systems (namely a few tens of atoms) by several groups. [START_REF] Pal | Nonadiabatic molecular dynamics for thousand atom systems: A tight-binding approach toward pyxaid[END_REF][START_REF] Humeniuk | Dftbaby: A software package for non-adiabatic molecular dynamics simulations based on long-range corrected tight-binding td-dft(b)[END_REF][START_REF] Stojanović | Nonadiabatic dynamics of cycloparaphenylenes with td-dftb surface hopping[END_REF][START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF][START_REF] Uratani | Non-adiabatic molecular dynamics with divide-and-conquer type large-scale excited-state calculations[END_REF]. Recently, we showed that this combined FSSH/TD-DFTB approach satisfactorily describes the relaxation of highly excited electronic states in polyacenes. [START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF][START_REF] Posenitskiy | Theoretical investigation of the electronic relaxation in highly excited chrysene and tetracene: The effect of armchair vs zigzag edge[END_REF] Yet, achieving statistically meaningful sampling of the initial conditions and/or propagating long-lasting trajectories beyond a few hundreds of femtoseconds proves to be numerically involved. Moreover, for large systems beyond a few hundreds or a few thousands atoms, the TD-DFTB scheme may no longer prove computationally feasible. In addition, the computational cost of non-adiabatic coupling calculations scales as powers of the the number of involved states and the system size, respectively. Thus, for large complexes, especially with high density of states, [START_REF] Smith | Hot electron cooling in silicon nanoclusters via landau-zener nonadiabatic molecular dynamics: Size dependence and role of surface termination[END_REF] the generation of training data for non-adiabatic couplings might be overly demanding. Also, ML potentials are not at present well transferable from small to large systems and there may be a need to avoid the explicit use of non-adiabatic couplings in the dynamics, so that consideration of simpler hopping algorithms is of particular interest. Historically, the Landau-Zener (LZ) approximation [START_REF] Zener | Non-adiabatic crossing of energy levels[END_REF]60] has been proposed to evaluate the hopping probability at diabatic state crossings. The LZ formulation has been adapted by Belayev and Lebedev [START_REF] Belyaev | Nonadiabatic nuclear dynamics of atomic collisions based on branching classical trajectories[END_REF] in order to deal with adiabatic states. Alternatively, Zhu and Nakamura [START_REF] Zhu | The two-state linear curve crossing problems revisited. ii. analytical approximations for the stokes constant and scattering matrix: The landau-zener case[END_REF][START_REF] Zhu | The two-state linear curve crossing problems revisited. iii. analytical approximations for stokes constant and scattering matrix: Nonadiabatic tunneling case[END_REF] have proposed an improved Landau-Zener formula for the switching probability. In recent years, multidimensional extensions of the Zhu-Nakamura (ZN) theory have appeared in the literature. [START_REF] Yu | Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: trans-cis photoisomerization in azobenzene[END_REF][START_REF] Hanasaki | An efficient approximate algorithm for nonadiabatic molecular dynamics[END_REF].

The present study is dedicated to the application of Deep Learning through the use of the SchNet architecture, to simplify quantum-classical dynamics in excited electronic states. The goals of the present study are to (i) obtain accurate machine-learned energies and forces for excited states based on TD-DFTB reference data; (ii) perform molecular dynamics and incorporate non-adiabatic effects via a simplified TSH approach without any a priori knowledge about conical intersections or non-adiabatic couplings; (iii) benchmark the simplified TSH scheme against the reference FSSH algorithm on the one hand and available experimental data on the other hand. The paper is organized as follows. In the next section, we briefly outline the basics of the TD-DFTB approach as well as some details regarding the SchNet architecture. Next, we present and discuss the results required for the benchmark study on phenanthrene, first assessing the accuracy of TD-DFTB and SchNet methods and then analyzing the electronic relaxation computed with different hopping schemes. Complementary on-the-fly analysis, i.e. along some selected trajectories, has been performed in order to shed light on the origin of the observed deviations between the simplified TSH and the FSSH approaches. Finally, the conclusions and perspectives are given.

Methods

In this section, we present a brief outline of the TD-DFTB and Deep Learning methods to calculate excited state properties and several TSH approaches to non-adiabatic dynamics, which are considered in this study.

Electronic structure calculations

The Self-Consistent Charge DFTB (SCC-DFTB) was developed by Elstner et al. [START_REF] Elstner | Selfconsistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] as an extension of the original DFTB framework for the ground state. [START_REF] Porezag | Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon[END_REF][START_REF] Seifert | Calculations of molecules, clusters, and solids with a simplified lcao-dft-lda scheme[END_REF] It is based on the second order expansion of the Kohn-Sham DFT total energy around a reference electronic density, so the final expression for the DFTB total energy of a system with M atoms reads

E SCC = Nocc i=1 n i M A=1 µ∈A M B=1 ν∈B c µi H 0 µν c νi + 1 2 M A=1 M B=1 ∆q A γ AB ∆q B + E rep , (1)
where n i is the occupation number of i-th molecular orbital (MO), µ and ν are the Kohn-Sham atomic orbital (AO) indices, ∆q A is the Mulliken charge of atom A, E rep is the atomic repulsive contribution, c µi are the Kohn-Sham MO coefficients and γ AB describes the Coulomb interaction between spherically symmetric charge distributions centered on atoms A and B with a short-range exchange-correlation contribution. The total energy E SCC is further minimized following the self-consistent procedure as proposed in the original work. [START_REF] Elstner | Selfconsistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] The remarkable computational efficiency of the DFTB approach comes from the fact that AO overlap and Hamiltonian matrix elements as well as repulsive potentials can be calculated only once and tabulated for a set of interatomic distances between different pairs of elements. In practice, they are actually precomputed based on some reference electronic structure method (most often DFT) and stored in an external file.

Linear response TD-DFTB was developed by Niehaus et al. [START_REF] Niehaus | Tight-binding approach to time-dependent density-functional response theory[END_REF] as a DFTB analogue of the conventional linear response TD-DFT. [START_REF] Casida | Time-Dependent Density Functional Response Theory for Molecules[END_REF] Excitation energies are given as eigenvalues Ω I of the following matrix equation:

A B B A X Y = Ω I I 0 0 -I X Y , ( 2
)
where I is the identity matrix, A and B are matrices with the elements given by

A ia,jb = (ε a -ε i )δ ij δ ab + 2K Σ ia,jb ; (3)
B ia,jb = 2K Σ ia,jb ; (4)
and indices i, j and a, b denoting the occupied and virtual MOs with energies ε i and ε a , respectively; Σ = S(T ) if singlet (triplet) excited states have to be computed. The coupling matrix elements K S ia,jb for singlet TD-DFTB transitions can be calculated using the generalized Mulliken approximation as follows [START_REF] Niehaus | Tight-binding approach to time-dependent density-functional response theory[END_REF]

K S ia,jb = M A=1 M B=1 q ia A γ AB q jb B , (5)
where q ia A are Mulliken atomic transition charges. It is worth mentioning that the TD-DFTB absorption spectrum neither involves doubly excited states, nor Rydberg states.

In order to perform molecular dynamics simulations coupled to TD-DFTB for electronic structure calculations, excited state energy gradients have to be developed. Their derivation for linear response TD-DFTB relies on the so-called Zvector method, which was initially applied by Furche and Ahlrichs [START_REF] Furche | Adiabatic time-dependent density functional methods for excited state properties[END_REF][START_REF] Furche | Erratum: "adiabatic time-dependent density functional methods for excited state properties[END_REF] to compute analytical forces within the TD-DFT approach. It was further adapted by Heringer et al. [START_REF] Heringer | Analytical excited state forces for the timedependent density-functional tight-binding method[END_REF][START_REF] Heringer | Erratum: "analytical excited state forces for the time-dependent density-functional tight-binding method[END_REF] to compute the gradients of the TD-DFTB excitation energy Ω I . Alternatively, one can follow derivations of TD-DFTB gradients with the long-range correction from Ref. [START_REF] Humeniuk | Dftbaby: A software package for non-adiabatic molecular dynamics simulations based on long-range corrected tight-binding td-dft(b)[END_REF].

Application of Deep Learning to molecular dynamics is essentially based on learning complex multidimensional PES, which has been reported for ground [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] and excited [START_REF] Westermayr | Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics[END_REF] states. We rely on a recently developed Deep Learning architecture called SchNet [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF], which is implemented in the open-source Python package SchNetPack. [START_REF] Schütt | Schnetpack: A deep learning toolbox for atomistic systems[END_REF] There are several reasons why this particular implementation has been chosen. First of all, SchNet is a continuous-filter convolutional neural network, which automatically generates filters that map one hidden layer to another one based on pairwise interactions between a given atom and the surrounding ones. [START_REF] Schütt | Schnet: A continuous-filter convolutional neural network for modeling quantum interactions[END_REF] It has been shown to yield accurate molecular representations and energies, especially when forces are included in the training. [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] Furthermore, it satisfies all the required symmetries and the resulting energies are rotationally and translationally invariant, while forces are rotationally equivariant, thus providing energy-conserving force models. [START_REF] Schütt | Schnetpack: A deep learning toolbox for atomistic systems[END_REF] A shifted softplus is used as an activation function:

ssp(x) = ln e x + 1 2 , (6)
due to an infinite order of continuity, which allows to obtain smooth PES, forces and derivative properties. SchNet is trained by minimizing the cost function J, which is computed as follows

J [E, F 1 , . . . , F M ], [ E, F 1 , . . . , F M ] = ρ||E -E|| 2 + 1 M M A=1 F A -- ∂ E ∂R A 2 , (7)
if both energies E and atomic forces F are included in the training. In the equation above, ρ is the trade-off between the energy and force loss, [START_REF] Pukrittayakamee | Simultaneous fitting of a potential-energy surface and its corresponding force fields using feedforward neural networks[END_REF] quantities with tilde denote the model (SchNet) predictions and plain symbols denote reference data. The second term in Eqn. ( 7) can be neglected if only energies are provided. The initial dataset is split into training, validation and test sets. The training set is used to optimize model parameters in order to minimize the aforementioned cost function. The validation set is used for an early stopping in order to prevent overfitting, which is a common issue in Machine Learning applications when model predictions are very accurate for the training data but significantly less accurate for new inputs. [START_REF] Behler | Constructing high-dimensional neural network potentials: A tutorial review[END_REF][START_REF] Westermayr | Machine learning for electronically excited states of molecules[END_REF] The test set is used to evaluate the final model accuracy. More details about the SchNet architecture and practical aspects can be found in Ref. [START_REF] Schütt | Schnet: A continuous-filter convolutional neural network for modeling quantum interactions[END_REF].

Non-adiabatic molecular dynamics

We now briefly go through the details of the quantum-classical TSH method for incorporating the non-adiabatic effects.

In the TSH approach, the nuclear wavepacket motion is simulated by an ensemble of independent classical trajectories.

Each trajectory evolves on a single electronic state at a given time with a probability to switch (hop) from the current state to another one.

In the Tully's Fewest-Switches Surface Hopping (FSSH) approach, the switching between excited states is controlled by the electron dynamics. Substituting the electronic wavefunction expanded in a basis of adiabatic electronic states into the time-dependent electronic Schrödinger equation, one derives the following equation for the propagation of the complex expansion coefficients C J (t): [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] i dC

J (t) dt = C J (t)E J (t) -i K =J C K (t)D JK (t), (8)
where E J is the adiabatic energy of state J, D JK is the non-adiabatic coupling between states J and K, which is calculated using a finite difference method as follows[73]

D JK (t + ∆t/2) ≈ 1 2∆t [ ψ J (t)|ψ K (t + ∆t) -ψ J (t + ∆t)|ψ K (t) ]. (9)
In the equation above and hereafter, ∆t denotes the nuclear time step that is used to propagate a trajectory. It is important to apply the decoherence correction on C J since the propagation of Eqn. ( 8) in FSSH is overcoherent, which means that electronic coherences C I C * J do not vanish after passing through the region of strong coupling between states I and J. Decoherence corrections have been shown to be crucial in a number of applications, [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF][START_REF] Subotnik | Can we derive tully's surface-hopping algorithm from the semiclassical quantum liouville equation? almost, but only with decoherence[END_REF] and we rely on the commonly used energy-based correction [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF] of the C J coefficients.

PREPRINT

The FSSH probability to switch from the active state I to another state K during the electronic time step ∆τ (used to propagate Eqn. 8) is estimated from the following equation [START_REF] Hammes | Proton transfer in solution: Molecular dynamics with quantum transitions[END_REF]

P FS IK (τ ) = max 0; -2∆τ Re (C * I (τ )C K (τ )) |C I | 2 D IK (τ ) , (10)
where |C I | 2 is the electronic population of a given excited state I. A uniform random number 0 < ξ < 1 is generated at each nuclear time step to determine whether the hop from the active state I to another state K is allowed from the quantum point of view. [START_REF] Hammes | Proton transfer in solution: Molecular dynamics with quantum transitions[END_REF] The hop is accepted if the following condition is fulfilled

K-1 J=1 P FS IJ < ξ K J=1 P FS IJ K = 1; (11)
ξ P FS I1 otherwise. ( 12
)
More details regarding the implementation of the Tully's FSSH scheme can be found in Refs. [START_REF] Humeniuk | Dftbaby: A software package for non-adiabatic molecular dynamics simulations based on long-range corrected tight-binding td-dft(b)[END_REF][START_REF] Stojanović | Nonadiabatic dynamics of cycloparaphenylenes with td-dftb surface hopping[END_REF][START_REF] Pal | Nonadiabatic molecular dynamics for thousand atom systems: A tight-binding approach toward pyxaid[END_REF] for the TD-DFTB approach in general and in Ref. [START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF] for the particular implementation in the deMon-Nano [76] code that has been used in this study.

Even though the FSSH is a commonly used approach to non-adiabatic dynamics, it has a number of limitations. For example, the wavefunctions overlap in Eqn. ( 9) can be computed differently within density-based methods like TD-DFT(B) [START_REF] Werner | Nonadiabatic dynamics within the time dependent density functional theory: Ultrafast photodynamics in pyrazine[END_REF][START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF], decoherence corrections are somewhat ad hoc [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF] and there are numerical instabilities in NACs that can require additional diabatization of Eqn. [START_REF] Pack Kaelbling | Reinforcement learning: A survey[END_REF]. [START_REF] Granucci | Direct semiclassical simulation of photochemical processes with semiempirical wave functions[END_REF][START_REF] Humeniuk | Dftbaby: A software package for non-adiabatic molecular dynamics simulations based on long-range corrected tight-binding td-dft(b)[END_REF] However, it is possible to avoid the aforementioned issues of the FSSH method if there is no electronic equation to propagate and no NACs to compute. This is one of the reasons why results of simplified methods based on the Landau-Zener approximation still appear in the literature. [START_REF] Belyaev | Nonadiabatic nuclear dynamics of atomic collisions based on branching classical trajectories[END_REF][START_REF] Xie | Accuracy of trajectory surface-hopping methods: Test for a twodimensional model of the photodissociation of phenol[END_REF][START_REF] Smith | Hot electron cooling in silicon nanoclusters via landau-zener nonadiabatic molecular dynamics: Size dependence and role of surface termination[END_REF][START_REF] Hanasaki | An efficient approximate algorithm for nonadiabatic molecular dynamics[END_REF][START_REF] Yu | Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: trans-cis photoisomerization in azobenzene[END_REF][START_REF] Chen | Deep learning for nonadiabatic excited-state dynamics[END_REF][START_REF] Hu | Inclusion of machine learning kernel ridge regression potential energy surfaces in on-the-fly nonadiabatic molecular dynamics simulation[END_REF] These schemes are essentially two-state models with one of the following switching probabilities: (i) based on the Belyaev-Lebedev (BL) approach [START_REF] Belyaev | Nonadiabatic nuclear dynamics of atomic collisions based on branching classical trajectories[END_REF] or (ii) based on the Zhu-Nakamura (ZN) theory [START_REF] Zhu | The two-state linear curve crossing problems revisited. ii. analytical approximations for the stokes constant and scattering matrix: The landau-zener case[END_REF][START_REF] Zhu | The two-state linear curve crossing problems revisited. iii. analytical approximations for stokes constant and scattering matrix: Nonadiabatic tunneling case[END_REF]. The BL hopping probability can be computed as follows

P BL IJ = exp - π 2 Z 3 IJ ZIJ , (13)
where Z IJ = |E I -E J | is the adiabatic energy gap between states I and J and ZIJ is its second time derivative at the crossing point t c , i.e. at the local minimum of Z IJ . It is worth mentioning that the equation above is applicable only at the diabatic curve-crossing point.

One of the key differences between the simplified TSH and the FSSH is that the hopping probability P BL IJ is evaluated only at the local minimum of the adiabatic energy gap Z IJ and only between neighboring states, which is a prerequisite of the underlying Landau-Zener approximation. A switch from an active state I to state J is performed if P BL IJ is larger than a random number 0 < ξ < 1. Notably, the aforementioned simplified TSH scheme naturally accounts for decoherence effects. [START_REF] Smith | Hot electron cooling in silicon nanoclusters via landau-zener nonadiabatic molecular dynamics: Size dependence and role of surface termination[END_REF] On the other hand, the electronic population |C I | 2 of excited state S I cannot be computed directly since coefficients C I are generally missing in the Landau-Zener approach. However, the population can be approximated based on the fractional occupation that is a ratio N I /N traj , where N I is the number of trajectories currently (at a given moment of time) running on a PES of S I and N traj is the total number of trajectories.

PREPRINT

To conserve the total energy after hopping, the nuclear velocities are rescaled uniformly by a factor β following the energy conservation law

E I + E kin = E J + β 2 E kin . (14)
Thus, the switch of two states can be still rejected if the energy gap E J -E I > E kin . Such hops are called "frustrated" or classically forbidden hops. Alternatively, one may perform a more accurate velocity rescaling based on energy gradients in order to conserve the total angular momentum. [START_REF] Bonhommeau | Dissociative ionization of neon clusters nen, n=3 to 14: A realistic multisurface dynamical study[END_REF][START_REF] Yu | Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: trans-cis photoisomerization in azobenzene[END_REF] 1

.3 Computational details

Dynamics in low-lying excited states of phenanthrene (C 14 H 10 ) has been chosen as a test case due to both its relatively large molecular size and the diversity of available experimental results, furthermore giving rise to some ambiguity [START_REF] Valerii | Ultrafast nonradiative transitions between higher excited states in organic molecules[END_REF] in their interpretation. Also, phenanthrene is among the most abundant species observed in meteorites. [START_REF] Brechignac | Photophysical studies of jet-cooled pahs: Emission spectra, lifetimes and van der waals clusters astrophysical implications[END_REF] The first part of the discussion in the next section is devoted to the comparison of TD-DFTB and TD-DFT absorption spectra. The TD-DFT absorption spectrum was computed with Gaussian 09 package [83] using BLYP functional and 6-31G(d,p) basis set.

All DFTB calculations presented in this study have been performed with the deMon-Nano [76] code using the mio-1-1 [START_REF] Elstner | Selfconsistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] set of parameters. The ability to reproduce energies and geometries corresponding to low-lying excited states of organic molecules has been critically assessed for different sets of DFTB parameters in Refs. [START_REF] Wanko | A global investigation of excited state surfaces within time-dependent density-functional response theory[END_REF][START_REF] Fihey | Performances of density functional tight-binding methods for describing ground and excited state geometries of organic molecules[END_REF].

Three TD-DFTB datasets have been generated independently for the S 2 , S 3 and S 4 adiabatic excited states of phenanthrene. These particular states have been chosen for the reasons given in the next Section. Each generated dataset contains 10000 points that have been sampled from a single N V T trajectory equilibrated at T = 500 K in the corresponding excited state during 50 ps using a chain of 5 Nosé-Hoover thermostats and ∆t = 0.5 fs. This temperature allows to sample a wider region of the PES, which is desirable due to the absence of any a priori knowledge about conical intersections between the states of interest. Excited state forces have been included in all generated datasets in order to benefit from a relatively low computational cost of reference TD-DFTB calculations and to enhance the training performance.

Next, each dataset is provided to SchNetPack with additional parameters that specify the network topology, i.e. the number of interaction blocks and the number of features, which are specified in Section 2.2. In fact, SchNet consists of a representation network (containing interaction blocks) and a prediction network. [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] We have only changed the topology of the former while the latter was kept fixed to the default pyramidal configuration, namely 3 dense layers with 24 input nodes (equal to the number of atoms) and 1 output node. Additionally, one may activate GPUs and/or include forces in the training. All SchNet models have been trained with a mini-batch stochastic gradient descent, using ADAM optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] with mini-batches of 100 samples, 1000 training epochs and the trade-off ρ = 1 in order to put more emphasis on the energy loss. [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] The trade-off parameter has not been fine-tuned since the default value was found to produce relatively accurate models. The learning rate is initially set to 10 -4 and can decay exponentially with a ratio 0.8 down to 10 -6 with a patience of 25 steps. The maximum number of epochs is set to 5000. The total number of model parameters varies between 238465 and 886529 depending on the SchNet topology for phenanthrene, thereby achieving a computational regime taking full advantage of the use of GPUs.

It is worth mentioning that by default SchNet training and evaluation is performed in single precision for floating point numbers. While this causes no issues during the training step, application of the resulting models to TSH simulations with small time steps might be affected. In particular, the BL approach requires smooth PES since hops are only possible at the local minima. We have performed TSH simulations with SchNet models evaluated in single or double precision. Convergence issues and numerical artefacts associated with the use of single-precision SchNet models are discussed in the Supplementary Information.

PREPRINT

Initial conditions for the TSH simulations were sampled from the thermal distribution of the ground state. A single trajectory was equilibrated at T = 300 K during 50 ps using a chain of 5 Nosé-Hoover thermostats and ∆t = 0.5 fs. Snapshots were taken every 50 fs to be further used as initial conditions.

Coupled FSSH/TD-DFTB simulations have been performed with the deMon-Nano code and serve as reference calculations for the simplified TSH scheme. Each classical trajectory has been propagated using the Velocity Verlet algorithm with ∆t = 0.25 fs during 1 ps and the conical intersection threshold [START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF] set to 10 meV for excited singlet states. Eqn. ( 8) is integrated using a 4-th order Runge-Kutta algorithm with an electronic time step ∆τ = 0.048 as.

In order to be consistent, the same set of initial conditions has been used for Tully's FSSH and for simplified TSH scheme.

The SchNetPack is natively interfaced with the Atomic Simulation Environment (ASE) [START_REF] Hjorth Larsen | The atomic simulation environment-a python library for working with atoms[END_REF] -an open-source set of Python tools for atomistic simulations. This is why the coupled TSH/SchNet simulations have been performed using the molecular dynamics driver of the ASE with the BL switching probability computed on the fly. The ASE driver calls the pre-trained SchNet model to compute adiabatic energies and forces at each nuclear time step, which are further used to run the nuclear dynamics (using Velocity Verlet algorithm) and to detect the local minimum of the energy gap. Due to the fact that the simplified TSH scheme requires data from three consecutive steps (to detect the local minimum of Z IJ and to compute ZIJ with finite differences), one has to actually come back to step t c and make the [t c ] → [t c + ∆t] step again on a new PES if the hop was accepted at t c . The second derivative of the energy gap ZIJ at the crossing point t c is computed using finite differences as follows

ZIJ | t=tc ≈ Z IJ (t c + ∆t) -2Z IJ (t c ) + Z IJ (t c -∆t) ∆t 2 . ( 15
)
We have benchmarked our implementation of the BL scheme for a one-dimensional model problem, namely Tully's simple avoided crossing model from Ref. [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF]. The computed transition probabilities (see Figure S2 in the Supplementary Information) are in good agreement with the work of Hanasaki et al. [START_REF] Hanasaki | An efficient approximate algorithm for nonadiabatic molecular dynamics[END_REF] and the observed deviations are associated with the use of different particle mass.

Finally, experimental results for phenanthrene have been extracted from Refs. [START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF][START_REF] Ohta | Electronic spectra and intramolecular dynamics of phenanthrene in a supersonic free jet[END_REF][START_REF] Brechignac | Photophysical studies of jet-cooled pahs: Emission spectra, lifetimes and van der waals clusters astrophysical implications[END_REF][START_REF] Blanchet | Towards disentangling coupled electronic-vibrational dynamics in ultrafast non-adiabatic processes[END_REF][START_REF] Dick | Determination of femtosecond lifetimes of higher excited singlet states by means of transient photophysical hole-burning. the s2 state of phenanthrene[END_REF]. We consider more specifically the full width at half maximum (FWHM) assessed from early studies performed in supersonic jets [START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF][START_REF] Ohta | Electronic spectra and intramolecular dynamics of phenanthrene in a supersonic free jet[END_REF][START_REF] Brechignac | Photophysical studies of jet-cooled pahs: Emission spectra, lifetimes and van der waals clusters astrophysical implications[END_REF] and the relaxation time from pump-probe experiments of Blanchet et al. [START_REF] Blanchet | Towards disentangling coupled electronic-vibrational dynamics in ultrafast non-adiabatic processes[END_REF] The approximate decay time τ app can be derived from the FWHM as follows [START_REF] Valerii | Ultrafast nonradiative transitions between higher excited states in organic molecules[END_REF]

τ app = 1 2πcΓ , ( 16
)
where Γ is the FWHM in cm -1 and c is the velocity of light.

PREPRINT 2 Results

This section is dedicated to the results of several TSH models applied to the electronic relaxation in phenanthrene. It includes: (i) the validation of the TD-DFTB excited states upon comparison of the computed absorption spectrum versus TD-DFT, CASPT2 and experimental results; (ii) the application of the SchNet model to learn energies and forces based on the TD-DFTB reference datasets; (iii) the analysis of the results of the FSSH/TD-DFTB simulations and their comparison with experimental findings; and finally (iv) the comparison of the simplified TSH simulations based on the Belyaev-Lebedev scheme, with the FSSH/TD-DFTB results, in order to assess the relevance and performance of non-adiabatic dynamics based upon machine-learned PESs with no a priori knowledge about couplings and conical intersections.

Absorption spectra of phenanthrene

First, the positions of several low-lying singlet states of neutral phenanthrene together with the corresponding FWHMs and decay times (where available) are deduced from the various experimental UV spectra and summarized in Table 1. We have computed the absorption spectrum of neutral phenanthrene at both the TD-DFTB and the TD-DFT levels of theory. These results together with selected experimental values are compiled in Figure 1. The agreement between TD-DFTB and TD-DFT spectra is reasonably good and both methods predict positions of two experimentally observed bright states around 284 nm and 266 nm with remarkable accuracy. It is important to note that the band at 284 nm (4.37 eV) corresponds to the S 0 → S 3 transition within both TD-DFT and TD-DFTB, whereas the available experimental studies attribute it to the S 0 → S 2 excitation. However, this might be due to the relatively low oscillator strength (≈10 -3 ) of S 1 compared to S 2 (0.06-0.07). For more details, see Table S1 in Supplementary Information. Moreover, the TD-DFTB ratio κ of oscillator strengths for S 3 /S 2 and S 3 /S 1 amounts to 1.5 and 79.4, respectively, the latter value being in very good agreement with κ = 84 reported by Amirav et al. [START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF] As for TD-DFT, κ for S 3 /S 2 and S 3 /S 1 is 1.4 and 113.7, respectively. The TD-DFTB results can also be compared with the ab initio CASPT2 calculations performed for neutral phenanthrene by González-Luque et al. [START_REF] González-Luque | Theoretical characterization of the absorption spectra of phenanthrene and its radical cation[END_REF] The first bright state in CASPT2 is S 3 with an excitation energy equal to 4.37 eV, which means that the TD-DFTB ordering is actually reasonable, while dark S 1 and S 2 states can be found at 3.42 eV and 4.26 eV, respectively. More recently, Nazari et al. [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF] have performed a detailed analysis of the PESs and of the ultrafast dynamics of monomeric phenanthrene and of some of its derivatives.

PREPRINT

Whereas both the TD-DFT and CASPT2 absorption spectra are in reasonably good agreement with previous studies, the ordering of the states in Ref. [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF] can be somewhat misleading since the higher-lying excited states (S 3 and above) have been discarded. Green and red vertical dashed lines indicate positions of some excited states from the experiments in Ne matrix [START_REF] Salama | Electronic absorption spectroscopy of matrix-isolated polycyclic aromatic hydrocarbon cations. ii. the phenanthrene cation (c14h10+) and its 1-methyl derivative[END_REF] and in the supersonic jet [START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF][START_REF] Ohta | Electronic spectra and intramolecular dynamics of phenanthrene in a supersonic free jet[END_REF], respectively (see Table 1). Inset in the left panel is the balls and sticks representation of phenanthrene (C 14 H 10 ) where black and cyan balls correspond to carbon and hydrogen atoms, respectively. Adapted from Ref. [START_REF] Posenitskiy | Non-adiabatic molecular dynamics of PAH-related complexes[END_REF].

SchNet training

The performance of several SchNet models is summarized in Table 2. Notably, all of them achieve chemical accuracy (errors smaller than 1 kcal/mol or 0.043 eV) on the TD-DFTB energies of phenanthrene. This is not surprising, taking into account the outstanding performance of SchNet for prediction of both ground state [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] and excited state [START_REF] Westermayr | Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics[END_REF] properties. It is worth mentioning that the real accuracy of the machine-learned quantities may vary depending on the underlying electronic structure method that has been used for training. For example, SchNet errors have exceeded 1 kcal/mol for wavefunction-based CASSCF data. [START_REF] Westermayr | Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics[END_REF] However and more generally, errors are larger for the SchNet model #6, trained without forces, compared to model #3, which is consistent with previous studies. [START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] Furthermore, including forces in the SchNet training allows to optimize the model much faster due to a significant amount of information provided by 3M energy gradients. The error does not vary significantly between models with 3 and 6 interactions, which has already been pointed out by Schütt et al. [START_REF] Schütt | Schnet -a deep learning architecture for molecules and materials[END_REF] for the QM9 [START_REF] Reymond | The chemical space project[END_REF][START_REF] Lorenz | 970 million druglike small molecules for virtual screening in the chemical universe database gdb-13[END_REF][START_REF] Ramakrishnan | Quantum chemistry structures and properties of 134 kilo molecules[END_REF] dataset. Yet, the difference in accuracy of SchNets with 1, 2 and 3 interactions can be more pronounced. [START_REF] Schütt | Quantumchemical insights from deep tensor neural networks[END_REF] Notably, using 3000 training and 1000 validation points, model #5 with 256 features achieves similar accuracy as model #4 that has been trained with twice as much data. Overall, 128 features and 3 interactions seem to be a reasonable choice for the SchNet architecture. All the aforementioned conclusions indicate an importance of the hyperparameter tuning for practical applications of neural networks. In fact, thoroughly crafted neural networks can achieve remarkable accuracy even on relatively small training sets. [START_REF] Schütt | Schnet: A continuous-filter convolutional neural network for modeling quantum interactions[END_REF] For all SchNet-based calculations presented below (Sections 2.4 and 2.5), we have used models #1, #3 and #7 for S 2 , S 3 and S 4 excited states, respectively.

FSSH/TD-DFTB simulations

In this subsection we assess the validity of the reference FSSH/TD-DFTB approach in the study of the radiationless electronic relaxation from the S 3 state of phenanthrene, by comparing the computed decay times with the experimental results from Table 1. An ensemble of 250 trajectories has been launched and analysed, which can be viewed in Figure 2. Notably, the shapes of the average population curves change slightly as the number of trajectories is increased from 100 to 250. Also, the computed decay times are 745 fs and 860 fs for 100 and 250 trajectories, respectively. Hence, convergence might not be fully achieved with respect to ensemble size but should at least be good enough. The S 3 population curve can be roughly divided into two parts. The first one is located in the time window [0;400] fs and corresponds to a relatively rapid decay of the initial population from S 3 towards higher-and lower-lying states. It can be approximated by an exponential decay rate of 550-600 fs. The population transfer after 400 fs appears to be slower than the initial one. These results are in qualitative agreement with the complex multichannel relaxation, which is expected for phenanthrene according to Ref. [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF]. As pointed out by Nazari et al. [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF], an internal conversion from the first bright singlet state consists of fast (100 fs) and slow (600 fs) contributions. The same group has also performed FSSH/TD-DFT simulations in the two lowest excited singlet states of a phenanthrene derivative. Their results are in PREPRINT qualitatively good agreement with experimental observations, despite neglecting the transitions towards higher-lying excited states. However, the major fraction (80%) of the population is rapidly (within less than 200 fs) transferred to the lowest excited state, in contrast with our FSSH/TD-DFTB simulations. Furthermore, the TD-DFT energy gaps (computed with ωB97XD and BHandHLYP functionals) [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF] between the two lowest excited states are smaller than the ones in the CASPT2 [START_REF] González-Luque | Theoretical characterization of the absorption spectra of phenanthrene and its radical cation[END_REF][START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF] or TD-DFTB (this work) spectra. This is consistent with faster and more complete (compared to the aforementioned FSSH/TD-DFTB results) relaxation reported by Nazari et al. in Ref. [START_REF] Nazari | Ultrafast dynamics in polycyclic aromatic hydrocarbons: the key case of conical intersections at higher excited states and their role in the photophysics of phenanthrene monomer[END_REF]. Overall, the FSSH/TD-DFTB approach applied to the electronic relaxation of neutral phenanthrene is in reasonably good agreement with the available experimental data. It is worth mentioning that comparison of dynamical simulations with experiments must be carried out with caution. Indeed, experimental access to excited-state dynamics is most often not direct, depending on the experimental probes and setups, involving either transient spectroscopy or time-resolved photo-electron spectroscopy. In the latter case, the probe signal stems from the ionization and realistic comparison with experiment may require additional ingredients in the simulation such as ionization probabilities from excited states. [START_REF] Tao | Ultrafast internal conversion in ethylene. i. the excited state lifetime[END_REF][START_REF] Ruckenbauer | Revealing deactivation pathways hidden in time-resolved photoelectron spectra[END_REF]

TSH/SchNet simulations

We start with the analysis of the TSH simulations with the BL hopping probability coupled to the SchNet-learned quantities (energies and forces) for two propagation steps, such as depicted in Figure 3. Clearly, there are no significant differences in the calculated electronic relaxation. The computed decay time is about 416 and 421 fs with ∆t = 0.25 and ∆t = 0.05 fs, respectively, which means that the TSH dynamics is converged with time step ∆t = 0.25 fs. The population transfer simulated with the BL/SchNet approach is slower than the one observed in the reference FSSH/TD-DFTB calculation, which is related to the fact that the number of BL-induced surface hops is significantly larger compared to that of FSSH. [START_REF] Xie | Accuracy of trajectory surface-hopping methods: Test for a twodimensional model of the photodissociation of phenol[END_REF]. It is worth mentioning another pitfall of the simplified TSH schemes that has been pointed out by Smith and Akimov [START_REF] Smith | Hot electron cooling in silicon nanoclusters via landau-zener nonadiabatic molecular dynamics: Size dependence and role of surface termination[END_REF], namely that BL-based approaches are well suited to study the dynamics in manifolds of states with energy gaps below 0.1-0.2 eV, but may fail for states separated by gaps 0.5 eV. Nonetheless in the present case, the BL/Schnet decay time seems to be in even closer agreement with the four experimental values in the range [0.42-0.52] ps (see Table 1), with respect to the reference FSSH calculation. Yet, a direct comparison with our theoretical results may be only semi-quantitative, depending on the experimental probe mechanism (see discussion in the last paragraph of Section 2.3). Additional TSH calculations have been performed with the TD-DFTB energies and forces instead of the machinelearned ones. The goal was to assess whether the aforementioned results are a direct consequence of using ML-based quantities or is due to the particular TSH implementation presented in this work. Figure 4b shows the occupation dynamics for an ensemble of 100 trajectories propagated using TD-DFTB energies and forces coupled to the BL hopping probability. The results are qualitatively similar to those obtained with the machine-learned quantities (see Figure 4a), yet with more rapid population transfer from S 3 to S 2 .

Discussion

To gain further insight into the performance of the TSH algorithms, we evaluated the dependence of the number of successful surface switches on the instantaneous adiabatic energy gap at the hopping event, which is compiled in Figure 5. The number of hops peaks at small energy gaps of about 0.02 eV and about 0.01 eV for BL and FSSH methods, respectively. Notably, the BL-induced switches are more frequent than the FSSH ones, despite the constraint that switches are only possible at the local minimum of the adiabatic energy gap. This is also consistent with the previous study on excited-state dynamics of a phenol molecule. [START_REF] Xie | Accuracy of trajectory surface-hopping methods: Test for a twodimensional model of the photodissociation of phenol[END_REF] Figure 5: Distribution of the number of successful surface hops computed with the Belayev-Lebedev (blue) and Tully's FSSH (red) approaches versus the corresponding adiabatic energy gap. Data collected from 200 trajectories.

We have performed additional analysis for one selected trajectory (see Figure 6) in order to demonstrate the numerical issues associated with the use of single-precision SchNet models. This trajectory has been propagated in S 3 with ∆t = 0.05 fs, coupled to TD-DFTB or SchNet (evaluated in single or double precision) for electronic structure calculations starting from the same set of initial conditions and with surface switches disabled. The adiabatic energy gap Z 23 between the S 3 and S 2 states and its second time derivative Z23 at each local minimum has been computed on the fly. As expected, the global energy gap evolution is not affected by the change of the SchNet precision. On the other hand, the distribution of Z23 values is very different for trajectory runs based on SchNet. In particular, the number of SchNet-visited local minima (each point in the bottom row corresponds to the finite difference calculation at the detected local minimum) is significantly larger for the single precision case compared to either TD-DFTB or double precision propagation. As can be seen from the inset in Figure 6, SchNet predicts several local minima when using single precision and zero minimum using double precision within the short considered period of time between 1.5 and 2.5 fs. Thus, we conclude that the TSH simulations should be performed with double-precision SchNet (see Supplementary Information for a more detailed discussion). 2). On the other hand, the training has to be done only once and the resulting SchNet model can be used in principal for an arbitrary number of trajectories. Thus, the TSH/SchNet simulations can be at least one order-of-magnitude faster than the FSSH/TD-DFTB ones while avoiding exponential scaling with respect to the number of excited states included in the propagation.

Conclusion

A detailed theoretical study dedicated to the application of Deep Learning to the non-adiabatic molecular dynamics of neutral phenanthrene, based on TD-DFTB calculations and simplified TSH switching probability has been presented. The results of the BL approach coupled to machine-learned PESs have been compared with those from Tully's FSSH/TD-DFTB approach and with available experimental data.

First of all, the accuracy of the TD-DFTB method has been assessed by calculating the absorption spectrum of neutral phenanthrene and upon comparison with higher-level electronic structure and available experimental data. The agreement for the low-lying excited singlet states is reasonably good. In a second step, the electronic relaxation from the bright S 3 state via the cascade of radiationless transitions has been investigated. The detailed analysis reveals good agreement between the simulated decay and the experimental results. Notably, the computed electronic relaxation in phenanthrene is significantly slower than what has been reported in our previous studies on polyacenes [START_REF] Posenitskiy | Nonadiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes[END_REF][START_REF] Posenitskiy | Theoretical investigation of the electronic relaxation in highly excited chrysene and tetracene: The effect of armchair vs zigzag edge[END_REF], which is well correlated with the large energy gap between the initially excited S 3 state and the one below in energy in the TD-DFTB spectrum. This study brings some insight in the photophysics of phenanthrene and has been considered as a reference for the evaluation of simplified approaches to non-adiabatic molecular dynamics. The present results are of interest for astrochemistry and laboratory experiments concerned with atto-or femto-second laser spectroscopy of carbonaceous molecules. Consequently and due to the diversity of available experimental and theoretical results, phenanthrene can be viewed as a benchmark system to either study photophysical processes or test advanced theoretical tools.

The performance of several Deep Learning (SchNet) models trained on the TD-DFTB datasets has been evaluated. Overall, the SchNet architecture has been found able to accurately reproduce the complex multidimensional PESs of phenanthrene, especially when forces are included in the training. Notably, all of them achieve chemical accuracy (errors smaller than 1 kcal/mol or 0.043 eV) on the TD-DFTB energies of phenanthrene. The mean absolute errors could be lowered to less than 0.02 eV for the energy and 0.08 eV/Å for the forces. One should precise the limits of the present conclusions. The real accuracy of the supervised Machine Learning methods is determined by the training data. Thus, the aforementioned SchNet models naturally inherit all advantages and disadvantages of the underlying TD-DFTB approach. The main drawback of some supervised Machine Learning algorithms is that there is no confidence in the predictions for inputs that lie beyond the space of training points. Nevertheless, this should not be an issue in this study since there is no fragmentation or isomerization involved.

Finally, a major focus has been given to the results computed with a simplified TSH scheme on SchNet-learned PESs.

The convergence of the Belyaev-Lebedev approach has been investigated. This TSH scheme coupled to SchNet for electronic structure calculations provides reliable results when SchNet is evaluated in double precision whereas the TSH dynamics propagated with small time steps fails to converge when the default single precision is used. The computed decay time for the S 3 state of phenanthrene is about 0.42 ps that is somewhat shorter than the decay in the reference FSSH calculations (0.55-0.60 ps). Nevertheless, both calculated decay times are in reasonably good agreement with the four experimental values in the range [0.42-0.52] ps (see Table 1). Indeed, a direct comparison with our theoretical results may be only semi-quantitative, depending on the experimental probe mechanism. As far as we know, this study is the first to investigate the electronic relaxation in excited states based on machine-learned quantities derived for a system with more than 20 atoms.

The extreme acceleration provided by Deep Learning will prove invaluable for investigating larger systems with numerous excited states and/or to allow for longer lasting trajectories or a much denser sampling of the initial conditions. There is no doubt that simple hopping schemes can be of valuable interest in parallel with the alternative route consisting in the application of supervised Machine Learning models to fit the non-adiabatic couplings thus enabling the use of the Tully's FSSH approach. [START_REF] Pavlo | Nonadiabatic excited-state dynamics with machine learning[END_REF][START_REF] Westermayr | Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics[END_REF]

Figure 1 :

 1 Figure 1: Absorption spectra (sticks in the left panel, convoluted in the right panel) of neutral phenanthrene computed with TD-DFT (orange) and TD-DFTB (blue) at the equilibrium geometry. All spectra have been truncated to 6.0 eV. Green and red vertical dashed lines indicate positions of some excited states from the experiments in Ne matrix[START_REF] Salama | Electronic absorption spectroscopy of matrix-isolated polycyclic aromatic hydrocarbon cations. ii. the phenanthrene cation (c14h10+) and its 1-methyl derivative[END_REF] and in the supersonic jet[START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF][START_REF] Ohta | Electronic spectra and intramolecular dynamics of phenanthrene in a supersonic free jet[END_REF], respectively (see Table1). Inset in the left panel is the balls and sticks representation of phenanthrene (C 14 H 10 ) where black and cyan balls correspond to carbon and hydrogen atoms, respectively. Adapted from Ref.[START_REF] Posenitskiy | Non-adiabatic molecular dynamics of PAH-related complexes[END_REF].

Figure 2 :

 2 Figure 2: Populations of the low-lying singlet excited states in phenanthrene averaged over an ensemble of 100 (left panel) or 250 (right panel) FSSH/TD-DFTB trajectories. The initial state is S 3 . Reproduced from Ref.[95].

Figure 3 :

 3 Figure 3: Occupations of the low-lying singlet excited states in phenanthrene averaged over an ensemble of 500 trajectories. The initial state is S 3 . Each trajectory has been propagated with an indicated time step ∆t on machinelearned PESs coupled to the Belyaev-Lebedev hopping probability.

Figure 4 :

 4 Figure 4: Occupations of the low-lying singlet excited states in phenanthrene averaged over an ensemble of 100 trajectories. Each trajectory has been propagated with the time step ∆t = 0.25 fs on PES provided by a) SchNet or b) TD-DFTB coupled to the Belyaev-Lebedev hopping probability.

Figure 6 :

 6 Figure 6: Adiabatic energy gap Z 23 between the S 3 and S 2 states (top row) and its second time derivative at the local minima (bottom row) along the selected trajectory propagated with ∆t = 0.05 fs computed with TD-DFTB (red), single-precision (SP) or double-precision (DP) SchNet models (green or blue, respectively). The inset corresponds to adiabatic energy gaps in the time window [1.5:2.5] fs computed with DP (solid blue line) or SP (green dashed line) SchNet models, see main text for discussion.

Table 1 :

 1 Selected results from the experiments with neutral phenanthrene in different environments.

	Environment	Attributed transition	Position, nm	FWHM, cm -1 Decay time, ps
	supersonic jet a	S 0 → S 1 S 0 → S 2	340.9 282.6	5.8 11.1	0.92 0.48
	supersonic jet b	S 0 → S 1 S 0 → S 2	340.9 282.6	-12.6	-0.42
	isopentane glass at 4 K c	S 0 → S 2	-	67±7	0.08
	supersonic jet d	S 0 → S 2	-	12.0	0.44
		S 0 → S 1	341.1	-	-
		S 0 → S 2	284.3	-	-
	Ne matrix at 4 K e	S 0 → S 3 S 0 → S 4	273.4 262.4	--	--
		S 0 → S 5	243.0	-	-
		S 0 → S 6	229.0	-	-
	gas phase f	S 0 → S 1 S 0 → S 2	341.6 283.7	--	-0.52

a Amirav et al.

[START_REF] Amirav | Statistical-limit line broadening in the s2 state of phenanthrene in supersonic jets[END_REF]

b Ohta and Baba

[START_REF] Ohta | Electronic spectra and intramolecular dynamics of phenanthrene in a supersonic free jet[END_REF]

c Dick and Nickel

[START_REF] Dick | Determination of femtosecond lifetimes of higher excited singlet states by means of transient photophysical hole-burning. the s2 state of phenanthrene[END_REF]

d Brechignac and Hermine

[START_REF] Brechignac | Photophysical studies of jet-cooled pahs: Emission spectra, lifetimes and van der waals clusters astrophysical implications[END_REF]

e Salama et al.

[START_REF] Salama | Electronic absorption spectroscopy of matrix-isolated polycyclic aromatic hydrocarbon cations. ii. the phenanthrene cation (c14h10+) and its 1-methyl derivative[END_REF]

f Blanchet et al.

[START_REF] Blanchet | Towards disentangling coupled electronic-vibrational dynamics in ultrafast non-adiabatic processes[END_REF]

Table 2 :

 2 Parameters of the considered SchNet models and the corresponding training times and errors. All models have used forces (except #6) and one NVIDIA®V100 GPU in the training. The MAEs (RMSEs) for energies and forces are given in eV and eV/Å, respectively. The MAE and RMSE have been computed for a test set (2000 structures for all models except #5 and #7 where the remaining 6000 and 3000 structures, respectively, have been used).

			Model parameters	Dataset sizes
	# State Features Interactions Training set Validation set Runtime, h Property	MAE (RMSE)
	1	S 2	128	6	6000	2000	14	energy force	0.014 (0.018) 0.075 (0.126)
	2	S 2	128	3	6000	2000	15	energy force	0.018 (0.023) 0.084 (0.138)
	3	S 3	128	6	6000	2000	13	energy force	0.016 (0.022) 0.068 (0.113)
	4	S 3	256	3	6000	2000	15	energy force	0.014 (0.018) 0.062 (0.106)
	5	S 3	256	3	3000	1000	5	energy force	0.017 (0.022) 0.075 (0.124)
	6	S 3	128	6	6000	2000	40	energy	0.029 (0.038)
								force	-
	7	S 4	128	6	6000	1000	13	energy force	0.023 (0.030) 0.137 (0.239)

Acknowledgements

We acknowledge the European Union (EU) and Horizon 2020 funding awarded under the Marie Skłodowska-Curie action to the EUROPAH consortium, grant number 722346. This work was performed using HPC resources from CALMIP (Grant 2020-P18019). E.P. acknowledges kind hospitality during his stay at the Aarhus University. The authors are grateful to the referees for several valuable comments that contributed to significantly improve the manuscript.

Data availability

All TD-DFTB datasets, pre-trained SchNet models (#1, #3, #7 from Table 2) and initial conditions used in this work are available free of charge on zenodo.org with the following DOI:10.5281/zenodo.4266393. Codes developed in this work are available on github.com/q-posev/ml-tsh.

Conflicts of interest

There are no conflicts to declare.

Supplementary Information

Table S1 contains additional information about several low-lying excited states of neutral phenanthrene extracted from TD-DFTB and TD-DFT absorption spectra. Figure S1 corresponds to the histogram of MAEs (between TD-DFTB energies and SchNet predictions at the same geometry) computed for 100 trajectories. Figure S2 corresponds to the BL benchmark calculations performed for a simple avoided crossing model. Figures S3,S4, S5 are similar to Figures 3,5,6, respectively, but correspond to TSH simulations based on single-precision SchNet models.