# Automated Selection of Nuclear Coordinates for Reduced Dimensionality Nonadiabatic Dynamics

Vincent Delmas, Alessandro Nicola Nardi, Isabella C. D. Merritt, Anthony Ferté, Ignacio Fdez. Galván, and Morgane Vacher*

<figure><img src="image_1.png" /></figure>

Cite This: J. Chem. Theory Comput. 2025, 21, 6611–6621

<figure><img src="image_2.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_3.png" /></figure>

Metrics & More

<figure><img src="image_4.png" /></figure>

Article Recommendations

<figure><img src="image_5.png" /></figure>

Supporting Information

**ABSTRACT:** Poor scaling of dynamics simulations with number of dimensions is currently a major limiting factor in the simulation of photochemical processes. In this work, we investigate ways to reduce the dimensionality of many-atom systems with a view toward enhancing computational efficiency while maintaining accuracy. Using mixed quantum-classical Trajectory Surface Hopping (TSH) simulations of three photoreactive molecules—trans-azomethane (tAZM), butyrolactone (Bulac), and furanone (Fur)—we explore two different dimensionality reduction techniques: Principal Component Analysis (PCA) and Normal Mode Variance (NMV). Dynamics simulations are run in full dimensionality and reduced dimensionality, employing either PCA or NMV, and the impact of the dimensionality reduction on selected electronic and geometric properties of the dynamics is evaluated. For all three molecules, both PCA and NMV can be used to select lower-dimensional spaces in which the full-dimensionality dynamics results are reproduced. PCA reduction outperforms NMV in all systems, allowing for a more significant dimensionality reduction without loss of accuracy. The improved accuracy of PCA is, for tAZM, mostly seen in the electronic properties while for both Fur and Bulac the advantage is clear in the ring-opening reaction itself as well. The present approach opens routes to simulation of larger photochemically relevant systems, through the use of automated dimensionality reduction, avoiding human bias.

<figure><img src="image_6.png" /></figure>

## 1. INTRODUCTION

Modeling excited-state dynamics allows us to understand how fundamental processes such as photochemical reactions, energy transfer, and transitions between electronic states occur. This is key for the understanding of natural processes and development of new technologies based on light-absorption, a field of significant current interest e.g., for green synthesis and optoelectronic devices.<sup>1,2</sup> Advances in computational chemistry have enabled detailed exploration of small molecular systems with remarkable precision.<sup>3–5</sup> However, accurate modeling of excited-state dynamics processes often proves challenging due to the inherent increase in computational cost as systems increase in size, with the dimensionality of the conformation space for a (nonlinear) molecule being $D = 3N_{at} - 6$, where $N_{at}$ is the number of atoms. The cost of dynamics simulations generally scales very poorly with dimensionality, with the severity of the scaling depending on the method used: in grid-based molecular quantum dynamics simulations, exponential scaling on $D$ is often labeled as “the curse of dimensionality”.<sup>6,7</sup> Practically, this limits full-dimensional dynamics studies to small molecules, with a different practical upper limit for $N_{at}$

depending on the excited state dynamics method employed along with any additional approximations within.

One approach to address the limitations imposed by increasing molecular size, often employed for fully quantum excited state dynamics methods such as Multi-Configuration Time-Dependent Hartree (MCTDH),<sup>8</sup> is to instead run the dynamics on a model of the molecule consisting of only a selected subset of dimensions, expected to be representative of the essential features of the dynamics. Determining the optimal subset of dimensions for such simulations is a critical aspect of model construction in MCTDH and other dynamics methodologies, with a compromise being necessary between accuracy and cost. Chemical knowledge of the system and process to be modeled, along with symmetry considerations,<sup>9</sup> often guides the selection of relevant degrees of freedom to be included in the model. Unfortunately, this lack of a systematic and

Received: January 21, 2025

Revised: June 11, 2025

Accepted: June 12, 2025

Published: June 24, 2025

<figure><img src="image_7.png" /></figure>
impartial approach implies a natural bias in the dynamics simulated, potentially limiting the accuracy of such models.

Prior literature has explored the idea of identifying the most important nuclear dimensions in a systematic manner, with the aid of lower-level dynamics such as the mixed-quantum-classical trajectory surface hopping (TSH) method.¹⁰ For instance, Gómez et al.¹¹ introduced an iterative approach using a feedback loop between MCTDH and linear vibronic coupling TSH,¹² to identify the minimum number of nuclear dimensions and electronic states to produce accurate dynamics. In their study, the least important nuclear modes are identified according to the amplitudes of key associated parameters (namely, the inter- and intrastate coupling to other modes). However, they note that this selection process could be carried out in a number of alternative ways, such as by using excited-state properties in their model, or identification of unimportant modes from the dynamics simulations.

In the present work, we evaluate the use of two different post-simulation analysis methods, namely, Principal Component Analysis (PCA)¹³,¹⁴ and Normal Mode Variance (NMV), to automatically and systematically reduce dimensionality. PCA and other linear and nonlinear methods have already been used quite extensively⁶,¹¹,¹⁵⁻¹⁹ in the field of molecular dynamics. However, this has been primarily with the aim of understanding multidimensional data created by dynamics simulations—i.e., for post-processing of high-dimensional data. In contrast, we propose to use PCA and NMV to extract meaningful information from dynamics performed in full dimensionality using for instance the mixed quantum-classical TSH method²⁰,²¹ and to identify the dimensions which are the most important when describing the reaction, based on this analysis. Selecting a relevant reduced dimensional space through such a protocol could be useful when aiming to run future simulations modeling photochemical reactions with higher accuracy (but higher cost) methods, for example using an improved electronic structure level, or in fully quantum dynamics simulations. Importantly, the suggested protocol assumes that the lower-level dynamics is able to capture the relevant features of the investigated process. Previous theoretical studies have benchmarked the surface hopping method against more accurate dynamics methods and demonstrated its suitability for simulating photochemical reactions,²²⁻²⁴ although it may be less suited for complex atto-photochemical reactions.²⁵ In the present work, a given combination of dynamics method (TSH) and electronic structure level (see Section 2.3) is used throughout. The dynamics simulated in the reduced dimensional space is compared to the full-dimensional results, to evaluate the quality of PCA/NMV as a dimensionality reduction technique.

