Chiral conversion and periodical decay in bridged-azobenzene photoisomerization: an ab initio on-the-fly nonadiabatic dynamics simulation

Abstract

With the implementation of global nonadiabatic switching probability algorithm based on the Zhu–Nakamura theory, we performed on-the-fly trajectory surface hopping simulations at the CASSCF level to analyze the cis ↔ trans photoisomerization mechanism of bridged-azobenzene upon S₁ excitation. The cis-to-trans isomerization process was initiated by the pedal-like twist of –N=N– moiety and observed hopping to S₀ mainly via the cis-form conical intersection region. While, the trans-to-cis process is induced by the nonsymmetrical torsion motion toward the central eight-membered ring and encountered surface hopping through another conical region lying in the middle of the isomerization pathway. In cis-to-trans photoisomerization, we also observed a particular chiral conversion pathway with surface hopping occurs at CI-trans conical region. The trajectories followed this pathway possess relative longer excited state lifetimes as the ∼40 fs initial stage oscillation in the cis Franck–Condon region. Moreover, the periodical population decay pattern of S₁ state in trans-to-cis process was observed with an interval of ∼16 fs.

---

1. Introduction

Azobenzene (Ab) derivatives have been widely studied experimentally and theoretically for their peculiar photochromic properties which enable them to act as potential candidates for molecular motors, optical data storage and photobiological devices etc.^1–3 The Ab molecule was proposed to photoisomerize via torsion, inversion or the mixed inversion-assisted-torsion pathways,^4–11 while in the Ab derivatives with geometrical hindrance, these isomerization pathways may be strengthened, weakened or even prohibited.^2,12–16

Many efforts have been devoted to enhance the photoresponse and photostability of azobenzene derivatives. The bridged-azobenzene (b-Ab) molecule which has been synthesized for decades^17,18 was recently found to exhibit superior photochromic properties than the parent Ab molecule and therefore it has become a topic of interests.^19–28 In the experimental investigations of Siewertsen et al.,^19,20 the well separated S₁(nπ*) absorption bands at λ = 404 nm (cis) and 490 nm (trans) together with large numbers of switching cycles was reported. They have also proposed that with an additional –CH₂–CH₂– bridge between the two benzene rings, the thermodynamics stability of cis and trans conformers has exchanged compare to parent Ab as the trans-b-Ab no longer adopts the planar equilibrium structure. Moreover, the severe conformational constraints led to a much accelerated trans → cis photoisomerization, however, comparing with the parent Ab, the cis → trans process is hardly affected. Meanwhile, relative higher quantum yields of Φ_cis→trans = 0.72 ± 0.04 and Φ_trans→cis = 0.50 ± 0.10 with corresponding time constants τ₁ = 70 (<50) and τ₂ = 270 (320) were derived.²⁰

Attracted by the high efficient photo-induced conformation conversion, many nonadiabatic trajectory surface hopping (TSH) molecular dynamics (MD) simulations have been performed to reveal the photoisomerization mechanisms of b-Ab. Böckmann and co-workers²¹ investigated the photoisomerization process of trans-b-Ab by using ab intio MD simulation at DFT level with the Tully's fewest-switches algorithm accounted for the nonadiabatic surface hopping. They have attributed the enhanced isomerization efficiency to favorable out-of-plane conformational preorientation of phenyl rings by the ethylenic bridge via ultrafast pedal motion of N atoms. Later, they have extended the investigation to both trans → cis and cis → trans processes with the same methodology.²⁶ Contrary to the accelerated trans → cis isomerization, the cis → trans process was found to be slightly hindered by a significant rearrangement of the bridge in the S₁ state. They proposed that both of the isomerization processes were proceeded by the NN twist mechanisms with equally or unequally contributing nitrogen sites, respectively, as those observed in parent Ab.^29,30 Furthermore, they have observed a nonreactive chirality switch induced by pseudorotation of the central 8-membered ring system in cis-b-Ab. In the dynamics simulation work of Carstensen et al.²² at the AM1 floating-occupation configuration interaction (FOCI-AM1) level, they proposed that the ethylenic bridge does not provide a steric hindrance in the cis → trans isomerization while the complex interplay of motions towards the central 8-membered ring blocks the forward reaction in the ground state and result in a reduced quantum yields. Based on potential energy surfaces (PESs) constructed at the CASSCF level with active space composed of 14 electrons in 12 orbitals, they have concluded that the cis → trans and the trans → cis photoisomerization proceed via different conical intersection regions. The two processes exhibit different features: barrierless versus ground state barrier and steep versus flat potential energy landscape, respectively. Density-functional based semiclassical electron-radiation-ion dynamics (SERID) simulations on the topic reaction system have been performed by Liu et al.²³ and Jiang and coworkers.^24,25 The former investigation²³ also provides PESs of S₁ and S₀ states constructed at CASPT2 level based on CASSCF optimized geometries and confirmed that only one conical intersection exists between S₁ and S₀ states. Moreover, they proposed that the internal energy redistribution controls the relaxation process from the S₁ Franck–Condon region to the conical region and therefore determines the lifetime. Jiang and coworkers^24,25 also confirmed that the ethylenic bridge does not hinder the rotation pathway and the S₁ → S₀ internal conversions occur in the same conical region in the cis ↔ trans processes. Surface hopping dynamics simulations based on the original Zhu–Nakamura theory has been carried on by Gao et al.,^27,28 which confirmed that the torsion around –N=N– moiety together with the twisting of phenyl rings in opposite direction is the dominant mechanism. They also reproduced the two types of CI regions and proposed that in the trans → cis and cis → trans isomerization, the azo-moiety initially approaches product region and then the two phenyl rings and ethylenic bridge adjust themselves to complete the process.

