First principles prediction of wavelength-dependent isomerization quantum yields of a second-generation molecular nanomotor

Abstract

Second-generation molecular nanomotors are gaining popularity within the biomedical field and intense research is being conducted to increase their efficiency for light-induced ultrafast photoisomerization. A key requirement for designing efficient molecular nanomotors is ensuring unidirectional rotation during isomerization and thermal helix inversion. Here, we used non-adiabatic trajectory surface hopping molecular dynamics based on TDDFT to study the excited state dynamics of the stable M- and metastable P-conformers of a second-generation Feringa-type molecular nanomotor. From the trajectories, we computed quantum yields for clockwise and anti-clockwise photoisomerization. Results show that the helicity of the initial structure dictates the direction of the isomerization: 52% of the trajectories starting from M-conformers isomerize clockwise, whereas 23% of the trajectories starting from P-conformers isomerize anti-clockwise. The quantum yield for clockwise isomerization can be maximized by excitation in the center of the absorption spectrum (350–400 nm). In this region, the M-conformer exhibits its maximum absorption and maximum (clockwise) isomerization quantum yield, whereas the P-conformer shows negligible excitation probability and (anti-clockwise) isomerization quantum yield. Moreover, we also find several trajectories, where thermal helix inversion occurs in the hot ground state, directly after excited state relaxation. While helix inversion was observed to proceed in both directions, its net effect favors unidirectional rotation. We further report excited state lifetimes and details about the structural dynamics.

---

1. Introduction

Light-driven molecular nanomotors (MNMs) hold significant potential for biomedical applications, including drug delivery, membrane permeation, and ion transport.^1,2 Feringa-type overcrowded alkene motors consist of two polycyclic alkene-based subunits, usually referred to as the stator and rotator, connected by a carbon–carbon double bond (Fig. 1). Ideally, light irradiation induces an E–Z photoisomerization³ (hν in Fig. 1) followed by a thermal helix inversion (Δ in Fig. 1).^4,5 When these two elemental steps are repeated and occur in the same direction, one subunit completes a full 360° rotation relative to the other. Substantial research efforts, both experimental^4,6–11 and computational,^5,12–17 have been focused on understanding the dynamics of MNMs and on optimizing MNMs to increase their efficiency. With regard to the latter, studies focused on: (i) increasing photoisomerization quantum yield,^6,12 (ii) increasing thermal helix inversion rates,^8,10,16 (iii) increasing rotational speed,^9–11 determined by the time taken for isomerization and helix inversion, and (iv) ensuring unidirectional rotation.^11,17,18 Computational studies focused on several aspects of the molecular mechanism and tuning factors of the process using a wide range of static and dynamic calculations,^15,19 at different levels of theory such as time-dependent density functional theory (TDDFT),¹³ its tight-binding version (TDDFTB),¹⁹ TDDFT restricted ensemble Kohn–Sham (REKS),^12,13 multiconfigurational methods,¹⁷ or ADC(2).¹⁵ This theoretical effort has led to increased understanding of the mechanism of MNMs and which factors influence its functioning. For instance, it has been shown that an increasing polarity of the solvent increases the excited state lifetime of the MNM, while increased viscosity produces the opposite effect.⁵ González and co-workers performed computational excited-state dynamics of a molecular rotor embedded in a solvent (DMSO), gauging the influence of the solvent on the MNM photoisomerization.¹⁵ Furthermore, it has been found that hydrogen bonds between MNM and solvent decrease the rotation efficiency.¹⁴ Also substituents play an important role in the efficiency of the photoinduced rotation: the gap between the heterolytic and homolytic partial bond breaking of the central carbon–carbon double bond should be minimal to lead to a faster photoisomerisation and enhanced quantum yield.¹³ Moreover, a stator that decreases steric repulsion in the fjord region also improves the photoisomerization rate. On the other hand, increasing the size of substituents at the stereogenic center accelerates the rotary motion because it reduces the energetic barrier for thermal helix inversion.¹¹ Despite the gained understanding of the working mechanism, there are still various aspects that require improvement of molecular properties. For instance, for application in biological tissue, activation by near-infrared instead of UV/vis light would be beneficial, as photons in this region exhibit deeper penetration depth in biological tissues.²⁰ Using two-photon excitation, near-infrared activation has recently been achieved,²¹ but it still requires high photon intensity, due to the low two-photon absorption cross sections. Another property to be increased is high photostability, to ensure the MNM can sustain as many cycles as possible without degradation.²²

Fig. 1 Graphical representation of the proposed stages of the unidirectional motion of CPNX. Helicity (M or P) is indicated. hν indicate photochemical steps, whereas Δ indicate thermal reactions. The stereogenic center is marked with *.

Feringa and co-workers designed a pioneering MNM that served as a blueprint for the development of first-generation MNMs, followed by second and third-generation versions.⁴ Unidirectional rotation is essentially a consequence of the chiral structure of the MNM.¹¹

While first-generation MNMs feature two stereogenic centers, second-generation MNMs contain only one, and third-generation MNMs are formed by merging two enantiomers from second-generation MNMs to create a pseudo-asymmetric center responsible for unidirectional rotation. The reduction in stereogenic centers has facilitated the synthesis of second- and third-generation MNMs, simplifying the process compared to their predecessors.²³