The dimensionality reduction protocol is tested on the three photochemical reactions illustrated in Figure 1, selected to cover a range of photochemical situations. The first reaction (Figure 1(a)) is the trans-to-cis photoisomerization of azomethane (tAZM), a prototypical double bond rotational isomerization which occurs on a 400 fs time scale according to prior simulations.²⁶ The second (Figure 1(b)) and third (Figure 1(c)) reactions are examples of bond breaking occurring during the ring-opening reactions of butyrolactone (Bulac) and furanone (Fur) respectively after excitation to the second excited state. These reactions differ in the mechanism of ring-opening, Bulac undergoing rapid ballistic bond dissociation within 50 fs, while bond dissociation in Fur

<figure><img src="image_8.png" /><figcaption>**Figure 1.** (a) trans-to-cis isomerization of azomethane, (b) ring opening of butyrolactone, and (c) ring opening of furanone. Reaction times and excitation wavelengths are taken from refs 26 and 27.</figcaption></figure>

occurs more slowly and gradually, ca. 140 fs, due to competition with ring-puckering.²⁷

This article is structured as follows: in Section 2 the methodology is introduced, including the dimensionality reduction procedure, the electronic and geometric properties used to evaluate the reactions, as well as the computational details. Following this, Section 3 presents the findings on the three studied molecules. Subsequently, Section 4 delves into the analysis, before the final conclusions are given in Section 5.

## 2. METHODOLOGY

### 2.1. Dimensionality Reduction Procedure.
The overall methodology used in this work is outlined in Figure 2. First the dynamics are performed in full dimensionality and visited geometries along the trajectories are extracted. Mixed quantum-classical simulations are typically performed in the $3N_{at}$ Cartesian coordinates space. The 6 additional dimensions compared to the $3N_{at} - 6$ conformation space of a nonlinear molecule arise from rotation and translation-quantum dynamics simulations are typically run in normal modes, without these motions. The Kabsch algorithm²⁸ is then subsequently applied to remove any rotational and translational differences from all geometries compared to the Franck-Condon geometry (assigning different weights to the different atoms depending on their mass). The Cartesian geometries are then converted to normal mode coordinates, using the ground-state minimum as the reference geometry and the normal mode vectors calculated at that reference geometry. One dimensionality reduction method is finally applied: PCA or NMV. The Principal Components (PCs) or the Normal Modes (NMs) with the greatest variance are selected, to form a new basis to simulate the dynamics in reduced dimensions.

The first dimensionality reduction approach used, PCA,¹³,¹⁴ is a linear transformation method, known to be a robust technique commonly used in machine learning, statistics, data analysis,²⁹,³⁰ and in dynamics.⁶,¹⁵⁻¹⁷ Its primary usage is to transform a data set with a high number of variables into a new set of uncorrelated variables, called principal components (PC), while capturing the maximum variance present in the data. The created lower number of PC retain the maximum amount of variance possible compared to the original data set, and are ordered by decreasing order of variance: the first PC explains the largest amount of variance in the data set, followed by the second, and so on. See *Supporting Information (SI)*, section 1 for more details.

The second dimensionality reduction approach used, NMV, is a very straightforward technique. It consists of computing
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Workflow of the dimensionality reduction procedure, starting from full dimensionality simulations.</figcaption></figure>

the variance of all geometries explored across the full-dimensional dynamics along each normal mode (as calculated at the initial geometry of the molecule). The normal modes are then ordered according to their explained variance, in decreasing order, and those containing the lowest variance values are then consecutively removed to reduce the dimensionality. See SI, section 2 for more details.

The nonadiabatic dynamics is then simulated in the reduced subspace selected by either PCA or NMV, using the same methodology, i.e., TSH, and computational settings as those in full dimensionality. It is important to realize that the electronic structure properties used in TSH, such as the energies and gradients, are not calculated directly in the lower dimensional subspace. Instead, at each time step, the gradient calculated analytically in the full-dimensional space is transformed—on the fly—into normal modes and projected onto the reduced subspace (using either the PC basis from PCA or the normal mode basis containing only the most relevant modes as determined by the NMV procedure). After this reduction step, the gradient is transformed back into Cartesian coordinates and used in the propagation of Newton's equation of motion via the velocity-Verlet algorithm. In this fashion, the geometry is naturally propagated within the reduced dimensional subspace. An alternative method to remove discarded

dimensions from the dynamics is to project, at each time step, the new geometry predicted using the full-dimensional gradient into the reduced space. Both the reduction on the gradients and the geometries were implemented and tested, leading to the same reduced dimension dynamics. In the present work, the reduction was applied on the gradients (forces). Importantly, initial conditions (geometries and velocities) need to be reduced before starting the dynamics. See SI, section 3 for a detailed step-by-step protocol for the dynamics in reduced dimensionality.

Since the electronic structure calculations are performed in full dimension, there is no increase in speed associated with running TSH in lower dimensions as implemented here; our purpose is solely to provide a proof-of-concept of high-quality dynamics run in the minimal number of dimensions, as well as to compare the two selection techniques (PCA/NMV) for dimensionality reduction protocols.

## 2.2. Dynamical Properties.
The time evolution of the electronic populations, as well as the time evolution of the key nuclear coordinates, for full and reduced dimensionality dynamics obtained with PCA and NMV methods are presented in SI. To quantitatively assess the accuracy of the dimensionality reduction in comparison to the full-dimensional ($3N_{at} - 6$) dynamics, two electronic and two geometric properties were selected for each molecule. These properties are presented in Table 1. For each pair of electronic and

<table><caption>**Table 1.** Electronic and Geometric Properties Chosen to Evaluate the Quality of the Dimensionality Reduction for Each Molecule</caption><thead><tr><th>properties</th><th>trans-azomethane</th></tr></thead><tbody><tr><td>electronic</td><td>S<sub>1</sub> and S<sub>0</sub> state population crossing time (in fs)<br>S<sub>0</sub> state population (in %) at the end of the simulation (280 fs)</td></tr><tr><td>geometric</td><td>time (in fs) to reach half of final cis yield<br>cis yield (in %) at the end of the simulation (280 fs)</td></tr><tr><td>properties</td><td>butyrolactone</td></tr><tr><td>electronic</td><td>time (in fs) for the S<sub>2</sub> population to reach the inflection point<br>S<sub>0</sub> state population (in %) at the end of the simulation (100 fs)</td></tr><tr><td>geometric</td><td>time (in fs) to reach half of the final C–O bond break yield<br>C–O bond break yield (in %) at the end of the simulation (100 fs)</td></tr><tr><td>properties</td><td>furanone</td></tr><tr><td>electronic</td><td>time (in fs) to reach half of the final S<sub>2</sub> population yield<br>S<sub>0</sub> state population (in %) at the end of the simulation (150 fs)</td></tr><tr><td>geometric</td><td>time (in fs) to reach half of the final C–O bond break yield<br>C–O bond break yield (in %) at the end of the simulation (150 fs)</td></tr></tbody></table>