Although such a lot dynamics simulations have been put forward on the photoisomerization of b-Ab, the mechanism for chirality conversion in the cis → trans photoisomerization and the properties of population decay on the excited state are still unclear. The main reasons are the insufficiency of the DFT method in dealing with conical intersections and the limited number of simulated trajectories in previous theoretical investigations. Thus a large scale full dimensional on-the-fly trajectory based nonadiabatic molecular dynamics simulation investigation employing high level ab initio calculation is required to resolve the existed problems. In this work, with implementation of recently developed Zhu–Nakamura nonadiabatic switching probability calculation formula which only using adiabatic potential energy surface and its gradients,⁵ we have performed an extensive on-the-fly trajectory surface-hopping dynamics simulations at the CASSCF ab initio quantum level to probe both the cis → trans and the trans → cis photoisomerization mechanism of the bridged-azobenzene. The quantum yields and lifetimes are estimated by analyzing 400 trajectories for cis-to-trans and 550 trajectories for trans-to-cis processes and the roles of the pathways involved in the processes will also be discussed. In the next section, we will give a brief description on the present trajectory surface-hopping method and ab initio method incorporated on-the-fly into the nonadiabatic molecular dynamic simulation algorithm. In Section 3, we will provide detailed analyses and explanations on the b-Ab photoisomerization mechanisms upon S₁ excitation.

2. Theoretical methods and computational details

Novel trajectory surface hopping method based on the improved Zhu–Nakamura theory (ZN switches) has been successfully applied to investigate nonadiabatic photodynamics of several reaction systems^4,5,31 and the algorithm is briefly described as follows. The on-the-fly trajectory is propagated by numerically integrating the Newtonian equation of motion with the velocity-Verlet method³² with the potential energy and its' gradients calculated by CASSCF method step by step. Along trajectory propagation, we detect avoided crossing point by computing minimum separation of two adjacent adiabatic potential energy surfaces within three consecutive time steps (this is usually executed in conical intersection region). At this avoided crossing point, we compute global nonadiabatic switching probability according to the improved Landau–Zener formula by Zhu and Nakamura,^33,34in which two unitless parameters, namely effective coupling and effective collision energy, are given byandin which, F₁ and F₂ are forces on two diabatic potential energy surfaces, V₁₂ is diabatic coupling, μ is reduced mass of diatomic molecule, E_X is energy at crossing point and E_t is potential energy plus kinetic energy component in direction of hopping vector (note + (−) in eqn (1) stands for F₂F₁ > 0 (F₂F₁ < 0)). Diabatic coupling V₁₂ in eqn (2) is directly computed by the energy separation between two adiabatic potential energy surfaces at avoided crossing. At this avoided crossing, we convert mass-scaled multidimensional forces into mass-scaled one-dimensional forces as given in eqn (4) and (5),⁵andwhere N is number of nuclei in molecule with the mass m_i (i = 1,2,…,N), α stands for x, y, and z component of Cartesian coordinates for the i-th nucleus, and multidimensional diabatic forces F₁^iα and F₂^iα are generalized from adiabatic forces within three consecutive time steps around detected avoided crossing along the trajectory. Detailed description along with definition of hopping direction and momentum change at an attempted trajectory hopping position is given in the ref. 4.

The simulated trajectories are started from Franck–Condon region of cis- or trans-b-Ab in S₁ state, and initial velocities are set up by the procedure summarized in ref. 4. The thermal kinetic energy with T = 300 K is added to all sampling trajectories with randomly distributing into initial Wigner velocities. Following the previous investigation in azobenzene,^4,5 the 0.5 fs time-step is selected in the entire simulation on b-Ab. For trajectories starting from both cis- and trans-b-Ab, the time limit is set up to 200 fs. In total, 400 trajectories for cis → trans photoisomerization process and 550 trajectories for the trans → cis case were run to get converged excited state lifetimes and quantum yields.

In the current work, the structures of cis- and trans-b-Ab in S₀, trans-b-Ab in S₁ and the minimum energy conical intersections (CI) between the S₀ and S₁ states were optimized by using the two-state averaged CASSCF method (SA2-CASSCF) with active spaces composed of 10 electrons in 8 orbitals (10e, 8o) and 6 electrons in 5 orbitals (6e, 5o), referred as CASSCF(10,8) and CASSCF(6,5), respectively. Based on the consideration of both the accuracy and efficiency, (which will be discussed in detail in next section) the CASSCF(6,5) was selected in the potential energy and energy gradients calculation embedded in molecular dynamics simulation. All the calculations were carried out with the 6-31G basis set. The CIs were optimized by using the MOLPRO 2009.1 package^35–37 and the optimization of other geometries as well as on-the-fly potential energies and their gradients' calculations were performed with the MOLCAS 7.5 program.³⁸ To account for the dynamical electron correlation, the corrected energy for optimized geometries at complete active space second order perturbation theory (CASPT2) and internally contracted multi-reference configuration interaction (ICMRCI) level were performed by employing the MOLCAS 7.5 and Xi'an CI program of our group,³⁹ respectively.

3. Results and discussions

3.1 Stationary structures and possible isomerization pathways

