Identification of Important Normal Modes in Nonadiabatic Dynamics Simulations by Coherence, Correlation, and Frequency Analyses

Nonadiabatic dynamics simulations of molecular systems with a large number of nuclear degrees of freedom become increasingly feasible, but there is still a need to extract from such simulations a small number of most important modes of nuclear motion, for example to obtain general insight or to construct low-dimensional model potentials for further simulations. Standard techniques for this dimensionality reduction employ statistical methods that identify the modes that account for the largest variance in nuclear positions. However, large-amplitude motion is not necessarily a good proxy for the inﬂuence of a mode on the excited-state wave function evolution. Hence, here we report a number of analysis techniques aimed at extracting from nonadiabatic dynamics simulations the vibrational modes that are most strongly affected by the electronic excitation process and that most signiﬁcantly affect the interaction of the electronic states. The ﬁrst technique identiﬁes coherent nuclear motion after excitation from the ratio between total variance and variance of the average trajectory. The second strategy employs linear regression to ﬁnd normal modes that have a statistically signiﬁcant effect on excitation energies, energy gaps, or wave function overlaps. The third approach uses time-frequency analysis to ﬁnd normal modes where the vibrational frequencies change in the course of the dynamics simulation. All three techniques are applied to the case of surface hopping trajectories of [Re(CO) 3 (Im)(Phen)] + (Im=imidazole; Phen=1,10-phenanthroline), showing that in this transition metal complex the nonadiabatic dynamics is dominated by a small number of carbonyl and phenanthroline in-plane stretch modes.

become increasingly feasible, but there is still a need to extract from such simulations a small number of most important modes of nuclear motion, for example to obtain general insight or to construct low-dimensional model potentials for further simulations. Standard techniques for this dimensionality reduction employ statistical methods that identify the modes that account for the largest variance in nuclear positions. However, large-amplitude motion is not necessarily a good proxy for the influence of a mode on the excited-state wave function evolution. Hence, here we report a number of analysis techniques aimed at extracting from nonadiabatic dynamics simulations the vibrational modes that are most strongly affected by the electronic excitation process and that most significantly affect the interaction of the electronic states. The first technique identifies coherent nuclear motion after excitation from the ratio between total variance and variance of the average trajectory. The second strategy employs linear regression to find normal modes that have a statistically significant effect on excitation energies, energy gaps, or wave function overlaps. The third approach uses time-frequency analysis to find normal modes where the vibrational frequencies change in the course of the dynamics simulation. All three techniques are applied to the case of surface hopping trajectories of [Re(CO) 3 (Im)(Phen)] + (Im=imidazole; Phen=1,10-phenanthroline), showing that in this transition metal complex the nonadiabatic dynamics is dominated by a small number of carbonyl and phenanthroline in-plane stretch modes.
File list (1) download file view on ChemRxiv all.pdf (30.83 MiB)

