Machine learning photodynamics reveal intersystem-crossing-driven ladderdiene ring opening

Abstract

Photochemical ring-opening reactions have become an essential tool for chemical syntheses under mild conditions with high atom economy. We propose a near-visible light-induced electrocyclic ring-opening reaction to afford cyclooctatetraene (COT) using carbonyl-functionalized tricyclooctadiene (1) based on our machine learning (ML) accelerated photodynamics simulations. Our CAM-B3LYP/cc-pVDZ calculations show that carbonyl group reduce the S₁-energy of 1 to 3.65 eV (340 nm) from 6.25 eV, approaching the visible light range. The multiconfigurational CASSCF(12,11)/ANO-RCC-VDZP calculations show small S₁ and T₁ energy gaps near an S₁-minimum region. Our ML-photodynamics simulations with 1000 FSSH trajectories revealed a stepwise ring-opening mechanism of 1 from the S₁, dominated by relatively fast S₁/T₁ intersystem crossings in 20 ps. The trajectories predict that the quantum yield to carbonyl-functionalized COT is 89%, suggesting the light-induced ring-opening reaction of 1 is highly efficient. This work demonstrates a predictive ML-photodynamics application for photochemical reaction design.

---

Introduction

Cyclooctatetraene (COT) is an essential molecular framework with an eight-membered ring with four π_CC-bonds. The unique chemical structure of COT can access versatile chemical products via ring-closing and cycloaddition reactions, which have been widely used in natural product synthesis,^1,2 drug design,³ and organic semiconductors.⁴COT features a low-lying triplet state (1.03 eV),⁵ notable triplet absorption (350 nm),⁶ and a short triplet lifetime (100 μs).⁶ Moreover, the triplet COT exhibits relatively low reactivity towards molecular oxygen. These features make COT an excellent triplet-state quencher for developing stable organic fluorophores⁷ and laser devices.⁸

The first synthesis of COT was reported in 1911 by Willstätter using pseudopelletierine (Fig. 1a).⁹ Later, many efforts have been made to improve the syntheses of COT and its derivatives with ethylene,¹⁰ butadiene,^11,12 and cycloocatdiene¹³ (Fig. 1b). However, most strategies require transition metal catalysts (e.g., Ni,^14–16 Cu,^11,12 and Pb¹⁷) and extensive thermal energy and reaction time. Thus, a cost-efficient metal-free synthesis of COT remains a challenging task. Photochemical reactions are gaining increasing importance in organic synthesis because of the mild reaction conditions,¹⁸ high atomic economy,¹⁹ and access to exotic molecular structures.²⁰ According to the literature, the photoexcitation of tricyclooctadiene (TOD) favors the electrocyclic ring-opening reactions toward COT.²¹ Our previous photodynamics studies on TOD showed the ring-opening of TOD occurs in the sub-picosecond timescale^22,23 (Fig. 1c). Our simulations reproduced the trend of experimentally observed cubane yields depending on the substituent effects and explained the reaction pathway toward COT formation. Nevertheless, the lowest light absorption (S₁) is far beyond the visible range (198 nm) due to the transitions to the ππ* states,²² which limit its applicability to produce COT.

Fig. 1 (a) Synthesis of COTvia elimination reactions of pseudopelletierine by Willstätter. (b) Synthesis of COT with transition metal-catalyzed thermal cycloaddition of alkene, dienes, and oxidation of cyclodiene. (c) Previously predicted photochemical 4π-disrotatory electrocyclic ring-opening of TOD toward COT. (d) Proposed photochemical ring-opening of the carbonyl-functionalized TOD (1).

Extending the π-conjugation of the vinyl groups with carbonyl groups to make α,β-unsaturated carbonyl compound can substantially reduce the energy S₁ state while maintaining the photochemical reactivity.^24,25 A recent computational and experimental study showed self-sensitized photochemical dimerization of isophorone under visible-light irradiation.²⁴ These findings prompted us to envision a theoretical system featuring a carbonyl-functionalized TOD (Fig. 1d), 1 where visible light could promote the electrocyclic ring-opening reaction to generate a COT derivative (4). We are unaware of any experimental or computational studies investigating the photochemical reaction of 1. It may result from time-resolved experiments that currently cannot resolve the structural information and inform the mechanism. Carbonyl compound are known to induce intersystem crossing during the excited-state dynamics,^26,27 which can extend the excited-state lifetime from tens to hundreds of picoseconds. The requisite multiconfigurational quantum mechanical calculations are computationally infeasible for the 10–100 s of picoseconds, even for a relatively small molecule like 1.

In this work, we combine the quantum mechanical calculations and machine learning (ML)-accelerated photodynamics approach to study the photochemical reactions of 1. We first computed the vertical excitation energy of 1 with the time-dependent density functional theory (TDDFT) and complete active space self-consistent field (CASSCF) methods. We then train neural networks (NN) based on the ground- and excited-state CASSCF energies and gradients to simulate the excited-state dynamics of 1 starting from the S₁-Franck Condon (FC) regions. The discussion section analyzes the trajectories of possible excited-state mechanistic pathways of 1. This manuscript aims to demonstrate a promising photochemical electrocyclic ring-opening reaction of 1 for the cost-efficient and metal-free synthesis of COT under mild conditions.

Results and discussion

Vertical excitation

The previous study reported the lowest absorption of the unsubstituted tricyclooctadiene (TOD) at 6.25 eV (198 nm), requiring ultraviolet light for photoexcitation.²² Towards a more chemoselective and sustainable methodology towards COT, we introduce a carbonyl group to modulate the photophysical properties of TOD. Here, we used the aldehyde group as smallest carbonyl functional group to reduce the computational cost. Table 1 summarizes the vertical excitation energies and oscillator strength of 1 computed by various methods.

Table 1 The vertical excitation energies oscillator strength of 1 using different computational methods^a

Methods	Energy (eV)						Oscillator strength
	S1	S2	S3	T1	T2	T3	S1	S2	S3