To explore the topology of potential energy surface along the photoisomerization channels starting from both isomers at S₁, we optimize the ground state minima, transition state and minimum-energy conical intersections between the S₁ and S₀ states at the SA2-CASSCF(10,8) and SA2-CASSCF(6,5), respectively. The optimized geometries at CASSCF level with different active spaces are listed in Table 1 and the relative energies with respect to S₀ cis and vertical excitation energies for these geometries are given in Table S1.† The atomic numbering is displayed in Scheme 1. Moreover, the optimized Cartesian coordinates for minima and CIs with different chirality at the CASSCF(10,8) and CASSCF(6,5) level are listed in Table S2.†

Table 1 Optimized geometry parameters for the ground state minima and S₀/S₁ conical intersections^a

Method	Geometry	NN	CNN^b	CNNC^c	CCCC^d	CCNN^e	CCCC^f
a Bond length is given in angstroms (Å), and bond angles and dihedral angles are given in degrees.b C7N2N1/C1N1N2.c C7N2N1C1.d C9C14C13C3.e C3C1N1N2/C9C7N2N1.f C7C9C14C13/C1C3C13C14.g This work.h Ref. 27.i Ref. 22.j Ref. 23.k Ref. 26.l Ref. 19.
CASSCF (6e,5o)/6-31G^g	cis	1.26	123.3/122.2	7.1	40.4	69.3/−72.3	34.3/−85.5
	trans	1.27	112.8/113.3	145.4	100.0	−65.3/−63.3	−44.8/−46.0
	CI-cis	1.33	127.3/117.4	94.7	85.9	−3.8/−87.0	−13.7/−76.4
	CI-trans	1.33	130.8/122.0	97.8	111.9	−26.5/−51.1	−48.8/−54.7
CASSCF (10e,8o)/6-31G^g	cis	1.26	122.6/122.0	7.4	35.7	69.8/−73.1	39.5/−84.1
	trans	1.27	113.0/113.0	145.7	100.7	−64.4/−64.4	−45.9/−45.9
	CI-cis	1.33	127.5/117.8	95.0	86.4	−4.5/−86.4	−14.4/−76.2
	CI-trans	1.30	135.2/123.4	97.2	110.0	−21.8/−56.0	−47.1/−57.7
CASSCF (6e,5o)/6-31G^h	cis	1.26	122.8/122.1	6.4	37.7	70.1/−71.8
	trans	1.27	113.3/112.8	145.4	100.0	−63.2/−65.3
CASSCF (10e,8o)/6-31G*^h	cis	1.25	121.3/120.9	6.3	30.7	72.7/−75.6
	trans	1.25	111.1/111.1	148.3	97.1	−66.1/−66.1
CASSCF (14e,12o)/6-31G*^h	cis	1.25	121.4/120.9	5.8	31.3	73.1/−74.9	42.5/−82.0
	trans	1.25	111.1/111.1	147.5	96.2	−65.6/−65.6	−43.6/−43.6
	CI-cis	1.28	134.7/121.6	95.3	103.4	−33.9/−95.3
	CI-trans	1.28	136.2/121.6	99.5	108.8	−34.3/−43.2
FOCI-AM1^i	cis	1.26	121.4/122.6	6.2	8.3	66.3/78.1
	trans	1.27	114.9/114.9	138.6	20.0	66.7/66.6
	CI-cis		129.7/111.4	90.6	74.2	5.6/105.3
	CI-trans		132.1/132.0	95.5	111.9	40.8/39.9
CASSCF (10e,8o)/6-31G*^j	cis	1.26	118.6/118.2	5.9	27.2	74.3/−80.4	47.5/−84.5
	trans	1.28	114.1/114.1	145.8	108.4	−64.0/−64.0	−50.3/−50.3
DFT/6-31G*^k	cis	1.26	121.7/121.6	6.6	35.3	72/113.4
	trans	1.27	112.0/112.1	146.3	98.6	125/127
Expt^l	cis	1.24	122.6/122.2	6.5	37.3	71.3/−72.0	37.5/−83.4
	trans	1.25	112.4/112.4	146.8	96.5	−64.3/−64.5	−42.9/−43.1

---

Scheme 1 Atomic numbering (geometry of trans-b-Ab).

As shown in Table 1, the internal coordinates of ground state cis and trans conformers optimized at CASSCF(6,5) are in very good agreement with those obtained by CASSCF(10,8) and also reproduces the experimental values¹⁹ very well. These results indicate that the dynamics simulation with CASSCF(6,5) could reasonably describe the geometrical changes in the isomerization process. The accuracy of energies calculated by CASSCF(6,5) was further verified by comparing with the relative energies to the S₀ cis isomer calculated by the CASPT2 and ICMRCI based on the CASSCF optimized geometry. Comparing with the experimental value of 70.8 and 58.4 kcal mol⁻¹,¹⁹ respectively, the calculated vertical excitation energies for S₀ cis and trans are overestimated with CASSCF(6,5), even with the ICMRCI correction (84.7 and 77.4 kcal mol⁻¹). It can be seen from Table 1, the calculations with (10e, 8o) active space gives better vertical excitation energies. However, the relative energies computed by CASSCF, CASPT2 and ICMRCI methods with (10e, 8o) and (6e, 5o) active spaces at the optimized geometries are on the whole quite close. With consideration of the computation cost and accuracy, the CASSCF(6,5) was used in the following dynamics simulations to obtain the potential energies and their gradients along the trajectories.

