Singlet machine learning photodynamics reveal competing inversion paths of methylated cyclooctatetrathiophene

Abstract

We used state-of-the-art machine-learning nonadiabatic molecular dynamics to investigate the stereochemical inversion reaction of a methylated thiophene-fused cyclooctatetraene derivative, MeCOTh. Minimum energy path calculations suggest that the steepest descent path on the S₁ surface of MeCOTh is towards a non-productive fluorescence decay pathway. Our machine learning photodynamics calculations revealed that relative stereochemical inversion occurs mainly on the S₁ surface (74% of trajectories), and we identified two competing inversion pathways. The first and main mechanistic pathway, seen in 62% of trajectories, showcases a “crown” structure with unidirectional sulfurs resulting from S–S closed-shell repulsions. The second pathway is the previously proposed inversion mechanism, which proceeds through a planar geometry of MeCOTh, and appeared in only 8% of trajectories. Our photodynamic simulations show that although excited-state Baird aromaticity contributes to the relative stereochemical inversion mechanism of MeCOTh, it is not the only electronic effect. Instead, the overall inversion mechanism is primarily governed by the interplay between Baird aromaticity and the S–S closed-shell repulsions.

---

Introduction

Photochemistry exemplifies the principles of green chemistry by utilizing renewable solar energy to facilitate reactions under mild conditions while minimizing resource use and waste production. Light enables the precise activation of chemical processes, such as [2 + 2] cycloaddition reactions¹ and gas evolutions towards strained compounds (e.g., cyclooctynes,² tetrahedrane,³ and bicyclo[1.1.0]butane⁴), that would otherwise be thermally or kinetically unfavorable. The specific stereochemistry,⁵ efficient atom economy, and easy implementations⁶ of these light-mediated reactions are of great interest in chemistry,^7–9 biology,^10–12 and material sciences.^13–15 Photochemical processes promote structural changes through isomerizations,^16–18 bond breakage^19–21/formation,^22,23 and or electronic changes, such as the inversion of aromaticity.^24–26 This precise spatiotemporal control has enabled the applications of photochemical reactions in energy storage materials^1,27,28 and drug delivery.^29–31

Many organic chromophores feature conjugated π-systems; some of which are aromatic in the ground state. Aromaticity is a fundamental concept in chemistry that has evolved since its introduction in the 19th century.^32,33 Aromatic compounds are exceptionally prevalent,^33–35 with approximately two-thirds of all identified molecules exhibiting some aromatic component.²⁴ The aromatic archetype, benzene, was first isolated in 1825 by Faraday.³⁶ Huckel presented his foundational work on aromaticity in 1931 using molecular orbital theory to describe cyclic, π-conjugated hydrocarbons with ‘4n + 2’-π electrons as ‘aromatic’.^34,37 In 1965, Breslow built upon Huckel's ‘4n + 2’ rule by introducing the counter-concept of ‘antiaromaticity,’ which states that cyclic, π-conjugated hydrocarbons with ‘4n’ π electrons are antiaromatic.^34,38 In 1966, Dewar built on aromaticity studies by using perturbation molecular orbital theory to calculate and analyze the heat of formation energies of pairs of open-chain polyenes and their cyclic counterparts (i.e., annulenes), offering qualitative evidence for the 4n + 2 and 4n rules.^32,39 In 1972, Baird conducted foundational research on excited-state aromaticity, identifying an inversion of Huckel's aromaticity rules in the triplet state (i.e., antiaromatic structure: 4n + 2 π and aromatic structure: 4n π).⁴⁰ Karadakov later extended this concept to the singlet excited state.^41,42 In the last decade, Sung et al. and Oh et al. have successfully experimentally observed the reversal of aromaticity in the excited states of hexaphyrins using time-resolved infrared spectroscopy, providing experimental evidence for the existence of Baird's aromaticity.^43,44

Baird's excited-state aromaticity and antiaromaticity have been leveraged to create innovative photochemical reactions^26,45–47 and photoresponsive materials,^48–50 utilizing the principles of excited-state aromaticity^51,52 or antiaromaticity relief.^53–55 Examples of these reactions include the photosolvolysis of 9-fluorenol⁵⁶ and cyclooctatetraene-fused acene dimers,^57,58 as it stabilizes essential mechanistic structures in the excited-state surfaces that would be unstable in the ground-state surface.²⁴ In both instances, the ground-state (S₀) antiaromatic planar geometry of an intermediate⁵⁹ (such as the 4π-cationic intermediate)⁵⁶ or transition state (like the 8π cyclooctatetraene core)^57,58,60 adopts an aromatic electronic configuration in the excited state, thereby facilitating its respective reaction. Scheme 1 illustrates the excited-state aromatic 4π-cationic intermediate, which facilitates the hydroxy extraction from the 9-fluorenol.^24,56

Scheme 1 Photosolvolysis of 9-fluorenol highlighting the stabilized aromatic 4π-cationic intermediate.^24,56

Cyclooctatetraene (3) was initially synthesized by Willstätter in 1911, sparking researchers' interest as it represented the next higher vinylogue of benzene.^61,62 Unlike benzene, cyclooctatetraene (COT) has polyolefinic properties, adopting a tub geometry in the S₀ state rather than a planar geometry.^62,63 As a 4n π system, its ground-state planar, bond-length equalized geometry (D_8h) is Huckel antiaromatic. Thus, the tub geometry is preferred to preclude the destabilizing conjugation of antiaromaticity.^42,61–63

The S₀ interconversion from 3 to 3b through a planar transition state (3-TS) was extensively studied and quantified to have a thermal barrier of approximately 10 kcal mol⁻¹.^24,62,64,65 Previous studies have verified the planar, aromatic geometry of 3-TS in the lowest triplet^66,67 and first and second singlet excited-states (S₁ and S₂)^62,67–69 excited states. Nonetheless, the impact of Baird aromaticity on the planarization energetics was unclear-3 planarizes without an activation barrier on the S₁-surface (Scheme 2).^25,70 To explore the energetics of inversion and the contributions of Baird aromaticity, Ueda et al. synthesized and examined the inversion of a methylated thiophene fused derivative of 3, MeCOTh.²⁵

Scheme 2 Proposed excited-state inversion mechanism for cyclooctatetraene (3) through a planar transition state, 3-TS, towards 3-b.^17,25,71

To provide experimental evidence for MeCOTh racemization, Ueda et al. used time-dependent circular dichroism to discover that MeCOTh could racemize thermally and photochemically.²⁵ The unirradiated samples showed no decrease in time-dependent circular dichroism signal, indicating that the reaction was light-mediated.²⁵ Time-resolved transient absorption spectroscopy identified only S₁ dynamics under λ = 365 nm irradiation, indicating selective S₁ excitation.²⁵ To evaluate the presence of Baird aromaticity in the S₁-state, MeCOTh was selectively promoted to the T₁-surface using a photosensitizer because Baird aromaticity is typically discussed in the triplet state.²⁵ This promotion facilitated a similar inversion, confirming the presence of Baird aromaticity on the S₁-surface.²⁵ The ground-state MeCOTh to MeCOTh-b inversion barrier was experimentally determined through circular dichroism decay profiles at temperatures of 40 °C, 50 °C, and 60 °C, with an activation barrier of 25.4 kcal mol⁻¹.²⁵ Applying the same method at 0 °C, 10 °C, and 20 °C to the irradiated samples resulted in the S₁ and T₁ inversion barrier values of 4.3 kcal mol⁻¹ and 4.0 kcal mol⁻¹, respectively.²⁵ This shows a decrease of 21.1 kcal mol⁻¹ and 21.4 kcal mol⁻¹ in these inversion barriers compared to their ground-state counterpart.²⁵ Ueda et al. proposed that the lower S₁ and T₁ inversion barriers stem from Baird aromatic stabilization of the planar transition state accessed during inversion.²⁵ They used quantum chemical calculations to investigate whether the photochemical planarization of MeCOTh, as seen in Scheme 3, is a consequence of Baird aromaticity.²⁵

Scheme 3 Proposed excited-state inversion mechanism for methylated thiophene fused derivative of 3, MeCOTh, through a planar transition state, MeCOTh-TS, towards the inverted geometry, MeCOTh-b.²⁵

Density functional theory calculations with B3LYP-D3(BJ)/6–311 + G(d,p)//B3LYP-D3(BJ)/6–31 G(d) determine that the S₀ inversion barrier was 29.9 kcal mol⁻¹, exceeding the experimental value by 4.5 kcal mol⁻¹.²⁵ In line with this trend, calculations using TD-B3LYP-D3(BJ)/6–311 + G(d,p)//B3LYP-D3(BJ)/6–31 G(d) found that the inversion barriers for S₁ and T₁ were 8.6 kcal mol⁻¹ and 9.0 kcal mol⁻¹, respectively; these values were over-predicted by 4.3 kcal mol⁻¹ and 5.0 kcal mol⁻¹ compared to the experimental values.²⁵ The calculated inversion barriers matched experimental findings, showing substantial decreases in activation free energies (21.3 kcal mol⁻¹ and 20.9 kcal mol⁻¹ on the S₁ and T₁ surfaces, respectively).²⁵ This decrease was attributed to Baird aromatic stabilization of the planar transition state.²⁵ These calculations produce static structures of limited value for a photochemical reaction, where subnanosecond timescales severely limit molecular equilibration. Thus, nonadiabatic molecular dynamics simulations are crucial for revealing the reaction mechanism and establishing structure-reactivity correlations. Such calculations investigate various mechanistic pathways that MeCOTh accesses as it transitions between the key mechanistic points identified by Ueda et al., enhancing their mechanistic understanding. We propose that Baird aromaticity promotes rapid inversions between MeCOTh and MeCOTh-b in the excited state and that the relative stereochemistries at the S₁/S₀ hopping points will remain in the final ground state geometries (i.e., no further S₀ inversions).

