From excited-state dynamics to product selectivity: theoretical insights into the photoinduced chemistry of aza-diarylethenes

Abstract

Aza-diarylethenes (aDAEs) are a novel class of photoswitches, which undergo reversible C–N bond formation upon irradiation. Previous experimental studies on aDAEs, differing in the heterocyclic aryl cores and substituents at the sp³-carbons of the cyclopentene, have revealed the formation of distinct types of photoproducts, on the one hand, DAE-typical closed ring isomers and on the other hand hitherto unknown zwitterionic structures. Here, we present a comprehensive comparative computational study of aDAE photocyclization mechanisms using nonadiabatic molecular dynamics simulations at the mixed-reference spin–flip time-dependent density functional theory level, to provide a unified mechanistic framework for understanding the divergent photoproduct formation. We show that photoexcitation drives ring closure via an isomerization conical intersection, but this crossing does not necessarily define the experimentally observed product: in one system it yields the closed isomer directly, whereas in another aDAE it accesses a ground-state pathway leading to the formation of a zwitterionic isomer.

---

1 Introduction

Photoswitches are key molecular tools in modern chemistry, enabling precise external control over molecular states through light irradiation. Due to their remarkable ability to undergo highly efficient and reversible switching, diarylethenes (DAEs) emerged as one key class of photoswitches over the past few decades.^1–5 They stand out due to their high thermal stability of both isomers, large photochromism and high photofatigue resistance.^2,4 These characteristics make diarylethenes particularly valuable for applications in optical data storage,^6,7 catalysis,⁸ host–guest interactions⁹ and photoresponsive materials.¹⁰ In typical DAEs, the open-ring isomer represents the thermodynamically favored state. Upon irradiation with ultraviolet light, it undergoes a photocyclization reaction, yielding the closed-ring isomer (see Fig. 1a). The closed form can subsequently be converted back to the open form through irradiation with lower-energy light, typically in the visible range. This process, occurring on the picoseconds timescale, is associated with little geometric but significant electronic rearrangements. In particular, there is an expansion of the π-system in the closed form, which leads to red-shifted absorption.

Fig. 1 Structural formulas of the diarylethene introduced by the group of Irie³³ (a) and the aza-diarylethenes presented by the group of Kobatake²⁶ (K, b) and the group of Dube²⁷ (D, c). All molecules are shown in their open form (left) and their experimentally observed photoproduct, i.e., either the closed or zwitter-ionic (ZI) form (right).

Since the pioneering work of Irie and Mohri¹ on thiophene-substituted diarylethenes, the field has expanded significantly and large scientific efforts have been made with respect to structural modifications to further improve the properties of DAEs. Isomerization mechanisms in DAEs have been subject to numerous experimental-spectroscopic^11–16 as well as computational studies.^17–21

The ring-closing reaction of DAEs is a 6π-electrocyclic reaction and in the simplest picture, when only considering the molecule core, the cyclization of all-cis-hexatriene (cHT) to 1,3-cyclohexadiene (CHD). This model system was revisited by Salazar and Faraji¹⁷ in 2020 using static spin–flip time-dependent density functional theory (SF-TDDFT). They showed that upon excitation of the open form, the bright ππ* S₁ state is populated, with S₂ being a dark ππ* state. An area of strong S₁/S₂ interactions is passed and, ultimately, a S₁/S₀ conical intersection (CI) is reached (nearly) barrier-less.¹⁷ The ring-opening and ring-closing reactions between cHT and CHD remain active topics of ongoing research.^22–25 In 2022, Salazar et al.¹⁸ investigated a non-symmetric DAE photoswitch within the same static theoretical framework. They found that upon light absorption, the open-ring isomer is excited to its bright S₂ state instead of the S₁ in the minimal model. Vibrational relaxation leads to C–C bond shortening and reaching a S₂/S₁ CI, from where two pathways are possible, which however both lead towards the same S₁/S₀ CI.¹⁸ Recently, the first full-size multireference dynamics study of a DAE switch was presented, where non-adiabatic molecular dynamics (NAMD) simulations were performed at ODM2/MRCI (multireference configuration interaction based on orthogonalization and dispersion corrected semiempirical methods) level.^19,20 The S₁ state was found to be the bright ππ* transition, and cyclization was described straightforwardly as a relaxation on the S₁ PES towards a S₁/S₀ CI. S₂ populations were found in several trajectories, but with short lifetimes only.

More recently, a novel subclass of diarylethenes, known as aza-diarylethenes (aDAEs), has emerged (see Fig. 1b and c). In these systems, traditional C–C bond connections are replaced with C–N bonds by substituting one of the thiophene rings in typical diarylethenes with nitrogen-containing heterocycles, introducing new electronic properties and reactivity patterns. This structural modification has led to the discovery of two different types of aDAEs:^26,27

On the one hand, aDAEs reported by the Kobatake group exhibit photochemical behavior closely resembling that of conventional diarylethenes, forming a similar closed-ring product upon irradiation.²⁶ Moreover, this group synthesized a novel aDAE that undergoes both photochemical and thermal ring-opening and ring-closing reactions.²⁸ Recently, they demonstrated effective control over the lifetimes of the closed-ring isomers by systematically varying the bridging units and heteroaryl cores.²⁹ Notably, in one representative system, K-open, which serves as the model compound in the present study, one thiophene ring is replaced by a thiazole, enabling formation of a closed-ring isomer through a C–N bond between the thiophene and thiazole units (see K-open and K-closed in Fig. 1b).

On the other hand, aDAEs reported by the Dube group represent a structurally distinct class, in which one thiophene unit is replaced by a pyrimidine moiety.²⁷ These systems further differ in the substitution patterns of the heterocycles and by the presence of hydrogen, rather than fluorine, at the central cyclopentene ring. These aDAEs have been shown to form unexpected, metastable zwitterionic isomers with an aromatic core rather than the conventional closed-ring structure (see D-open and D-ZI in Fig. 1c). The resulting zwitterionic species (D-ZI) features a negatively charged sulfur and a positively charged nitrogen atom.²⁷ This zwitterion-forming behavior endows D-ZI with distinctive reactivity and offers unique opportunities for manipulating chemical transformations with light. Specifically, these zwitterions enable double-bond isomerization,^27,30 act as photoinitiators for radical polymerization,³¹ and facilitate light-driven generation of thiol/thiolate moieties, which are key functionalities for organic synthesis, bioconjugation, and materials chemistry.³²

Despite the remarkable properties of these novel systems and their continued development in recent studies, no comprehensive study of the photoreaction mechanisms in aDAEs, nor of the factors governing the formation of different photoproducts, has been reported to date. Therefore, the present study aims to elucidate the mechanistic pathways leading to the formation of distinct photoproducts and to identify the structural features that determine the reaction outcome. To this end, we utilize mixed-reference spin–flip (MRSF) TDDFT for static and NAMD calculations to explore the complex photochemical pathways of different types of aDAEs. MRSF-TDDFT is a novel and promising method for excited-state studies providing a spin-adapted, multiconfigurational description that captures static and dynamical correlation and doubly excited configurations within a computationally efficient DFT framework.

Our results show that both D-open and K-open follow similar reaction pathways in the excited state. However, the experimentally observed photoproduct is not determined by the excited-state energy landscape alone but also by the topology of the ground-state surface, including relative stabilities and barrier heights. After passage through the photoisomerization conical intersection, relaxation on the ground state leads either to the closed-ring isomer (K-closed) or, in the case of D, to formation of the zwitterionic species (D-ZI).