With default conical intersection optimization procedure, we could only reproduce the four conical structures reported by Gao et al.²⁷ at CASSCF(10,8). (CI-trans type) By loosening convergence threshold in MOLPRO program,³⁷ we could also obtain similar conical geometries with CASSCF(6,5). Moreover, another type of conical regions, namely CI-cis, was roughly located by the dynamics simulation on the cis → trans isomerization with CASSCF(6,5), and then optimized by CASSCF(10,8) with default procedure and CASSCF(6,5) by loosening convergence threshold in MOLPRO.³⁷ This observations agree with previous reported CoIn-Z geometry by Carstensen et al.²² and CI₃/CI₄ by Gao et al.,^27,28 which were both optimized at the CASSCF level with larger active space involving 14 electrons in 12 orbitals. The four CI-trans type conical intersections are enantiomers mirroring to the perpendicular and parallel molecular plane, and can be distinguished by the signs of dihedral angles C7N2N1C1, C3C1N1N2 and C9C7N2N1. In which, the C7N2N1C1 angle is close to ∼95°, slightly smaller than CI-trans type CIs, but the other dihedral angles and bond angles relating to the central eight-membered ring are much different. In the process of approaching the CI-trans type conical regions, large amplitude benzene rings' open-close twist with respect to the central azo- and ethyl-bridge moiety and small amplitude N=N torsion are needed. While for the CI-cis case, the relative positions of the two benzene rings and ethyl-bridge stay close to cis-b-Ab, and the evident changes on C7N2N1C1 are mainly induced by the pedal-like twist motion of azo-bridge and related molecules. The two distinct types of conical regions would induce different photoisomerization mechanisms, and through which one the internal conversion occurs will be discussed in the next section in detail.

The trans-b-Ab molecule attracts more and more interests to researchers for its higher quantum yields and shorter excited state lifetimes than the parent trans-Ab. Here, we present the following explanations for its enhanced photoresponsibilities. Firstly, comparing with parent Ab, the S₁ state PES appears in the steeper slope from Franck–Condon region to CI-trans for the trans-b-Ab molecule (the energy difference is 6.9 kcal mol⁻¹ for trans-Ab⁴ and 39.2 kcal mol⁻¹ for trans-b-Ab). Thus, the trans-b-Ab molecule could obtain larger excited state driving force to arrive at the conical region faster on the excited state. Secondly, the trans-b-Ab no longer keeps flat, and thus do not need the long time vibrational redistribution to give enough kinetic energy into the out-of-plane motion as those in trans-Ab (∼350 fs in trans-Ab photisomerization⁵). The third reason is that the CNNC angle in trans-b-Ab (145.4) is smaller than that in b-Ab (180.0),⁴ while, in corresponding CI involved in the isomerization of these two trans form molecules are quite similar (97.8 and 97.4,⁴ respectively). Evidently, in trans-b-Ab, the isomerization pathway from trans side is shorter than the one in trans-Ab. With the above mentioned reasons, it is clear that the S₁ state lifetime for trans-b-Ab would be rather short and in the same time scales as the cis form conformer, which follow a steep downhill in both cis-b-Ab and parent cis-Ab. Moreover, the quantum yields of trans-b-Ab would also be enhanced by the larger driven force along isomerization path in S₁ state.

3.2 Photoisomerization reaction dynamics

To verify the reaction mechanism and explore the statistical properties for b-Ab photoisomerization upon S₁ excitation, we have performed ZN switches surface hopping dynamics simulations with 400 trajectories starting from cis and 550 trajectories starting from trans form. Based on the static electronic structure computations and dynamics simulations, we propose the following scenario for the photoisomerization processes of b-Ab as depicted in Fig. 1. Upon photoexcitation, the cis or trans conformers move along the respective isomerization pathways to visit conical regions and after surface hopping, the molecule may continue follow the isomerization path to product region or back to initial geometry. The trans → cis isomerization encounters surface hopping through CI-trans type conical regions, which are in agreement with previous investigations.²² While the cis → trans isomerization is more complicated, we have observed two competing processes corresponding to clockwise and counterclockwise torsion toward the N=N bridge. Moreover, the minority trajectories which experience helicity conversion (N=N clockwise motion) visit CI-trans type conical region, while the majority trajectories that keep helicity unchanged (N=N counterclockwise motion) hop to the S₀ state through CI-cis type conical region. The helicity conversion towards the CNNC dihedral angle in cis → trans photoisomerization has been reported for parent azobenzene^6–8 and also was proposed to be possible in the bridged-azobenzene.²⁶ In the current work, this chiral conversion pathway has been confirmed for the first time with implementation of the ZN switches surface hopping dynamics simulations.

Fig. 1 Potential energy profiles for the cis ↔ trans photoisomerization of bridged-azobenzene.