I. INTRODUCTION
Nonadiabatic dynamics simulations are nowadays an important pillar of computational chemistry, as it allows detailed investigations of many different photophysical and photochemical processes with a wide variety of applications. [1][2][3][4][5][6] The general goal of these simulations is to describe the simultaneous evolution of electronic and nuclear degrees of freedom in a molecular system after photoexcitation. This is a very challenging task, requiring a combination of reliable electronic structure theory and suitable nuclear dynamics techniques. If nuclear motion is described quantum-mechanically (e.g., using standard grid-based quantum dynamics methods 7,8 or multiconfigurational time-dependent Hartree (MCTDH) 9,10 ), then generally the electronic potential energy surfaces (PESs) need to be known in advance. A major bottleneck here for larger molecular systems is that the number of data points needed to describe these PESs scales roughly exponentially with the number of nuclear degrees of freedom, making quantum dynamics of large systems unfeasible. An approximate alternative is given by the ab initio molecular dynamics (AIMD) approach, where nuclear quantum effects are neglected but no a priori PES needs to be known because the electronic structure calculations are carried out on-the-fly at the current nuclear position. Thus, AIMD allows realistic nonadiabatic dynamics simulations in systems with many degrees of freedom.
Even though AIMD enables full dimensionality-a necessary condition for a realistic description of the system-there a) Electronic mail: sebastian.mai@univie.ac.at b) Electronic mail: leticia.gonzalez@univie.ac.at are still reasons why reduced-dimensional models are sometimes desirable. First, in some systems including all degrees of freedom is less important than describing nuclear quantum effects, but this can only be done in reduced models. Second, low-dimensional models that nevertheless capture the essence of the dynamics facilitate the interpretation of the results, in order to gain useful chemical insight. 11 One of the most convenient ways to generate such lowdimensional models is by applying some dimensionalityreduction techniques to full-dimensional AIMD trajectories. A significant body of literature was published on this topic, see, e.g., References 4,[11][12][13][14][15][16][17] One of the most common approaches is to take the coordinate data from the trajectories, compute the covariance matrix, and perform a principal component analysis (PCA). 4,15,17,18 In this way, the coordinate space is transformed into a vector basis where most of the statistical variance is contained in a subspace spanned by the first few vectors. The dimensionality reduction can then be accomplished by simply picking the desired number of vectors and discarding the rest. Related approaches are, e.g., multidimensional scaling (MDS), 16 ISOMAP, 16 or diffusion map, 11,14,19 which all construct different kinds of matrices and find the reduced space via their eigenvectors.
While the mentioned dimensionality reduction approaches have successfully been employed in molecular dynamics simulations from small organic molecules 15 to large proteins, 18 they all have in common that they exclusively consider the distribution of nuclear geometries. Here, we propose that for nonadiabatic dynamics this might not be the only important aspect of the simulations. In most nonadiabatic dynamics simulations, the evolution of the electronic wave function (e.g., the electronic populations) is of primary interest, but the influence of the nuclear degrees of freedom on this evolution is very complex. For example, it might be that some large-amplitude modes exhibit parallel PESs for all electronic states so that these modes do not promote state crossings where nonadiabatic population transfer might occur. On the contrary, there might be modes with fairly small amplitudes, but these modes induce vibronic couplings between states that would not interact without those modes (e.g., symmetry-breaking modes). Hence, the amplitude of a nuclear mode is not necessarily a good proxy for the importance of that mode for nonadiabatic dynamics.
In this contribution, we present a number of statistical analysis techniques that aim at identifying the nuclear degrees of freedom that have the largest importance with respect to the combined evolution of nuclei and electronic wave function. To this end, we investigate which modes show coherent motion or frequency shifts after excitation and which modes correlate the most with changes in energy gaps and wave function overlaps, as these quantities can be directly linked to nonadiabatic effects.
We apply the new analysis techniques to the recently published 20 full-dimensional ab initio excited-state dynamics simulations of [Re(CO) 3 (Im)(Phen)] + (Im=imidazole; Phen=1,10-phenanthroline)-shown in Figure 1-carried out with the surface hopping including arbitrary couplings (SHARC) approach. 21,22 [Re(CO) 3 (Im)(Phen)] + is a member of the class of Re(I) carbonyl diimine complexes that have gained large interest as bio-compatible photosensitizer for studies of electron transfer in proteins. [23][24][25] The complex' nonadiabatic dynamics was experimentally thoroughly characterized, [26][27][28] proposing that the first step after UV excitation is ultrafast intersystem crossing (ISC) from an initially excited singlet metal-to-ligand charge transfer ( 1 MLCT) state to a (triplet) 3 MLCT state, with minor involvement of a triplet intraligand ( 3 IL) state. For the same complex, in the last years a series of quantum dynamics simulations was published, [29][30][31] using MCTDH with linear vibronic coupling (LVC) models 29 with up to fifteen harmonic normal modes. The choice of these normal modes was made manually, based on inspection of the relevant gradient and vibronic coupling parameters employed LVC model. The availability of both a full-dimensional simulation and a low-dimensional simulation for comparison offers a nice opportunity to showcase the efficacy of our presented analysis techniques. Additionally, we present a detailed discussion of the underlying reasons that make the identified normal modes the most important ones.

II. METHODS
This section presents an overview over the different steps in our analysis: preparation of the data coming from the SHARC simulations, 20 transformation to normal mode coordinates, normal mode coherence analysis, normal mode correlation analysis, and frequency shift analysis.

A. Dynamics simulations
The SHARC simulations were recently published elsewhere. 20 Briefly, the SHARC method 21,22 was used to simulate 94 trajectories of [Re(CO) 3 (Im)(Phen)] + in water for 250 fs, including 7 singlet states (S 0 to S 6 ) and 8 triplet states (T 1 to T 8 ). The electronic structure level of theory was TDA-B3LYP with electrostatic embedding QM/MM 32 and a mixed-quality (triple ζ for Re, double-ζ for C, N, O, and H) basis set. The details of the level of theory are also described in Ref. 33.
Initial conditions were generated by first running a 10 ns MD trajectory and taking 500 equidistant snapshots. As described previously, 33 each snapshot was propagated in the ground state at QM/MM level for a random amount between 50 and 100 fs. At the end of these short trajectories, the excited states were computed and the initial excited state was selected.
Below, part of the analyses take into account the full trajectory data, i.e., the 500 ground state trajectories and the 94 excited-state trajectories. The ground state trajectories were shifted in time so that they all end at t = 0 fs, whereas the excited-state trajectories start at t = 0 fs. Since the 94 excitedstate trajectories are continuations of some of the ground state trajectories, this means that for the data analysis we employ 406 ground-state-only trajectories (starting between −50 and −100 fs and ending at 0 fs) and 94 ground + excited-state trajectories (starting between −50 and −100 fs and ending at 250 fs).

