Non-adiabatic excited state molecular dynamics of phenylene ethynylene dendrimer using a multiconfigurational Ehrenfest approach

Photoinduced dynamics of electronic and vibrational unidirectional energy transfer between meta-linked building blocks in a phenylene ethynylene dendrimer is simulated using a multiconfigurational Ehrenfest in time-dependent diabatic basis (MCE-TDDB) method, a new variant of the MCE approach developed by us for dynamics involving multiple electronic states with numerous abrupt crossings. Excited-state energies, gradients and non-adiabatic coupling terms needed for dynamics simulation are calculated on-the-fly using the Collective Electron Oscillator (CEO) approach. A comparative analysis of our results obtained using MCE-TDDB, the conventional Ehrenfest method and the surface-hopping approach with and without decoherence corrections is presented.


I. Introduction
Dendrimers are highly branched conjugated macromolecules that possess well-defined regular structures 1 with numerous peripheral groups, branched repeat units, and a core. Each of these components acts as an individual chromophore unit absorbing light in different characteristic ranges of the spectrum. 2 Their unique geometric and energetic properties 3 have made them the subject of a wide variety of technological 1,4 and academic studies. 5 Dendrimers have many useful properties, such as well-controlled synthesis, 6 exceptional light-harvesting capabilities over a broad region of the solar spectrum, and very efficient energy funneling 7 from the peripheral groups to the core. The efficient and controllable unidirectional energy transfer is associated with the p-conjugation of the regular arrays of coupled chromophore units in the dendritic macromolecule. 5b The time evolution of electronic excitation in organic conjugated materials, such as dendrimers, is determined by non-adiabatic dynamics involving multiple coupled electronic excited states. Following photoexcitation, the multiple photoinduced pathways to electronic and vibrational relaxation involve energy and/or charge transfer, internal conversion, and transition density localization/delocalization. Different theoretical approaches can deal with such photophysical processes.
A number of models, such as Pariser-Parr-Pople, 8 Frenkel exciton, 9 and Frenkel-Holstein, 10 have been used to describe quasi one-dimensional conjugated systems. The real-time Density Matrix Renormalization Group (DMRG) 11 technique is a powerful numerical method for studying these model systems.
In particular, PPE dendritic models can be described using Frenkel exciton Hamiltonians 12 with Coulombic coupling between excitations localized on the linear units. 13 In atomistic molecular dynamics simulations, surface hopping has been proved to be an efficient method to accurately deal with non-adiabatic processes in a large variety of systems. Nevertheless, it presents certain shortcomings. First, electronic decoherence effects are not accounted for in the standard surface hopping algorithm resulting in a number of computational inconsistencies that have been addressed by implementing a large variety of alternative decoherence corrections. 14 The inadequate treatment of decoherence can, in many cases, influence the accuracy of the simulated results. Moreover, there is no clear and direct strategy to improve surface hopping methods in order to introduce nuclear quantum effects such as delocalization, interferences and tunneling.
Among the most accurate fully quantum methods is timedependent Hartree (MCTDH), 15 which has been recently applied to simulate exciton dynamics in extended conjugated polymers. 16 However, this approach, in its original version, requires the parameterization of potential energy surfaces. A useful alternative is represented by a variety of methods relying on trajectoryguided Gaussian basis functions (TBF). The main advantage of these methods is that they can deal with direct or ab initio molecular dynamics 17 where excited state energies, gradients, and non-adiabatic coupling terms are evaluated on the fly, simultaneously with the nuclear propagation. In a variational multi-configuration Gaussian (vMCG) 18 approach, the trajectories for all TBFs are calculated together from the variational principle. vMCG is similar in spirit to MCTDH and it is among the most elegant methods of direct dynamics. However, because of its complexity, on the fly vMCG so far has been implemented only for a small number of basis Gaussians.
Another group of methods 17d,19 applies TBFs moving along classical trajectories. Despite that these methods are still fully quantum because these trajectories are employed only for propagating the basis, while the evolution of their amplitudes and, thus, of the total wave-function is determined by the time-dependent Schrödinger equation. The trajectories here are calculated independently allowing for a large and well-sampled basis. The methods involving branching, such as multiple spawning 17d,19a or multiple cloning 19g algorithms can hardly be used for dendrimers because the large number of electronic states and intersections between them would lead to the uncontrolled growth of the number of trajectories. But, on the other hand, the same makes dendrimers and other large conjugated molecules a nearly ideal object for Ehrenfest based approaches,20 where nuclear motion is defined by the average field, and amplitudes for all electronic states are propagated together with nuclear coordinates. In such systems, a wave packet undergoes frequent transition between many coupled electronic states, and as a result its motions can be frequently well described by the Ehrenfest trajectory even on long time scales.
In recent years we have been developing a group of methods based on multicofigurational Ehrenfest ansatz for the total wave function. 19d-g The central idea is to use Gaussian wavepackets guided by Ehrenfest average of the electronic states. Several variations of this idea have been presented. Our original version of ab initio MCE was using adiabatic electronic states generated by the electronic structure. Here we suggest a modification of the approach, which now has new features: each Gaussian trajectory carries its own electronic basis, which is adiabatic only in the center of the Gaussian. This is similar to the approach described in the appendix in ref. 19g, except that now we accurately take into account the overlaps between electronic states belonging to different trajectories, instead of assuming them to be Kronecker's delta as in ref. 19g. This modification makes the method applicable to large conjugated molecules where electronic states can change significantly on the length-scale of the Gaussian width. In particular, the wave function can change instantly at trivial unavoided crossings between two electronic states localized at remote parts of a large molecule. 21 To distinguish the new approach from the previous work we term it as multiconfigurational Ehrenfest in time-dependent diabatic basis (MCE-TDDB). We would like to emphasize that ''time-dependent diabatic basis'' used here should not be confused with ''diabatic basis'' commonly used in many studies: the trajectories here are still calculated in adiabatic basis, and the new representation affects only matrix elements between different TBFs.
We use the proposed MCE-TDDB method to simulate the excited state dynamics of the system composed of twoand three-ring linear polyphenylene ethynylene (PPE) chromophore units linked through meta-substitution as shown in Fig. 1(a). These chromophore units correspond to the building blocks of more complex phenylethynylperylene-terminated dendrimers, such as the nanostar. 5b,22 The meta-branching localizes the electronic excited states within each linear fragment, and the non-adiabatic coupling between linear fragments are responsible for quantum transitions resulting in the two-ringthree ring unidirectional electronic energy flow. Thus, this system is a good model for analyzing intramolecular electronic and vibrational energy transfer between chromophore units subsequent to the photoexcitation of dendrimer macromolecules.
We run the MCE-TDDB calculations in the basis of a swarm of interacting coherent state trains, 19g,23 which is the most advanced basis for a given number of calculated trajectories. The excited state energies, 13b,24 gradients 25 and non-adiabatic coupling terms 25a,26 needed to run the dynamics are calculated on the fly using the Collective Electron Oscillator (CEO) method. 5a,27 The results are compared with those given by the canonical Ehrenfest approach (EHR) and by the NA-ESMD (non-adiabatic excited-state molecular dynamics) 25a,28 method. The latter technique utilizes surface hopping guided by Tully's fewest-switches (FSSH) algorithm 29 and was specifically developed to describe photoinduced dynamics in large organic conjugated molecules. In particular, the NA-ESMD framework has been extensively applied before to study the intramolecular energy transfer in PPE dendrimers. 30 The paper is organized as follows. In Sections II and III we describe the MCE-TDDB method and computational details of simulations performed on our model system. In Section IV we present and discuss our results. Conclusions are given in Section V. II. Theory