To further reveal the origin of discrepancy between cis → trans and trans → cis isomerization mechanism, we have depicted the distribution of successful hopping points' in terms of CNNC dihedral angle versus energy gap, effective nonadiabatic coupling factors and effective collision energy in Fig. 2 and in terms of the other important internal coordinates in Fig. 3. As in Fig. 2a and b, the CNNC dihedral angle is distributed in the range [−80°, −90°] and [75°, 95°] for cis → trans and [95°, 110°] for trans → cis. It is clear that the surface hopping takes place in the range of absolute value of CNNC smaller (larger) than 97.8 and 94.7 at the conical intersection for cis → trans (trans → cis). This also agrees with the parent azobenzene⁵ and reflects the different motion direction when the reactant approaches the conical regions. In the cis-to-trans and trans-to-cis cases, the potential energy differences at hopping spots are distributed in 0.00 to 0.40 eV and 0.00 to 0.50 eV regions, respectively (Fig. 2 upper panel). However, the distribution of the effective coupling parameter a² and the effective collision energy b² are all in the region of a² ∼ [0.1,10⁴] (Fig. 2 middle panel) and b² ∼ [2200] (Fig. 2 bottom panel) for both cis-to-trans and trans-to-cis cases. According to the ZN-switches algorithm, the dimensionless parameters a² and b² determine the switching probabilities at the conical zones and the similar distribution of these parameters in both isomerization processes demonstrate the reliability of the current simulation. Moreover, for trans-to-cis isomerization, more reactive surface hopping points distributed in the region with larger energy gap, and consequently result in preference on smaller a² and b² values. For cis-to-trans isomerization, the hopping points for reactive and nonreactive trajectories cover similar regions. With only CNNC related distributions shown in Fig. 2, we could not distinguish how the CI-trans and CI-cis type conical regions contribute to the photoisomerization and therefore we display distributions of the other important internal coordinates at hopping spots in Fig. 3. Fig. 3a and b show the distribution of C9C7N2N1 versus C7N2N1C1 at surface hopping points for the cis-to-trans and trans-to-cis processes respectively. It can be seen that the counterclockwise channel in cis → trans process visits the CI-cis type conical region with C9C7N2N1 distributed in the range [−90°, −140°], while the clockwise channel in cis → trans process and all the trajectories of trans → cis visit chiral symmetric CI-trans type conical regions with C9C7N2N1 distributed in the range [0°, 30°] and [−20°, −65°], respectively. The variation of C7C9C14C13 and C9C7N2N1 dihedral angles indicates the twist motion of the central eight-membered ring in b-Ab. As shown in Fig. 3c and d, counterclockwise channel in cis → trans undergo surface hopping in the region which has inverse C7C9C14C13 value of [0°, 60°] and much smaller C9C7N2N1 value around −120° compared with the trans → cis pathway. The CNN angles reflect the inversion motion toward the azo-bridge, as shown in Fig. 3e, with C1N1N2 distributed in the range of [95°, 112°]. The clockwise ones possess larger C7N2N1 angle, which may be related to the flip of the azo-bridge. In Fig. 3f, more reactive surface hopping points were found in the region where C7N2N1 is larger than 117° in trans → cis process, indicating that the enhanced torsion motion would assist the reactive pathway.

Fig. 2 The distribution of successful hopping points in terms of the CNNC dihedral angle. (a) The first panel is for energy gaps, the second panel is for effective coupling parameters, and the third panel is for effective collision energy in the case of cis-to-trans. (b) The same as (a) but in the case of trans-to-cis.

Fig. 3 Geometrical parameter distribution for dihedral angles C7N2N1C1 versus C9C7N2N1 (upper panel), dihedral angles C7C9C14C13 versus C9C7N2N1 (middle panel), bond angles C1N1N2 versus C7N2N1 (bottom panel) at the surface hopping points involve in cis-to-trans (a, c and e) and trans-to-cis (b, d and f) processes respectively.

Within the 400 trajectories started at S₁ cis, 12 trajectories follow the clockwise pathway (8 isomerize to trans conformer and 4 back to cis form) and 388 trajectories follow counterclockwise pathway (217 isomerize to trans conformer and 171 back to cis form). While, for 550 trajectories started at S₁ trans, 244 isomerize to cis conformer and 306 back to cis form. As the C9C7N2N1 is the most important dihedral angle involving in the motion of photoisomerization, time evolutions of C9C7N2N1 dihedral angle for the simulated trajectories started at cis and trans are depicted in Fig. 4a and b, respectively. The trajectories follow clockwise and counterclockwise channels are clearly demonstrated in Fig. 4a, and the clockwise ones need more time to arrive the conical region with C9C7N2N1 around 95° as the flip of molecular plane.

Fig. 4 Time evolution of the C7N2N1C1 dihedral angle along simulated trajectories starting from (a) cis-b-Ab and (b) trans-b-Ab.

The quantum yields and the excited state lifetimes are estimated by the following procedure based on the simulated data. The quantum yield equal to the ratio of reactive ones in all trajectories with an standard error . The isomerization lifetime is obtained by fitting the time-dependent decay curve of S₁ state population to the exponential function f(t) = e^{−(t−t_d)/τ}, where t_d is the onset time of the S₁ population loss, and τ is the time constant for the exponential decay. The lifetime is estimated as the value of τ + t_d.

Simulated quantum yield and the lifetime for cis-to-trans (trans-to-cis) converge to 0.57 and 64 fs (0.44 and 35.5 fs) with 400 (550) sampling trajectories as summarized in Table 2. In current work, the estimated quantum yields for both the cis-to-trans and trans-to-cis processes show decrease trend with more trajectories taken into account. The simulated results are in good agreement with experimental values^19,20 and the other previous simulations based on relatively limited number of trajectories to some extent.^21–28

Table 2 Convergence of simulated quantum yields and lifetimes for the b-Ab photoisomerization in the S₁ state with respect to the number of trajectories in comparison with selected simulations and experimental results

	Number of trajectories	From cis-bridged-azobenzene		From trans-bridged-azobenzene
		Quantum yield	Lifetime (fs)	Quantum yield	Lifetime (fs)
a Ref. 22.b Ref. 23. The quantum yield in parenthesis is obtained by exclude the resonance trajectories.c Ref. 25.d Ref. 26.e Ref. 28.f Ref. 20.
SA2-CAS(6,5)	35	0.77 ± 0.092	55.5	0.54 ± 0.150	34.5
	100	0.63 ± 0.077	62.0	0.52 ± 0.096	34.0
	150	0.57 ± 0.071	62.5	0.50 ± 0.082	34.0
	200	0.58 ± 0.061	63.5	0.50 ± 0.071	35.0
	250	0.56 ± 0.056	64.5	0.44 ± 0.071	35.5
	300	0.56 ± 0.051	64.0	0.44 ± 0.065	35.5
	350	0.55 ± 0.048	63.5	0.44 ± 0.061	35.0
	400	0.57 ± 0.044	64.0	0.45 ± 0.056	35.0
	450			0.44 ± 0.053	35.5
	500			0.44 ± 0.051	35.5
	550			0.44 ± 0.048	35.5