B. Normal mode transformation
Even though the dynamics simulations were performed in Cartesian coordinates, for the statistical analysis we decided to work in harmonic normal mode coordinates, with the solvent coordinates neglected. This was mainly done in order to be compatible with the eventual generation of an LVC model based on the set of identified most important normal modes, in order to carry out MCTDH simulations in a way comparable to recent simulations. [29][30][31] Unlike the dimensionality-reduction techniques mentioned above (PCA, ISOMAP, diffusion map), 4,11-16 our goal here is not to find an optimal arbitrary coordinate basis, but rather to identify the most important modes within the employed normal mode basis.
In order to obtain clean normal mode coordinates from the transformation of the Cartesian coordinate data, the first step is to bring all geometries (all steps of ground and excited-state trajectories) into a standard orientation and to center them at the origin, because the normal mode coordinates are not capable of describing translation and rotation. Additionally, the normal mode basis is not suited to describe the free torsion of the imidazole ligand. 32 Hence, the following procedure was carried out with each geometry: (i) The center of mass of [Re(CO) 3 (Im)(Phen)] + is centered on the origin, and the Re(CO) 3 (Phen) moiety (disregarding the im ligand) is optimally overlaid with the reference geometry through the Kabsch algorithm. 34 (ii) The imidazole ligand is separately translated and rotated to overlap optimally with the imidazole of the reference geometry (again through the Kabsch algorithm). The latter step is necessary because an imidazole with very different rotation than in the reference will invalidate most normal mode coordinates. By removing the translation and rotation of the imidazole ligand, we obtain sensible results for most modes (e.g., the bond stretch modes of imidazole, or combined coordination sphere modes), although the information about a few modes (imidazole torsion, Re-N im bond length) will be lost in that way.
The aligned Cartesian coordinates {r α } are then transformed into normal mode coordinates {Q i }. This transformation requires the reference Cartesian coordinates {r ref α }, nuclear masses {M α }, Cartesian-to-normal mode transformation matrix K (with elements K αi ), and normal mode frequencies {ω i }: 35 Note that the resulting Q i are thus mass-and frequency-scaled, dimensionless displacements. The employed parameters for this transformation were taken from a frequency calculation of conformer B of [Re(CO) 3 (Im)(Phen)] + , using B3LYP/TZP and COSMO(water). 32 The geometry is given in Section S1 in the supporting information (SI). The computed normal mode displacements were employed in all statistical analysis described in the following.

C. Normal mode coherence analysis
Since there is no generally accepted definition of what constitutes the "most important" normal modes, we focus on three tasks: (i) identifying normal modes that behave differently in the ground state and excited state, (ii) identifying normal modes that have an effect on energy gaps, and (iii) finding normal modes that promote state crossings in the dynamics. The first of these three tasks we approach with the "normal mode coherence analysis" described in this subsection.
The general idea is based on earlier works by Plasser et al., 36 although we extend their idea to the combined analysis of the ground state and excited state data. The central concept is to compute on one side the total standard deviation over all trajectories and all time steps (total standard deviation) and on the other side the standard deviation of the average trajectory (coherent standard deviation). This allows distinguishing coherent motion (in-phase motion occurring synchronously in a majority of trajectories) from random motion, because the latter cancels out in the average trajectory. These standard deviations are then computed separately for different time windows, as exemplified in Figure 2.
The first statistical measure is the total average in some time window [t min ,t max ] for some set of trajectories: where ∆t is the simulation time step of 0.5 fs. We consider two sets of trajectories: (i) all 500 ground state trajectories (set "all") and (ii) the 94 trajectories that were continued in the excited state (set "exc"). For these two sets, we compute three total averages for each mode i:Q all i (−50, 0),Q exc i (−50, 0), Q exc i (0, 250). In Figure 2, these averages are denoted by black labels.
The second statistical measure is the corresponding standard deviation: Again, we compute three values for each mode: σ all i (−50, 0), σ exc i (−50, 0), and σ exc i (0, 250). These values are indicated in blue in Figure 2. These standard deviations measure the total motion that is occurring in mode i, with a large value indicating that the mode shows a large variance in displacements.
The third quantity is the coherent standard deviation: Note that the coherent standard deviation will never be larger than the total standard deviation.
With the above definitions, we can define some aggregate descriptors for the motion in a normal mode. The coherent ground state motion parameter  indicates whether there are coherent oscillations in the 500 ground state trajectories. Ideally, if the system was perfectly equilibrated at QM/MM level, this parameter should be vanish for all modes. In Figure 2, the value is 0.08 (maximum possible value of 1.0), indicating minor coherent motion in the ground state.
The selected coherence parameter indicates whether there is more coherent motion in the ground state for the "exc" set than for the "all" set. If this parameter is significantly larger than zero, then this is indication that the excitation probability depends on the value of Q i , such that only trajectories with a certain phase are excited. The excited-state coherence parameter indicates whether the normal mode shows stronger coherent oscillations in the excited state than in the ground state. A value of zero means that coherent oscillations are of the same magnitude in both cases. In Figure 2, motion is clearly more coherent than in the ground state, with a value of 0.13. For completeness, we also compute the ground state-excited state shift parameter which indicates whether the normal mode distribution shifts upon excitation. A large (positive or negative) value shows that upon absorption the minimum of the PES is shifted, leading to an excitation of the normal mode. In Figure 2, the value of shiftEX is +0.33, indicating a large shift that is easily visible in the figure. It can be identified as the difference between the yellow line at t = 0 fs and t = 250 fs in the right panel.
We note here that the normal-mode coherence analysis is based entirely on geometry data, neglecting the electronic states mentioned in the introduction. However, the special feature of this analysis is that the effect of the electronic excitation is considered. This allows identifying those normal modes that are important for describing the line shape of the absorption spectrum, as well as those that play a role in nuclear relaxation after excitation.