geometric properties, one is chosen to quantify the time scale of the process, and the other its yield. The accuracy of a dynamics simulation is assessed by calculating the absolute relative error (ARE) for each property P

$$ARE = \left| \frac{P_{\text{full}} - P_{\text{red}}}{P_{\text{full}}} \right| \qquad (1)$$

It involves comparing the value of the property calculated from the ensemble of trajectories, for both the full $P_{full}$ and reduced dimensionality $P_{red}$ dynamics simulations.

Finally, an overarching error (labeled as the “total” error) is defined as the average of the ARE over all properties
<figure><img src="image_10.png" /><figcaption>Figure 3. Simulation results for trans -azomethane. Distribution (%) of the trans , mix , and cis conformers at the end of the dynamics simulations as a function of number of dimensions, for both PCA (a) and NMV (b) reduction methods. Evolution of the electronic (c), geometric (d) and the total mean (e) errors, as a function of the number of dimensions, for each reduction method. In (c–e), the shaded area represents 10% of relative error compared to the dynamics performed in full dimensionality.</figcaption></figure>

$$
\text{ARE}_{\text{total}} = \frac{\text{ARE}_{\text{el1}} + \text{ARE}_{\text{el2}} + \text{ARE}_{\text{geom1}} + \text{ARE}_{\text{geom2}}}{4} \quad (2)
$$

With el and geom denoting the electronic and geometric
properties, respectively. The values of each error (electronic,
geometric and total), along with the average values of each
property (ARE) used to evaluate the reactions (Table 1), can
be found in SI.

**2.3. Computational Details.** For each selected set of dimensions for each molecule, 100 trajectory surface hopping<sup>20,21</sup> trajectories corresponding to 100 sampled initial conditions were simulated using the OpenMOLCAS code.<sup>31</sup> The initial conditions were generated using SHARC,<sup>32</sup> sampled from the Wigner distribution without temperature broadening at the Franck–Condon geometry, using harmonic frequencies calculated using OpenMOLCAS (see below for specifics of the electronic structure for each molecule). A set of 100 trajectories was found enough to converge the dynamics in full dimensionality, see ref 33. and its SI. *trans*-Azomethane was simulated for a total of 280 fs, butyrolactone for a total of 100 fs, and furanone for a total of 150 fs. A time step of 20 au (ca. 0.48 fs) for the integration of Newton’s equations of motion through the velocity-Verlet algorithm was used, with the integration of the electronic wave function at each time step split into 96 substeps. The energy based Persico–Granucci decoherence correction, with a decay factor of 0.1 was used throughout.<sup>34</sup> Energy conservation has been checked for all simulations, checking for active space stability, intruder states, and other simulation artifacts. A maximum of 10% of the trajectories were allowed to be rejected before an ensemble was considered to be no longer representative. Importantly, the dimensionality reduction does not lead to an increase in the number of trajectories not conserving total energy in comparison to simulation in full dimensionality.

The location of hops was determined using the Fewest-
Switches Surface Hopping algorithm<sup>20</sup> and after hopping

energy conservation was ensured by rescaling the kinetic
energy along the velocity direction. The Hammes-Schiffer–
Tully scheme³⁵ was used to approximate the nonadiabatic
couplings at each time step, using the biorthonormal wave
function overlap scheme as implemented in OpenMOL-
CAS.³¹,³³

For each molecule, electronic structure calculations were
carried out using the state-averaged complete active space self-
consistent field (SA-CASSCF) method, with active spaces and
state-averages defined as follows: SA2(6e, 4o) for trans-
azomethane, SA3(12e, 10o) for furanone, and SA3(10e, 8o)
for butyrolactone. The active spaces are identical to those
selected in refs 27 and 33, with full orbital descriptions found
in the SI of ref 33. For all molecules the 6-31G(d) basis set³⁶
with a single set of 3d polarization functions³⁷ on the C and N
centers was used. Resolution-of-identity based on the Cholesky
decomposition was used throughout to speed up electronic
structure calculations.³⁸

3. RESULTS

3.1. trans → cis Isomerization of Azomethane (tAZM).

The first photochemical reaction studied is the *trans*-to-*cis* photoisomerization of azomethane, which occurs after excitation to the first excited state. According to prior simulations,<sup>26,39</sup> this photochemical transformation takes place over a period of ca. 400 fs, through a rotation about the central C-N-N-C moiety. The excited state decay takes place through a conical intersection located roughly halfway between the two isomers, with a central dihedral around 90°.

The time evolutions of the electronic populations of S₀ and S₁ states, as well as the time evolutions of the dihedral torsion angle C-N-N-C, which describes the reaction coordinate, for the full and reduced dimension dynamics are presented in SI. The compositions of the principal components in terms of normal modes, as well as graphical representations of the first four components calculated from the full-dimensional
dynamics are reported in the SI. The first components consist of CH₃ rocking, NN stretching, CNN asymmetric bending and NN torsional motions.

In order to follow the progress of the reaction, the trajectories are classified into one of three conformers based on the value of this central dihedral angle: *cis* for angles under 45°, *mix* for those between 45° and 135°, and *trans* for angles of 135° or more. The fractions of conformers at the end of the simulations using this classification are shown in Figure 3 for both the (a) PCA and (b) NMV dimensionality reduction methods. A finer classification of the geometries based on the central dihedral angle leads to the same conclusions (see SI).

Both the NMV and PCA methods allow for some level of dimensionality reduction to be carried out without significant changes to the conformer distributions. In general, the distributions of *trans*, *mix*, *cis* at the end of the dynamics, shown in Figure 3(a,b), evolve quite similarly between the two reduction methods as the number of dimensions decreases. Differences to the reference final distribution obtained in full dimensionality (3N<sub>at</sub> − 6 = 24) begin to be seen using fewer than 20 dimensions, suggesting reduction from 24 to 20 dimensions is feasible. With 16 dimensions and fewer, the final balance of isomers shifts more toward the *trans* isomer: i.e., the *trans*-to-*cis* isomerization is reduced. While from these conformer distributions it seems like the PCA method may allow for slightly more dimensionality reduction to 18 dimensions, it is hard to conclude which method allows for more accurate dimensionality reduction based on this property analysis alone.

