Identification of important normal modes in nonadiabatic dynamics simulations by coherence, correlation, and frequency analyses

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–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 (nonadiabatic) on-the-fly ab initio molecular dynamics approaches, where nuclear quantum effects are neglected but no a priori PESs need to be known because the electronic structure calculations are carried out when needed at the current nuclear position. Thus, on-the-fly approaches are an important building block of many nonadiabatic dynamics simulations in systems with many degrees of freedom.

Even though on-the-fly techniques enable full dimensionality—a necessary condition for a realistic description of the system—there 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.¹¹ 

One of the most convenient ways to generate such low-dimensional models is by applying some dimensionality-reduction techniques to full-dimensional on-the-fly trajectories. A significant body of literature was published on this topic, see, e.g., Refs. 4 and 11–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),¹⁶ ISOMAP,¹⁶ 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¹⁵ to large proteins,¹⁸ 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²⁰ full-dimensional ab initio excited-state dynamics simulations of [Re(CO)₃(Im)(Phen)]⁺ (Im = imidazole; Phen = 1,10-phenanthroline)—shown in Fig. 1—carried out with the surface hopping including arbitrary couplings (SHARC) approach.^21,22 [Re(CO)₃(Im) (Phen)]⁺ is a member of the class of Re(I) carbonyl diimine complexes that have gained large interest as biocompatible photosensitizer for studies of electron transfer in proteins.^23–25 The complex’ nonadiabatic dynamics was experimentally thoroughly characterized,^26–28 proposing that the first step after UV excitation is ultrafast intersystem crossing (ISC) from an initially excited singlet metal-to-ligand charge transfer (¹MLCT) state to a (triplet) ³MLCT state, with minor involvement of a triplet intraligand (³IL) state. For the same complex, in the last few years a series of quantum dynamics simulations was published,^29–31 using MCTDH with linear vibronic coupling (LVC) models²⁹ with up to 15 harmonic normal modes. The choice of these normal modes was made manually, based on inspection of the relevant gradient and vibronic coupling parameters of the 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.

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

The nonadiabatic dynamics simulations were recently published elsewhere.²⁰ Briefly, the SHARC method^21,22 was used to simulate 94 trajectories of [Re(CO)₃(Im) (Phen)]⁺ in water for 250 fs, including 7 singlet states (S₀ to S₆) and 8 triplet states (T₁ to T₈). The main result was the observation of a two-step ISC process (see Sec. S1 in the supplementary material), with the first step purely electronic and the second step due to nuclear motion. The electronic structure level of theory was TDA-B3LYP with electrostatic embedding quantum mechanics/molecular mechanics (QM/MM)³² 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 based on a classical force field³² and taking 500 equidistant snapshots. As described previously,³³ each snapshot was propagated in the ground state at the QM/MM level for a random amount between 50 and 100 fs. This step served to equilibrate the bond lengths and angles from the less accurate classical force field to the QM/MM values, thereby avoiding that vibrations induced by the electronic excitation are confused with vibrations induced by the switch of forces. At the end of these short trajectories, the excited states were computed and the initial excited state was selected. Specifically, for each excited state within the excitation window of 2.8–3.2 eV, we decided stochastically³⁴ whether it is excited, based on probabilities proportional to its oscillator strength. See Sec. S2 in the supplementary material for a short description of the influence of the selection process.

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 excited-state 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).

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–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.

The solvent coordinates cannot be efficiently represented within a normal mode coordinate system. Besides the strong nonlinearity and anharmonicity of the solvent diffusion degrees of freedom, another problem is that one is usually not interested in the particular motion of each individual solvent molecule, but rather in their statistical arrangement around the solute. Hence, instead of explicit solvent coordinates, a better representation of this arrangement are radial distribution functions (RDFs), as the ones presented for [Re(CO)₃(Im) (Phen)]⁺ in Sec. S3 in the supplementary material. However, a detailed analysis of the solvent dynamics and its influence on the evolution of the solute is beyond the scope of the present manuscript. Hence, we will neglect the solvent coordinates in the following.