D. Normal mode correlation analysis
The four presented parameters cohGS, cohSE, cohEX, and shiftEX all serve to identify normal modes that exhibit significant nuclear motion related to the excitation process. However, as suggested in the introduction, the importance of a normal mode does not only depend on whether the mode undergoes significant motion. For example, there could be modes that show large coherent displacements, yet do not promote state crossings and therefore do not play a role in the evolution of the electronic wave function. Alternatively, modes with rather small displacements could affect energy gaps, leading to state crossings between the electronic states (see tasks (ii) and (iii) above).
For a normal mode to affect the electronic dynamics, one criterion would be that relevant state crossings occur along this mode. Such modes can be identified by investigating how the energy differences between the excited states correlate with the normal mode displacements. This is the basis for our second approach, "normal mode correlation analysis". Briefly, we use simple linear regression analysis, with the normal mode displacements Q traj i (t) as independent variables (vector X is the concatenation of Q quantities as dependent variables (vector Y). We then perform linear regression by computing the slope A as and intercept B as where N is the length of vectors X and Y, with a value of 47.000 (94 trajectories times 500 time steps). The coefficient of determination (Pearson correlation coefficient) R 2 is given by We use the R 2 value to judge whether a normal mode significantly affects the respective quantity. An example is shown in Figure 3, correlating the energy gaps ∆E(S 0 , T 1 ) with mode Q 94 . We note that for this analysis, we only employ the trajectory data at t > 0, since the excited states were not computed during the ground state trajectory simulations. We use three different kinds of dependent variables for the correlation analysis. The first kind are the energy gaps ∆E(S 0 , S 1 ) and ∆E(S 0 , T 1 ), taken from the spin-free energies from the SHARC output. The linear regression analysis of this data provides similar information to the coherence analysis above, i.e., which modes behave differently between ground state and excited states.
As a second set of dependent variables, we computed approximate diabatic energy differences. The diabatization procedure employs the state overlap matrices S(t,t + ∆t) containing the overlaps between the states of two consecutive time steps: where Ψ α (t) is an eigenstate of the molecular Coulomb Hamiltonian (MCH) 22 at time t. These overlap matrices were computed along all trajectories for use in the local diabatization method, but here we compute the product of all matrices along the trajectory to get the transformation matrix between the electronic basis at time t and time 0: We then use this transformation matrix to compute the diabatic energies based on the eigenenergies of the MCH: where the diagonal elements of H diab (t) are the sought-after diabatic energies. Note that this gives a different diabatic basis for each trajectory, due to the different initial geometries. We alleviate this issue by performing the linear regression for the energy differences for all pairs of states and taking the maximum R 2 value. This can be justified by the fact that a normal mode already becomes important if it modulates the energy gap between any one pair of states. To avoid that this parameter is redundant to the R 2 values from the energy gaps ∆E(S 0 , S 1 ) and ∆E(S 0 , T 1 ), we omit the S 0 in the set of pairs of states. Unlike the previous descriptors, the R 2 value from the diabatic energy differences does not measure if a normal mode behaves differently in ground and excited state, but rather measures whether a normal mode leads to excitedstate crossings.
As the third set of dependent variables, we use the matrix elements of the overlap matrices S(t,t + ∆t), as they are related to the nonadiabatic couplings and thus provide complementary information about state interactions. Due to its definition, the overlap matrix is a unit matrix if no state mixing occurs. Hence, deviations of the diagonal elements from 1 and deviations of the off-diagonal elements from 0 indicate state mixing. Therefore, we defined as dependent variables for the linear regression ln(1 − S αα ) and ln |S αβ |. These values were linearly fitted in the same way as the energies, and for each mode we obtained the maximum R 2 value among the different matrix elements. This maximum R 2 value measures thus whether the occurrence of state crossings correlates with the motion along the normal modes.

E. Frequency shift analysis
The third analysis that we performed-frequency shift analysis-is again related to task (i), i.e., finding modes that behave significantly different in ground and excited state. In the frequency shift analysis, we monitor the temporal change of vibrational frequency of the normal modes. This change can be obtained from (approximate) time-resolved vibrational spectra, which can be computed by means of a time-frequency analysis of the normal mode coordinates. For these simulations, we only employ the "exc" set of normal mode coordinate data, which was subjected to a wavelet transformation. 37 For the time-frequency analysis, we chose the complexvalued Morlet wavelet transformation. With this transforma- with the wave number ω, the dilution factor f = σ /2πcω∆t, and the mother wavelet: Here, the parameter σ allows choosing a trade-off between time and frequency resolution, with a larger σ giving better frequency resolution and worse time resolution. In the present case, we used σ = 8π. The wavelet transformations were carried out with the signal subpackage of SCIPY 1.0. 38 Using equation (15), we transformed each trajectory individually, and in the end added up the signals incoherently to obtain the total time-frequency power spectrum of normal mode i. From these spectra, we obtained the positions of the maxima at t = −20 fs and t = 150 fs by Gaussian fits. An example for the wavelet transformation, extraction of central frequencies, and obtaining of the frequency shift is given in Figure 4. We note here that the spectral resolution is relatively low (FWHMs of 150-200cm −1 ) due to the short simulation time and the desire to achieve a reasonable high temporal resolution. Despite this limitation, shifts in the maximum positions are visible when comparing the spectra before and after excitation, with a magnitude of +25 cm −1 to +30 cm −1 depending on the mode. Experimentally, the three C=O stretch frequencies shift by +39 cm −1 to +86 cm −1 for [Re(CO) 3 (Im)(Phen)] + , 39 although this shift was obtained after nanosecond delays, i.e., after full solvent relaxation. The instantaneous shift was only reported for the A'(1) mode (the highest-frequency mode of the three), with a value of only +9 cm −1 obtained through extrapolation. Hence, the frequency shifts that we obtain for modes Q 100 to Q 102 are in reasonable agreement with experiment.

A. Normal mode coherence analysis
In Figure 5, we show exemplary coherence analysis plots for four archetypical normal modes. Coherence analysis plots for all important modes (see below) can be found in section S2 in the SI. In panel (a), we show mode Q 41 , which is inactive. The left part of the panel shows full equilibration in the ground state. In the middle part, no coherent motion in the ground state can be seen, although slightly more noise is present than on the left side. Also after excitation (t > 0), no coherent motion can be seen. Consequently, all four parameters (cohGS, cohSE, cohEX, shiftEX) are close to zero.
In panel (b), mode Q 61 displays oscillations unaffected by excitation. The left part shows some some residual oscillations arising from incomplete equilibration during the ground state QM/MM simulations. However, middle and right parts show virtually identical oscillations that are not affected by the excitation selection or excited-state potentials. Thus, only cohGS is large.
Panel (c) shows mode Q 99 , which strongly affects the excitation probability. The left plot shows good equilibration in the ground state. The middle part shows significant oscillations, because only those trajectories were excited that have a negative value of Q 99 at t = 0. Interestingly, after excitation the vibrational motion is quite similar to the ground state, with no large shifts or changes in coherence. This behavior is identified by a large cohSE parameter.
Lastly, panel (d) presents mode Q 78 coherently oscillating after excitation. The left and middle plots show good equilibration in the ground state. After t = 0, all trajectories uniformly move to larger values of Q 78 , showing persistent in-phase oscillations for the entire simulation time. This leads to a large cohEX (and shiftEX) value.
The resulting parameters obtained from the coherence analysis are summarized in Figure 6. In Figure 6, we mark a number of modes by black color to indicate that we did not consider them for the discussion. This is for either of two reasons. (i) The normal mode involves rotation, torsion, or translation of the imidazole ligand; since for the normal mode transformation the imidazole was forced to a specific orientation, the extracted parameters for these modes are not meaningful. (ii) The normal mode has a frequency below 350 cm −1 (applying to Q 7 to Q 26 ); this means that within the 50 fs considered for the ground state statistics this mode has completed less than half an oscillation, making the statistics for this time interval biased and the extracted parameters meaningless. The other normal modes were considered in the analysis. The Figure also indicates the type of normal mode, marking carbonyl modes (stretches, bends) in red and in-plane vibrations of the Phen ligand in yellow (others in blue). These two classes of vibrations were found to be the most relevant for the dynamics of [Re(CO) 3 (Im)(Phen)] + . In panel (a) of Figure 6, we plot the value cohGS for each normal mode, describing the extent of coherent oscillations in the ground state. Here, especially some of the low-frequency modes show large values, while the high-frequency modes tend to have smaller values. These results show that the 50-100 fs simulation time used to equilibrate the ground state to the QM/MM forces is rather short, and longer ground state trajectories would have improved the equilibration of the lowfrequency modes. However, we are primarily interested in changes in the coherent motion between ground and excited state, which can easily be obtained by considering the cohGS parameter as a base line. This is done by substracting cohGS from the other coherence parameters, as defined above.
Panel (b) presents the selected coherence parameter cohSE, which is large if coherent oscillations are larger in the 94 ground state trajectories that were selected for excitation than in all 500 ground state trajectories. A relatively large number of normal modes at all frequencies shows significant values of this parameter, with the largest values obtained for Q 87 , Q 91 , Q 94 , and Q 99 . This indicates that the different normal modes significantly affect the excitation energies and oscillator strengths of the system, such that excitation tends to occur only at a specific position of the normal mode (e.g., at the outer turning point). One can thus assume that a correct description of the absorption spectrum requires the normal modes marked in Figure 6b, most importantly Q 87 , Q 91 , Q 94 , and Q 99 . These results could in principle originate simply from reducing the set of trajectories from 500 to 94, which naturally increases standard deviations through noise and therefore might lead to large cohSE values. However, we note that cohSE is obtained as the quotient of two standard deviations, which is expected to largely cancel this effect. Additionally, we used bootstrapping to verify that the 94 selected trajectories are not simply accidentally more coherent than the full set of 500 trajectories.
Panel (c) of Figure 6 shows the excited-state coherence parameter cohEX. A large value indicates that strong coherent oscillations are present in the excited-state trajectories. A negative value indicates that there is less coherent motion than in the ground state. Unlike in the first two panels, many modes show only very small values of cohEX, indicating that these modes are not significantly excited during the electronic excitation. The most strongly affected modes here are Q 33 , Q 40 , Q 52 , Q 70 , Q 78 , Q 79 , Q 83 , Q 87 , Q 97 , and Q 102 .
The norm of the ground-excited shift parameter |shiftEX| is shown in panel (d). It indicates whether the average value of the normal mode shifts during excitation, indicating that the ground and excited-state minima differ significantly. Modes with large |shiftEX| values will be strongly vibrationally excited during photon absorption. A notable correlation with panel (c) can be observed-modes that have a large shiftEX tend to also have a large cohEX value. The |shiftEX| values additionally identified modes Q 76 , Q 91 , and Q 94 , which have relatively small cohEX values.
B. Normal-mode correlation analysis Figure 7 shows the results of the linear regression analyses, where we tested how much the normal mode displacements are correlated with energy differences or wave function overlaps. These latter parameters control nonadiabatic interactions, hence the correlation analysis provides information about which normal mode contributes to the nonadiabatic population transfer. While Figure 7 only shows the Pearson correlation coefficients from the fits, in section S3 in the SI we present exemplary correlation analysis plots with explanations.
In panel (a), we present the Pearson correlation coefficient between normal mode displacement and the energy difference between S 0 and S 1 . Like the cohEX and shiftEX descriptors, this correlation coefficient provides information about which normal modes show very different forces between ground state and excited state. However, the correlation coefficient data is much less noisy than the cohEX data. Panel (b) of Figure 7 presents similar data, but for the energy difference between S 0 and T 1 . The results are very similar as in panel (a), which 19 33 is probably because S 1 and T 1 are both predominantly MLCT states and hence probably quite parallel on average. Panel (c) shows the maximal correlation coefficients between normal mode and any of the many excited state-to-excited state energy differences (diabatic). Finally, panel (d) gives the correlation coefficients between normal modes and deviations of the overlap matrix elements from unity. This parameter directly measures whether states change their wave function character, providing the most direct measure of nonadiabatic interactions.
Interestingly, all four correlation analyses identify the same set of normal modes, albeit with different weighting. The most prominent normal modes are Q 33 , Q 87 , Q 91 , Q 94 , Q 97 , Q 99 , and Q 102 . Additionally, modes Q 40 , Q 43 , Q 52 , Q 59 , Q 76 , and Q 79 show notable values. By comparing to the cohEX and shiftEX descriptors in Figure 6, we see that the correlation analysis identified several normal modes (Q 40 , Q 43 , Q 59 , Q 99 ) that are apparently very important for promoting state crossings, yet do not undergo a notable change in equilibrium position. The correlation analysis can therefore be used effectively to complement a purely displacement-based identification of important normal modes.  change of frequency indicates that the shape of the PES in the ground and excited states is significantly different, because the electronic excitation affects the involved bonds. Such frequency shifts should be kept in mind when setting up analytical model potentials (like the LVC models mentioned above).
As can be seen, most normal modes show relatively small frequency shifts of less than 10 cm −1 . Then main exception are the normal modes Q 96 to Q 102 . These are the highest-frequency Phen in-plane vibrations as well as the three carbonyl stretch modes. The Phen in-plane vibrations are systematically shifted to lower frequencies, whereas the carbonyl stretch modes are shifted to higher values. For the modes Q 97 , Q 99 , and Q 102 that were above identified as important for the dynamics, the frequency shifts indicate that a reduced-dimensional model might need to take care of describing the change in curvature of the corresponding PESs. For the other modes Q 96 , Q 98 , Q 100 and Q 101 , the other descriptors did not show any strong involvement in the dynamics. Hence, these modes might not be influential enough to warrant their inclusion in a reduceddimensionality model. Figures 6, 7, and 8 present a total of 9 descriptors pertaining to the relative importance of the normal modes in [Re(CO) 3 (Im)(Phen)] + . Out of these descriptors, cohGS and cohSE are mostly presented as baseline values, even though cohSE might be relevant for accurate computations of the spectral width. The frequency shift parameter is mostly relevant for choosing a good functional form of the PESs in analytical models. Regarding the choice of normal modes to include in a reduced-dimensional model, the remaining 6 descriptors (cohEX, shiftEX, R 2 S 1 , R 2 T 1 , R 2 diab , and R 2 ovl ) are of highest importance. Merging the sets of important modes suggested by these 6 descriptors leads to a list of 16 normal modes, which are presented in Table I. A combined importance measure (shiftX/10 plus cohEX/3 plus the four R 2 measures) is presented in Figure 9, together with a graphical depiction of these modes. As can be seen, the two descriptors that identified the largest number of modes as important are shiftEX and the correlation coefficient R 2 diab for the diabatic energy gaps. The other four descriptors did not identify any additional normal modes as important, but confirmed their assignment.

A. Character of important modes
Our different importance measures focus primarily on two aspects: (i) whether nuclear motion in a normal mode changes significantly (and coherently) after excitation, and (ii) whether the normal mode coordinates are correlated to relevant changes in energy gaps or overlaps. Based on our findings, there are two kinds of nuclear motion in [Re(CO) 3 (Im)(Phen)] + that fulfill these criteria. The first kind are the symmetric in-plane vibrations of the Phen ligand, as these most strongly affect the energies of the unoccupied π * orbitals, see Figure 10 for the example of Q 97 (the two lowest blue lines in the left plot). As these π * orbitals receive the excited electron in the MLCT and IL state, changing their orbital energy shifts the energy of these states. This induces many state crossings-like MLCT with IL, but also MLCT with MLCT-explaining the importance of these modes. The second kind of nuclear motion is related to carbonyl motion, as seen in the three modes Q 33 , Q 43 , and Q 102 . These modes are important because they strongly affect the energies of the three highest occupied orbitals, which are a mixture of metal d orbitals and carbonyl π * orbitals, as shown in Figure 10. The modification of these energies in turn affects the energies of the excited states, where the MLCT states are less affected (since they have one electron less in those orbitals) than the other states (ground state, IL). Consequently, excitation of MLCT states during absorption leads to significant forces along the carbonyl modes. Furthermore, these modes induce state crossings of the low-lying 3 IL state with the set of MLCT states. However, the carbonyl modes are less relevant in the population transfer among the different MLCT states.

B. Comparison to previous theoretical work
At this point, it is interesting to compare our set of important normal modes to the set employed in the 15-dimensional MCTDH calculations from Fumanal et al. 29 . Although they did not employ the same normal modes as in our analysis-because Ref. 29   Q 87 , Q 94 , Q 97 , Q 99 within the numbering used in the present manuscript). The 7 modes in the MCTDH model that we did not identify are Q 14 , Q 15 , Q 24 , Q 28 , Q 30 , Q 38 , and the combination of Q 42 and Q 44 . Out of these 7 modes, Q 14 , Q 30 , and Q 42+44 are non-symmetric coupling modes that are crucial in the MCTDH calculations because the model assumes C s symmetry and therefore requires some modes to mediate transitions from A to A states. This is different from our dynamics simulations, where the solvent cage and imidazole torsion breaks symmetry, so that no explicit coupling modes are required. However, as these modes still have approximately symmetric PESs with minima close to Q = 0, they do not show any strong coherent motion or shifts, and consequently are not identified as important by our analyses.
The other 4 modes present in Ref. 29 but not in our results are Q 15 , Q 24 , Q 28 , and Q 38 . Modes Q 15 and Q 24 were excluded from our analysis, as their frequencies are too low to unambiguously identify coherent motion, a shift in average position, or relevant correlations. Additionally, we expect that the harmonic, linear normal modes can not accurately represent these anharmonic modes. For mode Q 28 , we actually obtained some moderate descriptor values (notably in Figure 6c and 7c, d), although those were consistently smaller than for the selected 16 modes. Also for mode Q 38 , we obtained some modest indicators (Figure 6d). The normal modes Q 40 , Q 52 , Q 59 , Q 70 , Q 78 , Q 79 , Q 91 , and Q 102 were not included in the MCTDH simulations of Ref. 29. It might, however, be worthwhile to explore the influence of these modes in future MCTDH simulations.