FOCI-AM^a	250	0.23 ± 0.030		0.58 ± 0.030
DFTB^b	15	0.53 (0.73)		∼0.50
DFTB^c	65	0.49	∼200.0	0.45
DFT^d	30	0.32 ± 0.050		0.47 ± 0.100
CASSCF^e	35	0.66 ± 0.060	40.4	0.54 ± 0.080	51.3
Expt^f		0.72 ± 0.040		0.50 ± 0.100

---

Fig. 5 shows the distribution of the S₁ → S₀ hopping times and time-dependent decay curves of S₀ and S₁ population for the cis-to-trans and trans-to-cis processes. Periodical decay pattern in cis → trans isomerization of azobenzene has been presented in a recent letter,⁶ while in b-Ab photoisomerization, we observe similar feature in the trans → cis process. As illustrated in Fig. 5a, in cis-to-trans isomerization, most of the S₁ → S₀ hops happen in the region of 30–100 fs with maximum at 45 fs and therefore give birth to a smooth exponential S₁ population decay curve. In trans-to-cis isomerization, most of the S₁ → S₀ hops occur between 15 and 55 fs with maxima at 18 and 34 fs, respectively. The double peak feature in hopping time distribution result in zigzag S₁ population decay curve, in which the slope decreases in the 20–30 fs region. The periodic interval of S₁ → S₀ hopping in trans-to-cis process is estimated to be about ∼16 fs. In each period, the trajectories visit conical region on the S₁ PES and hop to S₀ or remain in the excited state until re-enter the conical region in next period to encounter surface hopping. As seen in Fig. 5b, more hops take place during the second interval, indicates the reverse direction motion (from cis side to trans side) toward C7N2N1C1 moiety is preferred in photoisomerization started from trans b-Ab. Based on the surface hopping algorithm, the motion direction of a trajectory keeps unchanged before and after hop event, and consequently the quantum yields are largely determined by the trajectories' motion direction toward reactive coordinates at surface hopping points. In trans-to-cis process, more trajectories hop to S₀ in the second interval, when the molecule motion direction is the same as that of the first period in cis-to-trans process. Considering the higher quantum yields in the cis-to-trans isomerization, the trans-to-cis isomerization would yield more trans products and consequently result in relative lower quantum yields. The maximum of hopping time distribution in cis-to-trans process (45 fs) is about 20 fs larger than in the trans-to-cis process, (26 fs, average of two peak at 18 and 34 fs) which agrees with the experimental lifetime difference (∼70 fs for cis-to-trans and <50 fs for trans-to-cis) of Siewertsen et al.²⁰

Fig. 5 Time dependent S₁ and S₀ states' populations (top panel) and distributions of the S₁/S₀ hopping time (bottom panel) for (a) cis-to-trans and (b) trans-to-cis isomerization.

In the current nonadiabatic dynamics simulation, we observe trajectories for the cis-to-trans isomerization via clockwise and counterclockwise route and trans-to-cis process. Time evolution of potential energies of the S₁ and S₀ states and five important angles involve in the photoisomerization, two NNC bond angles, two CCNN and CNNC dihedral angles, for reactive trajectories are displayed in Fig. 6, while the corresponding nonreactive ones are given in Fig. S1.† Fig. 6a shows a typical reactive trajectory for trans → cis process, and the most evident conformation change is in C7N2N1C1 dihedral angle. Within the process, the energy of S₁ state increases in the first 5 fs and then decreases to meet the increasing S₀ energy at 17.5 fs and then the surface hopping takes place. Simultaneously, the C7N2N1C1 dihedral angle decrease from 145° to ∼100° along with the other internal coordinates increase moderately. After surface hopping, the five angles continue to evolve as prior to hop, therefore, this isomerization process can be attributed to nonsymmetrical torsion toward central eight-membered ring moiety.

Fig. 6 Time evolution of potential energies of the S₁ and S₀ states and important angles involve in the photoisomerization for (a) trans to cis, (b) cis to trans (counterclockwise reactive), (c) cis to trans (clockwise reactive).

The potential energy profile in Fig. 6b for the counterclockwise cis → trans isomerization is similar to trans → cis process, where the S₁–S₀ potential energy gap decreases quickly to populate the S₁ → S₀ hop at 41 fs. In the stage prior to surface hopping, the C7N2N1C1 and C3C1N1N2 dihedral angles quickly increase and decrease respectively, while, the C9C7N2N1 decreases to −120° and oscillates around it. The features of these motions demonstrate the cis → trans isomerization initiated by pedal-like motion towards the –N=N– bridge and then the other parts that linked to the central eight-membered ring start to transform to trans form or back to cis form. Moreover, the values of the two NNC angles separated, with one increasing to ∼135° and the other one decreasing to ∼115°, indicating the participation of NNC inversion in the photoisomerization. Our studies reveal that the counterclockwise cis → trans isomerization passes through the CI-cis conical regions, and therefore confirms the two CI model presented in previous investigations.