a The dominant electronic configurations are labelled in the parenthesis.
CAM-B3LYP/cc-pVDZ	3.65 (nπ*)	5.14 (ππ*)	5.58 (σπ*)	2.87 (nπ*)	3.03 (ππ*)	3.74 (ππ*)	0.0006	0.0344	0.1171
CAM-B3LYP/aug-cc-pVDZ	3.68 (nπ*)	4.99 (ππ*)	5.40 (σπ*)	2.87 (nπ*)	3.09 (ππ*)	3.73 (ππ*)	0.0009	0.0384	0.1313
CAM-B3LYP/cc-pVTZ	3.70 (nπ*)	5.05 (ππ*)	5.47 (σπ*)	2.87 (nπ*)	3.10 (ππ*)	3.73 (ππ*)	0.0008	0.0355	0.1252
CASSCF(12,11)/ANO-RCC-VDZP	3.60 (nπ*)	6.92 (nπ* + ππ*)	7.45 (ππ*)	3.40 (nπ*)	3.55 (ππ*)	4.26 (ππ*)	0.00002	0.0020	0.0394
CASSCF(12,11)/ANO-RCC-VTZP	3.54 (nπ*)	6.86 (nπ* + ππ*)	7.40 (ππ*)	3.36 (nπ*)	3.52 (ππ*)	4.25 (ππ*)	0.00004	0.0017	0.0309
CASPT2 (12,11)/ANO-RCC-VDZP	3.78 (nπ*)	7.03 (ππ*^2)	7.48 (ππ*^2)	3.47 (nπ*)	3.74 (ππ*)	4.45 (ππ*)	0.00002	0.0021	0.0396
XMS-CASPT2 (12,11)/ANO-RCC-VDZP	4.01 (nπ*)	7.34 (ππ*^2)	7.61 (ππ*^2)	3.64 (nπ*)	3.92 (ππ*)	4.54 (ππ*)	0.00003	0.0065	0.0360

---

The CAM-B3LYP/cc-pVDZ calculations show the S₁, S₂, and S₃ energies of 1 are 3.65, 5.14 and 5.58 eV, respectively. These results indicate the carbonyl group substantially lowers the S₁ energy of 1 (3.65–3.70 eV) compared to TOD (6.25 eV).²² The S₁ of 1 is a nπ* state governed by the electronic transition from the oxygen lone pair of the carbonyl group to the lowest unoccupied π*-orbital of 1 (Fig. S1†). The computed oscillator strengths (0.0006–0.0008) imply a low-intensity absorption to the S₁ of 1 upon photoexcitation near the visible light (340 nm). The S₂ and S₃ are ππ* and σπ* states with larger absorption intensities and oscillator strengths than S₁, but photoexcitation to S₂ or S₃ requires UV light. Including diffuse functions in the double-zeta basis set (CAM-B3LYP/aug-cc-pVDZ) shows similar excitation energies, 3.68, 4.99, and 5.40 eV, respectively. Triple-zeta basis sets (i.e., CAM-B3LYP/cc-pVTZ) do not significantly affect the vertical excitation energies and characters of S₁, S₂, and S₃ (3.70, 5.05, and 5.47 eV, respectively).

The triplet calculations at the CAM-B3LYP/cc-pVDZ level show nearly degenerate T_1–3 states, where the T₁, T_2, and T₃ states are 2.87, 3.03, and 3.74 eV above S₀, respectively. The T₁ is a nπ* state, T2 and T3 are ππ*. The diffuse functions led to no significant effects on the triplet state energies and characters. The calculations with the cc-pVTZ basis set reproduced the similar results as cc-pVDZ. The small S₁–T₁ and S₁–T₂ gaps suggest that intersystem crossings are possible after direct photoexcitation to the FC-region.

We compared the CAM-B3LYP and ωB97X-D3 results to determine the influence of the chosen functionals on excited-state electronic structure calculations. Table S1† shows that ωB97X-D3 predicted systematically higher excited-state energies than CAM-B3LYP, which are 0.05–0.12 eV for S₁, T_1, and T₂ and 0.22–0.36 eV for S₂ and S₃. Despite the disagreement in excited-state energies and oscillator strengths, CAM-B3LYP and ωB97X-D3 produced the same electronic transitions in all singlet and triplet excited states (Table S2†). Thus, the CAM-B3LYP/cc-pVDZ results continue to guide our assembly of an active space for subsequent CASSCF calculations.

The possible electrocyclic ring-opening reactions of 1 after the photoexcitation to S₁ could lead to bond-breaking, which requires multiconfigurational descriptions of the potential energy surface and possible state-crossing regions. Thus, we employed the CASSCF calculations to study the excited-state relaxation of 1 from the S₁-FC point toward 4, including the singlet internal conversion and the singlet–triplet intersystem crossing mechanistic pathways. The active space includes 12 electrons and 11 orbitals that describe all electrons involved in the electronic excitation according to the CAM-B3LYP/cc-pVDZ calculations and the essential electrons and orbitals participating in the photochemical reactions, as shown in Fig. 2.

Fig. 2 The (12,11) active space of 1 with the average occupation numbers, computed at SA5-CASSCF(12,11)/ANO-RCC-VDZP level of theory. Isosurface value = 0.04.

The S₁ energy (3.60 eV) and electronic configuration of the 5-state-averaged (SA5)-CASSCF(12,11)/ANO-RCC-VDZP calculations agree with CAM-B3LYP/cc-pVDZ reference calculations (3.65 eV). The S₂ and S₃ energies are higher than the CAM-B3LYP/cc-pVDZ energies (6.92 and 7.45 eV, respectively) because the S₂ changed to a mix of nπ* and ππ* state and the S₃ includes double excitations in the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. The oscillator strengths of S_1–3 at the CASSCF(12,11)/ANO-RCC-VDZP level show the same trend as the CAM-B3LYP/cc-pVDZ results. The T_1–3 energies from the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations are 3.40, 3.55, and 4.26 eV, respectively, which are systematically higher than the results of the CAM-B3LYP/cc-pVDZ calculations while the electronic configurations are consistent with the CAM-B3LYP/cc-pVDZ results (Table 1). We tested a larger basis set (ANO-RCC-VTZP); it does not change the results but substantially increases the computation cost by 6 times.