Figure 3(d) shows the quantitative error obtained for the two chosen geometrical properties. We define accuracy in the dynamics results to be when ARE (compared to the full-dimensional dynamics) is below a threshold of 10%. This threshold is significantly larger than the variance expected from not fully converged simulations (see ref 33 and its SI). We see for the NMV a loss of accuracy at 18 dimensions. The PCA result is less clear, with an initial loss of accuracy at 17 dimensions followed by a decrease in ARE with 16 dimensions.

Figure 3(c) shows the absolute relative error in the electronic properties. Both reduction methods show a clear lowest dimensionality point beyond which the ARE dramatically increases. The PCA method maintains accuracy when decreasing to 15 dimensions, beyond which the ARE significantly increases, while the NMV method remains accurate only up to 17 dimensions. We note that the AREs are, in general for tAZM, lower for the geometric properties compared to the electronic ones.

For the NMV method, the total error dramatically increases when decreasing from 17 to 15 dimensions. The normal modes (mass weighted, with contributions on hydrogens hidden) removed for this dimensionality reduction are illustrated in Figure 4. It is noted that the normal modes are numbered according to their explained variance in decreasing order, consistently with the PCA modes. The first mode (a) appears to represent a symmetrical C≡N angle opening and reduction of the N≡N bond length, while the second mode (b) is less straightforward, possibly indicating the opening of a single C≡N = N angle. These two normal modes are likely to play significant roles in the *trans*-to-*cis* isomerization process. However, because the correlation between their variance and the main reaction pathway is not direct, these were removed, showing the limit of an approach only based on the variance. This is likely to be a notable problem for motions which vary

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Two normal modes (mass weighted, hydrogen contributions hidden for ease of interpretation) removed when decreasing from 17 to 16 (a) and 16 to 15 (b) dimensions for the azomethane molecule using the NMV method. The normal modes are numbered according to their explained variance (in decreasing order).</figcaption></figure>

only a small amount within the dynamics, but which are still important for the photochemical mechanism: previously suggested to be the case for the length of the N≡N bond.³⁹ Furthermore, it should be noted that the significance of all the preceding normal modes removed must also be considered, and the removal of these two modes may simply mark a pivotal moment where the dynamics undergoes a drastic change due to the overall set of modes removed.

Clearly, the total error as shown in Figure 3(e) is more influenced by the electronic errors and unambiguously shows that the PCA method exhibits a greater capacity to operate at lower dimensions when compared to the NMV method. Based on this total error, we predict that it would be possible to reduce the dimensionality required to describe this *trans*-to-*cis* photoisomerization of azomethane from 24 (3N<sub>at</sub> − 6) to only 16 (18) dimensions, using the PCA (NMV) reduction technique.

**3.2. Ring Opening of Butyrolactone (Bulac).** The second reaction studied is the photoinduced ring-opening of butyrolactone, triggered by excitation to the S₂ state. This reaction is characterized by a very rapid ballistic process, with the C=O bond predicted by previous simulations to break within the first 50 fs.²⁷

The time evolutions of the electronic populations of S₀, S₁ and S₂ states, as well as the time evolutions of the breaking C=O bond length, which describes the reaction coordinate, for the full and reduced dimension dynamics are presented in SI. The compositions of the principal components in terms of normal modes, as well as graphical representations of the first four components calculated from the full-dimensional dynamics are also reported in SI. The principal components consist of CC stretching, CO stretching, CO out of plane, CH₂ rocking and ring torsional motions.

Figure 5 illustrates the performance of PCA and NMV reduction methods for the Bulac molecule. Insets (a) and (b) in Figure 5 display the evolution of the percentage of C=O bonds broken over the course of the dynamics for each dimensional space, for PCA and NMV methods, respectively. The quantity of C=O bonds broken is already notably different to the full 30-dimensional dynamics with 24 dimensions selected by NMV (Figure 5(b)), while using PCA reduction notable differences in C=O breakage only appear with fewer than 7 dimensions (Figure 5(a)). PCA is again more accurate than NMV when reducing the dimensionality. The combination of geometric properties Figure 5(d) captures the difference between both methods clearly.
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Simulation results for butyrolactone. Top panel: evolution (%) of C–O bond break as a function of time (fs) for each dimensionality, using the PCA (a) and the NMV (b) method. Evolution of the electronic (c), geometric (d), and total mean errors (e), as a function of the number of dimensions for both applied reduction methods. In (c), (d), and (e), the shaded area represents 10% of relative error compared to the dynamics in full dimensions.</figcaption></figure>

Looking more closely, using the NMV method, the percentage of C–O bonds broken is no longer well-described decreasing to 22 dimensions, and then becomes significantly worse at 20 dimensions (Figure 5(b)). Not only is the amount of ring-opening decreased, but also some bond reformation occurs (see Figures 5b and SI). Figure 6 shows the 4 (mass weighted) normal modes removed by the NMV method for this decrease from 24 to 20 dimensions. Interestingly, the first three of these modes seem to exhibit only a small contribution toward the C–O bond breaking process which leads to ring-opening, with larger contributions on the hydrogens and other

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Four normal modes removed when decreasing from 24 to 22 dimensions (a, b) and from 22 to 20 dimensions (c, d) for the butyrolactone molecule. The normal modes are numbered according to their explained variance (in decreasing order).</figcaption></figure>

carbons in the ring. While these three modes do not seem directly involved in the bond breaking, they are probably important in the finalization of the ring opening and the prevention of ring reclosing. The final mode (d, normal mode 21) is the only one of these four seemingly important modes which contains a significant component on the C–O bond to be broken, and would have been intuitively expected to be necessary to describe the reaction if selecting modes manually.

This early failure of NMV contrasts to the much better performance of the PCA method, where each PC used as a dimension incorporates weighted contributions of all 30 normal modes. Figure 5(a) shows a slight deterioration of the ring-opening reaction upon removal of principal components 7 and 6. Interestingly, these components contain some CO out-of-plane and torsional motions (see SI).

Concerning the electronic properties as seen in Figure 5(c) both methods perform at an essentially equivalent level, with a steady and relatively low (less than 0.2) error throughout. Interestingly, for NMV reduction, the ARE increase to higher values for the geometric properties compared to the electronic ones, in contrast to the results for tAZM. As a result, the total error as shown in Figure 5(e) shows PCA to be more accurate. From this total error, we can conclude that using PCA we could reduce the full (in NM basis, 3N<sub>at</sub> − 6) 30-dimensional dynamics to only 7 dimensions, while NMV allows us to reduce to 24 dimensions. It is interesting to note that the error does not necessarily increase monotonically with the decrease of the number of dimensions. We do not believe that this is due to a poor convergence of the dynamics simulations but rather to the alternation of removals of coordinates that enhance and delay electronic relaxation, leading to some kind of cancellation of errors.