During 360° rotation, the MNM adopts several conformers (Fig. 1). Broadly, these can be classified into stable and metastable conformers, with opposite rotor helicity. While the absolute stereochemical configuration of the MNM is arbitrary, we use here the absolute configuration marked with * in Fig. 1, in which M-conformers (left-handed helicity) are stable, whereas P-conformers (right-handed helicity) are metastable. To ensure unidirectional rotation, both photoinduced isomerizations as well as thermal helix inversion should ideally proceed in the same direction; rotation in the opposite direction should not occur. While at thermal equilibrium, typically the stable conformer (the M-helical conformer in our example, Fig. 1) dominates the equilibrium of structures, this is not the case in non-equilibrium situations, such as the hot ground state after photorelaxation, as the primarily produced conformer is typically the metastable conformer. In this case, both left- and right-handed helicity conformers may coexist and both conformers are subject to photoexcitation in the next excitation cycle. The coexistence of several conformers could compromise unidirectional rotation, if both conformers, stable and metastable, undergo photoisomerization in opposite directions, as previously hypothesized.^12,16 Most computational studies that focussed on the photochemical isomerization have examined the primary isomerization of the stable conformer to the metastable conformer.^12–15,19 On the contrary, few research has been conducted to get insights into the excited-state relaxation of the metastable conformer (in our case P-conformer).^11,16 To address the different photochemical reactivity of the stable and metastable conformers, we investigate the excited-state decay pathways of both M- and P-conformers upon photoexcitation into the first excited singlet state (S₁). Specifically, we aim at verifying the hypothesis that conformers with opposite helicity undergo photoisomerization in different rotational directions.^12,16 Furthermore, we aim at determining the conditions at which the unidirectional rotation can be maximized. Using fewest switches surface hopping²⁴ based on TDDFT (TDDFT-SH), we study the dynamics of photoisomerization mechanism with a focus on predicting the quantum yield of clockwise and anti-clockwise photoisomerization, and its dependency on the excitation wavelength of the recently synthesized 9-(2′-methyl-2′,3′-dihydro-1′H-cyclopenta[a]naphthalen-1′-ylidene)-9H-xanthene (CPNX),^25–27 depicted in Fig. 1.

2. Computational details

2.1. Geometry optimization

Initially, we searched the conformational space of structures of CPNX using the software CREST,²⁸ employing metadynamics²⁹ simulations with extended tight-binding method (XTB)³⁰ in the gas phase. Structures within an energy window of 10 kcal mol⁻¹ with respect to the most stable structure found were further optimized using density functional theory (DFT)³¹ at the PBE0/def2-SVP level of theory^32–37 to locate the local minima of each species. Frequency calculations^38,39 were performed at the same level of theory to confirm that these structures corresponded to minima on the potential energy surface. The relative population of each conformer was assessed based on a Boltzmann distribution at 300 K. To verify the accuracy of the DFT results, the total energy using second-order approximate coupled cluster (CC2)/def2-TZVP was computed.^34–37,40 All DFT and CC2 calculations were performed with TURBOMOLE^41–43 and employed default convergence criteria as provided by the program.

2.2. Generation of ground state ensemble of structures

2.2.1. Replica exchange molecular dynamics. Five replica exchange molecular dynamics (REMD)⁴⁴ simulations using XTB were conducted in order to sample the conformational space of CPNX.³⁰ For each conformer, eight independent trajectories were generated at different temperatures (300, 350, 400, 450, 500, 550, 600, and 650 K). A time-step of 40 a.u. was used to propagate the nuclear degrees of freedom. Every 200 MD steps, the probability of switching between two trajectories at neighboring temperatures was evaluated. A comparison of this probability with a random number determined whether the current structures of the two trajectories would be exchanged. Further implementation details for REMD can be found elsewhere.^30,45,46 The total simulation time for each trajectory amounted to 9.5 ns, leading to a combined total of 47.5 ns, which yielded a ground-state ensemble of structures consisting of approximately 50 000 structures. REMD was carried out with a home-written script using TURBOMOLE.^41–43

2.2.2. Clustering of the structures. The structures from the REMD trajectories were clustered using the density-based spatial clustering of applications with noise (DBSCAN) algorithm, based on the structural differences between structures, as implemented in AmberTools 22.^47–49 For this algorithm, the user defines the minimum number of points that constitute a cluster and the neighborhood density (ε) around a point to determine whether it can form part of a cluster: it defines the maximum distance between two points to be considered neighbors. After selecting these parameters appropriately, the algorithm optimizes the number of clusters. In this study, ε was set to 0.5 Å, and the minimum number of points required to form a cluster was set to 1. Subsequently, the 372 starting structures for TDDFT trajectory surface hopping (TDDFT-TSH) (see below) were clustered using the same methodology to verify whether this reduced ensemble was representative of the complete REMD ensemble.

2.3. Computation of absorption spectra

The lowest 5 electronic excitation energies and oscillator strengths were computed for the ground state equilibrium structures and 372 selected structures. This was done using TDDFT in combination with the PBE0 exchange–correlation functional (TDPBE0) and the def2-SVP basis, and using CC2/def2-TZVP^34–37,40 with the resolution of identity approximation, as implemented in TURBOMOLE.^42,43,50 Absorption spectra were generated as a stick spectrum of the five lowest excitation energies and oscillator strengths broadened by Gaussian functions with a fixed full width at half maximum of 0.1 eV and converted to molar decadic extinction coefficients, as previously described.^51–54

2.4. Surface hopping dynamics

375 structures, selected at equal distances from the REMD ensemble, were used as starting structures for TDDFT-SH dynamics.^55–57 Three trajectories showed numerical errors along the excited-state dynamics and were discarded from the study. As a result, 372 trajectories are considered in this work. TDDFT-SH dynamics simulations were initiated in the first excited singlet state (S₁) and were propagated for 3 ps employing the hybrid PBE0 exchange–correlation functional^32,33 and the def2-SVP basis set,^34–37 as implemented⁵⁸ in TURBOMOLE V7.2.^41,42 TDDFT was used within the Tamm–Dancoff approximation (TDA) to minimize instabilities close to conical intersections.^59–61 Due to the life-time shorter than 20 ps, it is unlikely that triplet states play a major role in the dynamics of CPNX, as already indicated by Augulis and co-workers.²⁵ Therefore triplet states were not considered in our study. A total of two excited states were computed. A time-step of 40 a.u. was used to propagate the nuclear degrees of freedom. TDDFT-SH simulations employed analytic non-adiabatic coupling (NAC) vectors,^62–64 which allow proper velocity rescaling after surface hops to conserve the total energy, as proposed by Tully.²⁴ The hopping probability between excited states was evaluated at each time step, but in case of negative excitation energies, which occur at conical intersections within TDA,^60,61 a transition to the ground state was forced.

2.5. Wavelength-dependent product quantum yields