To verify the accuracy of the CASSCF results, we employed the state-specific and extended multistate (XMS)-CASPT2 (12,11)/ANO-RCC-VDZP calculations, which partially account for dynamical correlation. Both calculations show similar S₁ energies to the SA5-CASSCF(12,11)/ANO-RCC-VDZP results, confirming the nπ* character. However, the S₂ and S₃ states show (ππ*)² double excitations. Since TDDFT calculations cannot describe double excitations, it could result in discrepancies in the S₂ and S₃ energies compared to the CASPT2 and XMS-CASPT2 calculations. Nevertheless, the CASPT2(12,11)/ANO-RCC-VDZP and XMS-CASPT2(12,11)/ANO-RCC-VDZP calculations confirmed the relatively large S₁–S₂ and S₁–S₃ gaps (>3.0 eV) and the trends of their oscillator strengths observed in the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. These results are in line with the CAM-B3LYP/cc-pVDZ results (1.49 eV), which avoid S₁ state crossing with S₂ and S₃, especially in the FC region. The T_1–3 states computed with the CASPT2(12,11)/ANO-RCC-VDZP and XMS-CASPT2(12,11)/ANO-RCC-VDZP calculations show similar energies and consistent electronic configurations to the SA5-CASSCF(12,11)/ANO-RCC-VDZP and CAM-B3LYP/cc-pVDZ results (Table 1). Moreover, we benchmarked the spin–orbit couplings (SOC) between S₀, S₁, T_1, and T₂ (Table S4†). The SA5-CASSCF(12,11) calculations with the ANO-RCC-VDZP and ANO-RCC-VTZP basis sets show small differences of about 1 cm⁻¹ in the SOC norms. The SA5-CASSCF(12,11)/ANO-RCC-VDZP results are consistent with the CAM-B3LYP/ZORA-TZVP calculations, including the zero-order relativistic effect. Thus, the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations are suitable for studying the photochemical reactions from the low-lying S₁ state with the neighboring triplet states.

Ring-opening pathways

Our previous work showed the electrocyclic ring-opening reaction of TOD starts from a 4π-disrotatory ring-opening process of the cyclobutene ring.²² We anticipate two competing reaction pathways at the carbonyl-functionalized cyclobutene ring or the unsubstituted terminus of 1. Each ring-opening pathway goes through three mechanistic critical points from the S₁-FC point to an S₁/S₀ minimum energy conical intersection (MECI) towards the product. We computed the energy profiles of these critical structures in these two pathways using the SA5-CASSCF(12,11)/ANO-RCC-VDZP method to understand the ring-opening mechanism. The optimized structures and relative energies are collected in Fig. 3.

Fig. 3 (a) Potential energy profiles for the critical structures along with two electrocyclic ring-opening pathways at the SA5-CASSCF(12,11)/ANO-RCC-VDZP level. The red dotted lines show the breaking σ_CC bonds. (b) Optimized structures of MECIs. The red circles highlight the pyramidalized carbon, and the blue arches show the dihedral angles of over four connected carbon atoms. The lengths of breaking σ_CC bonds, substitution σ_CC bond, and π_CO bonds are labelled accordingly.

We define the reaction coordinate as the midpoint of the σ_CC bonds that undergo electrocyclic ring opening reactions. The electrocyclic ring-opening of 1 elongates along the reaction coordinate from 1.60 Å to 2.90 Å and 2.91 Å in two isomers (2 and 3) of the carbonyl-functionalized bicyclooctatriene (BOT) and to 3.02 Å in the product 4 (Fig. 3a). 3 is more stable than 2 because of its shorter σ_CC bond lengths to the carbonyl groups, which lowers the energy via stronger π-electron delocalization. The notable negative relative energies of 2, 3 and 4 indicate the ring-opening process of 1 is thermodynamically favored on the ground state, which provides the driving force for the ring-opening reactions.

We located two possible MECIs for each pathway. MECI-1 and MECI-2 show the breaking σ_CC bonds of 2.43 Å and 2.34 Å with dihedral angles of the π_CC bond with the carbonyl group, 48° and 39°, respectively; MECI-3 and MECI-4 show the breaking σ_CC bonds of 2.40 Å and 2.70 Å with dihedral angles of 47° and 61° on the other π_CC bond, respectively (Fig. 3b). The MECIs feature σ_CC bond lengths that lengthen with larger dihedral angles on the π_CC bonds. These dihedral angles resemble the pyramidalization of the carbon atoms in the S₁/S₀ conical intersection of the cis–trans isomerization of the ethylene.²⁸MECI-1 show a slightly longer σ_CC bond with the carbonyl group (1.48 Å) than MECI-2 (1.44 Å), while the MECI-3 and MECI-4 show similar bonds of 1.46 Å due to the similar cyclobutene moieties. The π_CO bond distance is nearly identical (1.20–1.21 Å) in all involved structures. The MECI energies show a trend of MECI-1 (2.89 eV) > MECI-3 (2.79 eV) > MECI-4 (2.70 eV) > MECI-2 (2.52 eV). We note that the energetic trends do not correlate to the geometrical changes of the π_CC bond but relate to the substitution σ_CC bond lengths. A possible explanation is that the shortened distance to the carbonyl groups could stabilize the diradicals via the extended delocalization, thus lowering the MECI energy. The energetic proximity of these MECIs suggests that multiple competing mechanistic pathways may exist.

The S₁-FC (3.60 eV) energy is nearly degenerate with T₂ (3.55 eV) and T₁ (3.40 eV) (Fig. 3a). These results prompted us to investigate S₁/T₁ and S₁/T₂ intersystem crossing mechanisms in 1. We interpolated the ring-opening pathways from the S₁-FC to the BOT intermediates to inform the energetic relationship between the singlet and triplet states. Fig. 4a illustrates a representative potential energy curve for the three lowest singlet and triplet states from 1 to 2viaMECI-2. The potential energy curves for other pathways show similar topology, as shown in Fig. S3.†

Fig. 4 (a) Interpolated reaction pathway for the first ring-opening step of 1, computed at SA5-CASSCF(12,11)/ANO-RCC-VDZP level. The blue arches show the dihedral angles of over four connected carbon atoms. (b) Optimized structures of the S₁/T₂ and T₂/T₁ minimum energy crossing points (MECP) with the relative energies computed at the SA5-CASSCF(12,11)/ANO-RCC-VDZP level.

