Research Article  Open Access
Constraint Trajectory SurfaceHopping Molecular Dynamics Simulation of the Photoisomerization of Stilbene
Abstract
Combining trajectory surface hopping (TSH) method with constraint molecular dynamics, we have extended TSH method from full to flexible dimensional potential energy surfaces. Classical trajectories are carried out in Cartesian coordinates with constraints in internal coordinates, while nonadiabatic switching probabilities are calculated separately in free internal coordinates by LandauZener and ZhuNakamura formulas along the seam. Twodimensional potential energy surfaces of ground and excited states are constructed analytically in terms of torsion angle and one dihedral angle around the central ethylenic C=C bond, and the other internal coordinates are all fixed at configuration of the conical intersection. At this conical intersection, the branching ratio from the present simulation is 48 : 52 (33 : 67) initially starting from trans(cis)Stilbene in comparison with experimental value 50 : 50. Quantum yield for transtocis isomerization is estimated as 49% in very good agreement with experimental value of 55%, while quantum yield for cistotrans isomerization is estimated as 47% in comparison with experimental value of 35%.
1. Introduction
Nonadiabatic dynamic based on onthefly trajectory surface hopping (TSH) approach is a powerful tool for investigation of the photochemical and photophysical processes involving electronically excited states with conical intersections. Trajectories are run on onthefly electronically adiabatic potential energy surfaces while nonadiabatic transitions are treated by mixing quantum/classical or semiclassical methods. In year 1971, the TSH method was already utilized for studying the nonadiabatic reactions by Tully and Pkeston [1] in which nonadiabatic transition probabilities were calculated by the LandauZener (LZ) formula along the seam. Later in year 1990, Tully [2] proposed the fewest switches (FS) algorithm in which nonadiabatic transition probabilities were calculated by the timedependent firstorder coupled equations along running trajectory. By extending LZ theory, Zhu and Nakamura (ZN) [3] developed the better nonadiabatic transition formula that was applied to the nonadiabatic reactions [4]. Nonadiabatic transition probability calculated by ZN formula is generally larger than that calculated by LZ formula. Extensive studies in comparison of LZ formula with Tully’s FS algorithm were carried out for photochemical ring opening in oxirane [5], in which FS transitions were taking place in a wide range of the conic zone while LZ transitions occurred locally around the conic zone. The other theoretical approaches such as Liouville dynamics [6], Bohmian dynamics [7], path integrals [8, 9], and multiple spawning [10] have been proposed. Nonadiabatic dynamic based on Tully’s FS algorithm has been widely applied for photochemistry of large molecular and biomolecular systems [11–19] and it has been well documented in the recent review paper [20].
Nonadiabatic transition probabilities calculated from Tully’s FS method require calculation of nonadiabatic coupling vectors while the transition probabilities calculated from LZ and ZN formulas require information only from two adiabatic potential energy surfaces around the crossing seam. The LZ and ZN methods are less expensive than the FS method for TSH. A nuclear motion can be treated separately from electronically nonadiabatic transitions in LZ or ZN method so that a step size of trajectory can be much longer than that of FS method where nuclear motion has to be integrated coupled with electronic motion. Trajectory switching in LZ or ZN method occurs at place very close to crossing seam while trajectory switching in FS method can happen at place far from crossing seam. As a result, number of onthefly trajectories can be quite demanding for convergent accuracy in FS method. However, LZ and ZN methods cannot take into account nonadiabatic transitions from noncrossing case. ZN method was successfully applied to nonadiabatic dynamics of two models of protonated Schiff base retinal up to 100 onthefly trajectories [21] and photoisomerization of bridged azobenzene [22, 23]. If nonadiabatic transitions are dominated by the crossing seam, LZ and ZN formulas can be applied.
It is a quite successful approach that the potential energies, energy gradients, and nonadiabatic coupling vectors are calculated onthefly for trajectories with full degree of freedoms for large systems. However, it can run limited number of onthefly trajectories with fulldimensional calculation if high level ab. initio quantum chemistry calculation is performed for potential energy surfaces and nonadiabatic coupling vectors. Although low level ab. initio calculation can increase number of running trajectories, it has limited accuracy for potential energy surfaces. An onthefly trajectory with reduced dimension for large systems provides an alternative way for studying nonadiabatic molecular dynamics. Especially when nonadiabatic transitions occur locally involving few degrees of freedom, trajectories can be run on reduceddimensional either onthefly or analytical potential energy surfaces depending on how many degrees of freedom are taken into account. As reduceddimensional degrees of freedom are usually happening on internal coordinates such as bond lengths, bond angles, and dihedral angles while trajectories are running on Cartesian coordinate systems, we need coordinate transformation between internal and Cartesian coordinates for solving constrained classical Hamiltonian equations. Thus, we need to solve coordinates, momentums, and the Lagrange multipliers simultaneously for which various numerical methods are available [24–30].
The purpose of the present study is to combine TSH method with constraint classical Hamiltonian equations to study nonadiabatic molecular dynamics. LZ or ZN methods are mostly suitable incorporated with TSH for constraint molecular dynamics because nonadiabatic transitions can be treated separately with nuclear motion [31, 32]. The present method provides an alternative way in onthefly TSH category for simulating large system with increasing number of trajectories and probes simulation in the nonadiabatic transition zone. It is useful complementary to fulldimensional onthefly TSH method. In order to apply the present method, we must investigate detail structure and degrees of freedom of conical intersections before deciding which degrees of freedom can be constrained for interesting photochemical processes. Nonadiabatic dynamics of Stilbene has only been studied by employing the ZN method in a short report [33] and the TSH method on FS algorithm within timedependent DFT local orbital framework [34, 35]. In contrast, there are many publications on onthefly TSH method for nonadiabatic dynamics of azobenzene [22, 23, 34–38]. Therefore, we apply the present method to study photoisomerization of Stilbene in detail.
The isomerization of cis and transStilbene has been extensively studied for more than forty years and the quantum yields of the photoisomerization were measured experimentally [39, 40]. Three mechanisms were characterized for the isomerization; they are the conventional onebond flip (OBF) [41], HulaTwist (HT) [42], and the aborted HT mechanisms [43]. The OBF mechanism involves a 180° rotation around the central ethylenic C=C bond. The HT mechanism is considered as the concerted torsion around the adjacent vinylphenyl bond and a remarkable bend of in plane C=C−C angle in accompanying the central ethylenic C=C bond rotation. The aborted HT mechanism is explained as that the rotation of one of two phenyl rings in Stilbene is aborted and turned back. The cistotrans and transtocis isomerizations are taken place by going through pathway of conical intersections which are named as OBFCI [44, 45] and HTCI [46], respectively.
The Stilbene has cis and transisomers that can be photoinduced and transferred from one to another in accompanying electronically nonadiabatic transition among various electronic states. Traditionally, the isomerization of cis and transconformations has been investigated by computing the ground state and the lowest excited state . The dominant excitation from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) corresponds to the transition which was confirmed to be a bright state [45]. In the recent years, the photoisomerization involving the conical intersections plus fluorescence spectra has been extensively investigated by both experimentalists and theoreticians [34, 35, 44–65]. The OBF mechanism describes isomerization mainly depending on phenyl rotation in which two phenyl rings twist around the central ethylenic C=C bond, but some other torsions might be also included [44–51]. The conical intersection OBFCI associated with the OBF mechanism was calculated to be the lowest in energy, so that the photoisomerization reactions were predicted to be in favor of this mechanism [45]. The further dynamical simulations were done by calculating evolution of HOMO and LUMO orbitals along with the nuclear motions including all nuclear coordinates [52, 53]. The vibrational normalmode analysis with specified PES involving the steep route indicated that the vibrational mode with frequency 240 cm^{−1} also strengthens the OBF mechanism [54]. Fuß and coworkers [43, 55, 56] suggested that the HT and the aborted HT mechanisms were deduced from the systematic features of the cis and transphotoisomerizations of nonpolar conjugated molecules and these were investigated by the experimental measurements as well [58–60].
Starting from the FranckCondon (FC) regions, PESs along the reaction coordinate (curved onedimension) were also investigated [54, 62, 65]. It was analyzed that onedimensional PESs with respect to the main twist of phenyl groups around the central C=C bond reflect the main dynamical effect for isomerization at OBFCI as mentioned above. The two torsion angles in terms phenyl rotation and the pyramidalization coordinate were utilized for constructing twodimensional PESs around the OBFCI [45]. These studies told that at the OBFCI the twodimensional PESs can well describe transcis nonadiabatic dynamics for Stilbene. The branching ratio is basically determined by potential energy surfaces near the OBFCI. The one torsion angle around the central ethylenic C=C bond plus one dihedral angle was a suitable choice for isomerization via the OBFCI (called as D1 and D2 in the present paper as shown in Figure 1). The potential energy surfaces of groundstate and the first excited state are constructed in terms of the D1 and D2 angles, and the other internal coordinates are fixed at the OBFCI configuration.
The rest of this paper is organized as follows. Section 2 first summarizes constrained molecular dynamics method and constructs twodimensional potential energy surfaces analytically for Stilbene. Trajectory surface hopping with onepassage nonadiabatic transition probability is carried out by LZ and ZN formula. Section 3 presents the implementations of the methods mentioned in Section 2, and results and discussions about branching ratio of isomerization and picture of evolution trajectory are also given therein. Concluding remarks are mentioned in Section 4 followed by three Appendices A, B, and C in which the detailed constraint relations, nonadiabatic transition probability, and hopping scheme are presented.
2. Constraint TSH Molecular Dynamics with LZ and ZN Formulas for Nonadiabatic Transitions
We first introduce constraint molecular dynamics for system with constraints, and then we construct twodimensional potential energy surfaces (PESs) for Stilbene molecule. Finally, we discuss how LZ and ZN formulas can be applied to TSH for calculating nonadiabatic transition probability along the seam. Of course, we can do onthefly nonadiabatic dynamics with four or more dimensional potential energy surfaces. However, it was well studied that twodimensional PESs are suitable for describing isomerization of Stilbene at OBFCI. Moreover, for two dimensions we can construct analytical PESs easily and can demonstrate how trajectories run around the seam.
2.1. Constraint Molecular Dynamics
We review constraint molecular dynamics for system with constraints in which the canonical Hamiltonian equations can be written as [24–30] where , , and are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respectively. The are constraint equations determined by internal coordinates that are fixed at the certain structures of nuclear configuration all the time. The and in (2) stand for unconstraint and constraint force acting on the corresponding atom, respectively. The Hamiltonian with Lagrange multipliers is given as where are to be determined. As in (3) are zero all the time, the energy conversation is satisfied.
The Lagrange multipliers are solved iteratively. With the given initial Lagrange multipliers, we solve and by the numerical integration of (1) and (2), and then the is inserted into the constraint equations (including bond lengths, bond angles, and dihedral angles) to get the Lagrange multipliers. This iterative process is eventually convergent. The technique to implement this iterative procedure is starting from the classical (Newtonian) equations [27]: Integrating both sides of (4) twice in time yields the constraint Cartesian coordinates at the time : where is solved by unconstrained force. The above can be inserted to every type of the constraint equations to derive nonlinear equations with respect to the corresponding Lagrange multipliers. Although there are a number of algorithms to solve , they only differ on how to solve these nonlinear equations. The SETTLE algorithm [28, 29] obtained the solutions of these nonlinear equations analytically for constraints, but it cannot be extended to the large system. The originally simple SHAKE algorithm [24, 25] was developed to solve bond length constraints which are converged linearly. Later advanced MSHAKE Newton method [24, 25, 30] and quasiNewton method [30] were straightforward and fast to solve the nonlinear equations. All of these methods involved the calculations of the Jacobian determinant about the constraint equations : where are then updated repeatedly by solving the system of linear equations iteratively: until the in which is small number with the prescribed tolerance of the constraints. The Jacobian determinant for the quasiNewton method is only computed once in the first iteration with the linear convergence, while the Jacobian determinant Newton’s method presents quadratic convergence. In the present work, the direct MSHAKE algorithm is employed and this method is not only to keep the quadratic convergence but also to simplify the derivation and implementation of the nonlinear equations. Detailed derivation of Jacobian determinant is given in Appendix A.
2.2. TwoDimensional Potential Energy Surfaces for Stilbene
Quenneville and Martínez [45] have done extensive calculation for the critical conical intersection OBFCI and potential energy surfaces in terms of two internal angles: phenyl rotation and the pyramidalization coordinates. They have done detailed comparison for the results calculated by complete active space selfconsistent field (CASSCF) and multireference perturbation theory (CASPT2) methods with basis sets of 631G., 6, and 63. Following their conclusion, we can safely utilize stateaveraged CASSCF method abbreviated as SANCAS(n/m), where N refers to the number of states included in the average, while N(=2), n(=2), and m(=2) are the number of active electrons and orbitals, respectively. This means that just and orbitals corresponding to HOMO and LUMO are involved in CASSCF calculation. Molcas 7.5 [66] program packages were employed for all calculations.
We have first used SA2CAS(2/2)/631G method to reproduce geometries of conical intersection OBFCI, the local minima on and potential energy surfaces given by [45]. Then, we concentrate on configuration of conical intersection at which the torsion angle D1 (C11C8C1C2) and dihedral angle D2 (H10C8C1H9) are −146.6° and −60.0°, respectively, as shown in Figure 1. We have carried out a preliminary test calculation for twodimensional potential energy surfaces with variables D1 and D2 (all bond lengths, all bond angles, and the other dihedral angles are fixed at OBFCI configuration). We found that and are better variables to describe the isomerization. If we set up interesting region for potential energy surface with energy not higher than 5 eV above OBFCI (where it is considered as zero point of potential energy), DD2 can be chosen as in [−20°, 90°] and DD1 is chosen as in [0°, −180°] that is sufficient to describe the present isomerization process. We use equal spacing 5° for grid points: DD1 is (37 grid points) and DD2 is (23 grid points). Total configurations are 851 plus OBFCI which all are computed by SA2CAS(2/2)/631G method, and then all ab. initio data are fitted into analytical function of potential energy surfaces by the least square method from Matlab 7.5 package [67]. Finally, we have adiabatic potential energy surfaces in the analytical form: where and , 3.3° and 3.3° are angles at OBFCI. The fitting parameters , , and to are all given in Table 1. Potential energy surfaces of and calculated from (8) are plotted in Figure 2 with zero point of potential energy at in a unit of electron volt. We computed the mean absolute error between calculated results from (8) and numerical data computed from method SA2CAS(2/2)/631G for about 100 randomly picked up configuration points and we found that this error is about 2.4 kcal/mol for the both potential energy surfaces, and less than 1.0 kcal/mol around OBFCI region. From analytical PES in (8), we have estimated local transisomer at 64.3° and 6.7° with energy 1.72 eV below OBFCI, and local cisisomer at 9.2° and 0.9° with energy 1.28 eV below OBFCI. Furthermore, the vertical excitation energies have been calculated to be 4.11 eV and 4.45 eV at local trans and cisminima, respectively, which are accidently in consistence with the experimental measurement 4.13 eV and 4.59 eV [68]. We confirmed that a global minimum of state is just right at conical intersection for the present twodimensional potential energy surface.