C. Limitations of Analysis
At this point, a brief discussion about the limitations of our different approaches to identify relevant normal modes for nonadiabatic dynamics is in order. First, our approach is limited by employing normal mode coordinates, based on a frequency calculation of the ground state of [Re(CO) 3 (Im)(Phen)] + . Within this coordinate system, some molecular motion of this molecule cannot be accurately represented, e.g., rotation or translation of the imidazole ligand with respect to the rest of the molecule or solvent motion. However, for other kinds of coordinate system our statistical identification techniques could also be applied.
Another limitation of the analyses is given by the rather short simulation time of the trajectories (50 fs in the ground state, 250 fs in the excited state). The central idea of our coherence analysis is to observe coherent nuclear motion within the different normal modes. However, if the frequency of a normal mode is very small (i.e., its period is very large), then a short simulation time does not allow observing coherent motion. This is the reason why we excluded all low-frequency modes (Q 7 -Q 26 ) from the coherence analysis. In order to employ the coherence analysis for such low-frequency modes, longer trajectories are required. The short trajectories also limit the resolution of the time-frequency transformation.
Furthermore, it is very difficult for our analysis approach to acknowledge normal modes that are (approximately) antisymmetric. While our dynamics simulations were carried out without any explicit symmetry, the effects of local symmetry effects can still be observed. Anti-symmetric modes produce symmetric PESs whose minimum is always located at zero for all states (or close to zero, if symmetry is weakly broken). Hence, anti-symmetric modes neither show a significant shift in average position or coherent motion nor do they show linear correlations in excitation energies, energy gaps, or overlaps. These modes might still play an important role in mediating state-to-state population transfer by acting as coupling modes. 29 Finally, a very difficult aspect of our analysis is to distinguish nuclear motion that is caused by changes in the electronic subsystem from nuclear motion that causes changes in the electronic subsystem. An example would be modes for which the ground state gradient is very different from the excited-state gradients, but all excited-state gradients are identical-such a mode will be strongly oscillating after excitation, but the excited-state populations would be the same with or without inclusion of this mode. Analyzing this kind of causality is difficult for all statistical methods similar to ours.