The interpolated structures in Fig. 4a show continuous elongations along the reaction coordinate from 1 (1.60 Å) to 2 (2.90 Å). The maximum dihedral angle of the π_CC bond at the MECI-2 structure is 39°. The S₁ state shows substantially smaller gaps to T₁ and T₂ than S₂ and T₃ in the first 5 steps, suggesting the T₁ and T₂ states are more mechanistically relevant to the ring-opening reaction of 1 than the S₂ and T₃ states. An S₁ minimum was located near the S₁-FC point and is connected to MECI-2 through an excited-state transition structure. The increasing S₁ energy could slow the S₁ → S₀ transition viaMECI-2. On the other hand, the depth of the S₁-minimum combined with the small S₁/T₁ and S₁/T₂ energy gaps could promote the intersystem crossing. The nearly degenerate T₁ and T₂ states in the S₁ minimum region could also generate the T₂/T₁ internal conversion to enhance the T₁ population. The geometry optimizations for the four interpolated pathways led to the same S₁/T₂ and T₂/T₁ minimum energy crossing points (MECPs). The S₁/T₂ MECP structure shows increased π_CO bond length (1.26 Å) and reduced length of the substitution σ_CC bond (1.42 Å), while the T₂/T₁ MECI shows a nearly identical structure to the S₀ minimum geometry (Fig. 4b). The relative energies of the S₁/T₂ and T₂/T₁ crossings are 3.49 eV and 3.39 eV above the S₀ minimum and under the S₁-FC point. These results suggest the intersystem crossing the subsequent triplet internal conversions is energetically favored. The T₁ state shows a relatively flat potential energy near the S₁ barrier, which requires less energy than S₁ and T₂ to afford the ring opening of 1. The interpolated T₁ barriers are 0.48 eV, 0.42 eV, 0.95 eV, and 0.84 eV in pathways 1–4, respectively (Fig. S3†), which are confirmed by the energy profiles at the XMS-CASPT2(12,11)/ANO-RCC-VDZP level. These results indicate pathway 2 could be the most favored reaction channel. Despite the mechanistic hints provided by this static approach, only non-adiabatic molecular dynamics can provide answers to the hypotheses stated here.

Machine-learning photodynamic

To quantify the light-induced electrocyclic ring-opening mechanism for 1, we enumerated all possible reaction pathways from the S₁-FC regions using nonadiabatic molecular dynamics simulations based on the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. The simulated absorption spectrum shows notable bands in 300–400 nm because the vibronic coupling enhances the transition dipole moment of the nπ* state, suggesting the S₁ of 1 is accessible for photoexcitation (Fig. S4†). As such, our simulations started from the S₁-FC regions, where the initial conditions are generated by Wigner sampling at the zero-point energy level. We included S₀, S₁, T₁, and T₂ because they are involved in the ring-opening reaction of 1, according to the potential energy profile (Fig. 4). We employed the generalized fewest switches surface hopping (FSSH) method to evaluate the probability of the internal conversions and intersystem crossings.²⁹ The nonadiabatic couplings (NACs) are computed by the curvature-driven time derivative coupling (kTDC) approach.^30,31 The SOCs are computed with the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. Due to the high computational cost (e.g., 4 days for a single trajectory), we propagated 300 trajectories in 400 fs with a time step of 0.5 fs.

We used our open-source ML-photodynamics code Python Rapid Artificial Intelligence Ab Initio Molecular Dynamics (PyRAI²MD) to train neural networks and accelerate the CASSCF calculations of the energies, gradients, and SOCs for the considered states. The ML-photodynamic simulations allowed us to propagate about 1000 trajectories in 20 ps with a timestep of 0.5 fs at substantially reduced computational costs (e.g., 16 minutes for a single trajectory instead of 200 days with SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations). The NN training employed adaptive sampling to iteratively collect under-sampled data, where the unphysical structures with large energy outliers were discarded. 9% of the NN trajectories early stopped in the excited-state due to large NN prediction uncertainty, and thus are not included in subsequent discussions. The CASSCF- and NN-trajectories show similar energy conservation behavior. They maintained a small average energy drift of 0.08 and 0.01 eV, and gradually increased to 0.31 and 0.15 eV after multiple surface hoppings between discontinuous potential energy surfaces, which could be improved with reduced simulation time step.

The CASSCF state populations plot in Fig. 5a shows that 69% of trajectories remained in S₁, while the S₀ and T₁ populations increased to 12% and 16%, respectively. The large remaining S₁ population suggests the reaction pathways from the S₁-FC region are still largely unexplored. The NN-trajectories show comparable state populations to the SA5-CASSCF(12,11)/ANO-RCC-VDZP reference results (Fig. 5b). The S₁ population reduced to 79% while the T₁ population increased to 15%. It fits an S₁/T₁ intersystem crossing time constant of 1.7 ps, which is close to 2.0 ps obtained from the CASSCF reference trajectories. We note the NN-predicted trajectories of the S₁ → S₀ transitions are 2%, which is lower than the reference trajectories (12%).

Fig. 5 (a) The state population dynamics of 1 in 400 fs nonadiabatic molecular dynamics with the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. (b and c) The state population dynamics of 1 in the first 400 fs and 20 ps of the ML-photodynamics simulations.

We performed an error analysis to further determine the origin of the discrepancy between the reference and NN-predicted trajectories and possible influence on the simulations results. Fig. S5† collects the NN-predicted energies and gradients for all MECIs. The S₀ and S₁ of MECI-1, MECI-2, and MECI-3 show notable errors larger than 0.10 eV. As such, the enlarged out-of-sample prediction errors could lead to inaccurate population transfer from S₁ to S₀, resulting in an underestimated S₀ population in the dynamics. To determine the influence of underestimated S₁ → S₀ transitions on the reaction channels, we turned off surface hopping (i.e., infinitely slow intersystem crossing) and still observed the ring-opening reaction of 1 (Table S7†). This test suggests the ring-opening reaction of 1 is not sensitive to the underestimation of S₁ → S₀ transitions.