Wavelength-dependent product quantum yields were calculated as previously described.^65,66 To this end, we clustered the trajectories into successful and unsuccessful trajectories, depending on whether they underwent complete Z–E isomerization around their central double bond (C9=C1′, defined in Fig. 3) or whether they decayed to the ground state unreactively. We further distinguished between clockwise and anti-clockwise isomerization. The wavelength-dependent product quantum yield (Φ_j(λ)) for the isomerization of conformer j was then obtained by dividing the absorption spectrum averaged over the initial structures of the corresponding successful trajectories by the average absorption spectrum of the initial structures of all trajectories, according towhere σ_i(λ) indicates the absorption spectrum of the initial structure with a conformation j of trajectory i, and N_succ,j and N_total,j indicate the number of successful trajectories and total number of trajectories, respectively. Spectra were calculated using CC2 as described above.

3. Results & discussion

3.1. Ground-state equilibrium structures

The XTB-level conformational search with CREST yields seven distinct conformers. Each conformer was further optimized using PBE0/def2-SVP (Cartesian coordinates of the structures, labeled I1–I7 are given in Section S1 of the ESI†). Based on the similarity of the dihedral angles (Table S1 in the ESI†) and total energies (Table S2 in the ESI†), we further reduced the number of conformers to four: anti-M, syn-M, anti-P, and syn-P (Fig. 2 and Table 1). The notation for these conformers is as follows: M-helical conformers are characterized by a left-handed helicity, while P-helical conformers exhibit a right-handed helicity. The anti-folded or syn-folded character is determined by the orientation of the stator and rotator relative to the olefinic plane (containing the C9=C1′ double bond). When both subunits are oriented in the same direction relative to the olefinic plane, the conformation is syn-folded; otherwise, it is anti-folded. I1 and I2 correspond to the same anti-M structure whereas conformers I5 and I6 are associated with the syn-M isomer. In the case of isomers I3 and I4, both represent a typical syn-P structure, although the methyl group points towards a slightly different direction and thus, their geometrical parameters and relative energies differ. The most stable structure for each conformer was used as the representative for anti-M (I1), syn-M (I5), anti-P (I7), and syn-P (I3).

Fig. 2 Graphical representation of the different equilibrium conformers of CPNX optimized using PBE0/def2-SVP. Color code: gray for carbon, white for hydrogen and red for oxygen.

Table 1 Statistical distribution of the CPNX conformers in terms of a Boltzmann distribution using PBE0/def2-SVP for the equilibrium structures, and from the ensemble (and representative sub-ensemble) of structures obtained from REMD simulations. Boltzmann populations arise from the addition of the individual populations of those structures that represent the same conformer (Table S2 in the ESI)

Conformer	Boltzmann population equilibrium structures	Population REMD (all structures)	Population REMD sub-ensemble (372 structures)
anti-M	0.9951	0.5460	0.5780
syn-M	3.630 × 10^−3	1.987 × 10^−3	1.882 × 10^−2
anti-P	1.940 × 10^−7	0.1770	0.1425
syn-P	1.310 × 10^−3	0.2750	0.2608

---

Populations obtained from total energies computed by CC2/def2-TZVP (Table S2 in the ESI†) and PBE0/def2-SVP (Table 1) are qualitatively similar. The anti-M structure is the most stable conformer, with a population of 99.51% at 300 K (Table 1). The syn-M and syn-P conformers have similar populations, 0.36% and 0.13%, respectively. Finally, anti-P shows negligible population, ∼0.00002%. These data are consistent with the relative stability for the conformers of the analogous sulfur-containing MNM studied by Durbeej et al.⁶⁷

3.2. Analysis and clustering of REMD trajectories

All structures from the five REMD simulations were clustered using the DBSCAN algorithm.⁴⁹ Four different clusters were identified and labeled in terms of the cluster centroid. These centroids coincide with the four conformers found within the conformational search (Table 1).

We notice differences in the populations of the REMD ensemble and the Boltzmann probabilities of the optimized structures (Table 1). The large population differences indicate that REMD did not converge because the number of transitions between structures of REMD trajectories at different temperatures was not enough to ensure proper capturing of the statistics of the different conformers, as conformers were trapped for long simulation periods in one local minimum.

From the REMD ensemble, 372 samples were taken as initial structures for the TDDFT-SH simulations (REMD sub-ensemble). The sub-ensemble has a similar conformer distribution as the full REMD ensemble (Table 1).

3.2.1. Geometric parameters. We analyzed the conformers in terms of five dihedral angles, defined in Fig. 3a. The dihedral angle distributions (Fig. 3b–f) are consistent with the values of the equilibrium structures (Table S1 in the ESI†) and can be used to identify and distinguish the different conformers. The distributions reveal that the helicity is associated with the sign of the dihedral angle β (Fig. 3c): negative β values indicate M-helicity, whereas positive β values correspond to P-helicity.

Fig. 3 (a) Lewis structure of CPNX and definition of characteristic dihedral angles. Distributions of dihedral angles of the structures of the REMD simulations color coded for each cluster are shown in panels (b)–(f), as labeled.

The conformers anti-M and syn-M can be distinguished by their distinct δ- and δ′-distributions, Fig. 3e and f, respectively. However, small overlaps between the two distributions exist; thus for a few cases, the angle distributions are not sufficient to distinguish the two conformers and more geometric parameters have to be considered.

The conformers syn-P and anti-P can also be distinguished based on their distinct δ- and δ′-distributions, but as in the case of the M-conformers, small overlapping regions between the distributions of the two P-conformers exist.

3.3. Absorption spectra and excitation energies

Analysis of the CC2 results of the lowest five electronic excited states of the equilibrium structures (Table 2) shows that, for all the conformers, S₁ is bright with oscillator strengths at least one magnitude higher than those of S₂–S₅. For all conformers, S₁ has a HOMO–LUMO character, corresponding to a π–π* transition (orbitals and transition densities shown in Fig. S1–S4 in the ESI†). S₁ excitation energies range between 2.74–3.75 eV, mainly located in the UV region of the electromagnetic spectrum, with syn-P being the only one that has its major band in the visible region. TDPBE0/def2-SVP excitation energies, oscillator strengths, and orbital characters (Table S3, ESI†) agree well with CC2 values. Compared to TDPBE0, CC2 excitation energies are slightly higher, with a maximum deviation of 0.319 eV for the first excitation of the syn-M conformer; all other excitations show deviations well below 0.3 eV.