We conducted a computational analysis of the inversion mechanism of MeCOTh through multiconfigurational (complete active space self-consistent field, CASSCF) and single-reference (time-dependent density functional theory, TD-DFT) quantum chemical calculations. In this study, we comprehensively outlined the reaction pathways for the inversion of MeCOTh to MeCOTh-b using our open-source machine learning code, Python rapid artificial intelligence ab initio molecular dynamics (PyRAI²MD),^72–74 which allowed us to perform machine-learned nonadiabatic molecular dynamics (ML-NAMD).

Active space and characterization of excitation on the optimized global minimum

The first step of the photochemical inversion reaction involves photon absorption to an electronically excited state. As such, we first computed the vertical excitation energies of MeCOTh with TD-DFT calculations using Gaussian16.⁷⁵ We compared our results to previous experimental studies that used a 365 nm light source (3.4 eV) to selectively promote MeCOTh to the S₁-state.²⁵ We performed vertical excitation energy calculations and characterized the nature of electronic transitions using range-separated hybrid functionals (i.e., CAM-B3LYP⁷⁶ and ωB97X-D77) with a double-zeta basis set (i.e., aug-cc-pVDZ^78,79). The predicted vertical excitation energy for the ground-state optimized MeCOTh, referred to as MeCOTh-S₀, was found to be 3.71 eV using CAM-B3LYP-D3(BJ)/aug-cc-pVDZ and 3.74 eV with ωB97X-D/aug-cc-pVDZ. Both approaches indicated that the S₀ → S₁ excitation corresponds to a π → π* transition in the COT-core, with oscillator strengths of 0.024 and 0.025, respectively. This consistent alignment with the experimental light source and TD-DFT calculations strongly validates our computational methods.

The multiconfigurational electronic structure of the MeCOTh during photoexcitation and photodynamic simulations requires multiconfigurational quantum mechanical techniques. Therefore, we employed CASSCF and complete active space second-order perturbation theory (CASPT2), both of which require an active space. We analyzed the orbital transitions suggested by TD-DFT and those relevant to the MeCOTh to MeCOTh-b isomerization (specifically, COT-core π orbitals). An active space comprising eight electrons was formulated across seven orbitals of the COT core; we excluded the highest-lying π*-orbital to prevent a double excitation to the S₂ state and enhance the efficiency of our computations (Fig. 1).^80,81

Fig. 1 CASSCF(8,7) active space for MeCOTh consisting of eight electrons distributed within four π-orbitals and three π*-orbitals of the COT-core. We noted the average electron occupancy of each orbital below the respective orbital. An isosurface value of 0.03 was used for all orbital images.