Fig. 6c shows a typical trajectory for clockwise cis → trans isomerization. Within the first 20 fs, the S₁–S₀ potential energy gap decreases fast and then oscillates around 40 kcal mol⁻¹ until 40 fs. The reason for the peculiar potential energy gap variation is that in the initial stage, the C7N2N1C1 dihedral angle oscillates in the region of [−20°, 10°] for about 40 fs to achieve the CNNC flip, and then the isomerization procedure starts with C7N2N1C1 quickly decreases to −100° and C9C7N2N1 increases to ∼20°. During this period of ∼30 fs, the S₁–S₀ potential energy gap decreases quickly again and surface hopping occurs at where the absolute value of those dihedral angles are similar to the ones of trans → cis isomerization shown in Fig. 6a. The reverse direction CNN change was also observed, which means the assistance of CNN inversion is also required in clockwise cis → trans isomerization. With the ∼40 fs trap in the semi-cis form, the clockwise cis → trans channel should possess a rather longer lifetime than the trajectories follow the other two isomerization pathway. Generally, the clockwise cis → trans isomerization process is also dominated by the nonsymmetrical torsion motion with respect to the azo and ethylenic bridges, while the additional CNN inversion is a promotion.

4. Concluding remarks

By employing the recently developed analytical switching probabilities calculation algorithm based on the Zhu–Nakamura theory, we have run trajectory-based nonadiabatic on-the-fly molecular dynamics simulations on the bridged-azobenzene photoisomerization for both cis-to-trans (400 sampling trajectories) and trans-to-cis (550 sampling trajectories) at the SA2-CASSCF(6,5)/6-31G quantum chemistry level. With all the simulated trajectories incorporated, the estimated quantum yields and lifetimes converge to 0.57 (0.44) and 64 (35.5) fs for the cis-to-trans (trans-to-cis) photoisomerization, respectively, which agree with transient absorption experiment observations and previous dynamics simulation investigations. The enhanced photoresponsibilities of trans-b-Ab, such as the higher quantum yields and shorter excited state lifetimes than the parent azobenzene are originated from the steeper PES in the S₁ state, shorter isomerization pathway and favorable non-planar initial structure. Furthermore, the most attractive feature in the azobenzene derivatives is the unusual quantum yields drop upon the ππ* excitation,^40–42 which will be the topic of our future investigations.

The cis ↔ trans photoisomerization of b-Ab was proposed to encounter surface hopping through one or two CI regions via twist motion toward central 8-membered ring. The current studies indicated that the cis-to-trans and trans-to-cis isomerization proceed via different conical regions with respective motion modes toward central N=N related moiety. Which agree with the previous investigations of Carstensen et al.²² and Gao et al.,^27,28 by using the CASSCF method. However, it is the main difference between the previous studies^21,23–26 and ours. The other simulation investigations which employing the DFT based methods prefer the model that the cis-to-trans and trans-to-cis isomerization are induced by the same conical region^21,23–26 with similar relaxation process in S₁ state as those for two CI model. In the current work, the two types of optimized CIs between S₀ and S₁ states are in cis-like and trans-like form, which possess rather different twisting angles within the central 8-membered rings except for the CNNC dihedral angles. Most of the trajectories starting from cis-b-Ab hop to S₀ via the CI-cis region, and this isomerization process was initiated by the pedal-like twist of –N=N– moiety and followed by the reorientation of the other neighboring dihedral angles. While, the trans-to-cis process is induced by the nonsymmetrical torsion motion toward the central eight-membered ring and encountered surface hopping through CI-trans region. We also observed the switch between these two isomerization pathways, which is accompanied with the chirality conversion staring from cis-b-Ab. This pathway confirms the existence of chirality change toward CNNC dihedral angle in photoisomerization of cis-b-Ab as those reported in parent cis-Ab.⁶ The trajectory follows this pathway stay in the Franck–Condon region for ∼40 fs to populate the CNNC flip and then converts to the torsion channel directing to the CI-trans region with minus CNNC value. Obviously, the lifetime in excited state for this isomerization channel would be prolonged. In the current simulation, the ratio of the chirality conversion channel is much smaller than that for the parent azobenzene (0.10–0.22).^6,8 One possible reason is that all trajectories in the current simulation start at the same geometry with different initial velocities. We will perform investigation with better initial condition which considering the influence of zero-point vibration energy in the near future. Moreover, in trans-to-cis process, the periodic S₁ population decay character was observed with an interval of ∼16 fs. The second period, in which the molecule moves along cis-to-trans direction, occupies more hopping events and therefore results in reduced quantum yields.

Acknowledgements

H. H. thanks support from the National Science Foundation of China under grant no. 21103136. Y. L. thanks support from the National Science Foundation of China under grant no. 21473134. C. Z. thanks support from the Ministry of Science and Technology of the Republic of China under grant no. 103-2113-M-009-007-MY3 and MOE-ATU project of the National Chiao Tung University. L. Y. thanks support from the Postdoctoral Fellowship by the Ministry of Science Technology of the Republic of China under grant no. 103-2811-M-009-048 and Scientific Foundation of Northwest University under grant no. 338050097.

References