Table 2 CC2 excitation energies in eV, oscillator strengths in a.u. (in parentheses), and dominant orbital contribution for the five lowest excited states of the ground state equilibrium structures of all four conformers optimized with PBE0/def2-SVP. H and L denote highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), respectively

	S1	S2	S3	S4	S5
anti-M	3.428	4.056	4.187	4.344	4.562
	(0.4463)	(0.002901)	(0.04247)	(0.02397)	(0.005384)
	H → L (95%)	H → L+1 (52%)	H−1 → L (33%)	H → L+1 (36%)	H−3 → L (38%)
		H−1 → L (34%)	H → L+1 (26%)	H−1 → L (28%)	H → L+1 (19%)
anti-P	3.755	4.140	4.299	4.515	4.717
	(0.3662)	(0.01518)	(0.04128)	(0.01377)	(0.03392)
	H → L (85%)	H → L+1 (31%)	H−1 → L (36%)	H → L+1 (44%)	H−3 → L (34%)
		H → L+2 (28%)	H → L+1 (19%)	H−1 → L (23%)	H → L+1 (20%)
syn-M	3.381	4.011	4.053	4.222	4.346
	(0.4517)	(0.005296)	(0.009919)	(0.07164)	(0.0009657)
	H → L (91%)	H → L+1 (30%)	H → L+1 (53%)	H−1 → L (54%)	H → L+3 (78%)
		H−3 → L (17%)	H → L+2 (15%)	H → L+1 (17%)	H-1 → L+3 (7%)
syn-P	2.745	3.778	3.854	3.919	3.973
	(0.4410)	(0.004166)	(0.008005)	(0.03583)	(0.002100)
	H → L (93%)	H → L+1 (63%)	H → L+2 (29%)	H → L+3 (33%)	H → L+3 (50%)
		H−3 → L (13%)	H−3 → L (19%)	H → L+2 (26%)	H−3 → L (12%)

---

In addition, the theoretical absorption spectrum of each conformer σ_conf(λ) was calculated as the average of the spectra of individual structures of the REMD sub-ensemble (Fig. 4) according to

where N_conf is the number of structures for a given conformer and σ_i(λ) denotes the absorption spectrum of one individual structure. The total absorption spectrum at thermal equilibrium was computed by weighting each conformer's spectrum according to its Boltzmann probability at 300 K. Since anti-M accounts for 99.51% of the structures (Table 1), the total spectrum is practically identical to that of anti-M.

Fig. 4 Absorption spectra for each conformer (solid lines, colors as indicated) computed from the REMD sub-ensemble (Table 3). For the total spectrum at thermal equilibrium (black), each conformer spectrum was weighted by the corresponding Boltzmann population (Table 1). Arrows indicate the experimental values of the wavelength at the maximum of the absorption band of the stable^26,27 and the meta-stable species.²⁵ Stick spectra show excitation energies and oscillator strengths (length representation) of the equilibrium structures (Fig. 2 and Table 2), colors as indicated.

The S₁ excitation energies of the equilibrium structures (stick spectrum in Fig. 4) align closely with the peak positions of the maximum absorption bands of each conformer. The spectra of anti-M and syn-M are very similar to each other. In contrast, the two P-conformers exhibit very different absorption spectra: anti-P has its main absorption band blue-shifted relative to the two M-conformer bands, whereas the main absorption band of the syn-P conformer appears on the red side of the spectrum. We find good agreement between our theoretical values and experimental lambda max values^25–27 previously reported for stable and meta-stable conformers, as indicated in Fig. 4. We conclude that in the window between 360–400 nm, M-conformers can be predominantly excited.

Besides the ground state equilibrium structures, we also generally find good agreement between CC2 and (TD)PBE0 potential energies along the torsional coordinate α, as further discussed in Section S2 of the ESI† and indicated by Fig. S5 in the ESI.†

3.4. TDDFT-SH dynamics

Within the first 3 ps of the TDDFT-SH simulations, 97.8% of initial M-conformers and 86.7% of initial P-conformers decayed to S₀. The results of the simulations of the 372 selected structures (Table 3) showed that out of the 222 trajectories starting from M-helical initial structures, 117 (52%) undergo clockwise photoisomerization (Fig. 5a and c), while 105 trajectories do not show any photochemical transformation within 3 ps (Fig. 5b and d). From the 150 trajectories starting from P-helical structures, 34 (23%) isomerize anti-clockwise (Fig. 6a and c), while 116 trajectories do not show any E-Z isomerization within 3 ps (Fig. 6b and d). None of the initial structures with M-helicity underwent anti-clockwise isomerization, and none of the initial structures with P-helicity underwent clockwise isomerization.

Table 3 Conformer distribution of TDDFT-SH the initial structures and number of trajectories undergoing clockwise and anti-clockwise photoisomerization. Unsuccessful trajectories did not undergo Z–E isomerization within 3 ps simulation time

Conformer	Initial structures	Clockwise	Anti-clockwise	Unsuccessful
anti-(M)	215	114	0	101
syn-(M)	7	3	0	4
anti-(P)	53	0	30	23
syn-(P)	97	0	4	93

---

Fig. 5 Evolution of dihedral angles α (a, b) and β (c, d) along the trajectories with initial M-conformers. Panels on the left (a, c) represent trajectories that isomerize, whereas panels on the right (b, d) account for trajectories without isomerization. Green and red represent the electronic state (S₁ or S₀, respectively) of the molecule along the trajectory. The black lines (left panels) belong to an example trajectory that undergoes thermal helix inversion in the hot ground state after E–Z isomerization; associated structures for this trajectory are labeled with numbers from 1 to 6; potential energies along the trajectory are shown in Fig. 8a. Yellow lines (right panels) correspond to the example trajectory displayed in Fig. 8b.

Fig. 6 Evolution of dihedral angles α ((a), (b)) and β ((c), (d)) along the trajectories with initial M-conformers. Panels on the left ((a) and (c)) represent trajectories that isomerize, whereas panels on the right ((b) and (d)) account for trajectories without isomerization. Orange and blue represent the electronic state (S₁ or S₀, respectively) of the molecule along the trajectory. The black lines (left panels) belong to an example trajectory that undergoes thermal helix inversion in the ground state after E–Z isomerization; associated structures for this trajectory are labeled with numbers from 1 to 6.