Working equations
Similar to the standard MCE approach, we represent the wavefunction on a trajectory-guided basis |c n (t)i: Here, basis functions |c n (t)i are composed of nuclear and electronic parts: where the electronic part is a superposition of several electronic eigenstates |f (n) I i, and the nuclear part |w n (t)i is a Gaussian wave-packet moving along an Ehrenfest trajectory: Parameter a determines the width of the Gaussians, % R n (t) and % P n (t) are the coordinate and momentum vectors of the n-th basis function center, and g n is a phase.
The basis |f (n) I ifor the electronic part of wave-function can be specified in different ways. In the ab initio multiple spawning method and in the original version of the ab initio MCE approach, the electronic part is represented in the basis of adiabatic eigenstates: |f (n) I i = |f I (r;R)i.
However for large p-conjugated molecules, such as our model system, the trajectories pass through a number of very sharp intersections and trivial unavoided crossings. 21 As a result, the adiabatic representation (4) is not optimal here because |f I (r;R)i can vary significantly at different points R on the length scale of a single nuclear Gaussian basis function. Thus, instead of using adiabatic eigenstates, in the MCE-TDDB approach the electronic part of each basis function is represented in its own diabatic basis composed of eigenstates for the center of this particular Gaussian In this representation, non-adiabatic coupling originates from the time-dependence of the electronic basis |f (n) I i, while for adiabatic basis it is due to the parametric dependence of |f (n) I i on R. It was shown in ref. 19g that both these representations lead to the same set of finale equations when the electronic wave function does not depend too strongly on nuclear coordinates. In this case the adiabatic electronic wave function (4) also needs to be calculated only at the center of the nuclear basis Gaussian, so in practice there is no difference between the two versions (4) and (5) if electronic wave function is smooth and changes slowly. However, this is not a case for the system under consideration, and a more rigorous approach must be developed. We will use electronic wave functions (5) and will take into account that the overlaps between the electronic eigenstates belonging to different Gaussians hf (n) I |f (m) J i can be very far from Kronecker's d IJ , even when these Gaussians are sufficiently close to each other and nuclear parts have a significant overlap.
The electronic overlaps hf (n) I |f (m) J i can be either calculated directly or propagated together with the basis. The latter approach is more convenient in this case, and the following equation for the time-dependence of the overlap integrals is used: Because summation in eqn (6) is limited to only a few lowest electronic states for which non-adiabatic coupling matrix elements d IJ are calculated, this method may slightly overestimate the electronic overlaps. Nevertheless, the accuracy of this approach is compatible with the accuracy of other approximations made in this work. On the other hand, eqn (6) cannot reproduce the state swaps for trivial unavoided crossings 21 because the area of high coupling is extremely localized in this case, much smaller than any reasonable computational step. In order to take trivial unavoided crossings into account, at every time step, we analyze the electronic overlaps calculated directly between electronic wave-functions at times t and t À Dt, and swap the states when such crossing takes place. As before, the motion of the centers of Gaussians is determined by the usual set of Hamilton's equations while phase g n evolves as The force % F n here is an Ehrenfest force that includes both usual gradient terms and additional terms related to non-adiabatic This journal is © the Owner Societies 2016 Phys. Chem. Chem. Phys., 2016, 18, 10028--10040 | 10031 derivative coupling: where V I ( % R n ) is the Ith potential energy surface and d IJ ( % R n ) = hf I |r R |f J i is the non-adiabatic coupling matrix element (NACME).
Also similar to our previous studies, 19f,g the evolution of the Ehrenfest amplitudes a (n) I is determined by the equation Similar to eqn (6) for electronic overlaps, eqn (10) cannot reproduce the amplitude swaps for trivial unavoided crossings, 21 and the amplitudes a (n) I are swapped directly when such crossing is identified. Eqn (7)-(10) form a complete set of expressions determining the evolution of the time-dependent basis |c n (t)i.
The above eqn (7)- (10) describe the evolution of the basis functions (2). Let us now consider the evolution of the whole wave-function (1) in the above-defined basis (2). Substituting eqn (1) into the time-dependent Schrodinger equation and taking into account non-orthogonality of the basis, we get The term hc m (t)|d/dt|c n (t)i in eqn (11), which originates from the time dependence of the basis, can be written as where w m dw n dt and Similar to the case for the potential energy matrix elements above, we write Finally, we substitute eqn (10), (19) and (20) into (17). It is easy to see that the terms of eqn (20) cancel out the off-diagonal terms of eqn (10). Thus, we obtain Eqn (6), (11)-(14), (16), (18) and (21) form a complete set for calculating time evolution of the amplitudes c n (t). This set of equations is similar to those used in ref. 19f; the difference is in overlaps and their time derivatives which include here not only nuclear but also electronic parts. The integrals for the Gaussian functions entering eqn (14) and (18) can be calculated analytically and are readily available, e.g., in ref. 19g.