In order to transform the Cartesian coordinate data into normal mode coordinates, one has to consider that due to diffusion the molecule might not be in the same orientation as the reference coordinates. Furthermore, the normal mode basis is not suited to describe the free torsion of the imidazole ligand,³² which means that an imidazole in a very different position than in the reference coordinates could lead to strongly distorted normal mode coordinates. Hence, the following procedure was carried out with each geometry prior to the transformation step: (i) the mass of [Re(CO)₃(Im) (Phen)]⁺ is centered on the origin, and the Re(CO)₃(Phen) moiety (disregarding the imidazole ligand) is optimally overlaid with the reference geometry through the Kabsch algorithm.³⁵ (ii) The imidazole ligand is separately translated and rotated to overlap optimally with the imidazole of the reference geometry (again through the Kabsch algorithm). In this way, the relative translation and rotation of imidazole is removed, and we obtain sensible normal mode coordinates for most modes, for example, bond stretches within the imidazole, but also combined coordination sphere modes. The improvement of the normal mode coordinates is presented in Sec. S4 in the supplementary material. The price to pay is that the information about a few modes (e.g., imidazole torsion, Re—N_im bond length) is lost in that way. As described below, these modes are neglected in the analyses.

The aligned Cartesian coordinates {r_A} are then transformed into normal mode coordinates {Q_i}. This transformation requires the reference Cartesian coordinates ${rAref}$⁠, nuclear masses {M_A}, Cartesian-to-normal mode transformation matrix K (with elements K_Ai), and normal mode frequencies {ω_i},³⁶ 

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)₃(Im) (Phen)]⁺, using B3LYP/TZP and COSMO(water).³² The geometry is given in Sec. S5 in the supplementary material. The computed normal mode displacements were employed in all statistical analyses described in the following.

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.,³⁷ 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 Fig. 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¯iall(−50,0)$⁠, $Q¯iexc(−50,0)$⁠, and $Q¯iexc(0,250)$⁠. In Fig. 2, these averages are denoted by red labels.

The second statistical measure is the corresponding total standard deviation σ,

Again, we compute three values for each mode: $σiall(−50,0)$⁠, $σiexc(−50,0)$⁠, and $σiexc(0,250)$⁠. These values are indicated in blue in Fig. 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 ρ,

defined as the standard deviation of the average trajectory. Again, we compute three values: $ρiall(−50,0)$⁠, $ρiexc(−50,0)$⁠, and $ρiexc(0,250)$⁠. These values are indicated in yellow in Fig. 2. Note that the coherent standard deviation will never be larger than the total standard deviation. Also note that, while the total and coherent standard deviations are different, this is not true for their respective averages—the “total average” is identical to the “temporal average of the average trajectory” (or “coherent average”). This can be seen in Fig. 2, where the thick blue and yellow lines exactly coincide.

With the above quantities, 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 the QM/MM level and a sufficient number of trajectories computed, this parameter should vanish for all modes. In Fig. 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 Fig. 2, motion is clearly more coherent in the excited state 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 Fig. 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 thick yellow lines at t = 0 fs and t = 250 fs in the middle and right panels.

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.

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 $Qitraj(t)$ as independent variables [vector X is the concatenation of $Qitraj(t)$ for all trajectories] and different 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² is given by

We use the R² value to judge whether a normal mode significantly affects the respective quantity. An example is shown in Fig. 3, correlating the energy gaps ΔE(S₀, T₁) with mode Q₉₄. While the slope of the linear regression (−0.135) is not large, i.e., the dependency of the energy gap on this normal mode is only weak, the correlation (R² = 0.104) is indeed statistically highly significant, i.e., very unlikely due to random chance, given the large number of data points (47 000) included in the linear regression. 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₀, S₁) and ΔE(S₀, T₁), taken from the spin-free energies from the SHARC output. The linear regression analysis of these data provides similar information to the coherence analysis above, i.e., which modes behave differently between the 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)²² 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² 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² values from the energy gaps ΔE(S₀, S₁) and ΔE(S₀, T₁), we omit the S₀ in the set of pairs of states. Unlike the previous descriptors, the R² value from the diabatic energy differences does not measure if a normal mode behaves differently in ground and excited states, but rather measures whether a normal mode leads to excited-state 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² value among the different matrix elements. This maximum R² value thus measures whether the occurrence of state crossings correlates with the motion along the normal modes.

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 states. 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.³⁸ 

For the time-frequency analysis, we chose the complex-valued Morlet wavelet transformation. With this transformation, the time-frequency distribution of $Qitraj(t)$ is

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.³⁹ 

Using Eq. (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 are given in Fig. 4. We note here that the spectral resolution is relatively low (FWHMs of 150–200 cm⁻¹) due to the short simulation time and the desire to achieve a reasonably 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⁻¹ to +30 cm⁻¹ depending on the mode. Experimentally, the three C=O stretch frequencies shift by +39 cm⁻¹ to +86 cm⁻¹ for [Re(CO)₃(Im) (Phen)]⁺,⁴⁰ 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⁻¹ obtained through extrapolation. Hence, the frequency shifts that we obtain for modes Q₁₀₀ to Q₁₀₂ are in reasonable agreement with experiment.

In Fig. 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 Sec. S6 in the supplementary material. In Fig. 5(a), we show mode Q₄₁, 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, and shiftEX) are close to zero.