**3.3. Ring Opening of Furanone (Fur).** The final reaction studied is the photoinduced ring-opening of the furanone
<figure><img src="image_14.png" /><figcaption>**Figure 7.** Simulation results for furanone. Evolution (%) of C–O bond break in function of time (fs) and dimensions for the PCA (a) and the NMV (b) methods. Evolution of the electronic (c), geometric (d) and the total mean (e) errors respectively, as a function of the number of dimensions for both reduction methods. In (c), (d), and (e), the shaded area represents 10% of relative error compared to the full dimension dynamics.</figcaption></figure>

molecule which, like the Bulac, is triggered by excitation to the S₂ state. This ring-opening is also characterized by a breaking of a C–O bond. However, unlike Bulac, this ring-opening process is more gradual, taking place on a ca. 180 fs time scale due to competition with ring-puckering.²⁷

The time evolutions of the electronic populations of S₀, S₁ and S₂ states, as well as the time evolutions of the breaking C–O bond length, which describes the reaction coordinate, for the full and reduced dimension dynamics are presented in SI. The compositions of the principal components in terms of normal modes, as well as graphical representations of the first four components calculated from the full-dimensional dynamics are also reported in SI. The principal components consist of CO stretching, CH out of plane, and COC bending motions.

**Figures 7**(a,b) illustrate the evolution of the percentage of C–O bonds broken over simulation time (150 fs), for both the PCA and the NMV methods, respectively. Using the PCA dimensionality reduction technique, two groups emerge toward the end of the dynamics. The lower percentage group contains the full-dimensionality dynamics, maintaining a similar percentage of C–O breakage by the end of the simulations (~57%) to as low as 14 dimensions. The higher percentage group is made up of the lower-dimensionality simulations, and is centered around 65% of C–O bond breakage by the end of the simulations. The most likely explanation for this result is that, when removing additional dimensions below the threshold of 14, part of the competing puckering mechanism is hindered, leading to a faster C–O bond breakage. **Figure 8**(b) depicts one of the PCA modes removed going from 14 to 12 dimensions and seems to support this hypothesis, with the carbons showing out-of-plane motion. Expressed in normal mode coordinates, the principal components 14 and 13 consist of CCC bending, CH out of plane, and CC stretching motions (see SI).

<figure><img src="image_15.png" /><figcaption>**Figure 8.** (a) Normal mode 10 removed when reducing from 10 to 9 dimensions using the NMV method, and (b) the PC 14 removed when reducing from 14 to 12 dimensions using the PCA method, for the furanone molecule. The normal modes are numbered according to their explained variance (in decreasing order), consistently with the principal components.</figcaption></figure>

In contrast, using the NMV dimensionality reduction method, the percentage of C–O bonds broken steadily rises as the dimensionality decreases. This suggests a gradual attenuation of the puckering mechanism with each removed normal mode. However, upon reducing to 9 dimensions, there is a notable decrease in the breakage of C–O bonds (from ~60%, to ~40%), indicating that the removed normal mode (shown in **Figure 8**(a)) instead principally limits the ring-opening process.

Upon examining the geometric properties (**Figure 7**(d)), it is clear that the PCA method allows for a more significant dimensionality reduction, with lower relative error than NMV at all dimensionalities, and accuracy until fewer than 16 dimensions are used. In contrast, the PCA and the NMV methods exhibit quite similar relative errors in electronic properties at all dimensionalities (**Figure 7**(c)). The total error, as a result, does indicate that PCA might allow for a reduction in dimensionality, remaining under the 0.1 threshold until reducing to 15 dimensions. However, the similar in scale relative errors in geometric and electronic properties, along
<figure><img src="image_16.png" /><figcaption>**Figure 9.** Top panel: evolution of the explained variance (%) by the PCA (blue) and NMV (orange) methods, as a function of the retained dimensions for (a) trans-azomethane, (b) butyrolactone, and (c) furanone. Bottom panel: evolution of the total error ARE<sub>total</sub> with the PCA (blue) and NMV (orange) methods as a function of the explained variance (%) for (d) trans-azomethane, (e) butyrolactone, and (f) furanone.</figcaption></figure>

with the similar performance for electronic properties, result in the total error giving a less clear indication for possible dimensionality reduction compared to using simply the geometric error. In fact, the quality of the dimension reduction for this reaction seems to be best evaluated by studying simply the percentage of broken C–O bonds. While the results are less clear for this molecule, there seems to be a possible dimensionality reduction from 24 to 16 dimensions using PCA, compared to no clear possible reduction using NMV.

## 4. DISCUSSION

In this Section, we start by discussing the differences in performances of the two reduction methods. The three molecules studied here all demonstrate that using PCA for dimensionality reduction yields higher quality results compared to the NMV method. This is, of course, expected since PCA is performed in normal mode representation. NMV consists of removing given normal modes while PCA eliminates linear combinations of them, allowing it to lose less variance per removed dimension. This is illustrated in Figures 9(a–c), which show the quantity of variance retained at each level of dimensionality for NMV and PCA. For instance, in the case of tAZM with only 15 dimensions selected with the NMV method, about 87% of the variance is kept, whereas with 15 dimensions selected with PCA, over 96% is still explained. This is also clearly seen for the other reactions in (b) and (c): PCA always leads to a slower loss of variance compared to NMV.

However, variance is not the whole story, as illustrated in Figures 9(d–f) which display the total ARE (%) as a function of the explained variance (%) for each molecule and method. For the tAZM molecule (d), the total ARE increases more rapidly with loss of variance under the PCA method compared to NMV, despite the similar accuracy of the two methods as already shown in Section 3.1. In the case of Bulac (e) PCA maintains a total ARE at 10% down to 95.5% of kept variance,

while NMV crosses the 10% error threshold at only 98% variance. For Fur (e), the evolution of error with variance is comparable between the two methods. We do not see a straightforward increase in error as a function of variance which follows a universal trend for both methods—the specific motions removed from the dimension space are more important than the simple value of remaining variance. The variance threshold below which error remains low is dependent on both the method of reduction, and the photochemical reaction, demonstrating the difficulty of establishing a generic variance threshold rule applicable to multiple molecular systems.