Population analysis
Population analysis is not trivial for the MCE-TDDB calculations when electronic states change very quickly because the order of electronic states can be different for different basis functions. Let us introduce the adiabatic population operator where, as before, |f (R) K i are electronic eigenfunctions in the point R. Then, the electronic state populations can be expressed as or, simplifying, A similar approach can be applied for calculations of any other electronic properties, such as electronic transition density matrices reflecting spatial extent of the respective many-body electronic wave-functions. If the quantity is described by a quantum operator N, then Using the same approximation as in (24)-(26) and assuming that N is real and depends only on the electronic degrees of freedom, we obtain where are matrix elements of N between eigenstates for the center of mth Gaussian and, thus, are readily available from the electronic structure data for the trajectories.

Basis improvement -coherent state trains
The finite size of the basis set is the most serious limitation of the MCE and related methods, such as MS, vMCG, etc. In order to increase the size of the basis, we apply coherent state train basis sets (TBS), 19g,23 where several Gaussian basis functions are moving along the same Ehrenfest trajectory but with a timeshift Dt. The most expensive part of ab initio molecular dynamics is the electronic structure calculations (energies, gradients, and NACMEs), and the computational cost can be greatly economized by using TBS because all basis functions in the train follow the same trajectory and, thus, repeatedly use the same electronic structure data. This trick allows us to increase the size of the basis by an order of magnitude practically without any additional computational cost. Fig. 2 summarizes the basis sampling ideas which were used in the current work.