In panel (b), mode Q₆₁ displays oscillations unaffected by excitation. The left part shows 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₉₉, which strongly affects the excitation probability. The left plot shows good equilibration in the ground state. The middle part shows significant oscillations, because those trajectories were more likely excited that have a negative value of Q₉₉ 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.

Finally, panel (d) presents mode Q₇₈ 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₇₈, 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 Fig. 6. In Fig. 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⁻¹ (applying to Q₇ to Q₂₆); 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)₃(Im) (Phen)]⁺.

In panel (a) of Fig. 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 low-frequency modes. However, we are primarily interested in changes in the coherent motion between ground and excited states, 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₈₇, Q₉₁, Q₉₄, and Q₉₉. This indicates that the different normal modes notably 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 Fig. 6(b), most importantly Q₈₇, Q₉₁, Q₉₄, and Q₉₉. 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 Fig. 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₃₃, Q₄₀, Q₅₂, Q₇₀, Q₇₈, Q₇₉, Q₈₃, Q₈₇, Q₉₇, and Q₁₀₂.

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₇₆, Q₉₁, and Q₉₄, which have relatively small cohEX values.

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 Fig. 7 only shows the Pearson correlation coefficients from the fits, in Sec. S7 in the supplementary material 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₀ and S₁. 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 are much less noisy than the cohEX data. Panel (b) of Fig. 7 presents similar data, but for the energy difference between S₀ and T₁. The results are very similar to those in panel (a), which is probably because S₁ and T₁ are both predominantly MLCT states and hence probably quite parallel on average. Panel (c) shows the maximal correlation coefficients between the 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₃₃, Q₈₇, Q₉₁, Q₉₄, Q₉₇, Q₉₉, and Q₁₀₂. Additionally, modes Q₄₀, Q₄₃, Q₅₂, Q₅₉, Q₇₆, and Q₇₉ show notable values. By comparing to the cohEX and shiftEX descriptors in Fig. 6, we see that the correlation analysis identified several normal modes (Q₄₀, Q₄₃, Q₅₉, and Q₉₉) 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.

Figure 8 presents the normal mode vibrational frequency shifts induced by the excitation. Exemplary frequency analysis plots with explanations can be found in Sec. S8 in the supplementary material. The vibrational frequencies are related to the curvature of the PESs and hence the strength of the involved bonds. A large 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⁻¹. The main exceptions are the normal modes Q₉₆ to Q₁₀₂. 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₉₇, Q₉₉, and Q₁₀₂ 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₉₆, Q₉₈, Q₁₀₀, and Q₁₀₁, 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 reduced-dimensionality model.

Figures 6–8 present a total of 9 descriptors pertaining to the relative importance of the normal modes in [Re(CO)₃(Im) (Phen)]⁺. Out of these descriptors, cohGS and cohSE are mostly intended 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, $RS12$⁠, $RT12$⁠, $Rdiab2$⁠, and $Rovl2$⁠) 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. As can be seen, the two descriptors that identified the largest number of modes as important are shiftEX and the correlation coefficient $Rdiab2$ for the diabatic energy gaps. The other four descriptors did not identify any additional normal modes as important, but confirmed their assignment.

We also computed a “combined importance measure” as the weighted sum of the 6 descriptors,