2 Methods

All density functional theory (DFT), linear-response time-dependent DFT (TDDFT) and equation-of-motion coupled cluster singles-doubles (EOM-CCSD) calculations were performed in Orca 5.0,^34,35 all mixed-reference spin–flip (MRSF) calculations were done in OpenQP.^36–39 For benchmarking purposes, second-order approximate coupled cluster singles and doubles (CC2) calculations, algebraic diagrammatic construction of second order (ADC(2)) calculations as well as Møller–Plesset perturbation theory up to second-order (MP2) calculations were performed in Turbomole 7.7.⁴⁰

2.1 Static calculations

First, we optimized ground-state (S₀) minima of the open (-open), closed (-closed) and zwitterionic (-ZI) configurations of D and K, namely D-open, D-closed, D-ZI and K-open, K-closed, K-ZI were obtained using DFT with the CAM-B3LYP⁴¹ functional and the def2-TZVP⁴² basis set. To account for solvent effects the Conductor-like Polarizable Continuum Model (CPCM)⁴³ of methanol was used to align with the experimental work for D,²⁷ most of which was performed in methanol (note that the experiments for K mentioned in this work were performed in hexane¹⁰).

Subsequently, the ground-state absorption properties were computed for all S₀ minima by means of TDDFT at the same level of theory as for the ground state calculations, using non-equilibrium linear-response (LR-)CPCM for implicit solvation (see Fig. 2). For benchmarking purposes, we further computed absorption properties for both open and closed isomers using MRSF TDDFT (hereafter referred to as only MRSF for simplicity), EOM-CCSD and ADC(2). For MRSF, we used the DTCAM-AEE (doubly-tuned Coulomb-attenuating method adaptive exact exchange) functional, which uses the Coulomb attenuating method and different amounts of exact exchange in the SCF and response part of the MRSF calculation. DTCAM-AEE has been designed to accurately describe conical intersections as well as NAMD simulations. Moreover, it has been shown to provide an improved description of charge-transfer states compared to commonly used functionals in linear-response TDDFT.⁴⁴ ADC(2) calculations were performed using the def2-TZVP basis set, with structures re-optimized at the MP2/def2-TZVP level of theory. For EOM-CCSD calculations, we used the domain-based local pair natural orbital (DLPNO) similarity transformed implementation with the more economical def2-SVP basis set, with structures re-optimized at the CC2/def2-TZVP level of theory. All MRSF, ADC(2), MP2, CC2 and EOM-CCSD calculations were performed in gas phase.

Fig. 2 Gaussian broadened (σ = 0.3 eV) absorption spectra for the open, closed and ZI isomers of D (a, top) and K (a, bottom) calculated with TD-CAM-B3LYP/def2-TZVP/CPCM(MeOH). Vertical excitation energies shown as sticks, and their lowest singlet excited states are highlighted. Visualization of the HOMO and LUMO orbitals of D-open (b, top) and K-open (b, bottom) indicating the main contribution of the bright S₁ state, calculated at the MRSF-TD-DTCAM-AEE/6-31G* level of theory.

Upon identifying a level of theory suitable for describing the Franck–Condon region of D-open and K-open, we calculated energy profiles along two reaction coordinates, the C–N distance connecting both heterocycles and the C–S distance in the thiophene ring (see yellow highlighted atoms in Fig. 3). Therefore, we performed relaxed scans at CAM-B3LYP/def2-TZVP/CPCM(MeOH) level of theory in both the S₀ and S₁ state. In the S₁ scans, LR-CPCM was used with equilibrium conditions. Moreover, we evaluated the single point energies at these S₁-relaxed geometries using DTCAM-AEE-MRSF-TDDFT/6-31G*, allowing for a comparison of the PESs. Based on the geometries obtained in the relaxed scans, we optimized S₁/S₀ minimum energy conical intersection points (MECI, hereafter simply referred to as CI) using DTCAM-AEE-MRSF-TDDFT/6-31G*.⁴⁵

Fig. 3 Potential energy curves of S₀ (solid lines, diamond symbols) and S₁ (solid lines, circle symbols) along the C–N (a) and C–S distance coordinates (b) for D (red) and K (blue), obtained from separate relaxed scans performed in these two singlet states using (TD-)CAM-B3LYP/def2-TZVP/CPCM(MeOH). Vertical S₀ energies at the S₁ relaxed geometries are indicated by the translucent lines and circle symbols. Therefore, vertical circle symbols correspond to identical geometries (within one molecule), whereas points with diamonds correspond to slightly different geometries. The location of D-CI_main and K-CI_main, optimized with MRSF-TDDFT using DTCAM-AEE/6-31G*, is indicated by star symbols in red and blue, respectively. Optimized transition state structures are highlighted with the cross symbols. Schematic representations of open, closed and ZI isomers are shown as insets, with the nitrogen, carbon and sulfur atoms of the scan coordinates highlighted in yellow. The orbital inset depicts the LUMO of D-open in the S₁ geometry as an example, revealing a favorable phase alignment of the π-contributions on the carbon and nitrogen atoms, consistent with a conrotatory ring closure in the S₁ state.

2.2 Non-adiabatic molecular dynamics simulations

For the non-adiabatic molecular dynamics (NAMD) simulations, we employed trajectory surface hopping (TSH) with the Zhu-Nakamura algorithm⁴⁶ as implemented in PyRAI2MD-hiam.^47–52 As electronic structure method, we chose MRSF, taking into account the five lowest singlet states using the DTCAM-AEE/6-31G* level of theory in gas phase. Among these, only the three lowest states (S_0,1,2) were included as active states in the surface hopping dynamics. We allowed hops between electronic states only at energy differences of 0.5 eV or below. The relatively new MRSF approach requires a restricted open-shell Hartree–Fock (ROHF) reference triplet state in order to access the response singlet states. As mentioned in earlier works, ROHF convergence is highly sensitive to small geometric variations.^53,54 Inspection of our trajectories revealed large energy jumps in the active state between consecutive steps in most trajectories, which we attribute to the ROHF reference converging to an incorrect solution at certain geometries. Following the implementation of routine checks in SHARC,^55,56 we consider any energy jump exceeding 0.7 eV to be problematic. Examining these jumps, we were unable to trace them back to specific geometric features, as geometries still appear reasonable at the problematic steps. The all-atom RMSDs between sequential geometries never exceeded 0.04 Å for D and 0.03 Å for K, and averaged below 0.02 Å for both. We refrained from completely discarding the affected trajectories, as the total number of trajectories required would otherwise be computationally very demanding. The fraction of problematic steps in all analyzed trajectories was around 0.5% for D and around 0.9% for K. Comparing productive and unproductive trajectories, they were affected to almost identical amounts in both cases (D: 0.6% for productive, 0.5% for unproductive trajectories; K: 1.0% for productive, 0.9% for unproductive trajectories). We will therefore assume in the following that isolated events throughout the trajectories and across the ensemble are unlikely to affect the overall outcome or interpretation.