NA-ESMD framework
In this work we compare the results of our approach with those of the NA-ESMD 25a,28 method. The NA-ESMD framework is based on the fewest switches surface hopping (FSSH) algorithm, 29 which utilizes the same mixed electronic states as in the Ehrenfest basis functions (2), and where amplitudes a (n) I are also propagated according to the same eqn (10). The difference is that the nuclear motion in the FSSH approach is guided not by the average force but by the gradient of one ''current'' potential energy surface at any given time. The hops from one ''current'' state to another are based on the non-adiabatic coupling strengths and a stochastic switching routine. The probability of a switch at a particular time-step is related to the transition flux between the ''current'' and other electronic states (see eqn (10)): An additional numerical procedure is used in the NA-ESMD algorithm to deal with trivial unavoided crossings, 21 as described in Section II-1.
The FSSH algorithm implemented in NA-ESMD simulations is designed to ensure the internal consistency between the so-called classical populations and the quantum populations. The classical populations are defined as a fraction of trajectories evolving on each potential energy surface: where N I (t) is the number of trajectories with ''current'' state I at any given time, and N traj is the total number of trajectories in the ensemble. The quantum populations are given by the share of each state in the mixed configuration averaged over the ensemble of trajectories: where a I (t) is the quantum amplitude of the Ith electronic state. However, in many cases a significant difference can be observed between P Cl I (t) and P Q I (t). 14a-c This disagreement is mainly caused by classically forbidden transitions (when nuclear kinetic energy is insufficient to conserve the total energy during the hops) and by the divergence of independent trajectories following passage through a region of strong coupling. Despite the fact that P Cl I (t) is the most frequently choice for electronic population analysis in surface hopping methods, in the present work we also analyze P Q I (t) since they are more directly related to the electronic populations calculated using MCE-TDDB and EHR methods.
Empirical corrections can be introduced to account for electronic decoherence. 14 Following extensive numerical tests, 14a we adopt the instantaneous decoherence (ID) approach, performed by resetting the quantum amplitude of the current state to unity after every attempted hop (regardless of whether hops are allowed or forbidden), thus removing the coherence of the quantum coefficients. This approach is based on the assumption that wavepackets traveling on different surfaces should immediately separate in the phase space and evolve independently. It has been previously shown 14a that ID improves FSSH internal consistency between classical and quantum systems while at the same time provides results that do not depend on external parameters and maintains physical relevance.