Interestingly, the superior performance of PCA compared to NMV is made clear by different types of property, depending on the molecule. For tAZM, the electronic properties most clearly differentiate between the two methods (Figure 3(b)), while for Bulac it is clearly the geometric properties which are the most helpful (Figure 5(d)). Finally, for furanone while the differences in geometric properties are clearer than the electronic ones, the advantage of PCA over NMV is best made obvious from the percentage of C–O broken bonds (Figures 7(a,b)). This variability in sensitivity demonstrates an interesting balance between the use of electronic and geometric properties to characterize performance in different reactions, showing a nuanced interplay between molecular and electronic structure depending on the reaction studied. In order to characterize the performance of NMV compared to PCA, properties thought to be representative of the reaction process of interest were selected. This process involves human input.

We now discuss more general aspects of the present approach. We note that it relies on an appropriate description of the dynamics by the lower-cost method. We have tested the robustness of the protocol with respect to the number of analyzed trajectories in the low-level dynamics in full dimensionality (see SI for details). Similar ARE were obtained
with fewer trajectories, indicating that the original ensemble of trajectories to be analyzed does not need to be converged quantitatively. A qualitative description of the process of interest may be sufficient. In principle, any relaxation pathway simulated by the low-level dynamics will be analyzed and taken into account in the dimension reduction protocol. However, the presence of a dominating pathway may steer the analysis and thus prevent the occurrence of minor pathways upon dimension reduction. If one encounters such a problem, a solution may be to separate the trajectories based on the followed pathway, to perform the dimension analysis separately for the different groups, and finally to merge the extracted subspaces.

We have also tested the impact of the total simulation time of the original ensemble of trajectories on the accuracy of the dimension reduction (see SI for details). The influence of the simulation time on the total ARE is in general more pronounced when a relatively low number of dimensions is used, whereas it is limited when a large number of dimensions is used. Our recommendation is to set the simulation time such that the whole process of interest is simulated, and to not expect a set of dimensions extracted from short-time dynamics to be relevant for longer time scale.

In addition, the current protocols are based on normal mode coordinates, since these are expected to perform better than Cartesian coordinates and would allow more direct applications to quantum dynamics simulations. However, normal modes are known to present limitations in case of large amplitude motions. Studying other systems of coordinates would be an interesting topic for a further study.

Finally, we also compared the relevant coordinates extracted with the NMV and PCA approaches to gradients and nonadiabatic coupling vectors calculated at the Franck–Condon geometry. The latter quantities are often used to construct “one shot” Linear Vibronic Coupling (LVC) models. The results show that the correlation is relatively low: see SI for more details. As previously noted,¹² if the dynamics explores regions far away from the Franck–Condon point, these directions are not enough to provide a reliable description of the dynamics of the photoexcited molecule.

## 5. CONCLUSIONS

In this work, we tested the performance of two dimensionality reduction methods, NMV and PCA, in order to identify the most important nuclear coordinates in nonadiabatic processes. The accuracy of both electronic and geometrical properties in reduced-dimensionality dynamics simulations have been assessed on three different photochemical reactions, including both double bond photoisomerization and ring-opening processes. In general, both reduction methods allow for the construction, in an automated manner, of lower dimensional spaces which well reproduce the full dimensional dynamics.

Using the PCA method, the dimensionality can be reduced from 24 to 16 for tAZM, from 30 to 7 for Bulac, and from 24 to 16 for Fur. Using NMV, the dimensionality can be reduced from 24 to 18 for tAZM, and from 30 to 24 for Bulac, with no clear reduction being possible for Fur. The PCA method always performs better than NMV, with accuracy in dynamics compared to the full-dimensional simulations being preserved when reducing to lower dimensional subspaces.

Depending on the molecule, the improvement offered by PCA vs NMV can be more or less significant. For instance, for the tAZM molecule we find that using PCA allows for a

reduction of only 1 or 2 additional dimensions compared to NMV, while for the Bulac molecule using PCA allows as many as 15 additional dimensions to be removed compared to NMV, and for the Fur molecule only PCA allows for any removal of dimensions.

The use of lower-dimensional spaces optimally constructed through methods such as PCA could provide a viable way to run expensive dynamics simulations, which are often limited due to the curse of dimensionality. One could imagine using the approach outlined in this work (PCA based on TSH variance) to select low-dimensional subspaces built through PCA in which higher cost dynamics could be run. For instance, one could use the selected reduced subspace to simulate the dynamics with more expensive on-the-fly dynamics methods like DD-vMCG or variants of quantum Ehrenfest. Alternatively, one can construct a Hamiltonian model based on the dimension reduction analysis for subsequent MCTDH simulations. We emphasize that building Hamiltonian models remains challenging, and this work only aims to provide guidelines for the first step, the selection of relevant nuclear coordinates. In order to generalize our findings further and learn more about which systems lend themselves better to more efficient dimensionality reduction through PCA, further investigations on more diverse and complex molecular systems are required.

## ASSOCIATED CONTENT

**SI Supporting Information**

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.jctc.5c00110.

Detailed explanation of both PCA and NMV dimensionality reduction techniques and detailed protocol of dynamics in reduced dimensionality; raw data for all molecules simulated at all dimensions; time evolutions of electronic populations and key nuclear coordinates, analysis of hopping geometries along key internal coordinates, finer-classication of *trans*-azomethane conformers, values of each descriptor used for evaluation, averaged over the trajectory ensemble, for each dimensionality in which the dynamics are run for both PCA and NMV reduction methods. Analysis of principal components in terms of normal modes and graphical representations for the former; test of the impact of the simulation time. Test of the robustness of the protocol with respect of the number of analyzed trajectories; comparison with gradient and nonadiabatic coupling vector coordinates (PDF)

## AUTHOR INFORMATION

**Corresponding Author**

Morgane Vacher – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; ● orcid.org/0000-0001-9418-6579; Email: morgane.vacher@univ-nantes.fr

**Authors**

Vincent Delmas – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France

Alessandro Nicola Nardi – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; ● orcid.org/0000-0002-0959-1690

Isabella C. D. Merritt – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; Univ. Bordeaux,
CNRS, Bordeaux INP, Institut des Sciences Moléculaires (ISM), F-33405 Talence, France; orcid.org/0000-0002-7965-7798

Anthony Ferté – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; orcid.org/0000-0001-8896-9219

Ignacio Fdez. Galván – Department of Chemistry—BMC, Uppsala University, SE-751 23 Uppsala, Sweden; Department of Chemistry—Ångström, Uppsala University, SE-751 20 Uppsala, Sweden; orcid.org/0000-0002-0684-7689

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jctc.5c00110