Initial conditions were generated based on the equilibrium geometry of D-open and K-open using Wigner sampling to obtain non-equilibrium geometries and we initialized all trajectories in S₁. We propagated the trajectories with a timestep of 0.5 fs up to a maximum of 2 ps for compound D and up to 1 ps for compound K due to their different photoreaction kinetics. For both D and K, 97 trajectories were successfully completed and are included in the analysis. For D, we had initialized three more, which were terminated prematurely due to convergence issues and unphysical behavior. For further analysis, we categorized our trajectories into productive trajectories describing the open → closed reaction and unproductive ones, including photophysical trajectories (i.e., open → open) and the formation of rare side products (ring-opened derivatives of open isomers, <2%).

To capture the ground-state dynamics of the closed isomers and investigate potential ZI isomer formation while including solvent effects, we performed ground-state Born–Oppenheimer molecular dynamics (BOMD) simulations. A similar approach has been used earlier to explain the formation of different photoinduced products, crossing transition states far from thermal equilibrium, of sulfine (H₂CSO).⁵⁷ We started these simulations with the geometries and velocities extracted from our surface hopping dynamics, using the snapshot at 10 fs after the hop to the ground state for each individual productive trajectory. We propagated the trajectories for 5 ps on S₀. For the BOMD, we used CAM-B3LYP/6-31G*/CPCM(MeOH) as electronic structure method and the dynamics framework provided by SHARC 3.^55,58

3 Results and discussion

3.1 Ground-state and absorption properties