III. Computational details
We have simulated the photoexcitation and intramolecular energy transfer between m-branched PPE units on the model molecule depicted in Fig. 1(a). Excited state energies, gradients and non-adiabatic couplings are calculated on the fly using the CEO approach. 5a,27 For all simulations presented here, we use the Austin Model 1 (AM1) 31 semiempirical level in combination with the configuration interaction singles (CIS) formalism to describe correlated excited states. This approach has worked well in our previous studies of similar systems. 30 The accuracy of excited-state energies in the AM1/CIS approximation has been benchmarked for a large variety of conjugated molecules, 32 such as PPV 33 and polyfluorenes. 34 In the particular case of our model PPE dendritic molecule (shown at the Fig. 1(a)), vertical excitation energies for the two lowest excited states calculated at the AM1/CIS level for an AM1-optimized structure have been compared with the results obtained using TDDFT (B3LYP/6-31G* model) as well as with the experimental results, 30,35 and a reasonable agreement has been achieved. More details about the benchmark calculations by SCS-AD(2), TDHF, TDDFT methods can be found in ref. 36, where the mechanism of intramolecular energy transfer suggested by NA-ESMD simulations at the AM1/CIS level 30 has been confirmed.
The evaluation of the accuracy of NACME calculations is a more difficult task for any electronic structure method, because direct experimental measurements of NACME are absent. NACME can only be evaluated indirectly from electron-phonon coupling elements and Frank-Condon overlaps. Such an evaluation suggests reasonable accuracy of NACME calculated by AM1/CIS for excited state potential energy surfaces in conjugated hydrocarbons. 33 Moreover, semi-quantitative direct comparisons of non-radiative relaxation rates calculated using our NA-ESMD approach and experimental time-resolved transient spectroscopy data 37 supports the conclusions of ref. 33. Also, both AM1/CIS calculations 38 and TDDFT 36 show that the major contribution to d 12 ( % R n ) comes from the ethynylene stretches, which represent the reactive coordinates of the energy transfer process.
Six singlet electronic states (S 1 -S 6 ) and their corresponding non-adiabatic couplings have been included in the simulations. The excited-state trajectories of 150 fs duration are propagated at constant energy in order to obtain a reasonable statistics. Time steps of 0.5 fs and 0.02 fs have been used in the ground and excited state dynamics, respectively. The Verlet finite difference algorithm 39 is used to integrate the nuclear equations of motion, while the electronic coefficients are propagated using the Runge-Kutta-Verner fifth-and sixth-order method based on the code designed by Hull,Enright,and Jackson. 40 In order to generate the initial conditions for excited-state dynamics, 1 ns of ground state molecular dynamics at 300 K was first performed using a Langevin friction coefficient g of 2.0 ps À1 . Snapshots of nuclei positions and momenta (conformational phase space) have been collected and used as initial conditions for subsequent photoexcitation dynamics modeling. The excited-state trajectories have been started from these initial configurations by instantaneously promoting the system to the state I selected according to a Frank-Condon window: where O laser is the frequency of the laser pulse, O I (R) and f I are the transition energy and normalized oscillator strength of the Ith excited state. The pulse is centered at 348 nm (the maximum of the absorption spectrum for the state S 2 ) and assumed to have a Gaussian shape f (t) p exp(Àt 2 /2T 2 ) with T = 42.5 fs corresponding to a FWHM (Full Width at Half Maximum) of 100 fs. In order to make the comparison consistent, the same set of initial positions, momenta and excited-state populations have been used both for principle trajectories in Ehrenfest dynamics and for NA-ESMD surface hopping. Equations of motions were integrated with the same time step, and at the same level of electronic structure theory.
We performed the MCE-TDDB simulations for 100 trajectory swarms. Each swarm consists of 10 trajectories: the principle one and 9 satellites. This procedure is similar to the one used previously in model calculation, 19d where first the initial wave packet created in the excited state after absorbing a photon was split into ''bits'' and then each ''bit'' was propagated separately on using a swarm of Ehrenfest configurations.
The initial conditions for the principle trajectories were generated as described above. The initial conditions for the satellites were taken 23 very close to those for the principle trajectory with a random shifts DR and DP generated according to the Gaussian distribution where b = 1000 is a compression parameter. At this value of b, the initial overlap between the principle and satellite trajectories in the swarm is about 0.93. Although this overlaps decays as dynamics run, Fig. 3 shows that this decay is sufficiently slow for our model system: at 150 fs time, the average overlap is still about 0.2, ensuring that quantum coupling between the trajectories in the swarm is present throughout our simulation time. The values of parameter a, which determines the width of the basis functions eqn (3), were taken according to the average tested parameters determined for hydrogen and carbon atoms by Martinez et al. 41 For each trajectory in a swarm, we use a coherent state train basis of 11 basis functions. Thus each ''bit'' was propagated on the basis of 10 Â 11 = 110 TBFs in total. A half of additional TBFs in the train are placed before and a half -after the main TBF with a time-shift of 0.6 fs, which provide the overlap of B0.7 between the nearest basis functions in the train. In order to avoid the necessity of backward propagation, the train basis was applied only after 3 fs of dynamics, when there was enough trajectory length to accommodate the tail half of the train; for first 3 fs, calculations were run in the basis of just 10 TBFs.
All the trajectories were calculated independently. Subsequently, the amplitudes c n (t) for all configurations in the basis were propagated using eqn (11) in a ''post-processing '' procedure. 19g This makes the numerical procedure particularly efficient because a very large number of trajectories can be calculated in parallel. Although the total CPU time can reach years, most of it is spent on calculating trajectories using our ab initio on the fly procedure described by eqn (7)-(10) above. Quantum equations for the amplitudes use approximations (16) and (20), which allow to calculate coupling matrix elements in eqn (11) on the basis of the information already available from the calculation of the Ehrenfest trajectories. This ''postprocessing'' technique is the central advantage of the group of methods based on our MCE methodology developed by us. The sampling techniques used in this work were tested on model systems, and it was demonstrated that the MCE can reproduce accurately fully quantum benchmark calculations. 42 Thus, MCE-TDDB and related methods 19f,g represent a very efficient ab initio direct dynamics methodology. The results of ab initio MCE-TDDB calculations were compared with those given by the standard Ehrenfest (EHR) and the NA-ESMD surface hopping methods. The EHR results were obtained from the same trajectories as in the MCE-TDDB, but only the principle trajectory in each swarm was taken into account, and the populations and transition densities for these principle trajectories were simply averaged over the ensemble. The NA-ESMD surface hopping calculations 28,29b,c were performed for a set of 100 excited-state trajectories with the same initial conditions as for principle trajectories in MCE-TDDB calculations. Simulations were also performed using the instantaneous decoherence corrections. More details about these kinds of simulations can be found elsewhere. 14a