The weights are chosen based on the range of values of these descriptors: the R² measures are approximately in the 0–0.1 range, whereas |cohEX| is in the 0–0.3 range and |shiftEX| in the 0–1.2 range; hence, |cohEX| and |shiftEX| are scaled down. The importance measure is presented in Fig. 9, together with a graphical depiction of the most 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 with relevant changes in energy gaps or overlaps. Based on our findings, there are two kinds of nuclear motion in [Re(CO)₃(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 Fig. 10 for the example of Q₉₇ (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₃₃, Q₄₃, and Q₁₀₂. 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 Fig. 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 ³IL state with the set of MLCT states. However, the carbonyl modes are less relevant in the population transfer among the different MLCT states.

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.²⁹ Although they did not employ the same normal modes as in our analysis—because Ref. 29 focuses on only one conformer of [Re(CO)₃(Im) (Phen)]⁺—the two sets of normal modes are mapped to each other in Table S1 in Sec. S9 in the supplementary material.

While our set of important degrees of freedom contains 16 normal modes, Ref. 29 employed 15 modes, out of which our analysis identified 8 as important (Q₃₃, Q₄₃, Q₇₆, Q₈₃, Q₈₇, Q₉₄, Q₉₇, and Q₉₉ within the numbering used in the present manuscript). The 7 modes in the MCTDH model that we did not identify are Q₁₄, Q₁₅, Q₂₄, Q₂₈, Q₃₀, Q₃₈, and the combination of Q₄₂ and Q₄₄. Out of these 7 modes, Q₁₄, Q₃₀, and Q₄₂₊₄₄ are nonsymmetric 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 break 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₁₅, Q₂₄, Q₂₈, and Q₃₈. Modes Q₁₅ and Q₂₄ 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 cannot accurately represent these anharmonic modes. For mode Q₂₈, we actually obtained some moderate descriptor values [notably in Figs. 6(c), 7(c), and 7(d)], although those were consistently smaller than for the selected 16 modes. Also for mode Q₃₈, we obtained some modest indicators [Fig. 6(d)].

The normal modes Q₄₀, Q₅₂, Q₅₉, Q₇₀, Q₇₈, Q₇₉, Q₉₁, and Q₁₀₂ 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.

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 as described here is limited by employing normal mode coordinates, based on a frequency calculation of the ground state of [Re(CO)₃(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. More generally, also molecular processes like isomerization or dissociation could not be represented in a normal mode coordinate system. However, in principle, the analyses do not require normal mode coordinates and any kind of coordinate system could be used as basis, since the analyses only require a coordinate vector for each trajectory time step. Nevertheless, we anticipate that not all coordinate systems would perform equally well. For example, the correlation analysis would probably not work properly if the coordinate basis is nonorthogonal. Additionally, the number of important modes might differ dramatically with the employed coordinate basis, as some bases are more suitable to compactly represent the motion of a molecule.

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₇–Q₂₆) 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 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 can still be observed. Approximately antisymmetric modes produce symmetric PESs, whose minima are always located at zero for all states (or close to zero, if symmetry is only weakly broken). Hence, antisymmetric 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.²⁹ 

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. What probably can be said about the causality in the dynamics of [Re(CO)₃(Im) (Phen)]⁺ is that the electronic state clearly causes changes in the nuclear motion, as evidenced by the normal modes that show coherent oscillations after electronic excitation. The inverse—nuclear motion affecting the electronic wave function—is more difficult to see, but there are at least two manifestations of it. First, the frozen-nuclei simulations in Ref. 20 showed that the populations of the spin-mixed states do only change when including nuclear motion; this is important because without changes in the electronic populations there is no nonadiabatic relaxation. Second, nuclear motion is needed to correctly describe the rise and decay of the triplet intraligand contribution to the wave function²⁰ that is found experimentally.²⁸ 

As a final remark, we will discuss how the presented analysis techniques could be applied to other nonadiabatic dynamics methods that are not surface hopping. The coherence analysis requires that for each normal mode a time-dependent density/distribution can be obtained. Hence, this analysis can be used easily with both trajectory-based methods (e.g., surface hopping, Ehrenfest dynamics,⁴¹ or multiple spawning dynamics⁴²) and with wave packet dynamics methods (variational multiconfigurational Gaussians,⁴³ multiconfigurational time-dependent Hartree,¹⁰ grid-based methods⁸). The correlation analysis does not even require explicit dynamics simulations to be available, since it is oblivious to the time evolution of the system and only investigates correlations between nuclear coordinates and energies/overlaps. Hence, a potential energy surface data base, with data coming from any on-the-fly dynamics method (surface hopping, Ehrenfest, or multiple spawning, direct-dynamics variational multiconfigurational Gaussians), would be all that is required. The frequency shift analysis in the way presented here—where trajectories are individually wavelet-transformed and then incoherently added—requires a trajectory-based method. However, one could alternatively wavelet-transform the average position along the modes to achieve essentially the same result. The latter approach would also be compatible with wave packet dynamics methods.

We have reported 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)₃(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)₃(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, and 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.

See the supplementary material for evolution of electronic populations; discussion of initial normal mode distribution; radial distribution functions; discussion of normal mode transformation; Cartesian coordinates of reference geometry; plots for all important modes for coherence, correlation, and frequency shift analysis; and comparison to literature normal mode selection.

The authors thank the Austrian Science Fund (FWF) through Project No. I2883 (Denetheor), the University of Vienna for financial support, and the Vienna Scientific Cluster (VSC) for very generous allocation of computational time. We thank A. Monari, C. Daniel, M. Fumanal, and E. Gindensperger for all the valuable discussions about the dynamics of the [Re(CO)₃(Im) (Phen)]⁺ complex.