To ensure accurate photophysical characterization using the (8,7) active space, we compared the vertical excitation energies and nature of transition predicted by TD-DFT, CASSCF, and CASPT2. The vertical excitation energy for MeCOTh-S₀ was calculated using a state average of the first five singlet states (i.e., S₀–S₄) for CASSCF (SA5-CASSCF(8,7)/ANO-S-VDZP) and applying the CASPT2 correction (CASPT2(8,7)/ANO-S-VDZP//SA5-CASSCF(8,7)/ANO-S-VDZP). We used a double-zeta basis set for all methods—aug-cc-pVDZ for TD-DFT and ANO-S-VDZP for CASSCF and CASPT2—to ensure that any observed differences between the methods stem from the treatment of electron correlation and not from basis set limitations. Both multiconfigurational methods overestimated the S₀ → S₁ vertical excitation energy, reporting values of 5.72 and 4.22 eV, respectively, compared to the TD-DFT and experimental findings. The discrepancies observed in the CAS methods arise from the relatively small active space, which does not include all conjugated π-orbitals. Both methods classify this excitation as a π → π* transition within the COT core and align with the TD-DFT results. Table 1 presents a summary of these findings.

Table 1 Benchmarked TD-DFT and CASSCF vertical excitation energies

Method	State	Energy (eV)	Wavelength (nm)	Oscillator strength	Nature
CAM-B3LYP-D3(BJ)/aug-cc-pVDZ	S1	3.71	334	0.024	π → π*
	S2	4.56	272	0.352	π → π*
	S3	4.82	257	0.139	π → π*
	S1	3.74	331	0.025	π → π*
	S2	4.60	270	0.376	π → π*
	S3	4.91	253	0.132	π → π*
SA5-CASSCF(8,7)/ANO-S-VDZP	S1	5.72	217	0.011	π → π*
	S2	6.62	187	0.112	π → π*
	S3	6.77	183	0.136	π → π*
SA5-CASPT2(8,7)/ANO-S-VDZP//SA5-CASSCF(8,7)/ANO-S-VDZP	S1	4.22	294		π → π*
	S2	5.64	220		π → π*
	S3	5.74	216		π → π*

---

After we compared the orbital transitions, vertical excitation energies, and oscillator strengths between CASSCF, CASPT2, and TD-DFT, we concluded that CAS(8,7), as shown in Fig. 1, appropriately captures the vertical excitation energies of MeCOTh at an equilibrium geometry on the ground-state (i.e., MeCOTh-S₀). We used SA5-CASSCF(8,7)/ANO-S-VDZP for all subsequent calculations.

Predicted absorption spectra and S₁ minimum energy path

To verify that the (8,7) active space captures the spectral properties of MeCOTh at nonequilibrium geometries, we compute the absorption spectrum of MeCOTh; the computed and experimental (Ueda et al.²⁵) spectra are compared. To that end, we generate 500 Wigner-sampled non-equilibrium structures based on the frequencies of MeCOTh-S₀ (Fig. 2a) and computed the vertical excitation energies [S₀ → S_n (n = 1–4)] and oscillator strengths for each structure with MS-CASPT2(8,7)/ANO-S-VDZP//SA(5)-CASSCF(8,7)/ANO-S-VDZP (Fig. 2b).⁸² The reported experimental spectrum²⁵ spanned a range of 200 to 500 nm, showing peaks at approximately 270 nm and 230 nm, with notable absorption at 300 nm that drops off to a low-absorption tail beyond 350 nm.²⁵ Therefore, we only plotted the peaks above 200 nm for direct comparison with the experimental absorption spectrum. The intensities of the calculated spectrum are normalized using the S₀ → S₂ transition. All experimental investigations of the MeCOTh → MeCOTh-b inversion dynamics utilized a light source of 365 nm to selectively excite MeCOTh to the S₁ surface using the longest wavelength absorption band.²⁵

Fig. 2 (a) The top view (top) and side view (bottom) overlay of 500 non-equilibrium structures from the Wigner ensemble of MeCOTh. (b) The predicted absorption spectrum was generated by aggregating the 500 vertical excitation energy calculations. The absorption intensity is normalized to the peak with the highest oscillator strength, S₃. All calculations utilized the MS-CASPT2(8,7)/ANO-S-VDZP//CASSCF(8,7)/ANO-S-VDZP method.⁸²

The computed spectrum (Fig. 2b) shows three absorption peaks centered at 317 nm, 253 nm, and 226 nm, corresponding to π → π* transitions. The S₀ → S₃ transition is represented by the yellow peak centered at 226 nm, which is the lowest-intensity peak. The blue peak, centered at 253 nm, corresponds to the S₀ → S₂ transition and has the highest intensity of all three peaks (i.e., the bright state). The red peak centered at 317 nm has the second-highest intensity and corresponds to the lowest-energy transition, S₀ → S₁; it features a low-absorption tail extending above 350 nm and diminishing completely by 400 nm. While the intensity of S₀ → S₂ is higher than the intensity of the S₀ → S₁ transition, the S₀ → S₂ transition cannot be accessed using the experimental light source (365 nm) because the excitation energy required exceeds the energy of the photon source.

To assess the accuracy of our computed spectrum, we matched the observed peaks in the computed spectrum to those observed in the experimental spectrum; matching peaks and spectral shapes between the experimental and computed spectra would indicate that our electronic structure method (MS-CASPT2(8,7)/ANO-S-VDZP//CASSCF(8,7)/ANO-S-VDZP) captures the spectral properties of MeCOTh correctly. The highest-intensity peak in the experimental spectrum is at approximately 230 nm, whereas the highest-intensity peak in the computed spectrum is at 253 nm (i.e., the S₀ → S₂ peak). This indicates that the computed spectrum is redshifted by 23 nm relative to the highest-intensity peak of the experimental spectrum. The second-highest-intensity peak in the experimental spectrum is at approximately 270 nm, whereas in the computed spectrum it is at 317 nm (i.e., the S₀ → S₁ peak). Like the S₀ → S₂ peak, the corresponding S₀ → S₁ peak is also red-shifted (47 nm) in the computed spectra relative to the peak observed experimentally.²⁵ The experimental spectrum shows significant absorption between the two peaks (i.e., between 230 and 270 nm).²⁵ This region matches the calculated overlap region that we observed between the S₀ → S₁ and S₀ → S₂ transition peaks (250 to 300 nm). After the peak at 270 nm, the experimental spectrum exhibits a shoulder region with significant absorption that drops into a low-absorption tail beyond 350 nm, which is used to selectively excite MeCOTh to the S₁ state.²⁵ Our calculated S₀ → S₁ peak at 317 nm has a low-absorption tail extending beyond 350 nm, matching the tail observed in the experiment.²⁵ The experimental spectrum has decreased absorption at higher wavelengths (i.e., < 230 nm). In our computed spectra, the S₀ → S₃ peak has the lowest intensity and shows significant overlap with the S₀ → S₂ peak; both S₀ → S₂ and S₀ → S₃ peaks decrease in intensity at higher wavelengths. Based on these results, we conclude that the calculated spectrum has peaks that are consistently redshifted (i.e., by +47 nm and +23 nm for S₀ → S_n (n = 1–2), respectively) compared to those present in the experimental spectrum reported by Ueda et al.²⁵ The discrepancy between λ_max and intensities can be attributed to the limitations of CASPT2 (e.g., a MAE of 0.11 eV),^83,84 basis set effect,^84,85 and to a relatively small active space that does not include all conjugated π-orbitals.

After assessing the spectral properties, we shifted our focus to elucidating the mechanism of the photochemical reaction. Given that the S₁ state is the only energetically accessible excited state under the experimental light source,²⁵ we calculated the minimum energy path (MEP) to identify the steepest descent path along the S₁ surface, starting from the Franck–Condon point of MeCOTh-S₀. We hypothesize that MeCOTh will adopt a more planar structure along S₁-MEP, with uniform bond lengths in the COT core (i.e., π_CC and σ_CC) because the core is Baird aromatic on the S₁-state, favoring a planar confirmation for maximal π-orbital overlap. Fig. 3 shows the S₁-MEP, the MeCOTh-S₀ geometry, the final MEP step geometry (MeCOTh-MEP-11), and an S₁ minimum geometry (MeCOTh-S₁).

Fig. 3 (a) The calculated MEP of MeCOTh along the S₁ surface at the SA(5)-CASPT2(8,7)/ANO-S-VDZP level. (b) The defined inversion angle, θ_inv, includes C₁–C₉–C₄–C₆, which are highlighted in green. We used θ_inv to monitor the planarization of the COT core. The bottom structures correspond to (c) the ground state, optimized SA(5)-CASPT2(8,7)/ANO-S-VDZP, MeCOTh-S₀, (d) the final step of MEP, MeCOTh-MEP-11, and (e) the optimized S₁ minimum, MeCOTh-S₁. All energies are relative to the ground-state-optimized geometry, MeCOTh-S₀.

Fig. 3a illustrates the MEP along the S₁ surface and the energies corresponding to the S₀–S₄ states; it includes 11 geometries. The S₁-MEP converges to a structure 3.11 eV (MeCOTh-MEP-11) above the MeCOTh-S₀. At MeCOTh-MEP-11, the S₁/S₀ energy gap is 2.74 eV; this large S₁–S₀ energy gap provides a direct path to a radiative decay channel. Fig. 3c–e show three geometries: MeCOTh-S₀, MeCOTh-MEP-11, and MeCOTh-S₁. We defined a planarity parameter, θ_inv, within the COT core to quantify the change in planarity within these three geometries and along the MEP (Fig. 3b); at a perfectly planar geometry, θ_inv would be 180°. MeCOTh-S₀ and MeCOTh-MEP-11 have θ_inv values of 132° and 139°, respectively (Fig. 3c and d). There is an increase of 7° in θ_inv from MeCOTh-S₀ → MeCOTh-MEP-11; the rise in θ_inv indicates increased planarity in MeCOTh along the S₁ surface as it approaches 180°. We assessed the bond lengths of π_CC and σ_CC in the COT core to investigate the bond length equalization that occurs concurrently with aromatic systems. The four π_CC measured 1.38 Å, 1.37 Å, 1.35 Å, and 1.37 Å while σ_CC bonds measured 1.47 Å, 1.48 Å, 1.47 Å, and 1.47 Å at MeCOTh-S₀. At MeCOTh-MEP-11, the π_CC measured 1.41 Å, 1.43 Å, 1.42 Å, and 1.41 Å, while the σ_CC measured 1.39 Å, 1.41 Å, 1.40 Å, and 1.41 Å. We find that all π_CC bond lengths have increased while all σ_CC bond lengths have decreased. This bond equalization coincides with an increase in planarization. Both geometric changes suggest that the electronic structure of MeCOTh-S₀ becomes Baird aromatic, consistent with our initial hypothesis.

We optimized an S₁ minimum, MeCOTh-S₁, using the final MEP geometry as an input (Fig. 3e). MeCOTh-S₁ is 2.65 eV above MeCOTh-S₀ and has a θ_inv value of 118°. This indicates a decrease of 14° and 21° from MeCOTh-S₀ and MeCOTh-MEP-11, respectively. The decreased θ_inv corresponds to a less planar structure than MeCOTh-S₀ or MeCOTh-MEP-11. We then measured the π_CC and σ_CC bond lengths to directly compare to MeCOTh-S₀ and MeCOTh-MEP-11. In MeCOTh-S₁, the four π_CC bonds measured 1.44 Å, 1.45 Å, 1.44 Å, and 1.45 Å, while the four σ_CC bonds measured 1.39 Å, 1.39 Å, 1.37 Å, and 1.40 Å. The elongated π_CC and shortened σ_CC indicate a double bond shift within the COT core. The decrease in θ_inv and the double bond shift resulted in a conformational change from a tub geometry in MeCOTh-S₀ and MeCOTh-MEP-11 to a boat–boat conformation⁸⁶ in MeCOTh-S₁. After S₁ optimization, the S₁/S₀ gap in MeCOTh-S₁ is 0.30 eV. The small S₁/S₀ gap suggests that the steepest path of descent is nonradiative, but it does not specify whether a relative stereochemical inversion occurs. While the MEP calculation provides detailed information on a single, steepest S₁ pathway, it omits dynamic information about the inversion mechanism.

We addressed the omitted dynamical effects by performing ML-NAMD simulations based on multiconfigurational quantum mechanical calculations. We generated 988 initial conditions via the Wigner sampling algorithm for production trajectories based on the frequencies of the MeCOTh-S₀. We captured nonadiabatic transitions (i.e., surface hops) using the fewest-switches surface hopping (FSSH) algorithm⁸⁷ for 5.5 ps with curvature-driven time-derivative coupling (kTDC) to evaluate nonadiabatic coupling based on the energy gaps predicted by the neural networks (NNs).^88,89 After 5.5 ps, 947 (96%) trajectories were on the S₀-state, while 41 (4%) were on the S₁-state. Based on the findings by Ueda et al., we hypothesize that longer excited-state lifetimes may lead to increased stereochemical inversions, owing to the reduced S₁-state inversion barrier. The increased inversions would directly affect the relative stereochemistries of the final geometries of the trajectories. In Fig. 4, we analyzed the relative stereochemistries of all trajectories by monitoring θ_inv. We tracked the θ_inv values until the final S₁/S₀ crossings, extending each trajectory by 500 fs on the S₀ surface (i.e., after the S₁/S₀ hop). This 500 fs window was set to ensure that the trajectories completely adopted their corresponding relative stereochemistries after the last S₁ → S₀ hop (indicated by black dots), preventing ground-state inversions. Thus, the total length of our trajectories varies and is directly related to their respective excited-state lifetime. After the 500 fs duration on the S₀ surface, we classified the trajectories according to their final geometry and plotted their respective θ_inv as a function of simulation time (ps) in Fig. 4.

Fig. 4 Trajectory map (θ_inv vs. simulation time) of MeCOTh → MeCOTh-b inversion. We separated the trajectories based on the geometry 500 fs after the S₁/S₀ crossing (a) retained, MeCOTh, and (b) inverted, MeCOTh-b. Black dots represent S₁/S₀ surface hopping points.

The optimized geometry, MeCOTh-S₀, has a θ_inv of 132°. Thus, we labeled structures with θ_inv < 180° as having retained relative stereochemistries (MeCOTh) when compared to MeCOTh-S₀. Structures with θ_inv > 180° are labeled as having inverted relative stereochemistries (MeCOTh-b) when compared to MeCOTh-S₀. We observed that 304 trajectories led to MeCOTh (Fig. 4a), while 643 led to MeCOTh-b (Fig. 4b). Trajectories leading to MeCOTh have an average simulation time of 1.7 ps (Fig. 4a), while trajectories leading to MeCOTh-b have an average simulation time of 1.2 ps (Fig. 4b). On average, trajectories leading to MeCOTh were on the S₁ surface 0.5 ps longer than those leading to MeCOTh-b. The difference in simulation time demonstrates that longer excited-state lifetimes increase the probability of inversions; at least two inversions are required to revert to MeCOTh.

We defined two torsional regions in Fig. 4 that represent the relative stereochemistries of MeCOTh (θ_inv < 180°) and MeCOTh-b (θ_inv > 180°). Trajectories in which the θ_inv values transition from less than 180° to over 180°, or vice versa, undergo a relative stereochemistry inversion. To quantify the relative stereochemical inversions in the excited state, we measured the average simulation time—directly linked to the excited-state lifetime—of the trajectories shown in Fig. 4a and b, along with the θ_inv transitions across these regions. We expected either no inversions or an even number of stereochemical inversions (i.e., θ_inv crossing from and to 180°) in Fig. 4a for trajectories leading to MeCOTh; a single or an odd number of relative stereochemical inversions is expected for trajectories that yield MeCOTh-b (Fig. 4b). In Fig. 4a, we observe that trajectories had an average simulation time of 1.7 ps and a maximum of 6 relative stereochemical inversions (i.e., θ_inv crossing from and to 180°). In Fig. 4b, we observe that trajectories had an average simulation time of 1.2 ps and a maximum of 3 relative stereochemical inversions (i.e., θ_inv crossing from and to 180°). From these findings, we concluded that a longer duration in the S₁ state promotes more inversions, as shown in Fig. 4a, where trajectories leading to MeCOTh require 2 or more inversions and remain, on average, 0.5 fs longer in the excited state. In contrast, trajectories that lead to MeCOTh-b remain on the excited state for shorter durations and have fewer inversions. These results are consistent with Ueda et al., who found that the barrier for relative stereochemical inversion in the excited state is 4.3 kcal mol⁻¹, due to Baird aromaticity, which facilitates the relative stereochemical inversion.²⁵

Hopping point analysis

We investigate the retention of relative stereochemistry and the prevalence of S₀-inversions by comparing the relative stereochemistries at S₁/S₀ hopping points and final geometries. We anticipate that the relative stereochemistries at the S₁/S₀ hopping points will align with the relative stereochemistries of the final geometries of the trajectories. (i.e., no S₀-state inversion). We plotted the S₁/S₀ hopping points against time in Fig. 5. A trajectory analysis reveals that MeCOTh can undergo one or more surface hops; thus, for clarity, we restricted the depicted hopping points to the final S₁ → S₀ hop, totaling 947 hopping points in Fig. 5. We use θ_inv to assess the relative stereochemistries of the S₁/S₀ hopping points.

Fig. 5 We implement color coding for the hopping points according to their relative stereochemistry 500 fs after S₁ → S₀ relaxation: orange denotes the relative stereochemistry of MeCOTh, and blue represents the relative stereochemistry of MeCOTh-b.

In Fig. 5, we identify three clusters of hopping points: 0–0.5 ps, 0.501–1 ps, and above 1.01 ps, which we have labeled based on their timing as “early,” “intermediate,” and “late.” The early cluster comprises 103 hopping points in the torsional region pertaining to the relative stereochemistry of MeCOTh (θ_inv < 180°). Within this cluster, 95 (92%) of the hopping points belong to trajectories leading to MeCOTh-b, while 8 (8%) are associated with trajectories that yield MeCOTh. The relative stereochemistry of the main product, MeCOTh-b (θ_inv = 228°), differs from the relative stereochemistry associated with the torsional region occupied by the early cluster (i.e., θ_inv < 180°). This finding contradicts our hypothesis; we expected that the relative stereochemistries of S₁/S₀ hopping points would align with those of the final geometries. The difference in relative stereochemistries between the S₁/S₀ hopping points and the final geometries within this cluster results from the S₁/S₀ crossings having increased momentum, which leads to the second relative stereochemical inversion of MeCOTh on the S₀ state (i.e., no retention of relative stereochemistries; Fig. S7).

The intermediate cluster includes 613 hopping points within the torsional region corresponding to MeCOTh-b (θ_inv > 180°). The intermediate cluster has 510 hopping points (83%) that belong to trajectories that lead to MeCOTh-b, while the remaining 103 hopping points (17%) are associated with trajectories yielding MeCOTh. In contrast to the early cluster, the relative stereochemistry of the major product from the intermediate cluster, MeCOTh-b (θ_inv = 228°), aligns with the relative stereochemistry associated with the conformational region occupied by that cluster (i.e., θ_inv > 180°). This finding indicates that the main pathway within the intermediate cluster leads to the relative stereochemistries at the S₁/S₀ hopping points remaining consistent in the final structures, supporting our initial hypothesis and previous work.²⁵

The late cluster contains 231 hopping points in the torsional region of MeCOTh (θ_inv < 180°). Within this cluster, 38 hopping points (16%) belong to trajectories leading to MeCOTh-b, and 193 hopping points (84%) belong to trajectories leading to MeCOTh. The main product of the late cluster, MeCOTh (θ_inv = 132°), is located within the torsional region occupied by this cluster. This finding indicates that the main pathway within the late cluster leads to conserving the relative stereochemistries of the S₁/S₀ hopping points. We conclude that the preserved relative stereochemistries between the S₁/S₀ hopping points of the intermediate and late clusters and the final geometries of their trajectories result from the MeCOTh fully adopting the relative stereochemistries of S₁/S₀ hopping points on the S₁ state and transitioning to the S₀ surface afterward. Consequently, the inversion barrier increases to 25.4 kcal mol⁻¹,²⁵ preventing ground-state reinversion. Relative stereochemistry conservation is absent in the trajectories that undergo S₁ → S₀ crossings within the early cluster, as trajectories that undergo early hopping have increased momentum relative to trajectories that undergo surface hopping at later timesteps (i.e., intermediate and late clusters) due to the initial relaxation from the Franck–Condon region (Fig. S8), affording the S₀ inversion from MeCOTh to MeCOTh-b.

To investigate the structural diversity within the clusters of surface-hopping points, we optimized three minimum-energy conical intersections (MECIs) for each relative stereochemistry (i.e., MeCOTh and MeCOTh-b), using the S₁/S₀ hopping points as initial guess structures. We hypothesize that the optimized MECIs of trajectories leading to MeCOTh and MeCOTh-b will be stereoisomers of one another, and that trajectories with S₁/S₀ hopping points exhibiting opposite relative stereochemistries compared to their final geometries might optimize to different MECI than those with matching relative stereochemistries. To test these hypotheses, we use θ_inv to assess the relative stereochemistries of the S₁/S₀ hopping points. Fig. 6 shows the differences in geometries and θ_inv for three representative MECIs- MeCOTh-MECI-early, MeCOTh-MECI-intermediate, and MeCOTh-MECI-late- and the three previously defined mechanistic critical points.

Fig. 6 We compared an SA5-CASSCF(8,7)/ANO-S-VDZP optimized minimum energy conical intersections from each cluster: MeCOTh-MECI-early, MeCOTh-MECI-intermediate, and MeCOTh-MECI-late to MeCOTh-S₀, MeCOTh-MEP-11, MeCOTh-S₁. Each structure has their inversion angle, θ_inv, below.

All optimized MECIs from S₁/S₀ hopping points in the early and late clusters have a θ_inv of 111°, regardless of the final geometry (MeCOTh or MeCOTh-b). This θ_inv value indicates that the relative stereochemistries of these MECIs correspond to MeCOTh. In contrast, all MECIs optimized from S₁/S₀ hopping points in the intermediate cluster had θ_inv = 249°, irrespective of the relative stereochemistries of the final geometry. This θ_inv value indicates that these MECIs have relative stereochemistries matching MeCOTh-b. Regardless of the cluster or the final relative stereochemistry, all MECIs remained at 3.25 eV above MeCOTh-S₀ (see Table S1 for all optimized MECIs).

To assess whether the optimized MECIs are structurally equivalent, we compared their bond lengths (π_CC and σ_CC) within the COT core and θ_inv to each other, as well as to other key mechanistic points (i.e., MeCOTh-S₀, MeCOTh-MEP-11, and MeCOTh-S₁; Fig. 6). In Fig. 6, MeCOTh-S₀ and MeCOTh-MEP-11 have tub geometries with θ_inv values of 132° and 139°, respectively. In contrast, MeCOTh-MECI-early, MeCOTh-MECI-intermediate, MeCOTh-MECI-late, and MeCOTh-S₁ have boat–boat geometries with θ_inv values of 111°, 249°, 111°, and 118°, respectively. All π_CC/σ_CC bond lengths remain consistent in all optimized MECIs (Table S1). The consistent π_CC/σ_CC bond lengths, boat–boat conformation, and comparable energies across all optimized MECIs imply that MECIs within each cluster are the same. Furthermore, the MECIs optimized from the early/late clusters are stereoisomers of those from the intermediate cluster. The similar π_CC/σ_CC bond lengths and shared conformation between MECIs and MeCOTh-S₁ indicate these structures occupy similar regions of the potential energy surface. The findings show no alternative MECI in the clusters for trajectories with opposite relative stereochemistries between the S₁/S₀ hopping points and their final geometries. Therefore, an S₀ inversion must take place in these trajectories.

Fig. 4 and 5 illustrate θ_inv against simulation time in picoseconds. We observed that all trajectories increase from θ_inv = 132° (MeCOTh) to θ_inv = 228° (MeCOTh-b) within 1 ps. This increase in θ_inv suggests that within the trajectories of MeCOTh → MeCOTh-b, a structure must exist where θ_inv = 180° (i.e., a planar cyclooctatetraene core). This illustration of the ensemble of ML-NAMD trajectories (Fig. 4) agrees with the MEP (Fig. 3a); the MeCOTh planarizes. We hypothesized that the dominant inversion mechanism is facilitated by a planar structure, which is supported by the findings of Ueda et al.,²⁵ our MEP (Fig. 3a), and the ML-NAMD trajectory maps (Fig. 4).

We defined new geometric parameters about the COT core to capture the planarity more accurately due to π_CC-torsions (Fig. 7). We define angles θ_1–4 to characterize the torsions resulting from the concomitant rotation of any two thiophenes, while angles θ_5–8 are assigned to capture the torsions within the thiophenes. In an ideal planar geometry, all angles (θ_1–8) should fall within ± 20° of the angles that represent planarity (i.e., 0°/360° or 180°).

Fig. 7 The geometrical parameters used to define an absolute planar geometry. We used (a) four dihedrals between the thiophenes (θ_1–4) and (b) four dihedrals within the thiophenes (θ_5–8) for eight dihedrals. If all the dihedrals are ± 20° from planarity (i.e., 0°/360° or 180°), we classify the geometry as planar.

Based on these geometrical thresholds (i.e., ± 20° from 0°/360° and 180°), we observe that 8% (71) of trajectories access planar geometries. In Fig. 8, we investigate a representative trajectory that undergoes the relative inversion mechanism via the planar geometry, as proposed by Ueda et al.²⁵ in more detail.

Fig. 8 A representative trajectory demonstrating the MeCOTh to MeCOTh-b inversion, driven by the previously hypothesized planar geometry. Snapshots are illustrated at 0 fs, 195 fs, 337 fs, and 883.5 fs during the simulation. The distance between two sulfur atoms from neighboring thiophene rings (d_SS) and the planarity angle (θ_inv) are emphasized below each structure.

We utilized angles θ_1–8 to capture the planar structures systematically; to ensure visual clarity in Fig. 8, we emphasize the distance between the two sulfur atoms, d_SS, and θ_inv. At the anticipated planar geometry, d_SS is expected to be a minimum, while θ_inv approaches 180°. At 0 ps, d_SS is 3.69 Å and θ_inv is 137°. By 195 fs, the trajectory has a d_SS of 3.24 Å and θ_inv of 157°. At 337 fs, d_SS = 3.08 Å and θ_inv = 175°. At 883.5 fs, MeCOTh shows a d_SS of 3.32 Å and θ_inv at 250°. Between 0 fs and 195 fs, d_SS decreased by 0.45 Å, and θ_inv increased by 20°, indicating a planarizing structure. This dynamical information aligns with the MEP calculation; the MeCOTh on the S₁-surface planarizes due to Baird aromaticity. Between 196 fs and 337 fs, d_SS decreased by 0.16 Å as θ_inv increased by 18°. This suggests that the planarization observed from 0 fs to 195 fs persists throughout the 196 fs to 337 fs interval; this resulted in a significantly planar structure at 337 fs. Between 338 fs and 883.5 fs, d_SS increased by 0.24 Å while θ_inv rose by 75°. This suggests that inversion occurs after the trajectory assumes a planar conformation, as d_SS rises concurrently with θ_inv approaching the torsional region that indicates the relative stereochemistry of MeCOTh-b. At the 883.5 fs geometry, the boat–boat conformation facilitates the S₁ → S₀ surface hopping event.⁷⁰ We conclude that these trajectories undergo a direct inversion, passing through a planar structure resulting from the Baird aromatic COT core. We expected to find these results replicated throughout our ML-accelerated NAMD simulations. However, our ML-NAMD trajectories indicate that other pathways outcompete it.

We identified a reaction pathway that involves a new geometry we labeled “crown,” characterized by the unidirectional orientation of all sulfur atoms (Fig. 9a). We used the defined structural parameters, including θ_inv, θ_1–8, and d_SS, to identify the crown-like geometries from the trajectories automatically. We illustrate the geometrical parameters in Fig. 9b–d.

Fig. 9 Geometrical parameters that define the (a) crown conformation. To systematically capture this geometry in our ML-NAMD trajectories, we tracked the following: (b) the inversion angle (θ_inv), sulfur–sulfur distance (d_SS), (c) four dihedral angles between the thiophenes (θ_1–4), and (d) four dihedral angles within the thiophenes (θ_5–8). A structure is classified as adopting the crown geometry when the following thresholds are met: θ_inv must deviate from planarity (outside ±20° from 180°), θ_1–4 should be near-planarity (within ±20° of 0°, 180°, or 360°), θ_5–8 should deviate from planarity (outside ±20° of 0°, 180°, or 360°), and at least one sulfur–sulfur distance must be less than 3.0 Å.

62% (588) of trajectories undergo the MeCOTh → MeCOTh-b inversion through a crown geometry. In Fig. 10, we investigate a representative trajectory that underwent the relative stereochemical inversion via the crown geometry in more detail.

Fig. 10 A representative trajectory demonstrating the MeCOTh to MeCOTh-b inversion, driven by the newly proposed crown geometry (t = 370 fs). Key geometries are illustrated at 0 fs, 195 fs, 370 fs, 486 fs, and 682.5 fs during the simulation. d_SS and θ_inv are emphasized below each structure.

We report the distance d_SS and θ_inv to ensure a direct comparison with Fig. 8. At 0 ps, d_SS is 3.64 Å, and θ_inv is 128°. By 195 fs, the trajectory has a d_SS of 3.18 Å and θ_inv of 140°. At 370 fs, d_SS measures 2.95 Å, while θ_inv is 151°. By 486 fs, d_SS is 3.26 Å, and θ_inv is 180°. Finally, at 682.5 fs, d_SS is 3.64 Å and θ_inv is 264°, respectively. Between 0 fs and 195 fs, d_SS decreased by 0.46 Å while θ_inv increased by 12°, indicating an increasingly planar structure within the first 200 fs of the simulation. These structural changes resemble those seen in the trajectory that inverts via the planar geometry (Fig. 10), demonstrating MeCOTh planarizing on the S₁-surface because of Baird aromaticity. From 196 fs to 370 fs, d_SS decreased by 0.23 Å, and θ_inv increased by 11°. The persistent decline in d_SS and increase in θ_inv indicate further planarity; however, by 370 fs, MeCOTh adopted a crown geometry. We attribute the crown geometry to a balance between the closed-shell repulsions between neighboring sulfurs and the stabilizing effects of Baird aromaticity. d_SS rises by 0.31 Å, while θ_inv increases by 29°, from 371 fs to 486 fs. Indeed, increased d_SS indicates that inversion occurs after the crown structure as the sulfur atoms (i.e., thiophene rings) separate. Within this interval, θ_inv increased to 180°, suggesting a planar structure, similar to the S₁ minimum energy structure reported by Ueda et al.²⁵ However, the non-planar characteristics of the 486 fs geometry imply that trajectories can invert their relative stereochemistries (i.e., pass θ_inv = 180°) without planarizing. Between 487 fs and 682.5 fs, d_SS increased by 0.38 Å and θ_inv by 84°, illustrating the completion of the relative stereochemical inversion of MeCOTh to MeCOTh-b. By 682.5 fs, MeCOTh exhibits a boat–boat conformation that facilitates the S₁ → S₀ surface hopping event.⁷⁰ We conclude that the mechanism illustrated in Fig. 10 is the preferred mechanistic pathway for the MeCOTh → MeCOTh-b inversion (62% of the trajectories).

Conclusion

This study provides a complete enumeration of the mechanistic pathways for the relative stereochemical inversion of MeCOTh to MeCOTh-b, while also assessing the role of Baird aromaticity in the mechanism. We employed multiconfigurational quantum mechanical computations, combined with ML-NAMD, to provide mechanistic insights into this relative stereochemical inversion in static and dynamic perspectives. We addressed the computational challenges associated with the SA5-CASSCF(8,7)/ANO-S-VDZP NAMD simulations by training multilayer perceptron NNs using our open-source software, PyRAI²MD, which enabled us to accelerate the photodynamic simulations by five orders of magnitude. We utilized our trained NNs to conduct ML-NAMD simulations and outlined the excited-state mechanistic pathways of MeCOTh → MeCOTh-b.

Our static, minimum energy path calculation led to a more planar structure (MeCOTh-MEP-11) relative to the optimized S₀-state geometry (MeCOTh-S₀); albeit, both had a tub conformation. We optimized MeCOTh-MEP-11 to an S₁-minimum with a boat–boat conformation (MeCOTh-S₁) featuring a small S₁/S₀ gap (0.30 eV). Our static calculation results show that the steepest-descent path on S₁ leads to a nonradiative decay channel, but it does not involve the relative stereochemical inversion. The discrepancy between experimental and computational results suggests that our static calculations lack essential dynamical effects. Our ML-NAMD simulations provide critical mechanistic insights and predict quantum yields of 68% for MeCOTh-b and 32% for MeCOTh. The ML-NAMD calculations demonstrate that 74% of trajectories undergo relative stereochemical inversion entirely on the S₁ surface, consistent with the lower S₁-state inversion barrier of 4.3 kcal mol⁻¹. In contrast, only 25% of trajectories undergo relative stereochemical inversions on the ground state, aligning with the higher S₀-state inversion barrier of 25.4 kcal mol⁻¹. All productive trajectories exhibit planarizing COT cores, consistent with MEP findings and previous studies that emphasize the role of Baird aromaticity in this reaction mechanism. We identified two competing pathways for inversion, influenced by differing electronic effects: Baird aromaticity and S–S closed-shell repulsions. In 62% of the trajectories, the inversion mechanism is mainly influenced by the balance between S–S closed-shell repulsions and Baird aromaticity, which directs the pathway towards relative stereochemical inversion via a crown geometry characterized by unidirectional sulfurs. Conversely, in 8% of the trajectories, Baird aromaticity overcomes S–S closed-shell repulsions, steering the trajectory towards relative stereochemical inversion through a planar geometry. Thus, the dominant inversion pathway passes through the newly labeled crown geometry on the S₁ surface rather than the previously proposed planar geometry.

While previous studies have focused on the aromaticity reversal and subsequent planarization of the COT-core of MeCOTh and other COT-derivatives as the main drivers of their aggregate-induced emission,⁴⁸ polymerization control,⁹⁰ and thermosalient effect,⁹¹ our findings suggest a competing driver for these macroscopic effects: the balance between S–S closed-shell repulsions and Baird aromaticity. This competing mechanistic pathway illustrates a potential competing driver for these macroscopic effects that may impede intramolecular motions,⁴⁸ disrupt hydrogen bonding networks,^90,91 and interfere with beneficial π-stacking,⁹¹ within these COT-based materials.

One of the main findings in the 2017 report by Ueda et al. was the temperature-dependent racemization of MeCOTh. Although exploring the temperature dependence of the reaction mechanism presented here is beyond the scope of this manuscript, we argue that our reported mechanism is temperature dependent. We hypothesize that an increase in temperature would provide sufficient thermal energy to facilitate the adoption of MeCOTh-crown, thus facilitating racemization. This remains to be explored in future studies.

Computational methods

Multiconfigurational methods

We performed multiconfigurational calculations using a state-averaged complete active space self-consistent field (SA-CASSCF) method in OpenMolcas 19.11.⁹² This approach utilizes the format SA(N)-CASSCF(m,n), where N signifies the number of singlet states averaged in the calculation. The variables m and n denote the number of electrons and orbitals in the active space. To be conscientious of the computational cost associated with CASSCF methodology, we aimed to construct the smallest possible active space that captured the photophysical properties of MeCOTh. To this end, we selected an active space consisting of eight electrons within seven orbitals (i.e., 4 π-orbital and 3 π*-orbital; Fig. 1). We optimized minima on the S₀ and S₁ surfaces, performed vibrational analysis to verify the stationary points, predicted the absorption spectra, and computed a minimum-energy path using SA5-CASSCF(8,7) with the ANO-S-VDZP basis set. We employed complete active space second-order perturbation theory to account for dynamic correlations in the minimum-energy path calculation. We applied multistate complete active space second-order perturbation theory to predict absorption spectra. To demonstrate the need for machine learning to study the inversion mechanism of MeCOTh, we benchmarked a single-step quantum-mechanical nonadiabatic molecular dynamics simulation of MeCOTh using SA5-CASSCF(8,7)/ANO-S-VDZP, which took 3.23 hours. Performing a single 1 picosecond QM-NAMD simulation, involving a series of 2000 single-point and gradient calculations, would take approximately 269 days.

Singe-reference method

We used density functional theory (DFT) to optimize the ground-state global minima of MeCOTh using the CAM-B3LYP-GD3(BJ)/6-31G(d,p). We used time-dependent density functional theory (TD-DFT) to calculate the excitation energies for the first 10 excited states of MeCOTh. The TD-DFT calculations were run using the CAM-B3LYP-GD3(BJ) and ωB97X-D functionals with the aug-cc-pVDZ basis set. All DFT calculations were performed using Gaussian 16.⁷⁵

Training data generation and NN training

We generated the initial training dataset in two sections: Wigner sampling and geometric interpolation. We performed Wigner sampling based on the frequencies of the ground-state, SA5-CASSCF(8,7)/ANO-S-VDZP optimized geometry of MeCOTh to generate 160 structures at the zero-point energy level. The geometric interpolation produced 8 structures, spanning from the Franck–Condon point of MeCOTh to an S₁/S₀ minimum energy conical intersection (MECI), and an additional 8 structures from the S₁/S₀ MECI to the inverted relative stereochemical product, MeCOTh-b; totaling 16 intermediate structures. We utilized the nuclear displacements from the 160 Wigner-sampled structures to adjust the interpolated geometries, resulting in 2560 initial data points. A SA5-CASSCF(8,7)/ANO-S-VDZP single-point calculation was performed on all data points to obtain their respective S₀ and S₁ energies and gradients. This data set was used to train our initial neural network (NN) potentials.

We employ the TensorFlow/Keras API in Python to construct and train the NN potential.⁹³ We constructed a fully connected feedforward multilayer perceptron NN that utilizes leaky softplus activation functions. The NN calculates the inverse distance of the input molecule to predict energies and forces. We simultaneously train the NN using energy and forces to maintain their physical relationship. The learning rate is set to 10⁻³ and incorporates a scheduler that reduces it to 10⁻⁴ and 10⁻⁵ when the validation loss plateaus. The dataset is split into training and validation sets with a 9 : 1 ratio. Additional information can be found in the SI.

To expand the initial training set with undersampled data, we used a committee model of two NNs (Table S2) to propagate 150 trajectories from the S₁ state for 10 ps with a step size of 0.5 fs. The standard deviation in the predicted energy and gradients by the NN committee served as a measure of uncertainty for the current prediction. The trajectories were halted when the standard deviation surpassed the energy and gradient thresholds. Based on our experience, we established the thresholds for energy and gradient at 0.04 hartree and 0.25 hartree·Bohr⁻¹, respectively, to effectively detect uncertain structures and enhance the NN potential. The final geometries of the halted trajectories were recalculated using SA5-CASSCF(8,7)/ANO-S-VDZP and added to the initial training set. Subsequently, adaptive sampling retrained the committee model using the updated dataset and restarted the trajectories. The process was repeated iteratively until no new structures emerged; this led to 52 iterations and a data set with 4207 total data points. The mean absolute errors in the predicted energies were 0.045 eVs and 0.046 eVs.

ML-photodynamics simulations

We used the NNs from the latest iteration of adaptive sampling (iteration 52) to propagate 988 ML-NAMD trajectories for 5.5 ps with a 0.5 fs time step. The probabilities of nonadiabatic electronic transitions were computed with Tully's FSSH method,⁸⁷ where we used the curvature-driven time derivative coupling (kTDC) to evaluate the nonadiabatic couplings based on the energy gaps predicted by the NNs.^88,89 The kTDC method showed accuracy when compared to ground-truth NAC values derived from quantum mechanical calculations.^88,89 It was especially effective for small energy gaps (e.g., <0.5 eV).⁸⁸

Author contributions

Conceptualization and methodology: SAL; analysis and investigation: CS and SAL; writing – original draft: CS; review and editing: SAL and CS; supervision, project administration, and funding acquisition: SAL.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The ML-photodynamics simulation code is open-sourced and released at: https://github.com/mlcclab/PyRAI2MD-hiam.

The NN models and initial conditions generated during this study are available at: https://doi.org/10.6084/m9.figshare.29294225.v1.

The complete training data and initial conditions can be requested from the corresponding author.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6sc00969g.

Acknowledgements

This work was supported by the National Science Foundation Center under NSF-CHE-2144556. All authors appreciate the assistance from the Northeastern Research Computing Team and the computing resources provided by the Massachusetts Life Science Center grant (G00006360).

References

1. S. Cho, J. Usuba, S. Chakraborty, X. Li and G. G. D. Han, Solid-state photon energy storage via reversible [2+2] cycloaddition of donor-acceptor styrylpyrylium system, Chem, 2023, 9(11), 3159–3171,  DOI:10.1016/j.chempr.2023.06.007.
2. A. A. Poloukhtine, N. E. Mbua, M. A. Wolfert, G.-J. Boons and V. V. Popik, Selective Labeling of Living Cells by a Photo-Triggered Click Reaction, J. Am. Chem. Soc., 2009, 131(43), 15769–15776,  DOI:10.1021/ja9054096.
3. G. Li, Q.-S. Li, Y. Xie, R. B. King and H. F. Schaefer, (Acetylene)dicobalt Carbonyl Derivatives: Decarbonylation of the H 2 C 2 Co 2 (CO) 6 Tetrahedrane, Organometallics, 2009, 28(12), 3390–3394,  DOI:10.1021/om801211g.
4. M. H. Chang and D. A. Dougherty, Photochemistry of 2,3-diazabicyclo[2.1.1]hex-2-ene. .beta. Carbon-carbon cleavage to a stereorandom triplet biradical, J. Am. Chem. Soc., 1982, 104(8), 2333–2334,  DOI:10.1021/ja00372a050.
5. S. Zhu, L. Jia, Q. Cheng, Q. Sun, X. Chen, H. Yu, Y. Han and H. Hou, Visible-Light-Induced Stereoselective Radical trans-Iodoalkylation of Terminal Alkyne with Iodoform, Org. Lett., 2024, 26(39), 8400–8404,  DOI:10.1021/acs.orglett.4c03241.
6. A. Noble, S. J. McCarver and D. W. MacMillan, Merging photoredox and nickel catalysis: decarboxylative cross-coupling of carboxylic acids with vinyl halides, J. Am. Chem. Soc., 2015, 137(2), 624–627,  DOI:10.1021/ja511913h.
7. S. D. Andrews and A. C. Day, The photolysis of 4-alkylidene-Δ1-pyrazolines: a route to trimethylenemethyls, Chem. Commun., 1966,(19), 667–669,  10.1039/C19660000667.
8. S. Kar, H. Sanderson, K. Roy, E. Benfenati and J. Leszczynski, Green Chemistry in the Synthesis of Pharmaceuticals, Chem. Rev., 2022, 122(3), 3637–3710,  DOI:10.1021/acs.chemrev.1c00631.
9. T. Long, S. Pan, S. Zhu and L. Chu, Catalyst-Free Intermolecular Sulfonyl/Fluoromethyl Heteroarylation of Vinyl Ethers via Visible-Light-Induced Charge Transfer, Chemistry, 2022, 28(15), e202104080,  DOI:10.1002/chem.202104080.
10. G. Dorman and G. D. Prestwich, Benzophenone Photophores in Biochemistry, Biochemistry, 1994, 33(19), 5661–5673,  DOI:10.1021/bi00185a001.
11. E. Smith and I. Collins, Photoaffinity labeling in target- and binding-site identification, Future Med. Chem., 2015, 7(2), 159–183,  DOI:10.4155/fmc.14.152.
12. A. B. Kumar, J. D. Tipton and R. Manetsch, 3-Trifluoromethyl-3-aryldiazirine photolabels with enhanced ambient light stability, Chem. Commun., 2016, 52(13), 2729–2732,  10.1039/c5cc09518b.
13. R. O. Al-Kaysi, F. Tong, M. Al-Haidar, L. Zhu and C. J. Bardeen, Highly branched photomechanical crystals, Chem. Commun., 2017, 53(17), 2622–2625,  10.1039/c6cc08999b.
14. H. Koshima, N. Ojima and H. Uchimoto, Mechanical Motion of Azobenzene Crystals upon Photoirradiation, J. Am. Chem. Soc., 2009, 131(20), 6890–6891,  DOI:10.1021/ja8098596.
15. A. Lennartson, A. Roffey and K. Moth-Poulsen, Designing photoswitches for molecular solar thermal energy storage, Tetrahedron Lett., 2015, 56(12), 1457–1465,  DOI:10.1016/j.tetlet.2015.01.187.
16. A. A. Beharry, L. Wong, V. Tropepe and G. A. Woolley, Fluorescence imaging of azobenzene photoswitching in vivo, Angew Chem. Int. Ed. Engl., 2011, 50(6), 1325–1327,  DOI:10.1002/anie.201006506.
17. Y. Yoneda, T. Konishi, K. Suga, S. Saito and H. Kuramochi, Excited-State Aromatization Drives Nonequilibrium Planarization Dynamics, J. Am. Chem. Soc., 2025, 147(14), 12051–12060,  DOI:10.1021/jacs.4c18623.
18. H. Lee, J. Tessarolo, D. Langbehn, A. Baksi, R. Herges and G. H. Clever, Light-Powered Dissipative Assembly of Diazocine Coordination Cages, J. Am. Chem. Soc., 2022, 144(7), 3099–3105,  DOI:10.1021/jacs.1c12011.
19. D. J. Shields, D. P. Karothu, K. Sambath, R. Ranaweera, S. Schramm, A. Duncan, B. Duncan, J. A. Krause, A. D. Gudmundsdottir and P. Naumov, Cracking under Internal Pressure: Photodynamic Behavior of Vinyl Azide Crystals through N(2) Release, J. Am. Chem. Soc., 2020, 142(43), 18565–18575,  DOI:10.1021/jacs.0c07830.
20. N. Yamamoto, M. Olivucci, P. Celani, F. Bernardi and M. A. Robb, An MC-SCF/MP2 Study of the Photochemistry of 2,3-Diazabicyclo[2.2.1]hept-2-ene: Production and Fate of Diazenyl and Hydrazonyl Biradicals, J. Am. Chem. Soc., 1998, 120(10), 2391–2407,  DOI:10.1021/ja971733v.
21. Y. Huang, X. Zheng, J. Wu, Y. Gao, Q. Ling and Z. Lin, Photoinduced pi-Bond breakage causing dynamic closing-opening shell transition of Z-type Diphenylmaleonitriles molecules, Nat. Commun., 2024, 15(1), 6514,  DOI:10.1038/s41467-024-50943-4.
22. S. Kobatake and D. Kitagawa, Photoresponsive Molecular Crystals for Light-Driven Photoactuators, Photosynergetic Responses in Molecules and Molecular Aggregates, H. Miyasaka, K. Matsuda, J. Abe and T. Kawai, Springer Singapore, 2020, pp 427–447.
23. E. Skolia and C. G. Kokotos, Photochemical [2 + 2] Cycloaddition of Alkenes with Maleimides: Highlighting the Differences between N-Alkyl vs. N-Aryl Maleimides, ACS Org. Inorg. Au, 2023, 3(2), 96–103,  DOI:10.1021/acsorginorgau.2c00053.
24. J. Yan, T. Slanina, J. Bergman and H. Ottosson, Photochemistry Driven by Excited-State Aromaticity Gain or Antiaromaticity Relief, Chemistry, 2023, 29(19), e202203748,  DOI:10.1002/chem.202203748.
25. M. Ueda, K. Jorner, Y. M. Sung, T. Mori, Q. Xiao, D. Kim, H. Ottosson, T. Aida and Y. Itoh, Energetics of Baird aromaticity supported by inversion of photoexcited chiral [4n]annulene derivatives, Nat. Commun., 2017, 8(1), 346,  DOI:10.1038/s41467-017-00382-1.
26. J. Drazenovic, C. J. Laconsay, N. Doslic, I. C. W. J and N. Basaric, Excited-state antiaromaticity relief drives facile photoprotonation of carbons in aminobiphenyls, Chem. Sci., 2024, 15(14), 5225–5237,  10.1039/d4sc00642a.
27. X. Li, S. Cho, J. Wan and G. G. D. Han, Photoswitches and photochemical reactions for optically controlled phase transition and energy storage, Chem, 2023, 9(9), 2378–2389,  DOI:10.1016/j.chempr.2023.05.029.
28. M. Manso, A. U. Petersen, Z. Wang, P. Erhart, M. B. Nielsen and K. Moth-Poulsen, Molecular solar thermal energy storage in photoswitch oligomers increases energy densities and storage times, Nat. Commun., 2018, 9(1), 1945,  DOI:10.1038/s41467-018-04230-8.
29. M. López-Cano, M. Scortichini, D. K. Tosh, V. Salmaso, T. Ko, G. Salort, I. Filgaira, C. Soler, D. Trauner and J. Hernando, et al., Photoswitchable Diazocine Derivative for Adenosine A3 Receptor Activation in Psoriasis, J. Am. Chem. Soc., 2025, 147(1), 874–879,  DOI:10.1021/jacs.4c13558.
30. B. Y. Wang, Y. C. Lin, Y. T. Lai, J. Y. Ou, W. W. Chang and C. C. Chu, Targeted photoresponsive carbazole-coumarin and drug conjugates for efficient combination therapy in leukemia cancer cells, Bioorg. Chem., 2020, 100, 103904,  DOI:10.1016/j.bioorg.2020.103904.
31. P. A. Shaw, M. Klausen and M. Bradley, A dual action coumarin-camptothecin polymer for light responsive drug release and photodynamic therapy, Polym. Chem., 2024, 15(2), 54–58,  10.1039/d3py01137b.
32. M. Rosenberg, C. Dahlstrand, K. Kilså and H. Ottosson, Excited State Aromaticity and Antiaromaticity: Opportunities for Photophysical and Photochemical Rationalizations, Chem. Rev., 2014, 114(10), 5379–5425,  DOI:10.1021/cr300471v.
33. G. Merino, M. Sola, I. Fernandez, C. Foroutan-Nejad, P. Lazzeretti, G. Frenking, H. L. Anderson, D. Sundholm, F. P. Cossio and M. A. Petrukhina, et al., Aromaticity: Quo Vadis, Chem. Sci., 2023, 14(21), 5569–5576,  10.1039/d2sc04998h.
34. H. Ottosson, A focus on aromaticity: fuzzier than ever before?, Chem. Sci., 2023, 14(21), 5542–5544,  10.1039/d3sc90075d.
35. P. v. R. Schleyer, Introduction: Aromaticity, Chem. Rev., 2001, 101(5), 1115–1118,  DOI:10.1021/cr0103221.
36. M. Faraday, On new compounds of carbon and hydrogen, and on certain other products obtained during the decomposition of oil by heat, Phil. Trans. Roy. Soc. Lond., 1825(115), 440–466,  DOI:10.1098/rstl.1825.0022.
37. E. Hückel, Quantentheoretische Beiträge zum Benzolproblem, Z. Phys., 1931, 70, 204–286,  DOI:10.1007/BF01339530.
38. Aromatic Character, Chem. Eng. News, 1965, 43(26), 90–100,  DOI:10.1021/cen-v043n026.p090.
39. M. J. S. Dewar, A molecular orbital theory of organic chemistry—VIII: romaticity and electrocyclic reactions, Tetrahedron, 1966, 22, 75–92,  DOI:10.1016/S0040-4020(01)82171-2.
40. N. C. Baird, Quantum organic photochemistry. II. Resonance and aromaticity in the lowest 3.pi.pi.* state of cyclic hydrocarbons, J. Am. Chem. Soc., 1972, 94(14), 4941–4948,  DOI:10.1021/ja00769a025.
41. P. B. Karadakov, Ground- and Excited-State Aromaticity and Antiaromaticity in Benzene and Cyclobutadiene, J. Phys. Chem., 2008, 112(31), 7303–7309,  DOI:10.1021/jp8037335.
42. P. B. Karadakov, Aromaticity and Antiaromaticity in the Low-Lying Electronic States of Cyclooctatetraene, J. Phys. Chem., 2008, 112(49), 12707–12713,  DOI:10.1021/jp8067365.
43. J. Oh, Y. M. Sung, H. Mori, S. Park, K. Jorner, H. Ottosson, M. Lim, A. Osuka and D. Kim, Unraveling Excited-Singlet-State Aromaticity via Vibrational Analysis, Chem, 2017, 3(5), 870–880,  DOI:10.1016/j.chempr.2017.09.005.
44. Y. M. Sung, J. Oh, K. Naoda, T. Lee, W. Kim, M. Lim, A. Osuka and D. Kim, A Description of Vibrational Modes in Hexaphyrins: Understanding the Aromaticity Reversal in the Lowest Triplet State, Angew Chem. Int. Ed. Engl., 2016, 55(39), 11930–11934,  DOI:10.1002/anie.201603631.
45. L. J. Karas, C. H. Wu and J. I. Wu, Barrier-Lowering Effects of Baird Antiaromaticity in Photoinduced Proton-Coupled Electron Transfer (PCET) Reactions, J. Am. Chem. Soc., 2021, 143(43), 17970–17974,  DOI:10.1021/jacs.1c09324.
46. J. Kim, J. Oh, S. Park, J. L. Zafra, J. R. DeFrancisco, D. Casanova, M. Lim, J. D. Tovar, J. Casado and D. Kim, Two-electron transfer stabilized by excited-state aromatization, Nat. Commun., 2019, 10(1), 4983,  DOI:10.1038/s41467-019-12986-w.
47. B. J. Lampkin, Y. H. Nguyen, P. B. Karadakov and B. VanVeller, Demonstration of Baird's rule complementarity in the singlet state with implications for excited-state intramolecular proton transfer, Phys. Chem. Chem. Phys., 2019, 21(22), 11608–11614,  10.1039/c9cp02050k.
48. Z. Zhao, X. Zheng, L. Du, Y. Xiong, W. He, X. Gao, C. Li, Y. Liu, B. Xu and J. Zhang, et al., Non-aromatic annulene-based aggregation-induced emission system via aromaticity reversal process, Nat. Commun., 2019, 10(1), 2952,  DOI:10.1038/s41467-019-10818-5.
49. S. Saito, S. Nobusue, E. Tsuzaka, C. Yuan, C. Mori, M. Hara, T. Seki, C. Camacho, S. Irle and S. Yamaguchi, Light-melt adhesive based on dynamic carbon frameworks in a columnar liquid-crystal phase, Nat. Commun., 2016, 7, 12094,  DOI:10.1038/ncomms12094.
50. P. P. Kalapos, P. J. Mayer, T. Gazdag, A. Demeter, B. Oruganti, B. Durbeej and G. London, Photoswitching of Local (Anti)Aromaticity in Biphenylene-Based Diarylethene Molecular Switches, J. Org. Chem., 2022, 87(15), 9532–9542,  DOI:10.1021/acs.joc.2c00504.
51. S. Paul, H. Kitakado, K. Suga, R. Kotani, N. Dey, R. Venkatramani, E. Matito, S. Saito and J. Dasgupta, Triplet conformation in chromophore-fused cyclooctatetraene dyes, J. Mater. Chem. C, 2023, 11(36), 12243–12253,  10.1039/d3tc02151c.
52. A. K. Pati, O. El Bakouri, S. Jockusch, Z. Zhou, R. B. Altman, G. A. Fitzgerald, W. B. Asher, D. S. Terry, A. Borgia and M. D. Holsey, et al., Tuning the Baird aromatic triplet-state energy of cyclooctatetraene to maximize the self-healing mechanism in organic fluorophores, Proc. Natl. Acad. Sci. U. S. A., 2020, 117(39), 24305–24315,  DOI:10.1073/pnas.2006517117.
53. T. Slanina, R. Ayub, J. Toldo, J. Sundell, W. Rabten, M. Nicaso, I. Alabugin, I. Fdez Galvan, A. K. Gupta and R. Lindh, et al., Impact of Excited-State Antiaromaticity Relief in a Fundamental Benzene Photoreaction Leading to Substituted Bicyclo[3.1.0]hexenes, J. Am. Chem. Soc., 2020, 142(25), 10942–10954,  DOI:10.1021/jacs.9b13769.
54. R. K. Mohamed, S. Mondal, K. Jorner, T. F. Delgado, V. V. Lobodin, H. Ottosson and I. V. Alabugin, The Missing C1-C5 Cycloaromatization Reaction: Triplet State Antiaromaticity Relief and Self-Terminating Photorelease of Formaldehyde for Synthesis of Fulvenes from Enynes, J. Am. Chem. Soc., 2015, 137(49), 15441–15450,  DOI:10.1021/jacs.5b07448.
55. R. Papadakis and H. Ottosson, The excited state antiaromatic benzene ring: a molecular Mr Hyde?, Chem. Soc. Rev., 2015, 44(18), 6472–6493,  10.1039/c5cs00057b.
56. P. Wan and E. Krogh, Evidence for the generation of aromatic cationic systems in the excited state. Photochemical solvolysis of fluoren-9-ol, J. Chem. Soc. Chem. Commun., 1985,(17), 1207–1208,  10.1039/C39850001207 10.1039/C39850001207.
57. R. Kotani, H. Sotome, H. Okajima, S. Yokoyama, Y. Nakaike, A. Kashiwagi, C. Mori, Y. Nakada, S. Yamaguchi and A. Osuka, et al., Flapping viscosity probe that shows polarity-independent ratiometric fluorescence, J. Mater. Chem. C, 2017, 5(21), 5248–5256,  10.1039/c7tc01533j.
58. T. Yamakado, S. Takahashi, K. Watanabe, Y. Matsumoto, A. Osuka and S. Saito, Conformational Planarization versus Singlet Fission: Distinct Excited-State Dynamics of Cyclooctatetraene-Fused Acene Dimers, Angew Chem. Int. Ed. Engl., 2018, 57(19), 5438–5443,  DOI:10.1002/anie.201802185.
59. D. Xing, F. Glocklhofer and F. Plasser, Proton transfer induced excited-state aromaticity gain for chromophores with maximal Stokes shifts, Chem. Sci., 2024, 15(43), 17918–17926,  10.1039/d4sc04692g.
60. D. Shukla and P. Wan, Evidence for a planar cyclically conjugated 8.pi. system in the excited state: large Stokes shift observed for dibenz[b,f]oxepin fluorescence, J. Am. Chem. Soc., 1993, 115(7), 2990–2991,  DOI:10.1021/ja00060a063.
61. F.-G. Klärner, About the Antiaromaticity of Planar Cyclooctatetraene, Angew. Chem., Int. Ed., 2001, 40(21), 3977–3981,  DOI:10.1002/1521-3773(20011105)40.
62. T. Nishinaga, T. Ohmae and M. Iyoda, Recent Studies on the Aromaticity and Antiaromaticity of Planar Cyclooctatetraene, Symmetry, 2010, 2(1), 76–97,  DOI:10.3390/sym2010076.
63. J. I. Wu, I. Fernandez, Y. Mo and P. Schleyer, Why Cyclooctatetraene Is Highly Stabilized: The Importance of “Two-Way” (Double) Hyperconjugation, J. Chem. Theory Comput., 2012, 8(4), 1280–1287,  DOI:10.1021/ct3000553.
64. A. Schild and B. Paulus, Multireference calculations for ring inversion and double bond shifting in cyclooctatetraene, J. Comput. Chem., 2013, 34(16), 1393–1397,  DOI:10.1002/jcc.23273.
65. C. Tokizaki, T. Yoshida and T. Takayanagi, Quantum transition state dynamics of the cyclooctatetraene unimolecular reaction on ab initio potential energy surfaces, Chem. Phys., 2016, 469–470, 97–104,  DOI:10.1016/j.chemphys.2016.02.005.
66. V. Gogonea, P. v. R. Schleyer and P. R. Schreiner, Consequences of Triplet Aromaticity in 4nπ-Electron Annulenes: Calculation of Magnetic Shieldings for Open-Shell Species, Angew. Chem., 1998, 37, 1945–1948.
67. F. Feixas, J. Vandenbussche, P. Bultinck, E. Matito and M. Sola, Electron delocalization and aromaticity in low-lying excited states of archetypal organic compounds, Phys. Chem. Chem. Phys., 2011, 13(46), 20690–20703,  10.1039/c1cp22239b.
68. H. Song, Y. Nam, D. Keefer, M. Garavelli, S. Mukamel and S. Tretiak, Nonadiabatic Molecular Dynamics Study of the Relaxation Pathways of Photoexcited Cyclooctatetraene, J. Phys. Chem. Lett., 2021, 12(24), 5716–5722,  DOI:10.1021/acs.jpclett.1c01397.
69. P. B. Karadakov and N. Preston, Aromaticity reversals and their effect on bonding in the low-lying electronic states of cyclooctatetraene, Phys. Chem. Chem. Phys., 2021, 23(43), 24750–24756,  10.1039/d1cp04394c.
70. M. Garavelli, F. Bernardi, A. Cembran, O. Castaño, L. M. Frutos, M. Merchán and M. Olivucci, Cyclooctatetraene Computational Photo- and Thermal Chemistry: A Reactivity Model for Conjugated Hydrocarbons, J. Am. Chem. Soc., 2002, 124(46), 13770–13789,  DOI:10.1021/ja020741v.
71. A. A. Kroeger and A. Karton, Graphene-induced planarization of cyclooctatetraene derivatives, J. Comput. Chem., 2022, 43(2), 96–105,  DOI:10.1002/jcc.26774.
72. J. Li, P. Reiser, B. R. Boswell, A. Eberhard, N. Z. Burns, P. Friederich and S. A. Lopez, Automatic discovery of photoisomerization mechanisms with nanosecond machine learning photodynamics simulations, Chem. Sci., 2021, 12(14), 5302–5314,  10.1039/d0sc05610c.
73. J. Li and S. A. Lopez, Machine learning accelerated photodynamics simulations, Chem. Phys. Rev., 2023, 4(3), 031309,  DOI:10.1063/5.0159247.
74. J. Li, R. Stein, D. M. Adrion and S. A. Lopez, Machine-Learning Photodynamics Simulations Uncover the Role of Substituent Effects on the Photochemical Formation of Cubanes, J. Am. Chem. Soc., 2021, 143(48), 20166–20175,  DOI:10.1021/jacs.1c07725.
75. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, et al., Gaussian 16 Rev. B.01, 2016.
76. T. Yanai, D. P. Tew and N. C. Handy, A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP), Chem. Phys. Lett., 2004, 393(1), 51–57,  DOI:10.1016/j.cplett.2004.06.011.
77. J.-D. Chai and M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections, Phys. Chem. Chem. Phys., 2008, 10(44), 6615–6620,  10.1039/B810189B 10.1039/B810189B.
78. T. H. Dunning Jr, Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen, J. Chem. Phys., 1989, 90(2), 1007–1023,  DOI:10.1063/1.456153.
79. R. A. Kendall, T. H. Dunning Jr. and R. J. Harrison, Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions, J. Chem. Phys., 1992, 96(9), 6796–6806,  DOI:10.1063/1.462569.
80. A. L. Thompson and T. J. Martinez, Time-resolved photoelectron spectroscopy from first principles: excited state dynamics of benzene, Faraday Discuss., 2011, 150, 293–311,  10.1039/c1fd00003a.
81. J. Li and S. A. Lopez, Multiconfigurational Calculations and Nonadiabatic Molecular Dynamics Explain Tricyclooctadiene Photochemical Chemoselectivity, J. Phys. Chem. A, 2020, 124(38), 7623–7632,  DOI:10.1021/acs.jpca.0c05280.
82. J. Finley, P.-Å. Malmqvist, B. O. Roos and L. Serrano-Andrés, The multi-state CASPT2 method, Chem. Phys. Lett., 1998, 288(2), 299–306,  DOI:10.1016/S0009-2614(98)00252-8.
83. M. Boggio-Pasqua, D. Jacquemin and P. F. Loos, Benchmarking CASPT3 vertical excitation energies, J. Chem. Phys., 2022, 157(1), 014103,  DOI:10.1063/5.0095887.
84. R. Sarkar, P. F. Loos, M. Boggio-Pasqua and D. Jacquemin, Assessing the Performances of CASPT2 and NEVPT2 for Vertical Excitation Energies, J. Chem. Theory Comput., 2022, 18(4), 2418–2436,  DOI:10.1021/acs.jctc.1c01197.
85. M. R. Silva-Junior, M. Schreiber, S. P. Sauer and W. Thiel, Benchmarks of electronically excited states: basis set effects on CASPT2 results, J. Chem. Phys., 2010, 133(17), 174318,  DOI:10.1063/1.3499598.
86. M. P. Freitas, Regio- and stereochemical behavior of cyclooctane and cyclodecane derivatives containing CF2 moieties, J. Fluorine Chem., 2025, 282, 110397,  DOI:10.1016/j.jfluchem.2025.110397.
87. J. C. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys., 1990, 93(2), 1061–1071,  DOI:10.1063/1.459170.
88. X. Zhao, I. C. D. Merritt, R. Lei, Y. Shu, D. Jacquemin, L. Zhang, X. Xu, M. Vacher and D. G. Truhlar, Nonadiabatic Coupling in Trajectory Surface Hopping: Accurate Time Derivative Couplings by the Curvature-Driven Approximation, J. Chem. Theory Comput., 2023, 19(19), 6577–6588,  DOI:10.1021/acs.jctc.3c00813.
89. Y. Shu, L. Zhang, X. Chen, S. Sun, Y. Huang and D. G. Truhlar, Nonadiabatic Dynamics Algorithms with Only Potential Energies and Gradients: Curvature-Driven Coherent Switching with Decay of Mixing and Curvature-Driven Trajectory Surface Hopping, J. Chem. Theory Comput., 2022, 18(3), 1320–1328,  DOI:10.1021/acs.jctc.1c01080.
90. T. Aoki, M. Ueda, T. Aida and Y. Itoh, Supramolecular Polymerization of a Photo-Fluttering Chiral Monomer: A Temporarily Suspendable Chain Growth by Light, J. Am. Chem. Soc., 2022, 144(16), 7080–7084,  DOI:10.1021/jacs.2c02176.
91. T. Takeda, M. Ozawa and T. Akutagawa, Jumping Crystal of a Hydrogen-Bonded Organic Framework Induced by the Collective Molecular Motion of a Twisted pi System, Angew Chem. Int. Ed. Engl., 2019, 58(30), 10345–10352,  DOI:10.1002/anie.201905075.
92. I. Fdez Galvan, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov and R. K. Carlson, et al., OpenMolcas: From Source Code to Insight, J. Chem. Theory Comput., 2019, 15(11), 5925–5964,  DOI:10.1021/acs.jctc.9b00532.
93. M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. T. Devin, Large-scale machine learning on heterogeneous distributed systems, arXiv, 2016, preprint, arXiv:1603.04467,  DOI:10.48550/arXiv.1603.04467.

---