2.3. TSH Scheme for LZ and ZN Methods
Section 2.1 implements classical trajectory simulation on a single potential energy within the Cartesian coordinate system plus constraints in internal coordinate systems. This not only simplified numerical calculation for integration of classical trajectory in molecular Hamiltonian but also avoided searching complicated form of kinetic operators in internal coordinates. However, nonadiabatic transition between two PESs is completely quantum effects, and it requires explicit form of the kinetic operators to calculate nonadiabatic coupling vector between two PESs. We focus on the torsion angle D1 and dihedral angles D2 (transferred into DD1 and DD2 for PESs depicted in Figure 2) and also present the following regular form of kinetic operators explained in Scheme 1, with where the dihedral angle D1 represents the relative torsion motion between two phenyl rings and the D2 represents the relative torsion motion between two hydrogen atoms. We assume that atoms 2 and 3 are fixed and two phenyl rings and two hydrogen atoms rotate on the circles around axis as shown in Scheme 1. For instance, we can think that and are rotation angles around axis for two phenyl rings, respectively. Thus, the moment of inertia is the reduced inertia from two phenyl rings where ( is atom mass in the one phenyl ring and is distance from atom to axis) and is estimated from the other phenyl ring. The same analysis can be applied for the moment of inertia with two hydrogen atoms: where ( is hydrogen atom mass and is distance from hydrogen atom to axis) and is estimated from the other hydrogen atom.
Onepassage semiclassical nonadiabatic transition probability (LZ and ZN formulas) can be calculated just from two adiabatic potential energy surfaces along the seam if we have the regular kinetic form of (9): where calculation of two parameters and is given in Appendix B.
We can simply look at the present twodimensional PESs in Figure 2 to find that the seam line is appearing at the . Seam line is defined as the trace of local minimum gap between two adiabatic PESs [69]. However, in the present work we propose to calculate the local maxima of the effective coupling parameter as introduced in Appendix C. Once the seam line is found, the nonadiabatic coupling vector is assumed to be perpendicular to the seam line and it appears in DD1 direction in the present case. Therefore, the effective collision energy parameter is calculated along the DD1 direction and it means that only angular momentum related to DD1 is changed during the trajectory surface hopping. After trajectory hopping from one potential energy surface to another, we need to redistribute this angular momentum change into six atoms (see Scheme 1) in the Cartesian coordinate systems and the detailed procedure is given in Appendix C.
The calculation shows that the seam line appears at with negligible dependence of DD2 angle. The effective coupling parameter along the seam line is calculated and depicted in Figure 3. Figure 3 shows that the diabatic region (where the onepassage transition probability is almost unity) is located in between 41° and 49° for angle DD2. The other region of DD2 belongs to effective nonadiabatic transition zone, especially for DD2 in between [9°, 39°] and [57°, 80°] where it is a strong nonadiabatic transition zone with .
3. Results and Discussions
Experimental measurement starts from global cisminimum and transminimum in which photoisomerization processes involve dynamics on fulldimensional potential energy surfaces. We construct the reduced twodimensional (2D) potential energy surfaces at OBFCI, and potential energy surface shows clear local cisminimum and transminimum separated by the seam line. We assume that trajectory starts from global cis (trans)minimum with photoexcitation from to PES and then reaches the certain area of the present local cis (trans)minimum on the 2D PES, from where trajectories start isomerization. Then, question is how to determine such area for both local cis and transisomers. Let us start with the test local cisarea and local transarea in rectangular region of [−19.2° < DD1 < −59.2°, 30.9° < DD2 < 70.9°] and [−180.0° < DD1 < −144.3°, 16.7° < DD2 < 56.7°], respectively, for both initial and final conditions of trajectories. Total energy for cistotrans isomerization is defined by the excitation energy 4.45 eV from to state estimated at local cisminimum (DD1 = −39.2°, DD2 = 50.9°). The initial conditions of DD1 and DD2 for trajectories are randomly picked up in the local cisarea vertically excited from to state; if this trajectory has vertical excitation energy which is larger than 4.45 eV, it is abandoned, and if this trajectory has vertical excitation energy which is smaller than 4.45 eV, it can be proceeded. Then, we specify initial kinetic energy for proceeding trajectory with Scheme 1 by choosing atoms 1, 4, 5, and 6 with nonzero velocities and all atoms in molecule with zero velocities. The proceeding trajectory which starts from local cisarea on state can have three ending ways; one is that it enters local cisarea on (called cistocis nonreactive trajectory), another is that it enters localtrans area on (called cistotrans reactive trajectory), and the third is that it enters outside of 2D PES boundary region depicted in Figure 2; for instance, DD1 is smaller than −180° or larger than 0° (called unreactive trajectory). An unreactive trajectory is due to that 2D PES has periodic properties in terms of DD1 and DD2 angles and we do count it as nonreactive trajectory.
Now we present how to distribute a given kinetic energy to atoms 1, 4, 5, and 6 with their initial velocities analyzed in the following. As shown in Scheme 1, the zcomponent of velocity is assumed to be zero so that for a given kinetic energy , we have the following relations: where is random number in the range of and () is mass of carbon (hydrogen) atom. The and components of velocities in (14) must satisfy because of constraints of bond lengths and bond angles. Moreover, we choose the angular moment along axis for atoms 1 and 4 and for atoms 5 and 6 as follows: This means that angular momentum along axis for the twist of atoms 1 and 4 is cancelled out from each other, so does for atoms 5 and 6. Equations (14), (15), and (16) all together have the eight related equations for the eight velocity components to be determined. However, the above procedure slightly violates conservation of total energy because two phenyl rings associated with atoms 1 and 4 in Scheme 1 can have induced velocities through the constraints for entire molecule. We can resolve this problem by choosing smaller kinetic energy in (14) than the previously specified one, and this can be obtained by subtracting induced kinetic energy from two phenyl rings. Of course, this has to be done iteratively with the constraint Hamiltonian equations (1) and (2) in the first step of classical trajectory integration.
We can do the same as mentioned above for transtocis isomerization, and the total energy is defined by the excitation energy 4.11 eV from to state estimated at local transminimum (DD1 = −164.3°, DD2 = 36.7°). The initial conditions of DD1 and DD2 for trajectories are randomly picked up in the local transarea vertically excited from to state. Three types of ending trajectories are also collected as mentioned above (transtotrans nonreactive, transtocis reactive, and unreactive trajectories which is classified as nonreactive). All numerical calculations are performed within the atomic unit. We use our own code for solving the constraint Hamiltonian equations and integration for classical trajectory in combination with the free MINPACK code [70] which contains Powell’s Dog Leg method [71] for solving the nonlinear equation solution.
For given initial coordinates and corresponding conjugate moments, we integrate (1) and (2) with initial the Lagrange multipliers . Then, outcome of coordinates is inserted into the Jacobian determinant in (7) to have the new that are inserted into (1) and (2) again. We do this iteratively until all Lagrange multipliers converge together with the update coordinates and momentums. For the second step of integration, the Lagrange multipliers are set up as the previous values, and in this way, we propagate trajectory until it ends with use of the fourthorder RungeKutta method (RK4) [72] by equal time step size 2.07 a.u. (about 0.05 fs).
Classical trajectory is propagated on a single adiabatic potential energy surface all the time. However, once trajectory is across the seam line, we compute effective coupling and effective collision energy according to the method given in Appendix C. Then, we can evaluate onepassage nonadiabatic transition probability by (12) and (13). We obtain the seam line by calculating local minima of effective coupling and found that they all are located at avoided crossing points. Comparing with random number in between 0 and 1, we decide to hop trajectory from one adiabatic potential energy surface to another if is larger than the random number; otherwise, it stays in the original potential energy surface. As trajectory hops to the other PES vertically, and it means that all coordinates are unchanged, but the moments are changed with method introduced in Appendix C.
3.1. Branching Ratios of CistoTrans and TranstoCis Isomerizations
By giving equal weight to all the sampling trajectories, we compute the branching ratios with initial condition starting from both local cisarea and local transarea. We carry out trajectory surface hopping simulation with number of trajectories 100, 500, 1000, and 2000, respectively, in order to check convergence of the branching ratios. The results in Table 2 show that the results are basically converged at 500 trajectories. The branching ratio for transtocis isomerization is 48 : 52 (experimental value is 50 : 50) in the present simulation. By considering a potential barrier in FrankCondon region of global transminimum which is taking 5% trajectories down to global transminimum by fluorescence [59, 61], we can estimate the quantum fields for transtocis isomerization as 0.95 × 52% = 49% which is in very good agreement with the experiment measurement 55.0%. The branching ratio for cistotrans isomerization is 33 : 67 (experimental value is 50 : 50) in the present simulation. There are two pathways when Stilbene is photoexcited from the global cisisomer to state (see [59, 61] and references therein); one pathway leads to the OBFCI with branching ratio 70% and another goes to the other conical intersection (branching ratio 30%) which is for side reaction of a ring closing to form DHP. If we take this branching ratio into consideration, we can estimate the quantum yield for cistotrans isomerization as 0.7 × 67% = 47% which is larger than experimental measurement with 35%. The present simulation shows isomerization in favor of reactive trajectories regardless of initial starting from local transarea or local cisarea and this means that branching ratio is not 50 : 50 at OBFCI.
 