V. CONCLUSIONS
We report a number of statistical analysis techniques to determine from a set of surface hopping trajectories the normal mode coordinates that are most important for the coupled electronic-nuclear dynamics. Three complementary strategies were employed. In the first one-"normal mode coherence analysis"-we compute for each normal mode the ratio between the total standard deviation and the standard deviation of the average trajectory to identify coherent motion along the normal modes. A comparison of the coherent motion in the ground state and in the excited states then identified modes where motion changed significantly after excitation. The second strategy-"normal mode correlation analysis"-employs linear regression techniques that test which normal modes have a statistically significant effect on excitation energies, energy gaps, and wave function overlaps. This allows identifying modes that will have a clear effect on the evolution of the electronic wave function. The third strategy-"frequency shift analysis"-uses wavelet transformation to obtain a time-resolved vibrational spectrum from which shifts in normal mode frequencies can be computed.
These techniques were applied to the dynamics simulations of the ultrafast ISC dynamics of [Re(CO) 3 (Im)(Phen)] + (Im=imidazole; Phen=1,10-phenanthroline) in solution. By working in linear, harmonic-oscillator normal modes, we neglect solvent motion and the imidazole torsion, but at the same time use a delocalized, uncorrelated set of coordinates that work well for the statistical analyses and can be readily interpreted. The results suggest a set of 16 normal modes that we deem necessary to comprehensively describe the dynamics in [Re(CO) 3 (Im)(Phen)] + . Interestingly, all these modes are either symmetric Phen in-plane deformation modes or stretch/bend modes of the Re-C and C=O bonds. Other modes, e.g., Phen out-of-plane modes, imidazole modes, modes involving hydrogen atoms, seem to play no significant role in the dynamics. We are confident that the presented techniques and results can provide new insight into the topic of identifying important degrees of freedom for nonadiabatic dynamics.

S1 Coordinates
The following coordinates correspond to the optimized (B3LYP/mixed basis) S 0 geometry of conformer B, which was obtained during the calculation of the normal modes used for the NMA.   Table S1 lists all normal modes that were either discussed in Ref. S1 (in their Figure 1 or in the 15-dimensional model of Figure 3d) or found to be important in the analyses of the present work. For ease of comparison, the table lists both the mode numbering for the present work and Ref. S1.

S-13
download file view on ChemRxiv all.pdf (30.83 MiB)