Moreover, we prepared an out-of-sample test set based on the S₁/T₁ surface hopping structures collected from the ML-photodynamics trajectories to inform the NN prediction accuracy for out-of-sample structures far from the equilibrium geometries. Table S8† shows that the test errors in energies, gradients, and SOC norms are 5, 2, and 6 times the validation errors. In contrast to the narrow error distributions of energies and gradients, the SOC norms show a wide error distribution from −60 to 60 cm⁻¹ with an MAE of 13.4160 cm⁻¹. To determine the sensitivity of ring-opening reaction to the ISC dynamics, we compared the ML-photodynamics with a constant bias by adding and subtracting the MAE to the predicted SOC norm. The 20 ps state population dynamics with reduced SOC norms (Fig. S7b†) show almost identical results to those without modifications on SOC norms (Fig. S7a†), suggesting the NN overestimations of the SOC norms do not significantly affect the S₁ → T₁ ISC rates. Increasing the predicted SOC norm mainly speeds up the T₁ → S₀ ISC after 10 ps but has little effect on the S₁ → T₁ ISC (Fig. S7c†). In all cases, the predicted yields of 4 are similar (Table S7†), suggesting a small influence of the NN prediction errors on the ring-opening reaction mechanism of 1. Thus, the NN errors will not significantly change our results and discussions on the photochemical ring-opening reaction mechanism of 1.

The NN-predicted state population dynamics in 20 ps show that 89% of the trajectories finished at the T₁ state; only 8% trajectories arrived to S₀ state via the S₁ → S₀ (3%), T₁ → S₀ (2%), T₂ → S₀ (3%) transitions (Fig. 5c). These results confirmed the S₁/T₁ intersystem crossing controls the S₁ non-radiative relaxation of 1. We measured the σ_CC bond lengths in 1–4 to understand the structural changes during the S₁ → T₁ transitions in the ML-photodynamics simulations (Fig. 6a).

Fig. 6 (a) Illustration of the σ_CC bonds R₁ and R₂ in 1–4. The scatter plots for the structural distributions of (b) the latest S₁/T₁ surface hopping points and (c) the final snapshots via the S₁ → T₁ transition in the ML-photodynamics simulations. The 2D space is defined by R₁ and R₂.

The structural distributions in Fig. 6b show four clustering regions at the latest S₁/T₁ surface hopping points. 41% of trajectories at region 1 correspond to the trajectories hop to T₁ without breaking the σ_CC bonds. 22% of the trajectories undergo a 4π-disrotatory ring-opening at the carbonyl-functionalized cyclobutene ring at region 2; 1% of the trajectories open the unsubstituted cyclobutene ring at region 3. It suggests the carbonyl group facilitates the ring-opening reaction. Moreover, 25% of the trajectories at region 4 indicate the ring-opening reaction of two σ_CC bonds could also occur during the S₁ → T₁ transitions. These results suggest that the ring-opening reactions are intrinsically favored in the S₁ and T₁ states of 1.

Fig. 6c shows the final structural distributions of the trajectories undergoing the S₁/T₁ surface hoppings. The number of trajectories at regions 1–3 substantially decreased to 2%, 6%, and 0%, respectively. As a result, 80% of the trajectories arrived in region 4. The trajectories forming 4 remained in T₁ after 20 ps due to slow T₁ → S₀ ISC. We considered them as final ground-state products since we did not observe any reverse reaction. The continuously increasing yield of 4 in our simulations suggests the second step ring-opening reaction of 1 is also thermodynamically favored.

The rare S₁ → S₀ conical intersection and T_n → S₀ intersystem crossings generally begin without breaking the σ_CC bonds of 1. The trajectories undergoing the S₁ → S₀ transitions directly afforded the ring-opening reaction. The trajectories in the T_1,2 → S₀ pathways first showed the above-mentioned S₁ → T₁ or S₁ → T₂ transitions and continued to hop to the ground state via the T₁/S₀ or T₂/S₀ intersections. These trajectories collectively contribute to another 7% of 4 on the ground state. In addition, we found 2% of 4 formed in the remaining S₁ trajectories. These results suggest even the minor reaction pathways favor the ring-opening reactions of 1. Our ML photodynamics simulations predicted an overall 89% quantum yield of 4.

We compared the NN-predicted trajectories in all possible reaction pathways to elucidate the ring-opening mechanisms of 1. The trajectory plots in Fig. 7a–c illustrate the changes in the σ_CC bond lengths, R₁ and R₂, as a function of simulation time.

Fig. 7 Plots for the randomly selected trajectories in the 20 ps ML-photodynamics simulations. Three top panels plot the trajectories undergoing the S₁/T₁ intersystem crossings at regions (a) 1, (b) 2, and (c) 4 in a space defined by intramolecular distance R₁ and R₂. Three bottom panels plot the trajectories undergoing the S₁/T₁ intersystem crossings at regions (d) 1, (e) 2, and (f) 4 in a space defined by the intramolecular distance R₁ and the π_CO bond distance R₃. The blue dots represent the latest S₁/T₁ surface hopping points in trajectories. The dark blue circles highlight the locations of the S₁/T₁ surface hopping regions.

Fig. 7a represents the trajectories undergoing the S₁/T₁ intersystem crossings near the S₁-FC region at R₁ = 1.62 Å and R₂ = 1.64 Å (region 1). After hopping to T₁, these trajectories proceeded to the intermediate 2 at R₁ = 2.92 Å and R₂ = 1.64 Å and then the final product 4 at R₁ = 3.34 Å, R₂ = 3.31 Å, showing a stepwise ring-opening mechanism. Fig. 7b shows the first ring-opening step near the intermediate 2 before the S₁ → T₁ transition. Then, the trajectories continued the second ring-opening step to form 4. Fig. 7c highlights the trajectories formed 4 before the S₁/T₁ intersystem crossings. These results suggest a stepwise mechanism of the photochemical ring-opening reactions of 1 that first breaks the σ_CC bond in the ring substituted with the carbonyl groups and then breaks the other σ_CC bond. The S₁/T₁ intersystem crossings can occur at the S₁-FC region, intermediate 2, and product 4 indicates the σ_CC bond-breaking does not correlate with the S₁/T₁ intersystem crossing structures. We plot the π_CO bond distance (R₃) in the carbonyl group with R₁ in the trajectories (Fig. 7d–f) to determine the origin of the S₁/T₁ intersystem crossings. The average value of R₃ is 1.22 Å at the S₁-FC region. It notably increases to 1.31 Å, 1.30 Å, and 1.46 Å at the S₁/T₁ intersystem crossing structures near the S₁-FC (Fig. 7d), intermediate 2 (Fig. 7e), and product 4 regions (Fig. 7f), consistent with a partial triplet electronic structure. Thus, the elongation of the π_CO bond is responsible for the S₁/T₁ intersystem crossings during the photochemical ring-opening reaction of 1.