Excited state lifetimes were determined using a monoexponential decay model. We find that M-conformers decay to the ground state faster than P-conformers, with excited state lifetimes of τ = 0.93 ps and τ = 1.31 ps, respectively (Fig. 7). Furthermore, the excited state lifetimes for molecules that underwent photoisomerization were τ = 1.21 ps for M-helicity and τ = 1.34 ps for P-helicity conformers. Thus, for both conformers, average lifetimes are larger for the reactive channel (E–Z isomerization) than for the unreactive channel (no E–Z isomerization). Inspecting the population changes (Fig. 7) more closely, we note that for both conformers, trajectories relaxing to S₀ before 0.5 ps are all unreactive. Thus, there seems to be an induction period of about 0.5 ps that is needed for trajectories to successfully isomerize. Inspecting the time evolution of α of reactive (Fig. 5a and 6a) unreactive (Fig. 5b and 6b) trajectories, we see that reactive trajectories start to reach values close to 90° (for M helicity) and −90° (for P helicity) at about 0.5 ps, while still being in the S₁ state. Unreactive trajectories, in contrast, have decayed to S₀ to a large fraction within the first 0.5 ps, before reaching 90° (for M) and −90° (for P). However, a few exceptions exist, where unreactive trajectories for both M and P conformers persist in the excited state beyond 0.5 ps. We further notice that when inspecting a density representation of the time evolution of α (Fig. S6, ESI†), that for M-conformers the average of α changes more drastically in the first 0.5 ps in reactive trajectories than in unreactive trajectories. For P-conformers, the density representation (Fig. S7, ESI†) reveals that starting structures of unreactive trajectories (panel b) exhibit less negative average α-values than starting structures for reactive trajectories (panel a).

Fig. 7 Population of each electronic state as a function of time for M- and P-like conformers. Color code: S₁ in green solid line, S₀ in red solid line, S₀ with successful isomerization in blue dashed line and S₀ with unsuccessful isomerization in orange dotted line.

These findings indicate that premature decay to S₀, before ≈0.5 ps and before reaching a value of ±90°, prevents trajectories from isomerizing. Since an α angle of ±90° is associated with a S₁–S₀ conical intersection, this is also reflected in the time evolution of the potential energies (Fig. 8) and the S₁–S₀ energy gap (Fig. 9): a large fraction of the unreactive trajectories (example in Fig. 8b) decay to the ground state before reaching a conical intersection, while reactive trajectories (example in Fig. 8a) persist long enough in S₁ to encounter a conical intersection, before they relax to S₀. This is further confirmed by inspecting the S₀–S₁ energy gap as a function of time (Fig. 9): for example, comparing reactive anti-M trajectories (panel a) with unreactive anti-M trajectories (panel b), we see that fewer unreactive trajectories reach energy gaps below 1 eV than reactive trajectories. We also note that about one half (52%) of all reactive trajectories (Fig. 9, panels a, c, e and g) decay to ground state due to negative excitation energies, indicating the proximity of a conical intersection; in unreactive trajectories this occurs in only 6% of all trajectories (Fig. 9, panels b, d, f and h).

Fig. 8 Potential energy of the S₂, S₁, and S₀ electronic states along (a) a reactive trajectory with helix inversion, (b) an unreactive trajectory. Structures displayed in Fig. 5 are marked and labeled for the reactive trajectory. The unreactive trajectory corresponds to the one depicted in Fig. 5, panels (b) and (d) in yellow, and structures are labeled.

Fig. 9 S₀–S₁ energy gap along the trajectories of each conformer. Left panels show reactive trajectories: (a) for anti-M, (c) for syn-M, (e) for anti-P and (g) for syn-P; right panels show unreactive trajectories: (b) for anti-M, (d) for syn-M, (f) for anti-P and (h) for syn-P. Green: trajectory is in S₁; red: trajectory is in S₀.

The time evolution of the NAC (Fig. S16, ESI†) shows that reactive trajectories exhibit larger magnitudes than unreactive trajectories, which is expected because they come close to the conical intersections, whereas most unreactive trajectories do not reach conical intersections. These trajectories decay to S₀, despite the low NAC magnitudes, which can occur based on the statistical nature of the TDDFT-SH approach.

Further inspecting the dynamics of the hot ground state after excited state decay, we observe several cases of thermal helix inversion, which can be identified by a sign change of β (Fig. 5 and 6, panels c and d). In the case of the reactive trajectories of the M-conformers, 9 trajectories exhibit helix inversion from the primarily formed P-conformer (positive β) to the more stable M-conformer (negative β), as exemplified by the black line in Fig. 5c and the snapshot structures ①–⑥: the primarily formed isomerized structure after isomerization (structure ② in Fig. 5) has P-helicity, whereas at later times (structures ④–⑥) the molecules adopts M-helicity. Also, the unreactive M-helical trajectories exhibit a few cases (7) of at least temporary thermal helix inversion, forming the less stable P-conformer, as indicated by a change from negative β to positive β (Fig. 5d). In the case of the reactive trajectories of P-conformers, we only see one trajectory (Fig. 6c, black line) that changes from the initially formed M-conformer (Fig. 6, structure ②) to the less stable P-conformer (Fig. 6, structure ⑤). For the unreactive trajectories of the P-conformers, however, we observe 10 cases where the stable M-conformer is formed after decay to the ground state, as indicated by a change of β from positive to negative (Fig. 6d). In summary, we see that thermal helix inversion can occur in the hot ground state within 3 ps upon photoexcitation in both M- and P-conformers. However, changes from P to M are generally more likely to occur, as expected from stability considerations. Thus, the net effect of thermal helix inversion in the hot ground state favors unidirectional rotation.

3.5. Wavelength-dependent product quantum yields