IV. Results
We have simulated the ultrafast intramolecular energy transfer that takes place after photoexcitation in the model branched dendritic molecule depicted in Fig. 1(a). This figure also shows the simulated absorption spectrum, in which the two largest contributions with the maxima at 378 nm (3.28 eV) and 348 nm (3.56 eV) come from S 1 and S 2 states corresponding to the excitation localized on two-and three-rings units respectively. There is no direct measurement of the absorption spectrum for this particular molecule. However, our calculated energies for S 1 and S 2 are in reasonable agreement with the values of 3.52 eV and 4.0 eV observed 35 in a larger molecule, the so-called nanostar, for the corresponding electronic states, which are also localized at the same two-and three-ring units.
The initial excited states are populated using a Frank-Condon window defined by a Gaussian shaped laser centered at 348 nm, which corresponds to the maximum of the contribution of the S 2 state. Fig. 4 presents the average populations of the lowest four electronic states as a function of time calculated using the NA-ESMD (a, b), MCE-TDDB (c), EHR (d) methods. For NA-ESMD, both classical P Cl I (t) and quantum P Q I (t) populations are presented. One can see that there is a significant difference between these two populations, which indicates that the so-called internal consistency requirement is violated here.
The dependencies for MCE-TDDB and EHR are very close to each other, except that MCE-TDDB, which relies on swarms and trains basis sets, produces slightly smother curves. On the one hand, this may be an indication that the Ehrenfest approach works really well for large complex systems with a singular photoproduct, such as PPE dendrimers. On the other hand, we know from our previous experience with model systems that using swarm/train basis in combination with trajectory branching (e.g., via cloning procedure 19g ) can further improve the convergence of the results and provide even quantitative agreement with multidimensional benchmark calculations. Cloning could introduce additional decoherence allowing basis functions for different electronic states to run away from each other. Such approaches, however, are very expensive for on the fly use, especially in the case of large molecules with many coupled electronic states, as a rapidly growing number of trajectories would require a huge amount of electronic structure calculations. We do not expect branching to change the result dramatically, but working on implementing even larger basis sets to prove that. Presented here MCE-TDDB technique provides a promising starting point for future developments of methods incorporating fully-quantum vibronic dynamics.
All population dependencies presented in Fig. 4 are similar qualitatively. In particular, one can see a noticeable oscillatory interchange of electronic populations between S 1 and S 2 states during the first 20 fs of dynamics after photoexcitation. The amplitude of the interchange between these states is revealed to be larger for Ehrenfest approaches with respect to surface hopping simulations. In all cases the ultrafast non-radiative decay of the electronic population on S 2 takes place, and by the time of 150 fs the S 1 population becomes significantly higher compared to its initial value. One can see that at the end of our excited-state dynamics, the S 1 population given by the MCE-TDDB and EHR methods lies in between the classical and quantum populations P Cl I and P Q I given by the NA-ESMD. In order to gain insight into the origin of the main differences between surface hopping and Ehrenfest approaches, we have performed the NA-ESMD simulations using decoherence corrections. The importance of incorporating decoherence schemes in the NA-ESMD simulations applied to large extended conjugated molecules have been previously demonstrated. 14a Decoherence has a large effect on the accuracy of the simulated dynamics both in terms of time scales for nonradiative relaxation and in promoting internal consistency of the FSSH algorithm.   5 shows the time-dependence of the average quantum populations of different electronic states for these NA-ESMD simulations using decoherence corrections. We observe that the oscillatory interchange of electronic populations between S 1 and S 2 (see Fig. 4) is damped in these simulations. This is expected because periodic collapses of the electronic wavefunction to a single state eliminate the interference between these two states. After that, the population of the S 1 state increases faster compared to quantum populations given by other approaches.
On the one hand, decoherence is not accounted for in the standard FSSH simulations, leading to internal inconsistencies of the method. 14a,c,n On the other hand, instantaneous decoherence assumes that the divergent wavepackets are immediately decoupled. According to the results shown in Fig. 4 and 5, the Ehrenfest approach data fall in between quantum populations given by the NA-ESMD simulations with and without decoherence corrections. The possible explanation is that the Ehrenfest force depends on the amplitudes of all electronic states taken into account, including sufficiently low populated higher-energy states. Then, continuous population changes due to numerous intersections between a fairly large number of electronic states lead to weak quasi-random fluctuations of Ehrenfest force, which provide natural decoherence. This does not happen in the surface hopping approach where the force is usually determined by a single occupied adiabatic state, and the hops are not frequent and mostly reflect only significant population changes.
The CEO approach 43 calculates the many-body eigenstates by diagonalizing the respective tetradic CIS operator in the molecular orbital (MO) basis. These CIS eigenstates are represented in the form of N occ Â N virt matrices (where occ and virt refer to occupied and virtual MOs), which can be transferred to the original atomic orbital (AO) basis. The resulting quantities are frequently denoted as transition density matrices (or electronic normal modes) and can be formally written as (r I ) mn hf g |c m + c n |f I i, where |f g i and |f I i are the many-electron wave-functions for the ground and Ith excited states, and c m + and c n are the electron creation and annihilation operators with indices n and m referring to the AO basis functions. Notably, owing to a unitary transformation between AO and MO basis and an orthogonality of all AO in the AM1/CIS approach, the normalization condition P n;m r I ð Þ nm 2 ¼ 1 24b,25a is maintained. 44 Matrix elements of the transition density matrices are subject to a straightforward interpretation: the diagonal elements of (r I ) nn reflect 45 light-induced changes in an orbital n due to promotion of the molecule to the Ith excited state defining its transition dipole moment (the expectation value of the dipole operator on the transition density matrix) and, subsequently, an oscillator strength f I . In contrast, the off-diagonal elements of (r I ) nm reflect probability amplitudes of electron transfer between orbitals n and m, describing the optically induced electronic coherences and charge-transfer phenomena. Overall, the transition density matrices directly and compactly reflect the properties of many-body wave-functions even beyond the CIS framework. Subsequently, these are convenient quantities to analyze photoinduced spatial electronic/excitonic delocalization, charge-transfer, etc. 13b In order to obtain the fraction of the transition density localized either on the two-ring or the three-ring linear PPE units, we sum up the atomic contributions belonging to each fragment as where X-ring (X = 2, 3) is either two-or three-ring linear PPE unit, and B-ring refers to the ring in between that is shared by both units. Fig. 6 shows localization of the electronic transition density of the four lowest electronic states at the minimum of the ground state potential energy surface. One can see that S 1 and S 2 states are mainly localized on the three-ring and tworing units, respectively. This means that the dendritic PPE molecule can be represented as an ensemble of weakly coupled linear chromophore units, 30c and the absorption spectra shown in Fig. 1 (b) are approximately a sum of the absorption of threering and two-ring units: 2,27a,b,46 the peak at 376 nm corresponds to the absorption of three-ring unit (S 1 state), and the shoulder around 348 nm is due to the additional contribution of the two-ring unit (S 2 state). Therefore, the excitation at our laser wavelength of 348 nm corresponds to an initial localization of the electronic transition density almost equally distributed between the two-ring and three-ring units. The electronic transition density of higher energy excited-states (S 3 and S 4 ) is qualitatively  different from that of S 1 and S 2 states: it is delocalized over both two-and three-ring fragments for S 3 , and strongly localized for the S 4 state.
Nuclear motion mixes electronic states providing ultrafast energy transfer between the photoactive units. Fig. 7 compares the time-dependencies of the average transition density fractions localized on the two-and three-ring linear PPE units ((r I ) X-ring 2 ) given by MCE-TDDB and NA-ESMD approaches. Note that the NA-ESMD density is calculated by averaging those of the ''current'' electronic state for all trajectories in the ensemble, similar to eqn (31) for classical electronic populations. The timedependencies of transition densities are consistent with our results for the electronic populations (see Fig. 4(a) and (c)): the NA-ESMD approach shows more efficient energy migration toward the three-ring unit, as well as more efficient relaxation from S 2 to the lowest state S 1 associated with this transfer.
It was shown previously that the electronic energy transfer is concomitant to the intramolecular vibrational energy redistribution, and the nuclear motions in the direction of the ethynylene bonds is playing a critical role in the process. 30b Fig. 8 compares the distributions of the two-ring ethynylene bond length in the ground state and at S 2 -S 1 crossing, defined as hop points for surface-hopping trajectories and as points with high NACMEs for Ehrenfest trajectories. Both types of calculations show that, in agreement with previous results, 30b S 2 -S 1 crossing is associated with ethynylene bond stretching in the two-ring part of the molecule.
In Fig. 9, we plot the ethynylene bond lengths for both twoand three-ring units as a function of time. Accompanying the initial excitation of S 2 , this bond of the two-ring unit is initially excited. After that, its excitation relaxes together with the ultrafast decay of the S 2 electronic population (Fig. 4). In contrast, the ethynylene bond of the three-ring unit becomes excited increasing its average amplitude of motion.
The difference in the amplitude of motion for the ethynylene bonds given by different approaches is consistent with our other results. At early times of the dynamics (B20 fs), the amplitude of motion of the ethynylene bond of the two-ring unit is larger for Ehrenfest than that for surface hopping simulation. This time corresponds to the time at which the oscillatory interchange of electronic populations between S 1 and S 2 states takes place (see Fig. 4). At longer times, we observe the reduction of the average vibrational amplitudes. This decay is faster for the Ehrenfest trajectories where the nuclear motion is coupled with amplitudes a (n) I for all electronic states providing better vibrational energy redistribution through continuous partial electronic transitions between many states. Fig. 8a indicates that the ethynylene bond of the two ring unit behaves differently for NA-ESMD and the Ehrenfest approach. The non-adiabatic contribution to the Ehrenfest force % F n (eqn (9)), not present in NA-ESMD simulations, introduces additional vibrational coherence at earlier times of simulations in the regions with large d 12 ( % R n ) values. At later times, the Ehrenfest force dependence on the continuous low population changes among higher-energy states seems to result in quasi-random fluctuations that lead to vibrational decoherence. The same additional vibrational coherence is responsible for higher amplitude of the oscillation of the electronic population, which both Ehrenfest and MCE-TDDB reveal as shown in Fig. 4. Fig. 7 Fraction of the transition density localized in the two-and threering linear PPE units as a function of time obtained from NA-ESMD (solid) and MCE-TDDB (dashed) simulations. Fig. 8 The probability distribution for the length of the ethynylene bond of the two-ring unit (labeled X in the inset) in the ground state (black), and at S 2 -S 1 crossing: the moments of effective hops in surface hopping dynamics (red), and moments of strong non-adiabatic coupling (greater than 0.02 Hartree) in Ehrenfest dynamics (green). This is consistent with the results of ref. 30, where a strong coupling between S 1 and S 2 states was shown to exist only during first 15-20 fs after photoexcitation.