1. H. M. D. Bandara and S. C. Burdette, Chem. Soc. Rev., 2012, 41, 1809–1825.
2. H. Dürr and H. Bouas-Laurent, Photochromism: molecules and systems, Gulf Professional Publishing, 2003.
3. B. L. Feringa, Molecular Switches, Wiley-VCH Verlag GmbH, Weinheim, FRG, 2001.
4. L. Yu, C. Xu and C. Zhu, Phys. Chem. Chem. Phys., 2015, 17, 17646–17660.
5. L. Yu, C. Xu, Y. Lei, C. Zhu and Z. Wen, Phys. Chem. Chem. Phys., 2014, 16, 25883–25895.
6. O. Weingart, Z. Lan, A. Koslowski and W. Thiel, J. Phys. Chem. Lett., 2011, 2, 1506–1509.
7. M. Pederzoli, J. Pittner, M. Barbatti and H. Lischka, J. Phys. Chem. A, 2011, 115, 11136–11143.
8. Y. Ootani, K. Satoh, A. Nakayama, T. Noro and T. Taketsugu, J. Chem. Phys., 2009, 131, 194306.
9. J. Shao, Y. Lei, Z. Wen, Y. Dou and Z. Wang, J. Chem. Phys., 2008, 129, 164111.
10. M. L. Tiago, S. Ismail-Beigi and S. G. Louie, J. Chem. Phys., 2005, 122, 094311.
11. T. Schultz, J. Quenneville, B. Levine, A. Toniolo, T. J. Martinez, S. Lochbrunner, M. Schmitt, J. P. Shaffer, M. Z. Zgierski and A. Stolow, J. Am. Chem. Soc., 2003, 125, 8098–8099.
12. C. R. Crecca and A. E. Roitberg, J. Phys. Chem. A, 2006, 110, 8188–8203.
13. T. Pancur, F. Renth, F. Temps, B. Harbaum, A. Krüger, R. Herges and C. Näther, Phys. Chem. Chem. Phys., 2005, 7, 1985–1989.
14. Y.-C. Lu, E. W.-G. Diao and H. Rau, J. Phys. Chem. A, 2005, 109, 2090–2099.
15. H. Rau and E. Lüddecke, J. Am. Chem. Soc., 1982, 104, 1616–1620.
16. F. Puntoriero, P. Ceroni, V. Balzani, G. Bergamini and F. Vögtle, J. Am. Chem. Soc., 2007, 129, 10714–10719.
17. W. W. Paudler and A. G. Zeiler, J. Org. Chem., 1969, 34, 3237–3239.
18. E. Tauer and R. Machinek, Liebigs Ann., 1996, 1996, 1213–1216.
19. R. Siewertsen, H. Neumann, B. Buchheim-Stehn, R. Herges, C. Näther, F. Renth and F. Temps, J. Am. Chem. Soc., 2009, 131, 15594–15595.
20. R. Siewertsen, J. B. Schönborn, B. Hartke, F. Renth and F. Temps, Phys. Chem. Chem. Phys., 2011, 13, 1054–1063.
21. M. Böckmann, N. L. Doltsinis and D. Marx, Angew. Chem., 2010, 49, 3382–3384.
22. O. Carstensen, J. Sielk, J. B. Schönborn, G. Granucci and B. Hartke, J. Chem. Phys., 2010, 133, 124305.
23. L. Liu, S. Yuan, W. H. Fang and Y. Zhang, J. Phys. Chem. A, 2011, 115, 10027–10034.
24. C. Jiang, R. Xie, F. Li and R. E. Allen, J. Phys. Chem. A, 2011, 115, 244–249.
25. C. Jiang, R. Xie, F. Li and R. E. Allen, Chem. Phys. Lett., 2012, 521, 107–112.
26. M. Böckmann, N. L. Doltsinis and D. Marx, J. Chem. Phys., 2012, 137, 22A505.
27. A. H. Gao, B. Li, P. Y. Zhang and K. L. Han, J. Chem. Phys., 2012, 137, 204305.
28. A. H. Gao, B. Li, P. Y. Zhang and J. Liu, Comput. Theor. Chem., 2014, 1031, 13–21.
29. M. Böckmann, N. L. Doltsinis and D. Marx, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 036101.
30. M. Böckmann, N. L. Doltsinis and D. Marx, J. Phys. Chem. A, 2010, 114, 745–754.
31. Y. Lei, H. Wu, X. Zheng, G. Zhai and C. Zhu, J. Photochem. Photobiol., A, 2016, 317, 39–49.
32. L. Verlet, Phys. Rev., 1967, 159, 98–103.
33. C. Zhu and H. Nakamura, J. Chem. Phys., 1992, 97, 8497.
34. C. Zhu and H. Nakamura, J. Chem. Phys., 1993, 98, 6208.
35. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053.
36. P. J. Knowles and H.-J. Werner, Chem. Phys. Lett., 1985, 115, 259.
37. G. Knizia, F. R. Manby and M. Schütz, et al., MOLPRO, a package of ab initio programs, Version 2009.1, ed. H.-J. Werner and P. J. Knowles.
38. G. KarlstrÖm, R. Lindh, P.-Å. Malmqvist, B. O. Roos, U. Ryde, V. Veryazov, P. O. Widmark, M. Cossi, B. Schimmelpfennig, P. Neogrady and L. Seijo, Comput. Mater. Sci., 2003, 28, 222–239.
39. Y. Wang, H. Han, Y. Lei, B. Suo, H. Zhu, Q. Song and Z. Wen, J. Chem. Phys., 2014, 141, 164114.
40. H. Satzger, C. Root and M. Braun, J. Phys. Chem. A, 2004, 108, 6265–6271.
41. M. Quick, A. L. Dobryakov, M. Gerecke, C. Richter, F. Berndt, I. N. Ioffe, A. A. Granovsky, R. Mahrwald, N. P. Ernsting and S. A. Kovalenko, J. Phys. Chem. B, 2014, 118, 8756–8771.
42. E. M. M. Tan, S. Amirjalayer, S. Smolarek, A. Vdovin, F. Zerbetto and W. J. Buma, Nat. Commun., 2015, 6, 5860.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra03788g

---