The trajectories in Fig. 8a show the S₁/S₀ surface hopping structures near the S₁-FC region (R₁ = 1.66 Å and R₂ = 1.66 Å). Like the trajectories hoped to T₁ in Region 1 (Fig. 7a), these trajectories proceeded to a stepwise ring-opening reaction of 1 to form 4. Fig. 8b illustrates the S₁ trajectories hoped to the ground state between the first and second steps of the ring-opening reaction of 1. The S₁/S₀ surface hopping structures resemble MECI-2 with R₁ = 3.14 Å, R₂ = 1.66 Å. The second ring-opening step was completed in the ground state at product 4, with R₁ = 3.14 Å and R₂ = 3.30 Å at the end of the simulations. The trajectories directly formed 4 in S₁ are shown in Fig. 8c, and the trajectories still follow a stepwise ring-opening pathway. Similar to the S₁/T₁ intersystem crossings, we found increased π_CO bond lengths from 1.22 Å at the S₁-FC regions to 1.58 Å, 1.54 Å, and 1.58 Å at the S₁/S₀ conical intersections near the S₁-FC (Fig. 8d), intermediate 2 (Fig. 8e), and product 4 regions (Fig. 8f). These findings indicate the same ring-opening mechanism of 1via the S₁/S₀ conical intersections and S₁/T₁ intersystem crossing but only differ in the time constants.

Fig. 8 Plots for the randomly selected trajectories in the 20 ps ML-photodynamics simulations. Three top panels plot the trajectories undergoing the S₁/S₀ surface hoppings at regions (a) 1, (b) 2, and (c) 4 in a space defined by the intramolecular distance R₁ and R₂. Three bottom panels plot the trajectories undergoing the S₁/S₀ surface hoppings at regions (d) 1, (e) 2, and (f) 4 in a space defined by intramolecular distance R₁ and the π_CO bond distance R₃. The red dots represent the latest S₁/S₀ surface hopping points in trajectories. The dark red circles highlight the locations of the S₁/S₀ surface hopping regions.

Combining the electronic state population distributions and structural classifications in the trajectories produces a comprehensive reaction network for the predicted photochemical electrocyclic ring-opening reaction of 1 toward 4 (Fig. 9). The detailed and quantitative connections over all reactant, intermediates, and products across multiple singlet and triplet states provide a deep understanding of the ring-opening mechanisms.

Fig. 9 Reaction network diagram for the photochemical electrocyclic ring-opening reaction of 1. The nodes with thin black circles correspond to the intermediate species, and the terminal nodes with bold-colored circles are products 1 (gray), 2 (pink), 3 (brown), and 4 (blue). The yields less than 1% are labelled with 0%. The blue arrows indicate the sources of product 4.

Conclusion

We combined the quantum mechanical calculations and machine learning photodynamics approach to enable long timescale photodynamics simulations for uncovering unexpected photochemical reaction pathways. Our showcase study on the stepwise ring-opening mechanisms of 1 show competing S₁/T₁ intersystem crossings and S₁/S₀ conical intersections, where the first ring-opening step starts at the carbonyl-functionalized cyclobutene ring. The CAM-B3LYP/cc-pVDZ calculations show that connecting to a carbonyl group substantially reduced the S₁ energy from 6.25 eV to 3.68 eV. The SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations show a close S₁/T₁ gap (0.2 eV) in the S₁-FC region along with the interpolated ring-opening reaction pathways. The 20 ps ML-photodynamics simulations show dominant S₁/T₁ intersystem crossings with a lifetime of 1.4 ps. The overall predicted quantum yield of 4 is 89%.

The trajectory analysis revealed a diverse reaction channel of the photochemical ring-opening reaction of 1 in S₁ and T₁ states. The ring-opening mechanism favors a stepwise σ_CC bond-breaking in 1, where the cyclobutene substituted with the carbonyl group opens before the ring on the other side. Further characterizations of the trajectories show that the elongation of the π_CO bond length of the carbonyl groups generates the S₁/T₁ intersystem crossing and S1/S0 conical intersections. Thus, the photochemical ring-opening reaction of 1 is independent of whether it is in the S₁ or T₁ states. As such, the S₁ excitation energy became the main driving force to break the σ_CC bonds of 1, which explains the high efficiency of the ring-opening reaction. Overall, our work demonstrates a feasible and light-controllable reaction tool for preparing the COT derivatives in a high yield.

Computational details

Quantum mechanical calculation

We optimized the gas-phase geometries of 1 using the density functional theory (DFT) at the PBE0/cc-pVDZ level. Frequency calculations confirmed the local minimum with no imaginary frequencies. The vertical excitation energies and oscillator strengths were computed at the CAM-B3LYP/cc-pVDZ, CAM-B3LYP/aug-cc-pVDZ and CAM-B3LYP/cc-pVTZ levels using the time-dependent DFT (TD-DFT) theory. The (TD)DFT calculations were carried out by the BDF program.^32–34 The SA5-CASSCF(12,11)/ANO-RCC-VDZP, CASPT2(12,11)/ANO-RCC-VDZP, and XMS-CASPT2(12,11)/ANO-RCC-VDZP calculations were performed with OpenMolcas.19.11.³⁵ The CAPST2 and XMS-CASPT2 used a level shift of 0.4 eV to avoid the intruder state problem.

Training data generation