^{a}Total energy of each classical trajectory is set up to the vertical excitation energy 4.45 eV at cisarea and 4.11 eV at transarea. ^{ b}References [59, 61]. 
Potential energy surface of state as shown in Figure 2 is down to hill much faster in DD1 direction than in DD2 direction, and thus in the average, speed of trajectory is much faster in DD1 direction than in DD2 direction. This means that effective collision energy is well above the seam and this makes trajectory have a big chance to hop down state in the first time to cross the seam line. This is why the branching ratio is in favor to reactive isomerization and both LZ and ZN formulas present the same results in Table 2. Meanwhile, initial vertical excitation energy is 3.17 eV above OBFCI at cisarea and is 2.39 eV above OBFCI at transarea, so that effective collision energy is higher for cistotrans than for transtocis. This explains why branching ratio is more in favor of reactive cistotrans than transtocis. We know that the conical intersection HulaTwist (HTCI) is higher in energy than in OBFCI, and these two CIs might be in close configuration. Therefore, we suspect that there is a big chance for cistotrans trajectory to access HTCI with high vertical excitation energy, but this is not the case for transtocis trajectory with low vertical excitation energy.
We have also tested how branching ratios are sensitive to choice of size of trans and cisarea, and we found that they are not so sensitive to the choice of the area size unless we define very small area for cisarea and transarea. This smallness is not reasonable because trajectory takes quite long time to reach local cis (trans)area from global cis (trans)isomer; it should have a chance to distribute in relatively large area as we tested above. We check sensitivity of branching ratios with respect to the slight change of seam line with DD1 = −103.25 and DD1 = −103.35, and we found that the results in Table 2 still stand. This is because effective coupling and collision energy parameters ( and ) are quite robust in the present formulation.
3.2. Detailed Analysis of Typical Classical Trajectories
It is interesting to look at a typical trajectory to analyze how isomerization is taking place. Figure 4 shows how a reactive trajectory for cistotrans isomerization varies with time by starting at local cisarea with photoexcitation energy 4.42 eV and total evolution time is 115 fs. Four snapshots of structure varying with time are shown in Figure 4(a); it can be seen that potential energy steeply decreases on state and increases on state with crossing seam line around 52 fs first and then crossing again shortly, but no hopping occurs due to the fact that nonadiabatic transition probability is smaller than random number there. Trajectory hopping from to state occurs at 80 fs where DD1 and DD2 are almost at conical configuration (where transition probability is almost a unit). After hopping, the trajectory ends at local transarea at 115 fs. Figure 4(b) shows that both DD1 and DD2 angles change with oscillation against time; DD1 drops considerably while DD2 keeps almost constant oscillation. It would be interesting to see original torsion angle D1 and dihedral angle D2 varying with time as shown in Figure 4(c). Figure 4(c) shows that the torsion angle D1 decreases monotonically from −100.0° (at local cisarea) to −185° (at local transarea), and this is exactly corresponding to what it is called as onebond flip process around the central C=C bond. Dihedral angle D2 decreases from the beginning at 7.5° to the end at −95° with threecircle oscillations. This is because light hydrogen atoms move much faster than heavy phenyl rings. The oscillation of D2 angle can be explained by its coupling with D1, when the terminal atoms of D2 are moving close to the terminal atoms of D1, leading to a large steric hindrance, and then far away from them periodically. These results provide the evidence that the electronically nonadiabatic transition would be dominated by the phenyl ring twist around the central double bond.
(a)
(b)
(c)
Figure 5 shows how a reactive trajectory for transtocis isomerization varies with time by starting at local transarea with photoexcitation energy 3.98 eV and evolution time is 133 fs. Potential energies on and states vary with time with four selected snapshots of structure as shown in Figure 5(a) for a reactive transtocis trajectory. There are threetime surface hopping from one potential energy surface to another during evolution of trajectory; the first hopping from to state happened at 50 fs, the second hopping from to state at 71 fs, and the last hopping from to state at 88 fs. After that, the trajectory is propagated on state until its end at 133 fs. Figure 5(b) shows that both DD1 and DD2 angles change with oscillation against time; DD1 rises considerably while DD2 keeps almost constant oscillation. Figure 5(c) shows that the torsion angle D1 increases monotonically from −188° (at local transarea) to −97° (at local cisarea), and this shows onebond flip picture around the central C=C bond again. Dihedral angle D2 increases from the beginning at −128° to the end at −5° with fourcircle oscillations. The reactive transtocis trajectory in Figure 5 evolutes longer time period than that of reactive cistotrans trajectory in Figure 4.
(a)
(b)
(c)
With randomly choosing 100 initial conditions as mentioned before in both transarea and cisarea, we have sampled 100 independent trajectories to analyze how these swarm of trajectories evolve in average. We plot DD1 and DD2 angle changes against time for these 100 trajectories in Figure 6; Figure 6(a) shows the most of trajectories started from cisarea end around 120 fs in average, while Figure 6(b) shows that the most of trajectories started from transarea end around 200 fs in average which is twice in time compared to trajectories starting from cisarea. This is because the vertical excitation energy for initial condition is 4.45 eV in cisarea and 4.11 eV in transarea, and thus speed of trajectory along DD1 direction on state potential energy surface is larger in cisarea than in transarea. Trajectories starting from cisarea have larger nonadiabatic transition probability to hop than trajectories starting from transarea. Figure 6(a) also shows that trajectories start from cisarea hop around 20 fs, 50 fs, and 100 fs, while trajectories start from transarea shown in Figure 6(b) hop around all times from 20 fs to 200 fs. The number of switching times from one potential energy surface to another is much greater in average for a trajectory staring from transarea than starting from cisarea. This also explains why the branching ratio at OBFCI is dependent on initial conditions; 33% (reactant) to 67% (product) for trajectory starting from cisarea and 48% (reactant) to 52% (product) for trajectory starting from transarea.
(a)
(b)
4. Concluding Remarks
Trajectory surface hopping with constraint molecular dynamics is proposed to investigate nonadiabatic molecular dynamics with reduceddimensional onthefly or analytical potential energy surfaces, and nonadiabatic transitions at the crossing seam are treated by LZ and ZN formulas. Since nonadiabatic transitions are taking place locally around the seam surfaces, we can probe trajectories around local region with reduced dimensional potential energy surfaces, and furthermore we can divide potential energy surfaces into different regions in which different constraints can be applied. In this way, we extend onthefly TSH method to deal with large system for nonadiabatic molecular dynamics. Of course, it goes back to fulldimensional onthefly TSH method if there is no constraint.
The present method is immediately applied to trans↔cis isomerization at OBFCI for Stilbene. Twodimensional PESs for the ground state and the lowest excited state can be well described by dihedral angles DD1 and DD2, in which the seam line is just right at DD1 = −103.3°. Dihedral angle D1 is corresponding to the twist of two phenyl rings around the central ethylenic C=C bond. The present simulation has shown an exact onebond flip behaviors of reactive trajectory, in which D1 changes monotonically from reactant to product for both cistotrans and transtocis isomerizations. Trajectory surface hopping simulation based on the present twodimensional potential energy surfaces can present quite quantitative information for isomerization of Stilbene. We have proposed to calculate two effective parameters and along trajectory moving on adiabatic potential energy surfaces. Effective coupling parameter is not only used for determining transition probability but also for determining the seam line; in the present twodimensional case this is similar to the previous work [69]. However, the present method is useful to deal with highdimensional nonadiabatic molecular dynamics as transition point is determined by local maxima of effective coupling parameter along evolution of trajectory. We do not need to calculate the seam line which is a very difficult task for highdimensional case.
The branching ratios of isomerization at OBFCI have been calculated by setting up initial condition at vertical excitation energy 4.45 eV and 4.11 eV for cisarea and transarea, respectively, and these two areas are local minima on ground state potential energy surface. Quantum yield for transtocis is estimated as 49% in compare with experimental results of 55%. This means that the OBFCI is really responsible for transtocis isomerization. On the other hand, quantum yield for cistotrans is estimated as 47% in compare with experimental results of 35%. This means that OBFCI may be just partly responsible for cistotrans isomerization. As vertical excitation energy at cisarea is 0.34 eV higher than that at transarea based on energy at OBFCI, we expect that there might be a great chance to access (HulaTwist) HTCI for cistotrans isomerization. In the near future, we should carry out trajectory surface hopping simulation based on potential energy surfaces connecting OBFCI with HTCI together. Nevertheless, we conclude that the present simulation is quite encouraged for further investigating on photoisomerization of Stilbene molecule and the present method is very useful in application to the large systems for nonadiabatic molecular dynamics with flexible dimensions of potential energy surfaces.
Appendices
A. Constraint Equations between Cartesian and Internal Coordinates
In order to calculate Jacobian determinant in (7), we need to build up relation equations between Cartesian and internal coordinates for constraint equations. For four atoms plotted in Figure 7, we have relation equations between Cartesian and internal coordinates for bond lengths , bond angles , and dihedral angles in the following equations, respectively: where Cartesian coordinates are , , , , , , , , and and the subscript of each variable is atom number. In terms of the Cartesian coordinate with the derivations of the constraint equations, equation (5) can be rewritten as in which with For dihedral angles, we have Finally, Jacobian determinants in (7) are given as where total atom number and .
B. Unified OnePassage Semiclassical Nonadiabatic Transition Probability
A unified onepassage semiclassical nonadiabatic transition probability was derived analytically as [73] where represents the ratio between adiabatic and diabatic potential energy gaps at transition point : where and correspond to the excited state (upper state) and the ground state (lower state), while and are diabatic potential energy surfaces. The adiabatic parameter that measures nonadiabatic transition strength is expressed in terms of the complex phase integrals: in which at the complex crossing point (its real part stands for transition point), is reduced mass of system, and is collision energy. In order to apply the above formula to multidimensional case, we generalize (B.3) in terms of adiabatic potential energy surfaces at the transition point . Applying the exponential model (see (2.20) of [41]), we can convert (B.3) into where It should be noted that (B.4) is exact for the exponential model and it becomes an approximation when it is generalized to the general twostate system. In comparison with the LandauZener transition formula in (12), we can decompose the adiabatic parameter into the following two diabatic parameters: Equations (B.6) and (B.7) are similar to the expressions derived by Zhu and Nakamura [3], but adiabatization form is much simpler. An effective coupling parameter in (B.6) is very useful to determine nonadiabatic transition zone and the seam surface for multidimensional potential energy surfaces, and an effective collision energy in (B.7) has an finite value at . The parameter defined in (B.2) which represents general type of nonadiabatic transition is only included in (B.7) for the effective collision parameter.
C. The Two Diabatic Parameters and Trajectory Surface Hopping along the Seam Surface
C.1. The Effective Coupling and the Seam Line
The effective coupling parameter in (B.6) is defined along the certain onedimensional direction, and the best direction is represented by the nonadiabatic coupling vector. However, we know that the direction perpendicular to the seam surface can be approximately considered as the direction of the nonadiabatic coupling vector. The seam surface can be calculated from two adiabatic potential energy surfaces without knowing any information of nonadiabatic coupling [69], in which the method requires to solve the certain partial differential equation. This is not easy to be generalized to the highmultidimensional case. We propose in the present work to calculate the local maxima of the effective coupling parameter , and then the collection of all the local maxima constructs the seam surface. The kinetic part of Hamiltonian in terms of two dihedral angles D1 and D2 as shown in (9) can be rewritten as where is the reduced moment of inertia that corresponds to the reduced mass required in (B.6). The variables and are given by At conical configuration point (where D10 = −146.6° and D20 = −60.0°), we estimate from (10) and (11) that 3.5 a.u. and 53.3 a.u. so that 8.6 a.u.
Let us start from conical configuration point (D10, D20) as reference point for the seam and they are transferred to variables and (see in Figure 8) by (C.2), and then we assume that and are a predefined point on the seam with direction normal to distance direction between and ; thus the point along coordinate can be expressed as Derivative of potential energy surfaces defined in (B.6) with respect to at can be estimated as This equation has to be transformed to the original variables D_{1} and D_{2} in terms of which the potential energy surfaces are given by where and are defined in (C.2), and is given as Finally, we can have an expression for effective coupling parameter in which When is rotating with respect to the conic point , the local maxima can be obtained. In this way, we can determine one point on seam surface as well as at that point. Then we can use this point as reference point to search next point on seam surface. Finally, we can determine whole seam surface and distribution on seam surface simultaneously.
C.2. The Effective Collision Energy and the Momentum Changes during Hopping
Classical trajectory is running in the Cartesian coordinate system, while the surface hopping is taking place in the system of internal dihedral angles D1 and D2. We need project internal coordinates to Cartesian coordinates after hopping; besides we cannot use all kinetic energy for hopping and only part of kinetic energy along hopping direction can be utilized. Let us write down classical kinetic energy in terms dihedral angles D1 and D2 as where dot stands for derivative with respect to time , and relation between and is given in (C.2). The last equality in (C.9) is defined as where is angular velocity along direction and is angular velocity along the direction normal to as shown in Figure 8. If we assume that change to after hopping, then they should satisfy where is unchanged and can be estimated from energy conservation as with , in which “+” (“−”) represents hopping from upper (lower) potential to lower (upper) potential. We discuss the classical allowed vertical hop in the present system, namely, since the collision energy is quite large. We now need to convert to that represent angular velocity change after hop along and direction shown in Figure 8. Equations (C.10) and (C.11) lead to the following relation: from which we can obtain This can be converted to original dihedral angles and system as where stand for angular velocity after hop in original dihedral angle system, and and are defined in (C.2). The effective collision energy for in (B.7) is then given by collision energy
Next, we derive relation between internal angular velocities and their corresponding velocities in Cartesian coordinate system. As we analyzed in the text that dihedral angle is defined by the Cartesian coordinates of four atoms; for example, the dihedral angle is given as where and (see Scheme 1 for definition of bond vectors , , and ). Taking derivative with respect to time in (C.17) with consideration that all bond angles and bond distances are fixed, we have where in which we utilized that all bond angles and bond distances are fixed again. Thus, we obtain where In the same way as we derived (C.20), we can have where (see Scheme 1 for definition of bond vectors and ) Assume that velocities after hop in Cartesian coordinate system are , , , , , and , and then they are related to angular velocities after hop as in which and are the same as in (C.20) and (C.22) since we consider the vertical hop. Inserting (C.24) and (C.20) and (C.22) into (C.15) leads to where is defined in (C.15). There are only two equations for six unknown Cartesian velocities after hop so that we need to make some physical intuitions to solve this problem. Let us assume that the relative motions after hop can in general be expressed as The first assumption is that all is equal to zero. Then, the relative motion does not change after hop () as the rotation is around bond 2–3 as shown in Scheme 1. Therefore, we can derive Finally, we can assume after hop that , , , , , and .
The moment of inertia is varying slowly with time as evolution of trajectory. For a convenience, we do surface hopping calculation with the constant moment of inertia at conic point. This introduces small error for energy conservation (this is actually small) because velocities in Cartesian coordinate after hopping might not satisfy constraint conditions exactly. We assume that the Cartesian velocities of atoms 2 and 3 are unchanged during hopping, so that we would do scaling for the Cartesian velocities of atoms 1, 4, 5, and 6 as , , , and to restore energy conservation in which , where is energy lose and is kinetic energy for atoms 1, 4, 5, and 6 before scaling. This can be done iteratively with constraint Hamiltonian equation. Actually, we confirm that is close to unity from numerical calculation during surface hopping ( is exact unity if hopping is taking place right at conic point).
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Yibo Lei would like to acknowledge the Postdoctoral Fellowship supported by National Science Council of China under Grant no. 1002811M009004. This work is supported by National Science Council of China under Grant no. 1002113M009005MY3. Yibo Lei would also like to acknowledge the support from National Science Foundation of China (Grant nos. 21033001, 21103136, and 21173166). Chaoyuan Zhu would like to thank the MOEATU Project of the National Chiao Tung University for support.
References
 J. C. Tully and R. K. Pkeston, “Trajectory surface hopping approach to nonadiabatic molecular collisions: the reaction of H^{+} with D_{2},” The Journal of Chemical Physics, vol. 55, no. 2, pp. 562–572, 1971. View at: Google Scholar
 J. C. Tully, “Molecular dynamics with electronic transitions,” The Journal of Chemical Physics, vol. 93, no. 2, pp. 1061–1071, 1990. View at: Google Scholar
 C. Zhu and H. Nakamura, “Theory of nonadiabatic transition for general twostate curve crossing problems. II. LandauZener case,” The Journal of Chemical Physics, vol. 102, no. 19, pp. 7448–7461, 1995. View at: Google Scholar
 C. Zhu, H. Kamisaka, and H. Nakamura, “New implementation of the trajectory surface hopping method with use of the ZhuNakamura theory. II. application to the charge transfer processes in the 3D ${\text{DH}}_{2}^{+}$ system,” Journal of Chemical Physics, vol. 116, no. 8, pp. 3234–3247, 2002. View at: Publisher Site  Google Scholar
 E. Tapavicza, I. Tavernelli, U. Rothlisberger, C. Filippi, and M. E. Casida, “Mixed timedependent densityfunctional theory/classical trajectory surface hopping study of oxirane photochemistry,” The Journal of Chemical Physics, vol. 129, no. 12, Article ID 124108, 2008. View at: Publisher Site  Google Scholar
 M. Sanier, U. Manthe, and G. Stock, “Quantumclassical Liouville description of multidimensional nonadiabatic molecular dynamics,” The Journal of Chemical Physics, vol. 114, no. 5, pp. 2001–2012, 2001. View at: Publisher Site  Google Scholar
 I. Burghardt and G. Parlant, “On the dynamics of coupled Bohmian and phasespace variables: a new hybrid quantumclassical approach,” The Journal of Chemical Physics, vol. 120, no. 7, pp. 3055–3058, 2004. View at: Publisher Site  Google Scholar
 P. Pechukas, “Timedependent semiclassical scattering theory. II. atomic collisions,” Physical Review, vol. 181, no. 1, pp. 174–185, 1969. View at: Publisher Site  Google Scholar
 J. S. Cao and G. A. Voth, “The formulation of quantum statistical mechanics based on the Feynman path centroid density. IV. algorithms for centroid molecular dynamics,” The Journal of Chemical Physics, vol. 101, no. 7, pp. 6168–6183, 1994. View at: Google Scholar
 M. BenNun, J. Quenneville, and T. J. Martínez, “Ab initio multiple spawning: photochemistry from first principles quantum molecular dynamics,” The Journal of Physical Chemistry A, vol. 104, no. 22, pp. 5172–5175, 2000. View at: Publisher Site  Google Scholar
 A. N. Alexandrova, J. C. Tully, and G. Granucci, “Photochemistry of DNA fragments via semiclassical nonadiabatic dynamics,” The Journal of Physical Chemistry B, vol. 114, no. 37, pp. 12116–12128, 2010. View at: Publisher Site  Google Scholar
 B. F. E. Curchod, T. J. Penfold, U. Rothlisberger, and I. Tavernelli, “Local control theory in trajectorybased nonadiabatic dynamics,” Physical Review A, vol. 84, no. 4, Article ID 042507, 2011. View at: Publisher Site  Google Scholar
 R. CrespoOtero and M. Barbatti, “Cr(CO)_{6} photochemistry: semiclassical study of UV absorption spectral intensities and dynamics of photodissociation,” The Journal of Chemical Physics, vol. 134, no. 16, Article ID 164305, 2011. View at: Publisher Site  Google Scholar
 E. Fabiano, T. W. Keal, and W. Thiel, “Implementation of surface hopping molecular dynamics using semiempirical methods,” Chemical Physics, vol. 349, no. 1–3, pp. 334–347, 2008. View at: Publisher Site  Google Scholar
 Z. Lan, E. Fabiano, and W. Thiel, “Photoinduced nonadiabatic dynamics of pyrimidine nucleobases: onthefly surfacehopping study with semiempirical methods,” The Journal of Physical Chemistry B, vol. 113, no. 11, pp. 3548–3555, 2009. View at: Publisher Site  Google Scholar
 O. Weingart, Z. Lan, A. Koslowski, and W. Thiel, “Chiral pathways and periodic decay in cisazobenzene photodynamics,” The Journal of Physical Chemistry Letters, vol. 2, no. 13, pp. 1506–1509, 2011. View at: Publisher Site  Google Scholar
 T. W. Keal, M. Wanko, and W. Thiel, “Assessment of semiempirical methods for the photoisomerisation of a protonated Schiff base,” Theoretical Chemistry Accounts, vol. 123, no. 12, pp. 145–156, 2009. View at: Publisher Site  Google Scholar
 V. Leyva, I. Corral, F. Feixas et al., “A nonadiabatic quantumclassical dynamics study of the intramolecular excited state hydrogen transfer in orthonitrobenzaldehyde,” Physical Chemistry Chemical Physics, vol. 13, no. 32, pp. 14685–14693, 2011. View at: Publisher Site  Google Scholar
 I. Tavernelli, B. F. E. Curchod, and U. Rothlisberger, “Mixed quantumclassical dynamics with timedependent external fields: a timedependent densityfunctionaltheory approach,” Physical Review A, vol. 81, no. 5, Article ID 052508, 2010. View at: Publisher Site  Google Scholar
 M. Barbatti, R. Shepard, and H. Lischka, “Computational and methodological elements for nonadiabatic trajectory dynamics simulations of molecules,” in Conical Intersections: Theory, Computation and Experiment, W. Domcke, D. R. Yarkony, and H. Köppel, Eds., World Scientific, Singapore, 2010. View at: Google Scholar
 T. Ishida, S. Nanbu, and H. Nakamura, “Nonadiabatic ab initio dynamics of two models of Schiff base retinal,” The Journal of Physical Chemistry A, vol. 113, no. 16, pp. 4356–4366, 2009. View at: Publisher Site  Google Scholar
 A. Gao, B. Li, P. Zhang, and K. Han, “Nonadiabatic ab initio molecular dynamics of photoisomerization in bridged azobenzene,” The Journal of Chemical Physics, vol. 137, no. 20, Article ID 204305, 2012. View at: Publisher Site  Google Scholar
 A. Gao, B. Li, P. Zhang, and J. Liu, “Photochemical dynamics simulations for transcis photoisomerizations of azobenzene and bridged azobenzen,” Computational and Theoretical Chemistry, vol. 1031, pp. 13–21, 2014. View at: Publisher Site  Google Scholar
 J.P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, “Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes,” Journal of Computational Physics, vol. 23, no. 3, pp. 327–341, 1977. View at: Google Scholar
 V. Kräutler, W. F. van Gunsteren, and P. H. Hünenberger, “A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations,” Journal of Computational Chemistry, vol. 22, no. 5, pp. 501–508, 2001. View at: Google Scholar
 G. Ciccotti and J.P. Ryckaert, “Molecular dynamics simulation of rigid molecules,” Computer Physics Reports, vol. 4, no. 6, pp. 346–392, 1986. View at: Google Scholar
 P. Gonnet, “PSHAKE: a quadratically convergent SHAKE in O(n^{2}),” Journal of Computational Physics, vol. 220, no. 2, pp. 740–750, 2007. View at: Publisher Site  Google Scholar
 S. Miyamoto and P. A. Kollman, “Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models,” Journal of Computational Chemistry, vol. 13, no. 8, pp. 952–962, 1992. View at: Google Scholar
 E. Barth, K. Kuczera, B. Leimkuhler, and R. D. Skeel, “Algorithms for constrained molecular dynamics,” Journal of Computational Chemistry, vol. 16, no. 10, pp. 1192–1209, 1995. View at: Publisher Site  Google Scholar
 P. E. Gill, W. Murray, and M. H. Wright, Numerical Linear Algebra and Optimization, AddisonWesley, Redwood City, Calif, USA, 1991.
 T. Chu, Y. Zhang, and K. Han, “The timedependent quantum wave packet approach to the electronically nonadiabatic processes in chemical reactions,” International Reviews in Physical Chemistry, vol. 25, no. 12, pp. 201–235, 2006. View at: Publisher Site  Google Scholar
 K. Han, J. Huang, S. Chai, S. Wen, and W. Deng, “Anisotropic mobilities in organic semiconductors,” Nature Protocol Exchange, 2013. View at: Publisher Site  Google Scholar
 Y. Lei, C. Zhu, Z. Wen, and S. Lin, “New implementation of semiclassical dynamic simulation on the photoisomerization of cis and trans isomers of free stilbene,” Acta Chimica Sinica, vol. 70, no. 17, pp. 1869–1876, 2012. View at: Publisher Site  Google Scholar
 A. J. Neukirch, L. C. Shamberger, E. Abad et al., “Nonadiabatic ensemble simulations of cisstilbene and cisazobenzene photoisomerization,” Journal of Chemical Theory and Computation, vol. 10, no. 1, pp. 14–23, 2014. View at: Publisher Site  Google Scholar
 Y. Ootani, K. Satoh, A. Nakayama, T. Noro, and T. Taketsugu, “Ab initio molecular dynamics simulation of photoisomerization in azobenzene in the n π state,” The Journal of Chemical Physics, vol. 131, no. 19, Article ID 194306, 2009. View at: Publisher Site  Google Scholar
 A. Toniolo, C. Ciminelli, M. Persico, and T. J. Martínez, “Simulation of the photodynamics of azobenzene on its first excited state: comparison of full multiple spawning and surface hopping treatments,” The Journal of Chemical Physics, vol. 123, no. 23, Article ID 234308, 2005. View at: Publisher Site  Google Scholar
 P. Cattaneo and M. Persico, “An abinitio study of the photochemistry of azobenzene,” Physical Chemistry Chemical Physics, vol. 1, pp. 4739–4743, 1999. View at: Publisher Site  Google Scholar
 M. Pederzoli, J. Pittner, M. Barbatti, and H. Lischka, “Nonadiabatic molecular dynamics study of the cistrans photoisomerization of azobenzene excited to the S_{1} state,” The Journal of Physical Chemistry A, vol. 115, no. 41, pp. 11136–11143, 2011. View at: Publisher Site  Google Scholar
 D. Gegiou, K. A. Muszkat, and E. Fischer, “Temperrature dependence of photoisomerization. VI. viscosity effect,” The Journal of the American Chemical Society, vol. 90, no. 1, pp. 12–18, 1968. View at: Publisher Site  Google Scholar
 D. Gegiou, K. A. Muszkat, and E. Fischer, “Temperature dependence of photoisomerization. V. effect of subsituents on the photoisomerization of stilbenes and azobenzenes,” Journal of the American Chemical Society, vol. 90, no. 15, pp. 3907–3917, 1968. View at: Publisher Site  Google Scholar
 A. Gilbert and J. Baggott, Essentials of Molecular Photochemistry, Blackwell Science, Oxford, UK, 1991.
 R. S. H. Liu, “Photoisomerization by HulaTwist: a fundamentalsupramolecular photochemical reaction,” Accounts of Chemical Research, vol. 34, no. 7, pp. 555–562, 2001. View at: Publisher Site  Google Scholar
 R. D. Sampedro, A. Cembran, M. Garavelli, M. Olivucci, and W. Fuβ, “Structure of the conical intersections driving the cistrans photoisomerization of conjugated molecules,” Photochemistry and Photobiology, vol. 76, no. 6, pp. 622–633, 2002. View at: Google Scholar
 Y. Amatatsu, “Ab initio study on the electronic structures of stilbene at the conical intersection,” Chemical Physics Letters, vol. 314, no. 34, pp. 364–368, 1999. View at: Publisher Site  Google Scholar
 J. Quenneville and T. J. Martínez, “Ab initio study of cistrans photoisomerization in stilbene and ethylene,” The Journal of Physical Chemistry A, vol. 107, no. 6, pp. 829–837, 2003. View at: Publisher Site  Google Scholar
 M. J. Bearpark, F. Bernardi, S. Clifford, M. Olivucci, M. A. Robb, and T. Vreven, “Cooperating rings in cisstilbene lead to an S_{0}/S_{1} conical intersection,” The Journal of Physical Chemistry A, vol. 101, no. 21, pp. 3841–3847, 1997. View at: Google Scholar
 J. H. Frederick, Y. Fujiwara, J. H. Penn, K. Yoshihara, and H. Petek, “Models for stilbene photoisomerization: experimental and theoretical studies of the excitedstate dynamics of 1,2diphenylcycloalkenes,” The Journal of Physical Chemistry, vol. 95, no. 7, pp. 2845–2858, 1991. View at: Google Scholar
 V. D. Vachev, J. H. Frederick, B. A. Grishanin, V. N. Zadkov, and N. I. Koroteev, “Stilbene isomerization dynamics on multidimensional potential energy surface. Molecular dynamics simulation,” Chemical Physics Letters, vol. 215, no. 4, pp. 306–314, 1993. View at: Publisher Site  Google Scholar
 V. D. Vachev, J. H. Frederick, B. A. Grishanin, V. N. Zadkov, and N. I. Koroteev, “Quasiclassical molecular dynamics simulation of the photoisomerization of stilbene,” The Journal of Physical Chemistry, vol. 99, no. 15, pp. 5247–5263. View at: Google Scholar
 G. Orlandi, L. Gagliardi, S. Melandri, and W. Caminati, “Torsional potential energy surface and vibrational levels in trans stilbene,” Journal of Molecular Structure, vol. 612, no. 23, pp. 383–391, 2002. View at: Publisher Site  Google Scholar
 W.V. Chiang and J. Laane, “Fluorescence spectra and torsional potential functions for transstilbene in its S_{0} and S_{1}(π,π*) electronic states,” The Journal of Chemical Physics, vol. 100, no. 12, pp. 8755–8767, 1994. View at: Google Scholar
 C. W. Jiang, R. H. Xie, F. L. Li, and R. E. Allen, “Transtocis isomerization of stilbene following an ultrafast laser pulse,” Chemical Physics Letters, vol. 474, no. 4–6, pp. 263–267, 2009. View at: Publisher Site  Google Scholar
 Y. Dou and R. E. Allen, “Detailed dynamics of a complex photochemical reaction: cistrans photoisomerization of stilbene,” The Journal of Chemical Physics, vol. 119, no. 20, pp. 10658–10666, 2003. View at: Publisher Site  Google Scholar
 S. Takeuchi, S. Ruhman, T. Tsuneda, M. Chiba, T. Taketsugu, and T. Tahara, “Spectroscopic tracking of structural evolution in ultrafast stilbene photoisomerization,” Science, vol. 322, no. 5904, pp. 1073–1077, 2008. View at: Publisher Site  Google Scholar
 W. Fuß, C. Kosmidis, W. E. Schmid, and S. A. Trushin, “The lifetime of the perpendicular minimum of cisstibene observed by dissociative intenselaser field ionization,” Chemical Physics Letters, vol. 385, pp. 423–430, 2004. View at: Publisher Site  Google Scholar
 W. Fuß, C. Kosmidis, W. E. Schmid, and S. A. Trushin, “The photochemical cistrans isomerization of free stilbene molecules follows a hulatwist pathway,” Angewandte Chemie International Edition, vol. 43, no. 32, pp. 4178–4182, 2004. View at: Publisher Site  Google Scholar
 D. C. Todd and G. R. Fleming, “Cisstilbene isomerization: temperature dependence and the role of mechanical friction,” The Journal of Chemical Physics, vol. 98, no. 1, pp. 269–279, 1993. View at: Google Scholar
 S. Abrash, S. Repinec, and R. M. Hochstrasser, “The viscosity dependence and reaction coordinate for isomerization of cisstilbene,” The Journal of Chemical Physics, vol. 93, no. 2, pp. 1041–1053, 1990. View at: Google Scholar
 R. J. Sension, S. T. Repinec, A. Z. Szarka, and R. M. Hochstrasser, “Femtosecond laser studies of the cisstilbene photoisomerization reactions,” The Journal of Chemical Physics, vol. 98, no. 8, pp. 6291–6315, 1993. View at: Google Scholar
 R. J. Sension, S. T. Repinec, and R. M. Hochstrasser, “Femtosecond laser study of energy disposal in the solution phase isomerization of stilbene,” The Journal of Chemical Physics, vol. 93, no. 12, pp. 9185–9188, 1990. View at: Google Scholar
 S. T. Repinec, R. J. Sension, A. Z. Szarka, and R. M. Hochstrasser, “Femtosecond laser studies of the cisstilbene photoisomerization reactions: the cisstilbene to dihydrophenanthrene reaction,” The Journal of Physical Chemistry, vol. 95, no. 25, pp. 10380–10385, 1991. View at: Publisher Site  Google Scholar
 R. Improta and F. Santoro, “Excitedstate behavior of trans and cis isomers of stilbene and stiff stilbene: a TDDFT study,” The Journal of Physical Chemistry A, vol. 109, no. 44, pp. 10058–10067, 2005. View at: Publisher Site  Google Scholar
 C. D. Berweger, W. F. van Gunsteren, and F. MüllerPlathe, “Molecular dynamics simulation with an ab initio potential energy function and finite element interpolation: the photoisomerization of cisstilbene in solution,” The Journal of Chemical Physics, vol. 108, no. 21, pp. 8773–8781, 1998. View at: Google Scholar
 J. Tatchen and E. Pollak, “Ab initio spectroscopy and photoinduced cooling of the transstilbene molecule,” The Journal of Chemical Physics, vol. 128, no. 16, Article ID 164303, 2008. View at: Publisher Site  Google Scholar
 A. Weigel and N. P. Ernsting, “Excited Stilbene: intramolecular vibrational redistribution and solvation studied by femtosecond stimulated raman spectroscopy,” The Journal of Physical Chemistry B, vol. 114, no. 23, pp. 7879–7893, 2010. View at: Publisher Site  Google Scholar
 G. Karlström, R. Lindh, P.Å. Malmqvist et al., “MOLCAS: a program package for computational chemistry,” Computational Materials Science, vol. 28, no. 2, pp. 222–239, 2003. View at: Publisher Site  Google Scholar
 Matlab, User Guide, The Mathworks, Torrance, Calif, USA, http://www.mathworks.com.
 J. K. Rice and A. P. Baronavski, “Ultrafast studies of solvent effects in the isomerization of cisstilbene,” The Journal of Physical Chemistry, vol. 96, no. 8, pp. 3359–3366, 1992. View at: Google Scholar
 C. Zhu, K. Nobusada, and H. Nakamura, “New implementation of the trajectory surface hopping method with use of the ZhuNakamura theory,” The Journal of Chemical Physics, vol. 115, no. 7, pp. 3031–3044, 2001. View at: Publisher Site  Google Scholar
 J. J. Moré, B. S. Garbow, and K. E. Hillstrom, “User guide for MINPACK1,” Argonne National Laboratory Report ANL8074, Argonne, Ill, USA, 1980. View at: Google Scholar
 M. J. D. Powell, Numerical Methods for NonLinear Algebraic Equations, Gordon and Breach, New York, NY, USA, 1970, Edited by P. Rabinowitz.
 L. F. Shampine and H. A. Watts, Mathematical Software III, Academic Press, New York, NY, USA, 1977, Edited by J. R. Rice.
 C. Zhu, “Unified semiclassical theory for the twostate system: analytical solutions for scattering matrices,” The Journal of Chemical Physics, vol. 105, no. 10, pp. 4159–4172, 1996. View at: Google Scholar
Copyright
Copyright © 2014 Yibo Lei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.