## Notes

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

The project is partly funded by the European Union (ERC, 101040356 - ATTOP, M.V., V.D. and A.N.N.). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. We also thank the Région des Pays de la Loire who provided postdoctoral funding for A.F., and Nantes University who provided doctoral support for I.C.D.M. The authors thank the CCIPL/Glicid mesocenter installed in Nantes and GENCI-IDRIS (Grant 2021-101353) for the generous allocation of computational time.

## REFERENCES

(1) Balzani, V.; Credi, A.; Venturi, M. Photochemical Conversion of Solar Energy. *ChemSusChem* **2008**, *1*, 26–58.

(2) Oelgemöller, M. Solar Photochemical Synthesis: From the Beginnings of Organic Photochemistry to the Solar Manufacturing of Commodity Chemicals. *Chem. Rev.* **2016**, *116*, 9664–9682.

(3) Pathak, S.; Ibele, L. M.; Boll, R.; Callegari, C.; Demidovich, A.; Erk, B.; Feifel, R.; Forbes, R.; Di Fraia, M.; Giannessi, L.; Hansen, C. S.; Holland, D. M. P.; Ingle, R. A.; Mason, R.; Plekan, O.; Prince, K. C.; Rouzée, A.; Squibb, R. J.; Tross, J.; Ashfold, M. N. R.; Curchod, B. F. E.; Rolles, D. Tracking the ultraviolet-induced photochemistry of thiophenone during and after ultrafast ring opening. *Nat. Chem.* **2020**, *12*, 795–800.

(4) Mai, S.; González, L. Molecular Photochemistry: Recent Developments in Theory. *Angew. Chem., Int. Ed.* **2020**, *59*, 16832–16846.

(5) Latka, T.; Shirvanyan, V.; Ossiander, M.; Razskazovskaya, O.; Guggenmos, A.; Jobst, M.; Fieß, M.; Holzner, S.; Sommer, A.; Schultze, M.; Jakubeit, C.; Riemensberger, J.; Bernhardt, B.; Helml, W.; Gatti, F.; Lasorne, B.; Lauvergnat, D.; Decleva, P.; Halász, G. J.; Vibók, A.; Kienberger, R. Femtosecond wave-packet revivals in ozone. *Phys. Rev. A* **2019**, *99*, No. 063405, DOI: 10.1103/PhysRevA.99.063405.

(6) Zauleck, J. P. P.; Thallmair, S.; Loipersberger, M.; de Vivieriedle, R. Two New Methods To Generate Internal Coordinates for Molecular Wave Packet Dynamics in Reduced Dimensions. *J. Chem. Theory Comput.* **2016**, *12*, 5698–5708.

(7) Liu, J.; He, X.; Wu, B. Unified Formulation of Phase Space Mapping Approaches for Nonadiabatic Quantum Dynamics. *Acc. Chem. Res.* **2021**, *54*, 4215–4228.

(8) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The multi-configurational time-dependent Hartree approach. *Chem. Phys. Lett.* **1990**, *165*, 73–78.

(9) Worth, G. A.; Meyer, H.-D.; Köppel, H.; Cederbaum, L. S.; Burghardt, I. Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics. *Int. Rev. Phys. Chem.* **2008**, *27*, 569–606.

(10) Mai, S.; González, L. Identification of important normal modes in nonadiabatic dynamics simulations by coherence, correlation, and frequency analyses. *J. Chem. Phys.* **2019**, *151*, No. 244115, DOI: 10.1063/1.5129335.

(11) Gómez, S.; Heindl, M.; Szabadi, A.; González, L. From Surface Hopping to Quantum Dynamics and Back. Finding Essential Electronic and Nuclear Degrees of Freedom and Optimal Surface Hopping Parameters. *J. Phys. Chem. A* **2019**, *123*, 8321–8332.

(12) Plasser, F.; Gómez, S.; Menger, M. F. S. J.; Mai, S.; González, L. Highly efficient surface hopping dynamics using a linear vibronic coupling model. *Phys. Chem. Chem. Phys.* **2019**, *21*, 57–69.

(13) Pearson, K. LIII. On lines and planes of closest fit to systems of points in space. London, Edinburgh Dublin. *London, Edinburgh Dublin Philos. Mag. J. Sci.* **1901**, *2*, 559–572.

(14) Wold, S.; Esbensen, K.; Geladi, P. Principal component analysis. *Chemom. Intell. Lab. Syst.* **1987**, *2*, 37–52.

(15) Capano, G.; Penfold, T. J.; Chergui, M.; Tavernelli, I. Photophysics of a copper phenanthroline elucidated by trajectory and wavepacket-based quantum dynamics: a synergetic approach. *Phys. Chem. Chem. Phys.* **2017**, *19*, 19590–19600.

(16) Tavadze, P.; Avendaño Franco, G.; Ren, P.; Wen, X.; Li, Y.; Lewis, J. P. A Machine-Driven Hunt for Global Reaction Coordinates of Azobenzene Photoisomerization. *J. Am. Chem. Soc.* **2018**, *140*, 285–290.

(17) Zhu, Y.; Peng, J.; Kang, X.; Xu, C.; Lan, Z. The principal component analysis of the ring deformation in the nonadiabatic surface hopping dynamics. *Phys. Chem. Chem. Phys.* **2022**, *24*, 24362–24382.

(18) Li, X.; Xie, Y.; Hu, D.; Lan, Z. Analysis of the Geometrical Evolution in On-the-Fly Surface-Hopping Nonadiabatic Dynamics with Machine Learning Dimensionality Reduction Approaches: Classical Multidimensional Scaling and Isometric Feature Mapping. *J. Chem. Theory Comput.* **2017**, *13*, 4611–4623.

(19) Li, X.; Hu, D.; Xie, Y.; Lan, Z. Analysis of trajectory similarity and configuration similarity in on-the-fly surface-hopping simulation on multi-channel nonadiabatic photoisomerization dynamics. *J. Chem. Phys.* **2018**, *149*, No. 244104, DOI: 10.1063/1.5048049.

(20) Tully, J. C. Molecular dynamics with electronic transitions. *J. Chem. Phys.* **1990**, *93*, 1061–1071.

(21) Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H⁺ with D₂. *J. Chem. Phys.* **1971**, *55*, 562–572.