The initial training data consists of two parts, generated by Wigner sampling and geometry interpolations. We first generate 20 non-equilibrium structures for 1 with the Wigner sampling at the zero-point energy level. Then, we collected 39 structures in the interpolated pathway from 1 to each MECI (Fig. S3†). Mixing the Wigner sampled and interpolated structures in all four pathways gave 3210 data points. The Wigner sampling uses the PyRAI²MD sampling tool;³⁶ the geometry interpolation uses the geodesic interpolation program.^37,38 The training data contain the S₀, S₁, T_1, and T₂ energies, gradients, and the norm of the spin–orbit couplings (SOCs) for S₁–T₂, S₁–T₁, S₀–T₂, and S₀–T₁, computed with the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculations. We performed adaptive sampling to collect the under-sampled data points in the initial training data. Two independently trained NNs were used as a committee model. We used the potential model (i.e., energies and gradients) together with the SOC model to propagate 100 trajectories from the S₁-FC points in 10 ps with a step size of 0.5 fs. We used the standard deviations of NN-predicted energies, gradients, and SOC norms to measure the prediction uncertainty, which early stopped the trajectories when the prediction uncertainty exceeded the empirical thresholds (energy: 0.05 hartree; gradient: 0.25 hartree bohr⁻¹; SOC norm: 60 cm⁻¹). The final data size is 6459. Details about the adaptive sampling are available in Table S5.† All training data are computed with the SA5-CASSCF(12,11)/ANO-RCC-VDZP calculation.

Neural networks

We implemented a fully connected feedforward NNs using TensorFlow/Keras API for Python. The NN computes the inverse distance matrix of the input molecule to predict the energies, gradients, and SOC norms. The predicted energy gap between two singlet or triplet states is used to compute the curvature-approximated time derivative coupling (kTDC),^30,31 derived from the Baeck–An approximation.³⁹ The NNs employ a leaky softplus activation function. The NNs used the first-order energy derivatives of atomic coordinates to fit the gradients, which ensures the physical relationship between energy and gradient. The energies and gradients of four electronic states are trained together with a combined loss function. Four NNs were trained simultaneously, where two NNs predict energy and gradient of four electronic states and two NNs predict SOC norms of four pairs of states. The training data are split into training and validation sets in a 9 : 1 ratio. The final validation mean absolute errors (MAEs) for NN-predicted energies are 0.0337–0.0340 eV. Additional information about the NN training and error evaluation is available in ESI.†

ML-photodynamics simulations

We prepared the initial conditions for 1 by Wigner sampling at the zero-point energy level. The ML-photodynamics simulations were set to 20 ps with a timestep of 0.5 fs. The probability of a nonadiabatic electronic transition was computed with Tully's fewest switches surface hopping (FSSH) method^40,41 using spin-diabatic representation, where the generalized formalism²⁹ was used to take the intersystem crossing into account. We used the curvature-approximated time-derivative coupling (kTDC) method^30,31 to evaluate the nonadiabatic couplings based on the NN-predicted energies when the S₁–S₀ and T₂–T₁ gaps are smaller than 0.5 eV. The ML-photodynamics approach is implemented in PyRAI²MD.³⁶ The nuclear timestep is 0.5 fs and the electronic timestep is 0.025 fs corresponding to 20 substeps in the surface hopping calculations.

Data availability

The code for PyRAI2MD can be found at https://github.com/mlcclab/PyRAI2MD-hiam. The data supporting this article have been included as part of the ESI† and can be found at https://github.com/mlcclab/PyRAI2MD_publications/tree/main/triplet_COT

Author contributions

Z. L. performed all calculations and prepared the first draft of the manuscript. Z. L. and J. L. analyzed the data. H. F. helped revise the manuscript. S. A. L. and J. L. conceived the idea of this work and revised the manuscript. J. L. oversaw the administration of this project.

Conflicts of interest

The authors delcare no conflicts of interests.

Acknowledgements

J. L. acknowledges the funding support by the National Natural Science Foundation Grants (22303053). Z. L. and J. L. thank the support of the high-performance computing resources at Hoffmann Institute of Advanced Materials at Shenzhen Polytechnic University. Z. L. and J. L. gratefully acknowledge HZWTECH for providing computation facilities. S. A. L. acknowledges the NSF (NSF-CHE-2144556) for funding. J. L. and S. A. L. appreciate the assistance from the Northeastern Research Computing Team and the computing resources provided by the Massachusetts Life Science Center Grant (G00006360).

References