We computed the wavelength-dependent product quantum yield for the isomerization of anti-M and syn-P conformers (Fig. 10a), based on eqn (1). The syn-M and anti-P conformers were ignored due to their negligible population at thermal equilibrium compared to anti-M and syn-P, respectively (Table 1). Focusing on the wavelength range between 320 to 520 nm, where sufficient sampling trajectories are present, we notice that with the exception of the long wavelength region above 460 nm, the anti-M conformer exhibits a significantly higher isomerization quantum yield than the syn-P conformer. For anti-M, the average absorption spectrum of initial structures of reactive trajectories (Fig. 10b, red, dashed) appears to be very similar to the total absorption spectrum (red, solid). Therefore, the maximum of clockwise isomerization quantum yield at ≈340 nm (Fig. 10a, red) is located close to the maximum absorption band at of the M-conformer at 360 nm (Fig. 10b, red). At this point, the probability of isomerization reaches approximately 60%. In contrast, the anti-clockwise isomerization quantum yield of syn-P exhibits two maxima red- and blue-shifted with respect to its maximum absorption band at 450 nm. The red-shifted peak reaches about 5%, whereas the blue-shifted peak reaches about 10%. Although the quantum yield of syn-P is nearly negligible at the peak of anti-M's major absorption band, excitation at around 375 nm results in a quantum yield of zero for syn-P, while for anti-M it remains around 50%. We conclude that excitation at 375 nm, would therefore maximize clockwise rotation by selectively targeting anti-M conformers. Nonetheless, at thermal equilibrium, the population of syn-P is approximately 300 times lower than that of anti-M (Table 1). As a result, excitation at the peak wavelength of anti-M's major absorption band will almost guarantee clockwise rotation.

Fig. 10 (a) Wavelength-dependent product quantum yield for anti-M and syn-P. (b) Normalized absorption spectra of all initials structures of the anti-M and syn-P conformers and absorption spectra of the initial structures of trajectories that isomerized (dashed). Vertical sticks represent the excitation energy and oscillator strength for the transition S₀ → S₁ of the initial structures.

4. Conclusion

In this work, we simulated the photodynamics of four conformers of a second-generation molecular rotor – anti-M, anti-P, syn-M, and syn-P – to investigate its ultrafast photoisomerization mechanism.

Our TDDFT-SH calculations confirm that the M- and P-conformers rotate in opposite directions, as previously hypothesized:¹⁸ M-helical conformers exclusively rotate clockwise, while P-helical conformers exclusively rotate anti-clockwise. Although our calculations suggest that the anti-M conformer is the only one significantly populated at thermal equilibrium, for non-equilibrium situations, which in practice could occur at high photon influx, all conformers are possibly populated. Therefore, it is advantageous if anti-clockwise rotation could be suppressed in order to maximize unidirectional rotation. Our results predict a significantly higher clockwise isomerization quantum yield for the anti-M than anti-clockwise isomerization quantum yield for the syn-P, which has the greatest stability among the P-conformers. Moreover, at around 375 nm there is no probability of isomerization of syn-P, while for anti-M the probability is around 50%, coinciding with the major absorption band of the latter. Thus we predict that excitation at 375 nm maximizes unidirectional rotation.

Furthermore, our simulations reveal that thermally induced helix inversion can occur in the hot ground state directly after excited state decay. While this has been observed in both directions, it generally favors unidirectional motion. Interestingly, non-reactive decay of P-conformers accelerates helix inversion to form the more stable M-conformer.

However, our simulations completely neglected solvent effects, which may influence both the thermal helix inversion rates and the isomerization quantum yields due to solvent viscosity. Since our results demonstrate that premature excited-state decay negatively impacts isomerization quantum yields, it is reasonable to expect that increased solvent viscosity would further reduce the quantum yield by hindering trajectories from reaching the conical intersection.

Regarding helix inversion in the hot ground state, solvent effects are also likely to decrease inversion rates. Acting as a heat bath, the solvent would dissipate nuclear kinetic energy generated during excited-state decay, making it more difficult to overcome the energetic barriers of the inversion process.

Nevertheless, our study provides a general strategy for evaluating the wavelength dependence of isomerization quantum yields, offering a means to optimize unidirectional rotation in CPNX.

Data availability

The data supporting this study are provided in the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health (NIH) under award number R16GM149410. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. We acknowledge technical support from the Division of Information Technology of CSULB. This work was partially supported by the Deutsche Forschungsgemeinschaft [DFG (German Science Foundation) grant TRR 325-444632635].

References