V. Conclusion
In this study we report further development and a new implementation of the ab initio MCE approach suggested previously, to model excited state dynamics of realistic molecular systems beyond the Born-Oppenheimer approximation. The method treats nuclear wave packet dynamics at a fully quantum level, while the electronic structure can be described by a variety of high accuracy ab initio methods. The approach is formally exact and physical observables can be calculated unambiguously without additional assumptions. Previously, the method has been validated only on the model and more simple systems reproducing benchmark results accurately.
The proposed ab initio MCE-TDDB approach is very efficient for large conjugated molecules, such as dendrimers, where the dynamics involving many continuously interacting electronic states, is well described by mean-field Ehrenfest trajectories. Although even larger basis sets with branching trajectories can be introduced, such calculations would be extremely expensive and unlikely to produce significantly different results.
As an application, here we simulate the exited state dynamics for a dendrimer molecule composed of two-and three-ring PPE chromophore units following laser pulse excitation at 348 nm. Our ab initio MCE-TDDB results are compared with those of surface hoping. Our calculations show that for dynamics involving a large number of electronic states, the approaches based on Ehrenfest trajectories naturally account for decoherence phenomena. Thus, our MCE-TDDB calculations validate decoherence corrections, previously introduced to the surface hoping technique in an ad hoc fashion.
Altogether, this makes Ehrenfest-based techniques, such as presented here MCE-TDDB, a method of choice for on the fly simulations of the excitation dynamics in large photoactive molecules.