Among the three studied conformations (open, closed, and ZI), the structures of D-open and K-open represent the global minimum geometries, consistent with experimental observations^10,27 and conventional DAEs.^20,26 The relative stability of the remaining two isomers differs between the compounds (with respect to the open form): K-closed (18.0 kcal mol⁻¹) is more stable than K-ZI (26.8 kcal mol⁻¹), whereas D-ZI (12.8 kcal mol⁻¹) is more stable than D-closed (21.0 kcal mol⁻¹). This energetic order of ground state minima provides first insights into why different photoproducts are observed in the experiments, namely, K-closed and D-ZI.^10,27 The same energetic order was qualitatively already proposed by the group of Dube,²⁷ but they reported smaller energy differences between the isomers (14.8 kcal mol⁻¹ for D-closed and 9.5 kcal mol⁻¹ for D-ZI.

To investigate the effect of the polar solvent methanol on relative stabilities, we optimized all isomers again in gas phase. The relative stability of the closed isomers with respect to the open forms remains largely unchanged, but the ZI-isomers are significantly destabilized (38.3 kcal mol⁻¹ for K in gas phase and 34.4 kcal mol⁻¹ in hexane, compared to 26.8 kcal mol⁻¹ in methanol, and 26.4 kcal mol⁻¹ for D in gas phase, compared to 12.8 kcal mol⁻¹ in methanol, see Fig. S3 and S4). This indicates a large effect of a polar solvent on the stability of the ZI-isomers, but not on the closed isomers.

Our excitation analysis shows that both, K-open and D-open have a bright ππ* S₁ state with charge-transfer (CT) character. It consists mostly of a HOMO → LUMO transition (see Fig. 2 and Tables S1–S4). Electron density is shifted from the π-orbital on the central double bond in the cyclopentene to the corresponding π*-orbital of that same double bond (cf. HOMO and LUMO orbitals in Fig. 2b and charge density differences in Tables S3 and S4). Moreover, upon excitation, electron density overall decreases on the thiophene and increases on the second heterocycle (pyrimidine in D and thiazole in K). These findings are in line with previous reports on different types of diarylethenes, which concluded that the C₂-symmetry at the central double bond in the LUMO together with strong contribution of frontier orbitals at the six pro-cycle atoms is a common feature of diarylethenes undergoing photocyclization reactions.^59,60

Absorption in K is red-shifted relative to D; for instance, the S₁ state of K-open (329 nm, f = 0.242) is red-shifted by 0.57 eV with respect to the S₁ state of D-open (286 nm, f = 0.201). Although the electron-withdrawing nature of the fluorine substituents in K (−I effect) stabilizes the HOMO and would typically lead to a larger HOMO–LUMO gap and hence a blue shift, the opposite trend is observed. This behavior can be rationalized by the different substituents on the heterocycles, most notably, the extended π-system introduced by the phenyl group in K, which most likely contributes to the observed shift.

In K-open, the lowest-energy absorption band also exhibits higher intensity than in D-open (note the different y-axis scalings in Fig. 2a). This enhancement arises from the presence of three bright excited states in K, in contrast to the predominantly single bright state in D. The additional contributions in K include ππ* excitations localized on the phenylthiophene unit (S₃, 278 nm, f = 0.498) and charge-transfer transitions involving the central cyclopentene π-bond as acceptor and both ring substituents (thiazole and phenylthiophene) as donors (S₂, 290 nm, f = 0.386). The maximum of the calculated lowest-energy absorption band of K-open agrees well with the experimentally obtained 297 nm, and the calculated 286 nm S₁ VEE for D-open is in good alignment with the experimentally observed shoulder around 300 nm.

Absorption bands of closed and ZI isomers are in both compounds significantly red-shifted (see Fig. 2a) compared to the open isomers, as observed experimentally,^26,27 which can be readily explained by the expansion of their π-systems. The lowest singlet states in D-closed and K-closed exhibit similar excitation energies (with K-closed absorbing +0.13 eV higher than D-closed). However, the S₁ state has a small oscillator strength in D-closed, such that again, absorption of K-closed appears red-shifted compared to D-closed (see Fig. 2 and Tables S1, S2). The calculated S₁ energy (470 nm) of K-closed aligns well with the experimentally measured absorption maximum at 527 nm,²⁶ whereas absorption of D-closed was not accessible with static absorption measurements.²⁷ Both S₁ states are ππ*-transitions with mainly HOMO–LUMO contribution, with HOMO and LUMO distributed over the whole π-system. Comparing the frontier orbitals of closed and open isomers, the character of HOMO and LUMO at the central bond in the five-membered ring flips (see Fig. S1), as can be derived for hexatriene to cyclohexadiene ring closure from molecular orbital theory¹⁷ and was observed for symmetrical as well as asymmetrical C–C diarylethenes in other works.^18–21,60 LUMOs in open isomers show a nodal plane at the cyclopentene double bond, and a similar contribution is observed in their HOMOs of closed isomers. K-ZI absorbs red-shifted compared to D-ZI (−0.43 eV) but appears dark, such that their first absorption bands are at similar energies. Again, their S₁ states are ππ*-transitions with mainly HOMO–LUMO contribution, with their HOMOs located at the opened thiophene and LUMOs distributed over the remaining π-system (see Tables S3 and S4 and Fig. S1). The qualitative features of our calculated absorption spectrum of D-ZI, an intense absorption maximum around 300 and a less intense maximum at 425 nm, agree well with experimental spectra.²⁷ Moreover, our calculated spectra for D-open and D-ZI agree well with calculated spectra from the group of Dube.²⁷

3.2 Potential energy curves

In typical DAEs, the distance of the ring-closing carbon atoms is considered the main reaction coordinate of photocyclization.^17–21 To gain initial insights into the potential energy landscape of aDAEs, we therefore performed relaxed scans along the corresponding C–N bond in both the S₀ and S₁ state, decreasing stepwise from the open- to the closed-form minima of K and D (cf. Fig. 3a). Relaxed scans of the C–S distance are employed to model a reaction between the closed and ZI isomers²⁷ (cf. Fig. 3b). Transition states were globally optimized based on the maxima along the reaction coordinates identified in these relaxed scans to obtain reaction barriers.

Open → closed. In the ground-state (S₀), the conversion from the open to the closed form of compounds D and K is characterized by an energy barrier of 36.1 kcal mol⁻¹ and 32.4 kcal mol⁻¹, respectively. For the thermal back-reaction from closed to open form, we obtain barriers of 15.1 kcal mol⁻¹ and 14.4 kcal mol⁻¹ for D and K, where the latter is in excellent agreement with the experimentally measured barrier of 14.1 kcal mol⁻¹ for K.²⁶ For compound D, only the barrier for D-ZI → D-open is experimentally accessible, which we will discuss below. According to the Woodward–Hoffmann rules, a 6π-electrocyclic reaction on the ground state is expected to proceed via a disrotatory pathway. However, inspection of the geometries along the reaction coordinate reveals that the individual motions of the two rings are at the same time subtle and too complex to be unambiguously classified as purely disrotatory or conrotatory. Analyzing the frontier orbitals, in this case, the HOMO, along with geometries, reveals that the phases of the π-contributions on the ring-closing nitrogen and carbon atoms are oriented in parallel, which points towards a disrotatory mechanism. However, the HOMO exhibits only small π-contributions at the nitrogen atom. Instead upon shortening of the C–N distance, the nitrogen lone pair seems to play a role in the ring-closing process (see supplementary information movies M1 and M2 for D and K, respectively). This is in line with an earlier study by the group of Kobatake, which concluded that the geometric rearrangements during ring-closing on the ground-state cannot be described as purely disrotatory and that the lone pair might contribute.²⁸

On the S₁ surface, the open-to-closed ring transformation proceeds without an energetic barrier for both compounds. Along this reaction coordinate, the S₁ and S₀ potential energy surfaces approach closely at d(C–N) of approximately 190 pm for K and 188 pm for D, indicating the presence of conical intersections (CIs) in these regions and efficient relaxation pathways to S₀. This suggests that both K-closed and D-closed are accessible photoproducts, although for D the closed form has not been experimentally reported²⁶ and may represent only a quasi-stationary intermediate. Analysis of the geometries and frontier orbitals along the d(C–N) coordinate reveals that ring closure on S₁ proceeds in a conrotatory manner, in accordance with the Woodward–Hoffmann rules for 6π electrocyclic reactions. In both D-open and K-open, the LUMOs display a favorable phase relationship between the π-orbitals at the nitrogen and thiophene carbon atoms (see Fig. 2b and Fig. 3a), which evolves into constructive orbital overlap as d(C–N) decreases (see SI M3 and M4 for D and K, respectively). This is in alignment with an earlier study by the group of Kobatake, which describes the photochemical ring-closure in an aDEA as conrotatory based on the molecular geometrical change observed in a relaxed scan.²⁸

For both compounds, the potential energy curves obtained from single-point calculations using MRSF-TD-DTCAM-AEE/6-31G* at the geometries along the relaxed scan agree well with the energies obtained from (TD-)DFT, verifying the use of (TD-)DFT for investigation of the potential energy surfaces (see Fig. S5) despite the well-known problems with single-reference methods.

Closed → ZI. The formation of the ZI species from the closed isomer is revealed in the potential energy curves along the d(C–S) coordinate (see Fig. 3b). They show that this transformation is for both compounds energetically accessible on S₀ and S₁ but that there are no other areas of strong S₀/S₁ interaction along this reaction coordinate. This suggests that the zwitterionic isomer arises from the closed-ring structure on the ground state after its photochemical formation from D-open or K-open.

On S₀, the closed → ZI conversion faces a barrier of 9.3 kcal mol⁻¹ and 17.0 kcal mol⁻¹ for D and K, respectively. We therefore hypothesize that D-closed could, in principle, be isolated at temperatures around −130 °C, provided it does not merely constitute a quasi-stationary species in the D-open → D-closed → D-ZI reaction cascade. In contrast, experiments conducted at −78 °C have only been able to isolate D-ZI, a finding that is consistent with both our and previously reported computed barriers.²⁷

For the thermal D-ZI → D-closed back-reaction, the barrier was not accessible experimentally, but instead, a barrier of 20.4 kcal mol⁻¹ was determined for D-ZI → D-open.²⁷ From our ground state PESs, we observe that with respect to D-ZI, the energetically highest point along our proposed reaction path is the transition state between D-closed and D-open. Therefore, we approximate the barrier for this two-step thermal reaction simply by the energetic difference between this transition state and D-ZI, resulting in 23.3 kcal mol⁻¹, which aligns well with the experimental result. With the same approximation, the group of Dube has suggested a barrier of 19.5 kcal mol⁻¹ for this process.

Solvent effects. To assess the effect of the polar solvent (methanol) on the reaction pathway, we further analyzed the potential energy curves in gas phase for D and K, and in hexane for K (see Fig. S3 and S4). The PES along the d(C–N) coordinate showed only minor changes, indicating little solvent influence. In contrast, the d(C–S) coordinate exhibited pronounced differences for gas phase and hexane compared to methanol, consistent with the strong solvation effect on ZI isomer stability. In gas phase and hexane, ZI isomers are less stable, and deviations from PESs in methanol appear already near the closed isomer at d(C–S) ≈ 230 pm. Thus, we neglected solvation when monitoring the open → closed reaction in S₁ (surface hopping dynamics) but included it in the ab initio S₀ dynamics immediately after hops to account for the correct description of the closed → ZI reaction.

3.3 Conical intersections

For D and K, our conical intersection search along the potential isomerization pathway resulted in MECI points, namely D-CI_main and K-CI_main (cf. structural insets in Fig. 4), that lie 2.0 eV and 1.4 eV below their corresponding Franck–Condon (FC) points, respectively (cf. star symbols in Fig. 3; we denote here that the first point on the S₁ curve is not identical to the FC geometry due to (partial) relaxation in the S₁ state). These MECI geometries show large structural similarities to their corresponding closed isomers. The primary distinction between them lies in the increased planarity of the six-membered ring core observed in the closed isomers (see Fig. S7) and the C–N distances: D-CI_main is located at d(C–N) of 177 pm (for comparison these distances are 148 pm in D-closed; 139 pm in D-ZI; and 331 pm in D-open); and K-CI_main at 189 pm (K-closed: 148 pm; K-open: 338 pm; K-ZI: 139 pm). This agrees well with previous studies of conventional diarylethenes, reporting similar CIs at C–C distances between 187 and 208 pm^18,20,21 and studies on the ring-closing reaction of cis-stilbene reporting 189 and 194 pm.^61,62 Moreover, studies on CHD/cHT and derivatives have reported similar CIs with C–C distances between 212 and 228 pm.^17,22–25 The differences between our C–N and reported C–C distances can be associated with the contraction of the orbitals at the N- compared to the C-atom, leading to generally shorter C–N vs. C–C bonds. We examined D-CI_main and K-CI_main in greater detail by analyzing the topology of the CI and the corresponding branching plane vectors, g̃ and h̃, derived from the gradient difference (g) and nonadiabatic coupling (h) vectors.⁶³ Both intersections are classified as sloped, single-path CIs (see Fig. S8 and Table S5).⁶⁴

Fig. 4 Projection of the productive, D-open → D-closed → D-ZI (a) and K-open → K-closed (b) MRSF-NAMD and DFT-AIMD trajectories onto the two main reaction coordinates (C–N and C–S distance). The temporal evolution is indicated by a color gradient (initial: white, final: dark red/blue) and the initial conditions (bright, unfilled symbols) and trajectory end points (dark, filled symbols) are highlighted. The inset shows an overlay of the hopping point geometries, whereas D-CI_main and K-CI_main are shown as black structures for comparison.

For D-CI_main, motion along the −g̃ vector corresponds primarily to the shortening of the C–N bond, defining the principal reaction coordinate. Displacement along −h̃ is dominated by the contraction of the C–C bond linking the cyclopentene and pyrimidine moieties, a deformation that likely promotes ring closure by converting a single bond in the open form into a double bond in the closed structure. While the reaction outcome cannot be rigorously predicted from the hopping-point geometries, their orientation in the branching plane suggests likely structural tendencies: displacements in the −g̃/−h̃ direction involve distortions that are compatible with ring closure, whereas motion towards +g̃/+h̃ correspond to motions leading back to the open form.

For K-CI_main, both g̃ and −h̃ are associated with the shortening of the reaction determining C–N bond. Additionally, −h̃ has contributions of the stretching of the central C–C bond, a geometric feature linked to the formation of closed isomers via the double-to-single bond change in the central bond. Therefore, the branching plane vectors indicate that hops in the −g̃/+h̃ quadrant lead back to the open isomer, whereas hops in the +g̃/−h̃ quadrant yield closed isomers.

Other S₁/S₀ CIs accessible from the FC region, representing different potential loss channels, were optimized from snapshots of our NAMD simulations. In the following, they will be denoted as D-CI_thio, D-CI_thio*, D-CI_twist and D-CI_twist* for D (0.77, 0.76, 1.26 and 1.45 eV below the FC point, see Table S6) and K-CI_thio, K-CI_thio*, K-CI_thiazole and K-CI_thiazole* for K (0.42, 0.55, 1.44 and 1.39 eV below the FC point, see Table S7).

3.4 Photoinduced dynamics

The static PES analysis provides a qualitative framework for understanding the potential reaction pathways following excitation of D and K. To complete this picture and capture the temporal evolution of the excited-state processes, including competing pathways and relaxation channels, we now turn to the discussion of the surface hopping dynamics results.

After S₁ excitation, both D and K rapidly exchange population with S₂ before slower S₁ → S₀ decay governing the reaction outcome. For a more detailed analysis, trajectories were categorized according to their photochemical outcome, determined by monitoring the C–N bond length: productive trajectories leading to the closed isomer, unproductive ones returning to the open isomer, and a few forming rare side products. Of the 97 trajectories simulated for each system, 11 for D and 5 for K remained in an excited state after 2 ps and 1 ps, respectively, and were excluded from the analysis, as their outcomes could not be classified. This exclusion introduces a slight bias towards shorter lifetimes, i.e., lead to an underestimation of the reaction rates relative to the complete ensemble.

In the following, we first focus on the geometric changes in the photoinduced trajectories, highlighting the structural evolution associated with productive and unproductive pathways. We then discuss the excited-state population dynamics and associated reaction rates, before examining the outcomes on the S₀ using ab initio ground-state molecular dynamics simulations, which provide insight into the closed → ZI reaction.

Structural evolution of D and K-open upon S₁ excitation. The structural evolution along the excited-state trajectories reveals the geometric changes that distinguish productive pathways, leading to ring closure, from unproductive pathways, which return to the open form. In all productive trajectories, ring closure proceeds via a conrotatory pathway, consistent with the S₁ relaxed scan. This is reflected in the evolution of the dihedral angles describing the orientation of the thiophene and N-heterocycle rings relative to the bridging unit (ϕ_thio and ϕ_pyr; see Fig. S9). In the Franck–Condon region of D-open, both dihedrals are positive and decrease toward zero or slightly negative values as the system approaches the hopping points. This indicates rotation of both heterocycles in the same direction and supports a conrotatory mechanism. Although, in case of D, in principle either nitrogen atom of the pyrimidine ring could participate in the electrocyclic ring closure, all productive trajectories involve the nitrogen already closer to the thiophene ring, indicating that our initial conditions precondition the system for this specific ring closure.

While a broader distribution of ring orientations is expected under experimental solution conditions, the present study focuses on photochemical processes, i.e. trajectories involving chemical transformation rather than photophysical relaxation pathways, and therefore restricts the analysis to this productive orientation to obtain statistically meaningful sampling. Consequently, the computed lifetimes and yields reflect the dynamics of a selected subset of trajectories defined by the chosen initial ensemble and should not be interpreted as quantitative predictions of experimental values.

In the following, we discuss the respective yields, geometric pathways and hopping points in more detail, starting with D and subsequently turning to K and a comparison between both molecules.

D-open. For D, 51% of the trajectories that reach S₀ follow productive pathways towards D-closed, while the remaining 49% return to the open form or produce minor side products. In all productive trajectories, S₁ → S₀ hops occur close to the optimized D-CI_main, with an average d(C–N) of 179 pm at the hopping points (vs. 177 pm in the optimized geometry). After hopping to the ground state, d(C–N) shortens further, and the four-ring core flattens, leading to formation of D-closed. The structural distribution of these hops is illustrated as a geometry overlay in Fig. 4a and Table S6.

Based on the geometries of all hops in the vicinity of D-CI_main, labeled by reaction outcome (productive/unproductive), the differentiation between productive and unproductive pathways can be described by three key geometric features of the central six-membered ring. They were identified via supervised random forest classification using bond lengths, bond angles, and torsion angles of the respective ring as selected and computed with the help of shnitsel-tools⁶⁵ (cf. SI, Section S3). Among these, the reaction-determining C–N bond length is the most important, accompanied by two dihedrals describing the planarity of the newly formed six-membered ring in productive pathways: the dihedral around the pyrimidine C–N bond (φ₁) and the dihedral around the newly formed C–N bond (φ₂) in D-closed (see structural insets in Fig. S11).

Analysis of the distribution of hopping points along these coordinates reveals three distinct geometric groups (cf. Fig. S9 and S10). A subgroup of productive hops clusters at short C–N distances, with any hop at d(C–N) ≤ 170 pm being productive; these also share dihedral values within −20° < φ₁ < +30° and −20° < φ₂ < +50°, reflecting a highly planar orientation of the six pro-cycle atoms, which form the new six-membered ring. The remaining productive hops, at d(C–N) > 170 pm, and unproductive hops at similar C–N distances, display larger dihedral values, corresponding to less flattened ring geometries. These two groups cannot be further separated based on geometric descriptors, hopping times, or the energy gaps between S₀ and S₁ at the hop. Overall, 15 of 42 unproductive hops (36%) occur near D-CI_main, with an average d(C–N) of 190 pm, only slightly longer than in productive hops (see gray symbols in Fig. S10 and S11).

The remaining unproductive hops outside the vicinity of D-CI_main occur at geometries with longer C–N distances, with an average d(C–N) of 281 pm, and different structural features: One subgroup of these unproductive hops (40% of unproductive hops, 18% of total hops) shows a large twist of the two heterocycles at the bridging cyclopentene due to pyramidalization of both carbon atoms of the central double bond. These hops occur close to two other CIs, D-CI_twist and D-CI_twist* (see top right columns in Table S6). Furthermore, eight hops (9% of all hops) occur at structures bearing a ring-opened thiophene-unit, reflecting the optimized CIs D-CI_thio and D-CI_thio* (see bottom row in Table S6).

K-open. Of the 97 trajectories reaching S₀ in K, only 30% follow productive pathways to form K-closed, resulting in a lower photochemical reaction yield compared to D (51%). Again, all productive trajectories follow a conrotatory pathway similar to the trajectory of the relaxed scan on the S₁ state. As in D, this behavior is evident from the evolution of the dihedral angles describing orientation of the heterocycles (ϕ_thio and ϕ_thiazole). Both dihedral angles exhibit positive values at their initial points and decrease towards zero or small negative values at their final points (see Fig. S9).

Analysis of the hopping points shows that productive hops, like for D, occur close to the optimized K-CI_main (with average C–N distances of 199 pm vs. 189 pm in the optimized K-CI_main cf. Table S7). Most unproductive hops (59; 92% of unproductive and 64% of all trajectories) also cluster near K-CI_main, with an average C–N distance of 221 pm.

Analysis along the same three key coordinates used for D (C–N distance and the two dihedrals φ₁ and φ₂ describing ring planarity) shows that, unlike in D, productive and unproductive hops do not form distinct clusters (see Fig. S12 and S13). Pairwise correlations reveal that both types of hops near K-CI_main occupy largely overlapping regions, analogous to the mixed hopping region around D-CI_main, indicating that these internal coordinates do not discriminate reaction outcomes. As observed for D-CI_main, neither the energy gap at the hop (ΔE₀₁) nor the hopping time (t_hop) correlate clearly with the reaction outcome (Fig. S11).

The remaining five unproductive hops (8% of unproductive, 5% of all hops) show geometries which structurally differ more from K-CI_main and show either an open thiophene or an open thiazole. Those hops occur in proximity to CIs K-CI_thio, K-CI_thio* and K-CI_thiazole (see structures of these hopping points in the last two columns of Table S6). In contrast to D, no hops at geometries exhibiting a twisted cyclopentene moiety were observed, which indicates that the loss channel via a twisted CI observed in D is missing in K and results in more trajectories reaching the K-CI_main region.

Comparison. Comparing D and K highlights distinct branching behavior near their respective main CIs. At K-CI_main, all hops occur at C–N distances close to or larger than that of the optimized CI, with dihedrals associated with core planarity in the ranges of −60° < φ₁ < −20° and +40° < φ₂ < +70°, respectively. These values are comparable to the mixed cluster near D-CI_main (−60° < φ₁ < −20° and +50° < φ₂ < +80°; cf. Fig. S11 and S13), indicating similarly non-discriminating geometries. In contrast, D exhibits an additional cluster of 19 productive-only hops (32% of trajectories passing near D-CI_main) at shorter C–N distances and substantially smaller absolute dihedrals, reflecting a more planar ring-closing segment. Consequently, hops near K-CI_main arise from geometries that do not favor a particular outcome, whereas a significant fraction of hops near D-CI_main originate from structures already closer to the closed isomer, leading to consistently productive outcomes. Overall, 73% of trajectories hopping near D-CI_main are productive, compared to only 32% near K-CI_main. Despite more trajectories reaching K-CI_main than D-CI_main, the overall fraction of unproductive trajectories remains higher in K than in D (70% vs. 49%) due to the different branching near the CI. We note again that the quantum yields reported here, obtained from our NAMD simulations, are not necessarily directly comparable to experimental values. Additional nonproductive decay channels may exist that are not captured within the present dynamical framework. Furthermore, the neglect of explicit solvent effects in the surface-hopping simulations may lead to an increase of the computed quantum yields, as previously demonstrated in studies comparing light-driven double isomerizations in molecular rotors in methanol and in the gas phase.^66,67

Photochemical kinetics of D- and K-open. Upon excitation to the S₁ state, both D and K exhibit qualitatively similar early-time dynamics. Within the first few femtoseconds, significant population transfer from S₁ to S₂ occurs, yielding nearly equal populations of the two states. Subsequently, the S₂ population decays back to S₁, from which internal conversion to the ground state proceeds on a slower (sub-picosecond) timescale (see Fig. 5 and Fig. S15).

Fig. 5 Population plots from MRSF-NAMD simulations for all trajectories of D (a) and K (b). For comparison, the population plots of the productive, i.e., open → closed form trajectories, are shown (thin lines). To highlight their differences, the area between the curves of all trajectories and the productive ones is highlighted by colored areas.

To compare the dynamical behavior of productive and unproductive trajectories, we analyzed these subsets separately using a sequential kinetic model: S₂ → S^hot₁ → S^cold₁ → S₀. We divided S₁ in S^hot₁ and S^cold₁ to account for the observed bi-exponential kinetics in the S₀ and S₁ population curves upon back-population of S₁ from S₂. This behavior can be associated with the substantial geometric rearrangements within S₁ that are required to reach any S₁/S₀ CI. In the models S₁ lifetimes are obtained by the sum of the lifetimes of S^hot₁ and S^cold₁.

Comparing D and K, the overall dynamics in K are faster. This is evident from the S₁→ S₀ decay times: 500 fs (productive) and 648 fs (unproductive) for D, vs. 393 fs (productive) and 306 fs (unproductive) for K. In D, productive trajectories decay more rapidly than unproductive ones, consistent with the hopping point distribution, which shows unproductive hops extending beyond 1.3 ps (see right column in Fig. S10). In contrast, in K unproductive trajectories are slightly faster than productive ones, as also reflected in the temporal distribution of hopping points (see right column in Fig. S12). We denote here again that due to the exclusion of trajectories still in an excited state at their final points, the lifetimes given here might be slightly biased towards shorter lifetimes.

Fate of D- and K-closed on S₀. We investigated the fate of all productive isomers (44 trajectories for D and 28 trajectories for K) on the ground state with Born–Oppenheimer molecular dynamics (BOMD) simulations, initialized with geometries and velocities taken 10 fs after the final S₁→ S₀ surface hop. As the surface hopping simulations are performed in the gas phase, the molecules retain substantial excess kinetic energy following the nonadiabatic transition. Under experimental conditions, this excess energy would be dissipated into the surrounding medium through vibrational cooling, leading to an overall slower dynamics. Accordingly, even though the BOMD simulations include implicit solvent effects, the resulting reaction times are still likely to be faster than those observed experimentally. Thus, the BOMD results should be viewed primarily as a source of mechanistic understanding rather than as providing precise quantitative kinetic predictions.

Our results show that 80% of D-closed form D-ZI within 5 ps simulation time. A single trajectory (2%) returns to D-open within the same simulation time, whereas the remaining 18% stayed in D-closed. The picosecond timescale of this ground-state reaction obtained from our BOMD simulations despite a hypothetical micro- to nanosecond lifetime of D-closed obtained from static calculations of relative energies and barrier heights can be attributed to the elevated kinetic energy from irradiation. It enables D-closed to overcome the ground-state barrier more rapidly than would occur in a reaction starting from the ground state equilibrium geometry. On the other hand, neither K-closed → K-ZI nor K-closed → K-open back reaction was observed in the same simulation time. Instead, all trajectories remained in K-closed.

Following the evolution of the two main reaction coordinates in all productive trajectories reveals that shortening of the C–N bond to form the central six-membered ring and elongation of the C–S distance to open the thiophene ring in D are in fact two separate steps. It is only after reaching the closed isomer that the C–S bond starts to elongate, as illustrated by the orthogonal evolution of both reaction coordinates along the trajectories (see Fig. 4).

In conclusion, our dynamic results rationalize the experimental observation of different products following irradiation of D and K. Both compounds form their closed isomers on sub-picosecond timescales via funneling through a diarylethene-typical CI. However, D subsequently converts on the ground state to D-ZI, whereas K remains as K-closed. The lack of experimental detection of D-closed, despite its initial formation, can be attributed to its short lifetime, further decreased by excess kinetic energy from irradiation.

Conclusions

This study investigated the photoinduced ring-closure dynamics of the open aza-diarylethenes (aDAEs) D-open and K-open. Static and nonadiabatic dynamical calculations show that excitation to the bright first singlet excited state (S₁) directs both systems towards conical intersections closely resembling those known for diarylethenes. While several minor deactivation pathways are accessible, the overall quantum yield in our simulations is predominantly determined by the branching behavior in the vicinity of these intersections. Short C–N distances and near-planar core geometries at the hopping points favor productive ring closure, consistent with a 6π-azacyclic mechanism. From this region, the excited-state population either relaxes back to the open form or proceeds to the closed isomer on the ground state, with surface-hopping simulations indicating sub-picosecond ring closure for both compounds.

Despite their similar excited-state dynamics, the two systems diverge markedly after closure. In D, the closed isomer rapidly undergoes a thermal ground-state transformation to the zwitterionic product D-ZI, which according to our ground-state dynamics takes place within a few picoseconds, whereas in K the zwitterionic form is energetically inaccessible and K-closed remains the terminal product. These findings identify ground-state isomer stability as the key factor governing the distinct photochemical outcomes.

Overall, this work provides a detailed mechanistic picture of aDAE photoactivation and highlights how subtle structural and energetic differences control reaction fate.

Author contributions

M. H. performed, analyzed, and visualized all simulations reported in this work under the supervision of C. M. who conceived and oversaw the study. H. D. contributed as expert on aDAEs through initializing mechanistic investigations and discussions of their excited state behavior. All authors participated in the writing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study, including XYZ geometries of minima and minimum-energy conical intersections, absorption properties, relaxed scans and trajectories, are openly available in Zenodo, at https://doi.org/10.5281/zenodo.17601075, reference number 17601075.

The supplementary information (SI) contains benchmarks of electronic structure theory methods for absorption properties and scans along reaction coordinates, characterization of minimum energy conical intersections, analysis of hopping points, and kinetic modeling. See DOI: https://doi.org/10.1039/d6cp01330a.

Acknowledgements

We gratefully acknowledge financial support from the Fonds der Chemischen Industrie (Kekulé Fellowship, M. H.), the Competence Center Engineering of Advanced Materials of the Friedrich-Alexander-Universität Erlangen-Nürnberg (EAM Starting Grant, C. M.) and the Bayerische Akademie der Wissenschaften (C. M.). Furthermore, C.M. acknowledges funding by the Hector Fellow Academy. This work was funded by Deutsche Forschungsgemeinschaft (DFG) within the Collaborative Research Center ‘ChemPrint’ (project 538767711, CRC 1719). All high-performance computational work was performed at the Erlangen National High Performance Computing Center (NHR@FAU).

Notes and references

1. M. Irie and M. Mohri, J. Org. Chem., 1988, 53, 803–808.
2. M. Irie, Chem. Rev., 2000, 100, 1685–1716.
3. K. Matsuda and M. Irie, J. Photochem. Photobiol., C, 2004, 5, 169–182.
4. M. Irie, T. Fukaminato, K. Matsuda and S. Kobatake, Chem. Rev., 2014, 114, 12174–12277.
5. H.-B. Cheng, S. Zhang, E. Bai, X. Cao, J. Wang, J. Qi, J. Liu, J. Zhao, L. Zhang and J. Yoon, Adv. Mater., 2022, 34, 2108289.
6. K. Uchida, M. Saito, A. Murakami, S. Nakamura and M. Irie, Adv. Mater., 2003, 15, 121–125.
7. T. Fukaminato, T. Doi, N. Tamaoki, K. Okuno, Y. Ishibashi, H. Miyasaka and M. Irie, J. Am. Chem. Soc., 2011, 133, 4984–4990.
8. B. M. Neilson and C. W. Bielawski, J. Am. Chem. Soc., 2012, 134, 12693–12699.
9. M. Takeshita and M. Irie, Tetrahedron Lett., 1998, 39, 613–616.
10. S. Kobatake, H. Hasegawa and K. Miyamura, Cryst. Growth Des., 2011, 11, 1223–1229.
11. P. R. Hania, A. Pugzlys, L. N. Lucas, J. J. D. de Jong, B. L. Feringa, J. H. van Esch, H. T. Jonkman and K. Duppen, J. Phys. Chem. A, 2005, 109, 9437–9442.
12. H. Jean-Ruel, R. R. Cooney, M. Gao, C. Lu, M. A. Kochman, C. A. Morrison and R. J. D. Miller, J. Phys. Chem. A, 2011, 115, 13158–13168.
13. S. Shim, I. Eom, T. Joo, E. Kim and K. S. Kim, J. Phys. Chem. A, 2007, 111, 8910–8917.
14. H. Miyasaka, S. Araki, A. Tabata, T. Nobuto, N. Malaga and M. Irie, Chem. Phys. Lett., 1994, 230, 249–254.
15. H. Miyasaka, T. Nobuto, A. Itaya, N. Tamai and M. Irie, Chem. Phys. Lett., 1997, 269, 281–285.
16. N. Ohtaka, Y. Hase, K. Uchida, M. Irie and N. Tamai, Mol. Cryst. Liq. Cryst., 2000, 344, 83–88.
17. E. Salazar and S. Faraji, Mol. Phys., 2020, 118, e1764120.
18. E. Salazar, S. Reinink and S. Faraji, Phys. Chem. Chem. Phys., 2022, 24, 11592–11602.
19. M. Martyka and J. Jankowska, J. Photochem. Photobiol., A, 2023, 438, 114513.
20. M. Martyka and J. Jankowska, Phys. Chem. Chem. Phys., 2024, 26, 13383–13394.
21. J. Jankowska, M. Martyka and M. Michalski, J. Chem. Phys., 2021, 154, 204305.
22. Y. Liu, R. Xu, D. M. Sanchez, T. J. Martìnez and T. J. A. Wolf, Annu. Rev. Phys. Chem., 2025, 76, 615–638.
23. E. G. Champenois, D. M. Sanchez, J. Yang, J. P. Figueira Nunes, A. Attar, M. Centurion, R. Forbes, M. Guhr, K. Hegazy, F. Ji, S. K. Saha, Y. Liu, M.-F. Lin, D. Luo, B. Moore, X. Shen, M. R. Ware, X. J. Wang, T. J. Martìnez and T. J. A. Wolf, Science, 2021, 374, 178–182.
24. Y. Liu, D. M. Sanchez, M. R. Ware, E. G. Champenois, J. Yang, J. P. F. Nunes, A. Attar, M. Centurion, J. P. Cryan, R. Forbes, K. Hegazy, M. C. Hoffmann, F. Ji, M.-F. Lin, D. Luo, S. K. Saha, X. Shen, X. J. Wang, T. J. Martìnez and T. J. A. Wolf, Nat. Commun., 2023, 14, 2795.
25. T. J. A. Wolf, D. M. Sanchez, J. Yang, R. M. Parrish, J. P. F. Nunes, M. Centurion, R. Coffee, J. P. Cryan, M. Guhr, K. Hegazy, A. Kirrander, R. K. Li, J. Ruddock, X. Shen, T. Vecchione, S. P. Weathersby, P. M. Weber, K. Wilkin, H. Yong, Q. Zheng, X. J. Wang, M. P. Minitti and T. J. Martìnez, Nat. Chem., 2019, 11, 504–509.
26. S. Hamatani, D. Kitagawa and S. Kobatake, J. Phys. Chem. Lett., 2023, 14, 8277–8280.
27. M. Sacherer, S. Gracheva, H. Maid, C. Placht, F. Hampel and H. Dube, J. Am. Chem. Soc., 2024, 146, 9575–9582.
28. S. Hamatani, D. Kitagawa and S. Kobatake, Angew. Chem., Int. Ed., 2024, 63, e202414121.
29. D. Kitagawa, T. Mori, H. Sawa, S. Hamatani and S. Kobatake, Chem. Eur. J., 2025, 31, e202501077.
30. M. Sacherer and H. Dube, Angew. Chem., Int. Ed., 2025, 64, e202415961.
31. M. Sacherer and H. Dube, Adv. Mater., 2025, 37, 2411223.
32. A. V. Stepanov, A. A. Shatrova, D. V. Berdnikova and A. G. Lvov, Org. Biomol. Chem., 2025, 23, 5908–5928.
33. K. Uchida and M. Irie, Chem. Lett., 2006, 24, 969–970.
34. F. Neese, WIREs Comput. Mol. Sci., 2012, 2, 73–78.
35. F. Neese, WIREs Comput. Mol. Sci., 2022, 12, e1606.
36. W. Park, K. Komarov, S. Lee and C. H. Choi, J. Phys. Chem. Lett., 2023, 14, 8896–8908.
37. S. Lee, M. Filatov, S. Lee and C. H. Choi, J. Chem. Phys., 2018, 149, 104101.
38. S. Lee, E. E. Kim, H. Nakata, S. Lee and C. H. Choi, J. Chem. Phys., 2019, 150, 184111.
39. V. Mironov, K. Komarov, J. Li, I. Gerasimov, H. Nakata, M. Mazaherifar, K. Ishimura, W. Park, A. Lashkaripour, M. Oh, M. Huix-Rotllant, S. Lee and C. H. Choi, J. Chem. Theory Comput., 2024, 20, 9464–9477.
40. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begusic, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov, S. Fürst, S. Gillhuber, R. Grotjahn, M. Kaupp, M. Kehry, M. Krstic, F. Mack, S. Majumdar, B. D. Nguyen, S. M. Parker, F. Pauly, A. Pausch, E. Perlt, G. S. Phun, A. Rajabi, D. Rappoport, B. Samal, T. Schrader, M. Sharma, E. Tapavicza, R. S. Treβ, V. Voora, A. Wodynski, J. M. Yu, B. Zerulla, F. Furche, C. Hattig, M. Sierka, D. P. Tew and F. Weigend, J. Chem. Theory Comput., 2023, 19, 6859–6890.
41. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
42. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
43. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001.
44. K. Komarov, W. Park, S. Lee, M. Huix-Rotllant and C. H. Choi, J. Chem. Theory Comput., 2023, 19, 7671–7684.
45. B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson and T. L. Windus, J. Chem. Inf. Model., 2019, 59, 4814–4820.
46. L. Yu, C. Xu, Y. Lei, C. Zhu and Z. Wen, Phys. Chem. Chem. Phys., 2014, 16, 25883–25895.
47. J. Li, P. Reiser, B. R. Boswell, A. Eberhard, N. Z. Burns, P. Friederich and S. A. Lopez, Chem. Sci., 2021, 12, 5302–5314.
48. J. Li, R. Stein, D. M. Adrion and S. A. Lopez, J. Am. Chem. Soc., 2021, 143, 20166–20175.
49. J. Li and S. A. Lopez, Chem. Eur. J., 2022, 28, e202200651.
50. J. Li and S. A. Lopez, Acc. Chem. Res., 2022, 55, 1972–1984.
51. L. Wang, C. Salguero, S. A. Lopez and J. Li, Chem, 2024, 10, 2295–2310.
52. Z. Li, F. J. Hernandez, C. Salguero, S. A. Lopez, R. Crespo-Otero and J. Li, Nat. Commun., 2025, 16, 1194.
53. M. Brady and R. Crespo-Otero, J. Chem. Phys., 2025, 163, 234118.
54. E. G. F. de Miranda, R. Souza Mattos, S. Mukherjee, J. M. Toldo, C. H. Choi, M. T. D. N. Varella and M. Barbatti, J. Chem. Theory Comput., 2026, 22, 1–19.
55. S. Mai, D. Avagliano, M. Heindl, P. Marquetand, M. F. S. J. Menger, M. Oppel, F. Plasser, S. Polonius, M. Ruckenbauer, Y. Shu, D. G. Truhlar, L. Zhang, P. Zobel and L. Gonzalez, SHARC3.0: Surface Hopping Including Arbitrary Couplings Program Package for Non-Adiabatic Dynamics, 2023, https://sharc-md.org/.
56. S. Mai, B. Bachmair, L. Gagliardi, H. G. Gallmetzer, L. Grünewald, M. R. Hennefarth, N. M. Høyer, F. A. Korsaye, S. Mausenberger, M. Oppel, T. Piteša, S. Polonius, E. S. Gil, Y. Shu, N. K. Singer, M. X. Tiefenbacher, D. G. Truhlar, D. Vörös, L. Zhang and L. González, SHARC4.0: Surface Hopping Including Arbitrary Couplings – Program Package for Non-Adiabatic Dynamics, https://sharc-md.org/, 2025.
57. B. Mignolet, B. F. E. Curchod and T. J. Martìnez, Angew. Chem., Int. Ed., 2016, 55, 14993–14996.
58. S. Mai, P. Marquetand and L. González, WIREs Comput. Mol. Sci., 2018, 8, e1370.
59. C. Climent, Z. Xu, M. O. Wolf and D. Casanova, J. Phys. Chem. Lett., 2024, 15, 8042–8048.
60. D.-H. Xu, L. Li, X.-Y. Liu and G. Cui, Molecules, 2021, 26(9), 2724.
61. S. Karashima, X. Miao, A. Kanayama, Y.-I. Yamamoto, J. Nishitani, N. Kavka, R. Mitric and T. Suzuki, J. Am. Chem. Soc., 2023, 145, 3283–3288.
62. M. Farmani, W. Park, S. Lee and C. H. Choi, Photochem. Photobiol. Sci., 2025, 24, 955–962.
63. B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem., 2007, 58, 613–634.
64. I. Fdez Galvan, M. G. Delcey, T. B. Pedersen, F. Aquilante and R. Lindh, J. Chem. Theory Comput., 2016, 12, 3636–3653.
65. K. Höllring, T. E. Röhrkasten and C. Müller, Digit. Discov., 2026 10.1039/D5DD00299K.
66. A. Blanco-Gonzalez, M. Manathunga, X. Yang and M. Olivucci, Nat. Commun., 2024, 15, 3499.
67. M. Filatov, M. Paolino, R. Pierron, A. Cappelli, G. Giorgi, J. Leoonard, M. Huix-Rotllant, N. Ferre, X. Yang, D. Kaliakin, A. Blanco-Gonzalez and M. Olivucci, Nat. Commun., 2022, 13, 6433.

---