1. V. Garca-López, D. Liu and J. M. Tour, Light-activated organic molecular motors and their applications, Chem. Rev., 2019, 120, 79–124.
2. V. Garca-López, F. Chen, L. G. Nilewski, G. Duret, A. Aliyan, A. B. Kolomeisky, J. T. Robinson, G. Wang, R. Pal and J. M. Tour, Molecular machines open cell membranes, Nature, 2017, 548, 567–572.
3. T. Stauch and A. Dreuw, Predicting the Efficiency of Photoswitches Using Force Analysis, J. Phys. Chem. Lett., 2016, 7, 1298–1302.
4. N. Koumura, R. Zijlstra, R. van Delden, N. Harada and B. Feringa, Light-driven monodirectional molecular rotor, Nature, 1999, 401, 152–155.
5. M. Klok, L. P. B. M. Janssen, W. R. Browne and B. L. Feringa, The influence of viscosity on the functioning of molecular motors, Faraday Discuss., 2009, 143, 319–334.
6. P. Roy, A. S. Sardjan, W. R. Browne, B. L. Feringa and S. R. Meech, Excited State Dynamics in Unidirectional Photochemical Molecular Motors, J. Am. Chem. Soc., 2024, 146, 12255–12270.
7. C. L. F. van Beek and B. L. Feringa, Coupled Rotary Motion in Molecular Motors, J. Am. Chem. Soc., 2024, 146, 5634–5642.
8. N. Koumura, E. M. Geertsema, M. B. van Gelder, A. Meetsma and B. L. Feringa, Second Generation Light-Driven Molecular Motors. Unidirectional Rotation Controlled by a Single Stereogenic Center with Near-Perfect Photoequilibria and Acceleration of the Speed of Rotation by Structural Modification, J. Am. Chem. Soc., 2002, 124, 5037–5051.
9. D. Pijper, R. A. van Delden, A. Meetsma and B. L. Feringa, Acceleration of a Nanomotor: Electronic Control of the Rotary Speed of a Light-Driven Molecular Rotor, J. Am. Chem. Soc., 2005, 127, 17612–17613.
10. M. M. Pollard, A. Meetsma and B. L. Feringa, A redesign of light-driven rotary molecular motors, Org. Biomol. Chem., 2008, 6, 507–512.
11. J. Vicario, M. Walko, A. Meetsma and B. L. Feringa, Fine Tuning of the Rotary Motion by Structural Modification in Light-Driven Unidirectional Molecular Motors, J. Am. Chem. Soc., 2006, 128, 5127–5135.
12. A. Nikiforov, J. A. Gamez, W. Thiel and M. Filatov, Computational Design of a Family of Light-Driven Rotary Molecular Motors with Improved Quantum Efficiency, J. Phys. Chem. Lett., 2016, 7, 105–110.
13. M. Filatov, M. Paolino, R. Pierron, A. Cappelli, G. Giorgi, J. Léonard, M. Huix-Rotllant, N. Ferré, X. Yang, D. Kaliakin, A. Blanco González and M. Olivucci, Towards the engineering of a photon-only two-stroke rotary molecular motor, Nat. Commun., 2022, 13, 6433.
14. M. Filatov, M. Paolino, D. Kaliakin, M. Olivucci, E. Kraka and S. K. Min, Impact of solvation on the photoisomerisation dynamics of a photon-only rotary molecular motor, Commun. Phys., 2024, 7, 219.
15. J. Wen, S. Mai and L. González, Excited-State Dynamics Simulations of a Light-Driven Molecular Motor in Solution, J. Phys. Chem. A, 2023, 127, 9520–9529.
16. A. Kazaryan, Z. Lan, L. V. Schäfer, W. Thiel and M. Filatov, Surface Hopping Excited-State Dynamics Study of the Photoisomerization of a Light-Driven Fluorene Molecular Rotary Motor, J. Chem. Theory Comput., 2011, 7, 2189–2199.
17. C. García-Iriepa, M. Marazzi, F. Zapata, A. Valentini, D. Sampedro and L. M. Frutos, Chiral Hydrogen Bond Environment Providing Unidirectional Rotation in Photoactive Molecular Motors, J. Phys. Chem. Lett., 2013, 4, 1389–1396.
18. N. Koumura, E. M. Geertsema, A. Meetsma and B. L. Feringa, Light-Driven Molecular Rotor: Unidirectional Rotation Controlled by a Single Stereogenic Center, J. Am. Chem. Soc., 2000, 122, 12005–12006.
19. G. Díaz Mirón, C. R. Lien-Medrano, D. Banerjee, M. Monti, B. Aradi, M. A. Sentef, T. A. Niehaus and A. Hassanali, Non-adiabatic Couplings in Surface Hopping with Tight Binding Density Functional Theory: The Case of Molecular Motors, J. Chem. Theory Comput., 2024, 20, 10602–10614.
20. C. Ash, M. Dubec, K. Donne and T. Bashford, Effect of wavelength and beam width on penetration in light-tissue interaction using computational methods, Lasers Med. Sci., 2017, 32, 1909–1918.
21. D. Liu, V. Garca-López, R. S. Gunasekera, L. Greer Nilewski, L. B. Alemany, A. Aliyan, T. Jin, G. Wang, J. M. Tour and R. Pal, Near-infrared light activates molecular nanomachines to drill into and kill cells, ACS Nano, 2019, 13, 6813–6823.
22. V. García-López, D. Liu and J. M. Tour, Light-Activated Organic Molecular Motors and Their Applications, Chem. Rev., 2020, 120, 79–124.
23. D. R. Pooler, A. S. Lubbe, S. Crespi and B. L. Feringa, Designing light-driven rotary molecular motors, Chem. Sci., 2021, 12, 14964–14986.
24. J. C. Tully, Molecular Dynamics with Electronic Transitions, J. Chem. Phys., 1990, 93, 1061.
25. R. Augulis, M. Klok, B. L. Feringa and P. H. V. Loosdrecht, Light-driven rotary molecular motors: an ultrafast optical study, Phys. Status Solidi C, 2009, 6, 181–184.
26. M. Klok, W. R. Browne and B. L. Feringa, Kinetic analysis of the rotation rate of light-driven unidirectional molecular motors, Phys. Chem. Chem. Phys., 2009, 11, 9124–9131.
27. M. Klok, N. Boyle, M. T. Pryce, A. Meetsma, W. R. Browne and B. L. Feringa, MHz unidirectional rotation of molecular rotary motors, J. Am. Chem. Soc., 2008, 130, 10484–10485.
28. P. Pracht, F. Bohle and S. Grimme, Automated exploration of the low-energy chemical space with fast quantum chemical methods, Phys. Chem. Chem. Phys., 2020, 22, 7169–7192.
29. A. Laio and M. Parrinello, Escaping free-energy minima, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
30. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB–An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671.
31. M. Häser and R. Ahlrichs, Improvements on the direct SCF method, J. Comput. Chem., 1989, 10, 104.
32. J. P. Perdew, M. Ernzerhof and K. Burke, Rationale for mixing exact exchange with density functional approximations, J. Chem. Phys., 1996, 105, 9982–9985.
33. C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys., 1999, 110, 6158–6170.
34. A. Schäfer, H. Horn and R. Ahlrichs, Fully optimized contracted Gaussian basis sets for atoms Li to Kr, J. Chem. Phys., 1992, 97, 2571–2577.
35. A. Schäfer, C. Huber and R. Ahlrichs, Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr, J. Chem. Phys., 1994, 100, 5829–5835.
36. F. Weigend, Accurate Coulomb-fitting basis sets for H to Rn, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
37. F. Weigend and R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
38. P. Deglmann, F. Furche and R. Ahlrichs, An efficient implementation of second analytical derivatives for density functional methods, Chem. Phys. Lett., 2002, 362, 511.
39. P. Deglmann and F. Furche, Efficient characterization of stationary points on potential energy surfaces, J. Chem. Phys., 2002, 117, 9535.
40. O. Christiansen, H. Koch and P. Jørgensen, The second-order approximate coupled cluster singles and doubles model CC2, Chem. Phys. Lett., 1995, 243, 409–418.
41. TURBOMOLE V6.3, TURBOMOLE GmbH, Karlsruhe, 2011; available from https://www.turbomole.com.
42. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding and C. Hättig, et al., TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations, J. Chem. Phys., 2020, 152, 184107.
43. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begusic, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov and S. Fürst, et al., TURBOMOLE: Today and tomorrow, J. Chem. Theory Comput., 2023, 19, 6859–6890.
44. Y. Sugita and Y. Okamoto, Replica-exchange molecular dynamics method for protein folding, Chem. Phys. Lett., 1999, 314, 141.
45. C. Cisneros, T. Thompson, N. Baluyot, A. C. Smith and E. Tapavicza, The role of tachysterol in vitamin D photosynthesis - a non-adiabatic molecular dynamics study, Phys. Chem. Chem. Phys., 2017, 19, 5763–5777.
46. E. Tapavicza, Time-Dependent Density Functional Theory, Jenny Stanford Publishing, 2022, pp. 141–197.
47. D. A. Case, T. E. Cheatham III, T. Darden, H. Gohlke, R. Luo, K. M. Merz Jr., A. Onufriev, C. Simmerling, B. Wang and R. J. Woods, The Amber biomolecular simulation programs, J. Comput. Chem., 2005, 26, 1668–1688.
48. J. Shao, S. W. Tanner, N. Thompson and T. E. Cheatham, Clustering Molecular Dynamics Trajectories: 1. Characterizing the Performance of Different Clustering Algorithms, J. Chem. Theory Comput., 2007, 3, 2312–2334.
49. M. Ester; H.-P. Kriegel; J. Sander and X. Xu, A density-based algorithm for discovering clusters in large spatial databases with noise. Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, 1996, p 226–231.
50. C. Hättig and A. Köhn, Transition moments and excited-state first-order properties in the coupled-cluster model CC2 using the resolution-of-the-identity approximation, J. Chem. Phys., 2002, 117, 6939–6951.
51. S. A. Epstein, E. Tapavicza, F. Furche and S. A. Nizkorodov, Direct photolysis of carbonyl compounds dissolved in cloud and fog droplets, Atmos. Chem. Phys., 2013, 13, 9461–9477.
52. C. W. Grathwol, N. Wössner, S. Swyter, A. C. Smith, E. Tapavicza, R. K. Hofstetter, A. Bodtke, M. Jung and A. Link, Azologization and repurposing of a hetero-stilbene-based kinase inhibitor: towards the design of photoswitchable sirtuin inhibitors, Beilstein J. Org. Chem., 2019, 15, 2170–2183.
53. E. Tapavicza, G. F. von Rudorff, D. O. De Haan, M. Contin, C. George, M. Riva and O. A. von Lilienfeld, Elucidating an Atmospheric Brown Carbon Species–Toward Supplanting Chemical Intuition with Exhaustive Enumeration and Machine Learning, Environ. Sci. Technol., 2021, 55, 8447–8457.
54. D. L. Sofferman, A. Konar, J. N. Mastron, K. G. Spears, C. Cisneros, A. C. Smith, E. Tapavicza and R. J. Sension, Probing the Formation and Conformational Relaxation of Previtamin D3 and Analogues in Solution and in Lipid Bilayers, J. Phys. Chem. B, 2021, 125, 10085–10096.
55. E. Tapavicza, A. M. Meyer and F. Furche, Unravelling the details of vitamin D photosynthesis by non-adiabatic molecular dynamics simulations, Phys. Chem. Chem. Phys., 2011, 13, 20986.
56. E. Tapavicza, G. D. Bellchambers, J. C. Vincent and F. Furche, Ab initio non-adiabatic molecular dynamics, Phys. Chem. Chem. Phys., 2013, 15, 18336–18348.
57. E. Tapavicza, T. Reutershan and T. Thompson, Ab initio simulation of the ultrafast circular dichroism spectrum of provitamin D ring-opening, J. Phys. Chem. Lett., 2023, 14, 5061–5068.
58. F. Furche and R. Ahlrichs, Adiabatic time-dependent density functional methods for excited state properties, J. Chem. Phys., 2002, 117, 7433.
59. S. Hirata and M. Head-Gordon, Time-dependent density functional theory within the Tamm-Dancoff approximation, Chem. Phys. Lett., 1999, 314, 291–299.
60. B. G. Levine, C. Ko, J. Quenneville and T. J. Martinez, Conical intersections and double excitations in density functional theory, Mol. Phys., 2006, 104, 1039.
61. E. Tapavicza, I. Tavernelli, U. Rothlisberger, C. Filippi and M. E. Casida, Mixed time-dependent density-functional theory/classical trajectory surface hopping study of oxirane photochemistry, J. Chem. Phys., 2008, 129, 124108.
62. I. Tavernelli, E. Tapavicza and U. Rothlisberger, Nonadiabatic coupling vectors within linear response time-dependent density functional theory, J. Chem. Phys., 2009, 130, 124107.
63. I. Tavernelli, E. Tapavicza and U. Rothlisberger, Non-adiabatic dynamics using time-dependent density functional theory: Assessing the coupling strengths, J. Mol. Struct. THEOCHEM, 2009, 914, 22.
64. R. Send and F. Furche, First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance, J. Phys. Chem., 2010, 132, 044107.
65. T. Thompson and E. Tapavicza, First-Principles Prediction of Wavelength-Dependent Product Quantum Yields, J. Phys. Chem. Lett., 2018, 9, 4758–4764.
66. E. Tapavicza, T. Thompson, K. Redd and D. Kim, Tuning the photoreactivity of Z-hexatriene photoswitches by substituents - a non-adiabatic molecular dynamics study, Phys. Chem. Chem. Phys., 2018, 20, 24807–24820.
67. B. Oruganti, C. Fang and B. Durbeej, Computational design of faster rotating second-generation light-driven molecular motors by control of steric effects, Phys. Chem. Chem. Phys., 2015, 17, 21740–21751.

---

Footnotes

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

‡ Equal contribution.

---