Constraint Trajectory Surface-Hopping Molecular Dynamics Simulation of the Photoisomerization of Stilbene

1 Key Laboratory of Synthetic and Natural Functional Molecule Chemistry of Ministry of Education, College of Chemistry & Materials Science, Shaanxi Key Laboratory of Physico-Inorganic Chemistry, Northwest University, Xi’an 710069, China 2Department of Applied Chemistry, Institute of Molecular Science and Center for Interdisciplinary Molecular Science, National Chiao-Tung University, Hsinchu 300, Taiwan


Introduction
Nonadiabatic dynamic based on on-the-fly 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 on-the-fly 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 DH 2 + nonadiabatic reactions by Tully and Pkeston [1] in which nonadiabatic transition probabilities were calculated by the Landau-Zener (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 time-dependent first-order 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 DH 2 + 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][12][13][14][15][16][17][18][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 2 International Journal of Photoenergy 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 on-the-fly 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 on-the-fly 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 on-the-fly for trajectories with full degree of freedoms for large systems.However, it can run limited number of on-the-fly trajectories with full-dimensional 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 on-the-fly 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 reduced-dimensional either on-the-fly or analytical potential energy surfaces depending on how many degrees of freedom are taken into account.As reduced-dimensional 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][25][26][27][28][29][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 on-the-fly 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 on-the-fly 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 time-dependent DFT local orbital framework [34,35].In contrast, there are many publications on on-the-fly TSH method for nonadiabatic dynamics of azobenzene [22,23,[34][35][36][37][38].Therefore, we apply the present method to study photoisomerization of Stilbene in detail.
The isomerization of cis-and trans-Stilbene 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 one-bond flip (OBF) [41], Hula-Twist (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 vinyl-phenyl 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 cis-to-trans and trans-to-cis isomerizations are taken place by going through pathway of conical intersections which are named as OBF-CI [44,45] and HT-CI [46], respectively.
The Stilbene has cis-and trans-isomers 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 trans-conformations has been investigated by computing the ground state  0 and the lowest excited state  1 .The dominant excitation from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) corresponds to the  1  * 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,.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][45][46][47][48][49][50][51].The conical intersection OBF-CI 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 normal-mode 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 trans-photoisomerizations of nonpolar conjugated molecules and these were investigated by the experimental measurements as well [58][59][60].
Starting from the Franck-Condon (FC) regions, PESs along the reaction coordinate (curved one-dimension) 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 OBF-CI as mentioned above.The two torsion angles in terms phenyl rotation and the pyramidalization coordinate were utilized for constructing twodimensional PESs around the OBF-CI [45].These studies told that at the OBF-CI the two-dimensional PESs can well describe trans-cis nonadiabatic dynamics for Stilbene.The branching ratio is basically determined by potential energy surfaces near the OBF-CI.The one torsion angle around the central ethylenic C=C bond plus one dihedral angle was a suitable choice for isomerization via the OBF-CI (called as D1 and D2 in the present paper as shown in Figure 1).The potential energy surfaces of ground-state  0 and the first excited state  1 are constructed in terms of the D1 and D2 angles, and the other internal coordinates are fixed at the OBF-CI configuration.
The rest of this paper is organized as follows.Section 2 first summarizes constrained molecular dynamics method and constructs two-dimensional 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.

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 two-dimensional 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 on-the-fly nonadiabatic dynamics with four or more dimensional potential energy surfaces.However, it was well studied that two-dimensional PESs are suitable for describing isomerization of Stilbene at OBF-CI.Moreover, for two dimensions we can construct analytical PESs easily and can demonstrate how trajectories run around the seam.

Constraint Molecular Dynamics.
We review constraint molecular dynamics for system with   constraints in which the canonical Hamiltonian equations can be written as [24][25][26][27][28][29][30] where   ,   , and   are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respectively.The   () = ( = 1, 2, . . .,   ) are constraint equations determined by   internal coordinates that are fixed at the certain structures of nuclear configuration all the time.The  uc

𝑖
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  uc  ( + Δ) 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   = 3 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 M-SHAKE Newton method [24,25,30] and quasi-Newton method [30] were straight-forward and fast to solve the nonlinear equations.All of these methods involved the calculations of the Jacobian determinant about the constraint equations   = 0: where   are then updated repeatedly by solving the system of linear equations iteratively: until the max |  | <  ( = 1, . . .,   ) in which  is small number with the prescribed tolerance of the constraints.The Jacobian determinant for the quasi-Newton 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 M-SHAKE 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.

Two-Dimensional Potential Energy Surfaces for Stilbene.
Quenneville and Martínez [45] have done extensive calculation for the critical  1 / 0 conical intersection OBF-CI 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 self-consistent field (CASSCF) and multireference perturbation theory (CASPT2) methods with basis sets of 6-31G., 6-31G * , and 6-31G * * .Following their conclusion, we can safely utilize state-averaged CASSCF method abbreviated as SA-N-CAS(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 SA-2-CAS(2/2)/6-31G method to reproduce geometries of conical intersection OBF-CI, the local minima on  0 and  1 potential energy surfaces given by [45].Then, we concentrate on configuration of conical intersection at which the torsion angle D1 (C11-C8-C1-C2) and dihedral angle D2 (H10-C8-C1-H9) are −146.6∘ and −60.0 ∘ , respectively, as shown in Figure 1.We have carried out a preliminary test calculation for two-dimensional potential energy surfaces with variables D1 and D2 (all bond lengths, all bond angles, and the other dihedral angles are fixed at OBF-CI configuration).We found that DD1 = (D1 + D2)/2 and DD2 = (D2 − D1)/2 are better variables to describe the isomerization.If we set up interesting region for potential energy surface  1 with energy not higher than 5 eV above OBF-CI (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 0, −5, . . ., −180 (37 grid points) and DD2 is −20, −15, . . ., 90 (23 grid points).Total configurations are 851 plus OBF-CI which all are computed by SA-2-CAS(2/2)/6-31G 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: +   +  24 cos ( + ) + where  = DD1 and  = DD2,  0 = −103.3∘ and  0 = 43.3∘ are angles at OBF-CI.The fitting parameters  1 ,  1 , and  0 to  32 are all given in Table 1.Potential energy surfaces of  0 and  1 calculated from ( 8) are plotted in Figure 2 with zero point of potential energy at ( 0 ,  0 ) = 0 in a unit of electron volt.We computed the mean absolute error between calculated results from (8) and numerical data computed from method SA-2-CAS(2/2)/6-31G 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 OBF-CI region.From analytical PES in (8), we have estimated local trans-isomer at DD1 = −164.3∘ and DD2 = 36.7 ∘ with energy 1.72 eV below OBF-CI, and local cis-isomer at DD1 = −39.2∘ and DD2 = 50.9∘ with energy 1.28 eV below OBF-CI.Furthermore, the vertical excitation energies have been calculated to be 4.11 eV and 4.45 eV at local trans-and cis-minima, 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  1 state is just right at conical intersection for the present two-dimensional potential energy surface.

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 Φ1 and Φ2 are rotation angles around -axis for two phenyl rings, respectively.Thus, the moment of inertia  D1 is the reduced inertia from two phenyl rings where  Φ1 = ∑     2  (  is atom mass in the one phenyl ring and   is distance from atom  to -axis) and  Φ2 is estimated International Journal of Photoenergy from the other phenyl ring.The same analysis can be applied for the moment of inertia  D2 with two hydrogen atoms: where  1 =    2  (  is hydrogen atom mass and   is distance from hydrogen atom  to -axis) and  2 is estimated from the other hydrogen atom.
One-passage 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  2 and  2 is given in Appendix B. We can simply look at the present two-dimensional PESs in Figure 2 to find that the seam line is appearing at the DD1 = −103.3∘ .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  2 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  2 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 DD1 = −103.3∘ with negligible dependence of DD2 angle.The effective coupling parameter  2 along the seam line is calculated and depicted in Figure 3. Figure 3 shows that the diabatic region (where the one-passage 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 0.1 ≤  2 ≤ 10.

Results and Discussions
Experimental measurement starts from global cis-minimum and trans-minimum in which photoisomerization processes involve dynamics on full-dimensional potential energy surfaces.We construct the reduced two-dimensional (2D) potential energy surfaces at OBF-CI, and  0 potential energy surface shows clear local cis-minimum and trans-minimum separated by the seam line.We assume that trajectory starts from global cis (trans)-minimum with photo-excitation from  0 to  1 PES and then reaches the certain area of the present local cis (trans)-minimum on the 2D  1 PES, from where trajectories start isomerization.Then, question is how to determine such area for both local cis-and trans-isomers.Let us start with the test local cis-area and local trans-area 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 cis-to-trans isomerization is defined by the excitation energy 4.45 eV from  0 to  1 state estimated at local cis-minimum (DD1 = −39.2∘ , DD2 = 50.9∘ ).The initial conditions of DD1 and DD2 for trajectories are randomly picked up in the local cis-area vertically excited from  0 to  1 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 cis-area on  1 state can have three ending ways; one is that it enters local cis-area on  0 (called cis-to-cis nonreactive trajectory), another is that it enters local-trans area on  0 (called cis-to-trans 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  init to atoms 1, 4, 5, and 6 with their initial velocities analyzed in the following.As shown in Scheme 1, the z-component of velocity is assumed to be zero so that for a given kinetic energy  init , we have the following relations: International Journal of Photoenergy 7 where  is random number in the range of [0, 1] and  1 ( 5 ) 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  init 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 trans-tocis isomerization, and the total energy is defined by the excitation energy 4.11 eV from  0 to  1 state estimated at local trans-minimum (DD1 = −164.3∘ , DD2 = 36.7 ∘ ).The initial conditions of DD1 and DD2 for trajectories are randomly picked up in the local trans-area vertically excited from  0 to  1 state.Three types of ending trajectories are also collected as mentioned above (trans-to-trans nonreactive, trans-tocis 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   () = 0.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 fourth-order Runge-Kutta 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  2 and effective collision energy  2 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  2 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.

Branching Ratios of Cis-to-Trans and Trans-to-Cis Isomerizations.
By giving equal weight to all the sampling trajectories, we compute the branching ratios with initial condition starting from both local cis-area 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 transto-cis isomerization is 48 : 52 (experimental value is 50 : 50) in the present simulation.By considering a potential barrier in Frank-Condon region of global trans-minimum which is taking 5% trajectories down to global trans-minimum by fluorescence [59,61], we can estimate the quantum fields for trans-to-cis isomerization as 0.95 × 52% = 49% which is in very good agreement with the experiment measurement 55.0%.The branching ratio for cis-to-trans isomerization is 33 : 67 (experimental value is 50 : 50) in the present simulation.There are two pathways when Stilbene is photo-excited from the global cis-isomer  0 to  1 state (see [59,61] and references therein); one pathway leads to the OBF-CI 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 cis-to-trans 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 trans-area or local cis-area and this means that branching ratio is not 50 : 50 at OBF-CI.
Potential energy surface of  1 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  2 is well above the seam and this makes trajectory have a big chance to hop down  0 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 OBF-CI at cis-area and is 2.39 eV above OBF-CI at trans-area, so that effective collision energy  2 is higher for cis-to-trans than for trans-to-cis.This explains why branching ratio is more in favor of reactive cis-to-trans than trans-to-cis.We know that the conical intersection Hula-Twist (HT-CI) is higher in energy than in OBF-CI, and these two CIs might be in close configuration.Therefore, we suspect that there is a big chance for cis-to-trans trajectory to access HT-CI with high vertical excitation energy, but this is not the case for trans-tocis trajectory with low vertical excitation energy.
We have also tested how branching ratios are sensitive to choice of size of trans-and cis-area, and we found that they are not so sensitive to the choice of the area size unless we define very small area for cis-area and trans-area.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 ( 2 and  2 ) are quite robust in the present formulation.

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 cis-to-trans isomerization varies with time by starting at local cis-area with photo-excitation 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  1 state and increases on  0 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  1 to  0 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 trans-area 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 cis-area) to −185 ∘ (at local trans-area), and this is exactly corresponding to what it is called as one-bond flip process around the central C=C bond.Dihedral angle D2 decreases from the beginning at 7.5 ∘ to the end at −95 ∘ with three-circle 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.
Figure 5 shows how a reactive trajectory for trans-tocis isomerization varies with time by starting at local transarea with photo-excitation energy 3.98 eV and evolution time is 133 fs.Potential energies on  0 and  1 states vary with time with four selected snapshots of structure as shown in Figure 5(a) for a reactive trans-to-cis trajectory.There are three-time surface hopping from one potential energy surface to another during evolution of trajectory; the first hopping from  1 to  0 state happened at 50 fs, the second hopping from  0 to  1 state at 71 fs, and the last hopping from  1 to  0 state at 88 fs.After that, the trajectory is propagated on  0 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 trans-area) to −97 ∘ (at local cis-area), and this shows one-bond flip picture around the central C=C bond again.Dihedral angle D2 increases from the beginning at −128 ∘ to the end at −5 ∘ with four-circle oscillations.The reactive trans-to-cis trajectory in Figure 5 evolutes longer time period than that of reactive cisto-trans trajectory in Figure 4.
With randomly choosing 100 initial conditions as mentioned before in both trans-area and cis-area, 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 trans-area end around 200 fs in average which is twice in time compared to trajectories starting from cis-area.This is because the vertical excitation energy for initial condition is 4.45 eV in cis-area and 4.11 eV in trans-area, and thus speed of trajectory along DD1 direction on  1 state potential energy surface is larger in cis-area than in trans-area.Trajectories starting from cis-area have larger nonadiabatic transition probability to hop than trajectories starting from trans-area. Figure 6(a) also shows that trajectories start from cis-area hop around 20 fs, 50 fs, and 100 fs, while trajectories start from trans-area 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 trans-area than starting from cis-area.This also explains why the branching ratio at OBF-CI is dependent on initial conditions; 33% (reactant) to 67% (product) for trajectory starting from cis-area and 48% (reactant) to 52% (product) for trajectory starting from trans-area.

Concluding Remarks
Trajectory surface hopping with constraint molecular dynamics is proposed to investigate nonadiabatic molecular dynamics with reduced-dimensional on-the-fly 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 on-the-fly TSH method to deal with large system for nonadiabatic molecular dynamics.Of course, it goes back to full-dimensional on-the-fly TSH method if there is no constraint.
The present method is immediately applied to trans↔cis isomerization at OBF-CI for Stilbene.Two-dimensional 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 one-bond flip behaviors of reactive trajectory, in which D1 changes monotonically from reactant to product for both cis-to-trans and trans-to-cis isomerizations.Trajectory surface hopping simulation based on the present two-dimensional potential energy surfaces can present quite quantitative information for isomerization of Stilbene.We have proposed to calculate two effective parameters  2 and  2 along trajectory moving on adiabatic potential energy surfaces.Effective coupling parameter  2 is not only used for determining transition probability but also for determining the seam line; in the present two-dimensional case this is International Journal of Photoenergy similar to the previous work [69].However, the present method is useful to deal with high-dimensional nonadiabatic molecular dynamics as transition point is determined by local maxima of effective coupling parameter  2 along evolution of trajectory.We do not need to calculate the seam line which is a very difficult task for high-dimensional case.
The branching ratios of isomerization at OBF-CI have been calculated by setting up initial condition at vertical excitation energy 4.45 eV and 4.11 eV for cis-area and trans-area, respectively, and these two areas are local minima on ground state potential energy surface.Quantum yield for trans-tocis is estimated as 49% in compare with experimental results of 55%.This means that the OBF-CI is really responsible for trans-to-cis isomerization.On the other hand, quantum yield for cis-to-trans is estimated as 47% in compare with experimental results of 35%.This means that OBF-CI may be just partly responsible for cis-to-trans isomerization.As vertical excitation energy at cis-area is 0.34 eV higher than that at trans-area based on energy at OBF-CI, we expect that there might be a great chance to access (Hula-Twist) HT-CI for cis-to-trans isomerization.In the near future, we should carry out trajectory surface hopping simulation based on potential energy surfaces connecting OBF-CI with HT-CI 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.

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   1  2 , bond angles   1  2  3 , and dihedral angles   1  2  3  4 in the following equations, respectively: , and ⃗  = ⃗  4 − ⃗  3 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 (A.7) Finally, Jacobian determinants in (7) are given as where  = total atom number and  = , , .

B. Unified One-Passage Semiclassical Nonadiabatic Transition Probability
A unified one-passage semiclassical nonadiabatic transition probability was derived analytically as [73] where  represents the ratio between adiabatic and diabatic potential energy gaps at transition point  0 :  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 two-state system.In comparison with the Landau-Zener transition formula in (12), we can decompose the adiabatic parameter  into the following two diabatic parameters:

C. The Two Diabatic Parameters and Trajectory Surface Hopping along the Seam Surface
C.1.The Effective Coupling  2 and the Seam Line.The effective coupling parameter  2 in (B.6) is defined along the certain one-dimensional 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 high-multidimensional case.We propose in the present work to calculate the local maxima of the effective coupling parameter  2 , 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  = √ D1  D2 is the reduced moment of inertia that corresponds to the reduced mass required in (B.6).The variables  and  are given by Let us start from conical configuration point (D10, D20) as reference point for the seam and they are transferred to variables  0 and  0 (see in Figure 8) by (C.2), and then we assume that  1 and  1 are a predefined point on the seam with  direction normal to distance direction between ( 0 ,  0 ) and ( 1 ,  1 ); thus the point (, ) Finally, we can have an expression for effective coupling parameter in which When  is rotating with respect to the conic point ( 0 ,  0 ), the local maxima  2 can be obtained.In this way, we can determine one point on seam surface as well as  2 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  2 distribution on seam surface simultaneously.

C.2. The Effective Collision Energy 𝑏 2 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 ( Ḋ1, Ḋ2) and ( ẋ , ẏ ) is given in (C.2).The last equality in (C.9) is defined as where V 1 is angular velocity along  direction and V 2 is angular velocity along the direction normal to  as shown in Figure 8.

If we assume that
after hopping, then they should satisfy where V 2 is unchanged and  can be estimated from energy conservation as (C.12) with Δ 0 = ( + −  − )/2, in which "+" ("−") represents hopping from upper (lower) potential to lower (upper) potential.We discuss the classical allowed vertical hop in the present system, namely,  ≥ 0 since the collision energy is quite large.We now need to convert (V  1 , V  2 ) to ( ẋ  , ẏ  ) that represent angular velocity change after hop along  and  direction shown in Figure 8 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 D1 is given as .The moment of inertia  = √ D1  D2 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 ̇⃗  =  2 ̇⃗  6 to restore energy conservation in which  = √1 − Δ/, 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).

Figure 2 :
Figure 2: Analytical two-dimensional PESs for the ground state  0 and the first excited state  1 with respect to the combined internal angles DD1 and DD2 [33].

Figure 4 :
Figure 4: A typical reactive trajectory for cis-to-trans isomerization starts from local cis-area with photo-excitation energy 4.42 eV.(a) Snapshots of potential energy surfaces and structures varying with time, (b) DD1 and DD2, and (c) D1 and D2 varying with time.

Figure 5 :Figure 6 :
Figure 5: A typical reactive trajectory for trans-to-cis isomerization starts from local trans-area with photo-excitation energy 3.98 eV.(a) Snapshots of potential energy surfaces and structures varying with time, (b) DD1 and DD2, and (c) D1 and D2 varying with time.

Figure 7 :
Figure 7: Schematic diagram for relations between Cartesian and internal coordinates for four atoms.

2 . (B. 7 )
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  2 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  2 in (B.7) has an finite value at  =  0 .The parameter  defined in (B.2) which represents general type of nonadiabatic transition is only included in (B.7) for the effective collision parameter.

2 Figure 8 :
Figure 8: Schematic diagram in which  is assumed as direction of nonadiabatic coupling vector.

Table 2 :
Branching ratios (numbers in parentheses are calculated from Landau-Zener formula) for cis-to-trans and trans-to-cis isomerizations at one-bond flip conical intersection (OBF-CI).Total energy of each classical trajectory is set up to the vertical excitation energy 4.45 eV at cis-area and 4.11 eV at trans-area. [59,61]erences[59,61].
[41]where  + and  − correspond to the excited state (upper state) and the ground state (lower state), while  2 and  1 are diabatic potential energy surfaces.The adiabatic parameter  that measures nonadiabatic transition strength is expressed in terms of the complex phase integrals: ) =  − ( * ) at the complex crossing point  * (its real part  0 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  0 .Applying the exponential model (see(2.20)of[41]),we can convert (B.3) into * Derivative of potential energy surfaces defined in (B.6) with respect to  at ( 1 ,  1 ) can be estimated as where  1 and  2 are defined in (C.2), and  is given as . Equations (C.10) and (C.11) lead to the following relation: ẏ  cos  − ẋ  sin  =  ( ẏ cos  − ẋ sin ) ,  +  sin 2  (1 − ) sin  cos  (1 − ) sin  cos  cos 2  + sin 2   , Ḋ2  ) stand for angular velocity after hop in original dihedral angle system, and  1 and  2 are defined in (C.2).The effective collision energy for  2 in (B.7) is then given by collision energy ẋ  cos  + ẏ  sin  = ẋ cos  + ẏ sin , Taking derivative with respect to time  in (C.17) with consideration that all bond angles and bond distances are fixed, we have  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   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 as the rotation is around bond 2-3 as shown in Scheme 1.Therefore, we can derive (C.26)The first assumption is that all   is equal to zero.Then, the relative motion ( ̇⃗   3 − ̇⃗   2 ) does not change after hop ( 3 = 1)