1. L. Kelebekli, M. Çelik, E. Şahin, Y. Kara and M. Balci, Stereospecific synthesis of a new class of aminocyclitol with the conduramine D-2 configuration, Tetrahedron Lett., 2006, 47(39), 7031–7035.
2. L. Kelebekli, Y. Kara and M. Balci, Stereospecific synthesis of a new class of compounds: bis-homoconduritol-A,-D, and-F, Carbohydr. Res., 2005, 340(12), 1940–1948.
3. A. Xiong, G. Ma, K. Maeda, T. Takata, T. Hisatomi, T. Setoyama, J. Kubota and K. Domen, Fabrication of photocatalyst panels and the factors determining their activity for water splitting, Catal. Sci. Technol., 2014, 4(2), 325–328.
4. T. Nishinaga, T. Ohmae, K. Aita, M. Takase, M. Iyoda, T. Arai and Y. Kunugi, Antiaromatic planar cyclooctatetraene: a strategy for developing ambipolar semiconductors for field effect transistors, Chem. Commun., 2013, 49(47), 5354–5356.
5. 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, 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.
6. T. NatháDas and K. IndiraáPriyadarshini, Triplet of cyclooctatetraene: Reactivity and properties, J. Chem. Soc., Faraday Trans., 1994, 90(7), 963–968.
7. Q. Zheng, S. Jockusch, Z. Zhou, R. B. Altman, J. D. Warren, N. J. Turro and S. C. Blanchard, On the mechanisms of cyanine fluorophore photostabilization, J. Phys. Chem. Lett., 2012, 3(16), 2200–2203.
8. V. T. Mai, V. Ahmad, M. Mamada, T. Fukunaga, A. Shukla, J. Sobus, G. Krishnan, E. G. Moore, G. G. Andersson and C. Adachi, Solid cyclooctatetraene-based triplet quencher demonstrating excellent suppression of singlet–triplet annihilation in optical and electrical excitation, Nat. Commun., 2020, 11(1), 5623.
9. R. Willstätter and E. Waser, Über Cyclo-octatetraen, Ber. Dtsch. Chem. Ges., 1911, 44(3), 3423–3445.
10. J. H. Rigby and N. C. Warshakoon, A convenient synthesis of 1, 2-disubstituted cyclooctatetraenes, Tetrahedron Lett., 1997, 38(12), 2049–2052.
11. C. Wang, J. Yuan, G. Li, Z. Wang, S. Zhang and Z. Xi, Metal-mediated efficient synthesis, structural characterization, and skeletal rearrangement of octasubstituted semibullvalenes, J. Am. Chem. Soc., 2006, 128(14), 4564–4565.
12. Y. Yamamoto, T. Ohno and K. Itoh, Selective synthesis of fused cyclooctatetraenes by [4 + 4] coupling between two different diene units, Chem.–Eur. J., 2002, 8(20), 4734–4741.
13. S. Majumder and A. L. Odom, Simple and convenient one-pot synthesis of cyclooctatetraene, Tetrahedron Lett., 2008, 49(11), 1771–1772.
14. W. Reppe, O. Schlichting, K. Klager and T. Toepel, Cyclisierende polymerisation von acetylen I über cyclooctatetraen, Justus Liebigs Ann. Chem., 1948, 560(1), 1–92.
15. P. A. Wender, J. P. Christy, A. B. Lesser and M. T. Gieseler, The Synthesis of Highly Substituted Cyclooctatetraene Scaffolds by Metal-Catalyzed [2 + 2 + 2 + 2] Cycloadditions: Studies on Regioselectivity, Dynamic Properties, and Metal Chelation, Angew. Chem., Int. Ed., 2009, 48(41), 7687–7690.
16. P. A. Wender and J. P. Christy, Nickel (0)-catalyzed [2 + 2 + 2 + 2] cycloadditions of terminal diynes for the synthesis of substituted cyclooctatetraenes, J. Am. Chem. Soc., 2007, 129(44), 13402–13403.
17. S. Blouin, V. Gandon, G. Blond and J. Suffert, Synthesis of Cyclooctatetraenes through a Palladium-Catalyzed Cascade Reaction, Angew. Chem., Int. Ed., 2016, 55(25), 7208–7211.
18. J. Xuan, X.-K. He and W.-J. Xiao, Visible light-promoted ring-opening functionalization of three-membered carbo-and heterocycles, Chem. Soc. Rev., 2020, 49(9), 2546–2556.
19. Y. Zheng, W. Huang, R. K. Dhungana, A. Granados, S. Keess, M. Makvandi and G. A. Molander, Photochemical intermolecular [3σ+ 2σ]-cycloaddition for the construction of aminobicyclo [3.1. 1] heptanes, J. Am. Chem. Soc., 2022, 144(51), 23685–23690.
20. D. P. Poudel, A. Pokhrel, R. K. Tak, M. Shankar and R. Giri, Photosensitized O2 enables intermolecular alkene cyclopropanation by active methylene compounds, Science, 2023, 381(6657), 545–553.
21. R. Gleiter and S. Brand, Photochemistry of Bridged and Unbridged Octaalkyl-Substituted syn-Tricyclo [4.2. 0.02, 5] octa-3, 7-diene Derivatives, Chem.–Eur. J., 1998, 4(12), 2532–2538.
22. 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.
23. 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.
24. J. Xie, X. Zhang, C. Shi, L. Pan, F. Hou, G. Nie, J. Xie, Q. Liu and J.-J. Zou, Self-photosensitized [2+ 2] cycloaddition for synthesis of high-energy-density fuels, Sustainable Energy Fuels, 2020, 4(2), 911–920.
25. M. A. Theodoropoulou, N. F. Nikitas and C. G. Kokotos, Aldehydes as powerful initiators for photochemical transformations, Beilstein J. Org. Chem., 2020, 16(1), 833–857.
26. S. Maeda, K. Ohno and K. Morokuma, A theoretical study on the photodissociation of acetone: Insight into the slow intersystem crossing and exploration of nonadiabatic pathways to the ground state, J. Phys. Chem. Lett., 2010, 1(12), 1841–1845.
27. L. Favero, G. Granucci and M. Persico, Dynamics of acetone photodissociation: a surface hopping study, Phys. Chem. Chem. Phys., 2013, 15(47), 20651–20661.
28. B. G. Levine and T. J. Martínez, Isomerization through conical intersections, Annu. Rev. Phys. Chem., 2007, 58(1), 613–634.
29. G. Cui and W. Thiel, Generalized trajectory surface-hopping method for internal conversion and intersystem crossing, J. Chem. Phys., 2014, 141(12), 124101.
30. 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.
31. X. Zhao, I. C. 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.
32. W. Liu, F. Wang and L. Li, Relativistic density functional theory: The BDF program package, Recent advances in relativistic molecular theory, 2004, pp. 257–282.
33. Y. Zhang, B. Suo, Z. Wang, N. Zhang, Z. Li, Y. Lei, W. Zou, J. Gao, D. Peng and Z. Pu, BDF: A relativistic electronic structure program package, J. Chem. Phys., 2020, 152(6), 064113.
34. W. Liu, F. Wang and L. Li, The Beijing density functional (BDF) program package: methodologies and applications, J. Theor. Comput. Chem., 2003, 2(02), 257–272.
35. I. Fdez Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov and R. K. Carlson, OpenMolcas: From source code to insight, J. Chem. Theory Comput., 2019, 15(11), 5925–5964.
36. 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.
37. X. Zhu, K. C. Thompson and T. J. Martínez, Geodesic interpolation for reaction pathways, J. Chem. Phys., 2019, 150(16), 164103.
38. M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean and M. Devin, Tensorflow: Large-scale machine learning on heterogeneous distributed systems, arXiv, 2016, preprint, arXiv:1603.04467,  DOI:10.48550/arXiv.1603.04467.
39. K. K. Baeck and H. An, Practical approximation of the non-adiabatic coupling terms for same-symmetry interstate crossings by using adiabatic potential energies only, J. Chem. Phys., 2017, 146(6), 064107.
40. J. C. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys., 1990, 93(2), 1061–1071.
41. S. Hammes-Schiffer and J. C. Tully, Proton transfer in solution: Molecular dynamics with quantum transitions, J. Chem. Phys., 1994, 101(6), 4657–4667.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sc07395a

---