(22) Gómez, S.; Spinlove, E.; Worth, G. Benchmarking non-adiabatic quantum dynamics using the molecular Tully models. *Phys. Chem. Chem. Phys.* **2024**, *26*, 1829–1844.

(23) Ibele, L. M.; Curchod, B. F. E. A molecular perspective on Tully models for nonadiabatic dynamics. *Phys. Chem. Chem. Phys.* **2020**, *22*, 15183–15196.

(24) Janoš, J.; Slaviček, P. What Controls the Quality of Photodynamical Simulations? Electronic Structure Versus Non-adiabatic Algorithm. *J. Chem. Theory Comput.* **2023**, *19*, 8273–8284.

(25) Tran, T.; Ferté, A.; Vacher, M. Simulating Attochemistry: Which Dynamics Method to Use? *J. Phys. Chem. Lett.* **2024**, *15*, 3646–3652.

(26) Ruckenbauer, M.; Barbatti, M.; Sellner, B.; Muller, T.; Lischka, H. Nonadiabatic Photodynamical Simulations in Solution. *J. Phys. Chem. A* **2010**, *114*, 12585–12590.

(27) Schalk, O.; Galiana, J.; Geng, T.; Larsson, T. L.; Thomas, R. D.; Fdez. Galván, I.; Hansson, T.; Vacher, M. Competition between ring-puckering and ring-opening excited state reactions exemplified on 5H-furan-2-one and derivatives. *J. Chem. Phys.* **2020**, *152*, No. 064301, DOI: 10.1063/1.5129366.

(28) Kabsch, W. Automatic indexing of rotation diffraction patterns. *J. Appl. Crystallogr.* **1988**, *21*, 67–72.
(29) Bro, R.; Smilde, A. K. Principal component analysis. *Anal. Methods* **2014**, *6*, 2812–2831.

(30) Karamizadeh, S.; Abdullah, S. M.; Manaf, A. A.; Zamani, M.; Hooman, A. An overview of principal component analysis. *J. Signal Inf. Process.* **2013**, *04*, 173–175.

(31) Li Manni, G.; Fdez. Galván, I.; Alavi, A.; Aleotti, F.; Aquilante, F.; Autschbach, J.; Avagliano, D.; Baiardi, A.; Bao, J. J.; Battaglia, S.; Birnoschi, L.; Blanco-González, A.; Bokarev, S. I.; Broer, R.; Cacciari, R.; Calio, P. B.; Carlson, R. K.; Carvalho Couto, R.; Cerdán, L.; Chibotaru, L. F.; Chilton, N. F.; Church, J. R.; Conti, I.; Coriani, S.; Cuéllar-Zuquin, J.; Daoud, R. E.; Dattani, N.; Decleva, P.; de Graaf, C.; Delcey, M. G.; De Vico, L.; Dobrautz, W.; Dong, S. S.; Feng, R.; Ferré, N.; Filatov(Gulak), M.; Gagliardi, L.; Garavelli, M.; González, L.; Guan, Y.; Guo, M.; Hennefarth, M. R.; Hermes, M. R.; Hoyer, C. E.; Huix-Rotllant, M.; Jaiswal, V. K.; Kaiser, A.; Kaliakin, D. S.; Khamesian, M.; King, D. S.; Kochetov, V.; Krošnicki, M.; Kumaar, A. A.; Larsson, E. D.; Lehtola, S.; Lepetit, M.-B.; Lischka, H.; López Ríos, P.; Lundberg, M.; Ma, D.; Mai, S.; Marquetand, P.; Merritt, I. C. D.; Montorsi, F.; Mörchen, M.; Nenov, A.; Nguyen, V. H. A.; Nishimoto, Y.; Oakley, M. S.; Olivucci, M.; Oppel, M.; Padula, D.; Pandharkar, R.; Phung, Q. M.; Plasser, F.; Raggi, G.; Rebolini, E.; Reiher, M.; Rivalta, I.; Roca-Sanjuán, D.; Romig, T.; Safari, A. A.; Sánchez-Mansilla, A.; Sand, A. M.; Schapiro, I.; Scott, T. R.; Segarra-Martí, J.; Segatta, F.; Sergentu, D.-C.; Sharma, P.; Shepard, R.; Shu, Y.; Staab, J. K.; Straatsma, T. P.; Sørensen, L. K.; Tenorio, B. N. C.; Truhlar, D. G.; Ungur, L.; Vacher, M.; Veryazov, V.; Voß, T. A.; Weser, O.; Wu, D.; Yang, X.; Yarkony, D.; Zhou, C.; Zobel, J. P.; Lindh, R. The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry. *J. Chem. Theory Comput.* **2023**, *19*, 6933–6991.

(32) Mai, S.; Avagliano, D.; Heindl, M.; Marquetand, P.; Menger, M. F. S. J.; Oppel, M.; Plasser, F.; Polonius, S.; Ruckenbauer, M.; Shu, Y.; Truhlar, D. G.; Zhang, L.; Zobel, P.; González, L. SHARC3.0: Surface Hopping Including Arbitrary Couplings – Program Package for Non-Adiabatic Dynamics 2023 https://sharc-md.org/.

(33) Merritt, I. C. D.; Jacquemin, D.; Vacher, M. Nonadiabatic Coupling in Trajectory Surface Hopping: How Approximations Impact Excited-State Reaction Dynamics. *J. Chem. Theory Comput.* **2023**, *19*, 1827–1842.

(34) Granucci, G.; Persico, M. Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping. *J. Chem. Phys.* **2007**, *126*, No. 134114.

(35) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. *J. Chem. Phys.* **1994**, *101*, 4657–4667.

(36) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. *J. Chem. Phys.* **1972**, *56*, 2257–2261.

(37) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. *Theor. Chim. Acta* **1973**, *28*, 213–222.

(38) Aquilante, F.; Lindh, R.; Pedersen, T. B. Unbiased Auxiliary Basis Sets for Accurate Two-Electron Integral Approximations. *J. Chem. Phys.* **2007**, *127*, No. 114107.

(39) Sellner, B.; Ruckenbauer, M.; Stambolić, I.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. Photodynamics of Azomethane: A Nonadiabatic Surface-Hopping Study. *J. Phys. Chem. A* **2010**, *114*, 8778–8785.

<figure><img src="image_17.png" /></figure>

CAS BIOFINDER DISCOVERY PLATFORM™

STOP DIGGING
THROUGH DATA
—START MAKING
DISCOVERIES

CAS BioFinder helps you find the
right biological insights in seconds

Start your search

<figure><img src="image_18.png" /></figure>