Hindawi Publishing Corporation International Journal of Photoenergy Volume 2014, Article ID 132149, 18 pages http://dx.doi.org/10.1155/2014/132149

Research Article

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

Yibo Lei,1,2 Shaomei Wu,1 Chaoyuan Zhu,2 Zhenyi Wen,1 and Sheng-Hsien Lin2

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, Xian 710069, China

2 Department of Applied Chemistry, Institute of Molecular Science and Center for Interdisciplinary Molecular Science, National Chiao-Tung University, Hsinchu 300, Taiwan

Correspondence should be addressed to Yibo Lei; leiyb@nwu.edu.cn and Chaoyuan Zhu; cyzhu@mail.nctu.edu.tw Received 3 March 2014; Revised 15 April 2014; Accepted 16 April 2014; Published 7 July 2014 Academic Editor: Yusheng Dou

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.

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 Landau-Zener and Zhu-Nakamura formulas along the seam. Two-dimensional potential energy surfaces of ground S0 and excited Sj 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 trans-to-cis isomerization is estimated as 49% in very good agreement with experimental value of 55%, while quantum yield for cis-to-trans isomerization is estimated as 47% in comparison with experimental value of 35%.

1. 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 DH2+ 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 DH2+ 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 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-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 nona-diabatic transition zone. It is useful complementary to full-dimensional 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-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=Cbond 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 photo-induced 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 S0 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 nn* 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 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-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 one-dimensional PESs with respect to the main twist of phenyl

D1: C11-C8-C1-C2 D2: H10-C8-C1-H9

Figure 1: Atom numbering of Stilbene at the Sl/S0 conical intersection (OBF-CI) configuration in terms of two dihedral angles D1 and D2 at which D10 = -146.6° and D20 = -60.0°.

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 pyra-midalization coordinate were utilized for constructing two-dimensional 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 S0 and the first excited state Sj 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 one-passage 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 Nc 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.

2.1. Constraint Molecular Dynamics. We review constraint molecular dynamics for system with Nc constraints in which the canonical Hamiltonian equations can be written as [2430]

dH (q, p, X)

i = 1,2,..., 3N,

dH (q, p, X) dqi

= fr + ft = -

dV(q) % dgk (q) d1i k dcli

where m;, qt, and pt are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respectively. The gk(q) = (k= 1,2,... ,NC) are constraint equations determined by Nc 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 Nc Lagrange multipliers Xk(t) is given as

H(q,p,V = T +V(^) + 1 Xk9k (q),

where Xk(t) are to be determined. As gk(q) 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 qt and pt by the numerical integration of (1) and (2), and then the qt 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]:

miqi (t) = -T~

V(q) + TXk (t)9k (q)

Integrating both sides of (4) twice in time yields the constraint Cartesian coordinates at the time t + At:

qt (t + At) = quc (t + At) + mi1(At)2 f. (t)

= qW (t + At)-m:l(At)2 Yh

d9k (q)

where ^¡c(t + At) is solved by unconstrained force. The above qt(t + At) can be inserted to every type of the constraint equations to derive Nc nonlinear equations with respect to the corresponding Nc Lagrange multipliers. Although there are a number of algorithms to solve Xk(t), they only differ on how to solve these nonlinear equations. The SETTLE algorithm [28, 29] obtained the solutions of these nonlinear equations analytically for Nc = 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 gk = 0:

' d01 d01 d01

dX 1 dX 2 dX Nc

d02 d02 d02

dX 1 9X 2 ' dX Nc

d0Nr d0Nr

dX 1 dX 2

where Xk are then updated repeatedly by solving the system of linear equations iteratively:

until the max\gk\ < t (k = 1,..., Nc) in which t 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.

2.2. Two-Dimensional Potential Energy Surfaces for Stilbene. Quenneville and Martinez [45] have done extensive calculation for the critical SXIS0 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 n and n* 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 S0 and 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 DDI = (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 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 DDI is chosen as in [0°, -180°] that is sufficient to describe the present isomerization process. We use equal spacing 5° for grid points: DDI is 0, -5,..., -180 (37gridpoints) andDD2 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:

W(X' y) = c0 exp [-a1{x - x0)2 -b1{y- y0f ]

+ q cos + c2 cos + c3 cos (x)

+ c4 cos (y) + c5 cos (2x) + c6 cos (2y)

t \ \x + y

+ Cj cos (3x) + cs cos (3y) + c9 sin

'2x + y " x + 2y

+ c11 sin

+ c10 sin

+ c12 sin (x + y) + c13 sin (2x + y) + c14 sin (x + 2y) + c15 sin (2x + 2y) + c16 sin (3x + y) + c17 sin (x + 3y) + c18 sin (3x + 2y) + c19 sin (2x + 3y)

. , [X + y

+ c20 sin (3x + 3y) + c21 cos —2—

'2x + y] \x + 2y - + c23 cos -

. 2 . 23 -2

+ Cn cos

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

+ C24 cos (x + y) + C25 cos (2x + y) + C26 cos (x + 2y) + C27 cos (2x + 2y) + C28 cos (3x + y) + C29 cos (x + 3y) + c30 cos (3x + 2y) + c31 cos (2x + 3y) + c32 cos (3x + 3y),

where x = DDI and y = DD2, x0 = -103.3° and y0 = 43.3° are angles at OBF-CI. The fitting parameters ax, fy, and c0 to c32 are all given in Table 1. Potential energy surfaces of S0 and Si calculated from (8) are plotted in Figure 2 with zero point of potential energy at W(x0, y0) = 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 DDI = -164.3° and DD2 = 36.7° with energy 1.72 eV below OBF-CI, and local cis-isomer at DDI = -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 Si state is just right at conical intersection for the present two-dimensional 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

Table 1: The coefficients in (1) for two-dimensional S0 and Sj PESs around one-bond flip conical intersection (OBF-CI) (co ~ c32 in units of eV) [33].

Index a1 61 Ç) Cj ¿2

S0 35 3.5 0.495231 364.2436 166.2537

Si 39 3.9 -0.52219 256.2847 -18.8715

Index C3 C4 C5 C6 C7

S0 -63.2233 -228.305 -12.7773 25.79321 -3.78713

Si -74.255 -85.6398 -4.51131 16.61385 -1.65933

Index C8 G, C10 C11 C12

S0 -3.97288 9.467412 106.3829 -36.167 -72.9472

Si -2.44485 -219.158 171.3529 98.32809 -110.789

Index C13 C14 C15 C16 C17

S0 5.205437 -0.08421 -4.58729 5.134825 3.857498

Si 9.156298 -4.96308 -1.1299 0.059523 4.654567

Index C18 C19 C20 C21 C22

S0 -6.20916 0.171133 2.090849 -538.156 -56.4412

Si -1.20414 -1.73863 0.658395 -277.529 87.56788

Index C23 C24 ¿25 C26 C27

S0 254.7531 157.3086 17.19558 -58.19 -8.161

Si 112.7934 47.53832 -7.42468 -41.5302 7.965267

Index C28 C30 C31 C32

S0 7.267131 4.132838 -3.08607 1.734321 0.088946

Si 1.25399 3.375719 -0.42815 -0.71119 -0.09718

H 5 6 D2 H . a = r2 - r

b= r 3 - r 2

c= r4- r3 —> _^ _^

d= r 2 — r 5

e= r6- r3

Scheme 1

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

ft2 Ô2

ft2 Ô2

2ID1 3D12 2ID2 ÔD22'

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 z-axis as shown in Scheme 1. For instance, we can think that 01 and 02 are rotation angles around z-axis for two phenyl rings, respectively. Thus, the moment of inertia ID1 is the reduced inertia from two phenyl rings

J®1 + I

where J01 = m,- p2 (m,- is atom mass in the one phenyl ring and p{ is distance from atom i to z-axis) and I02 is estimated

from the other phenyl ring. The same analysis can be applied for the moment of inertia ID2 with two hydrogen atoms:

IH\IH2 Ih1 + IH2

where IH1 = mHp'^ (mH is hydrogen atom mass and pH is distance from hydrogen atom H to z-axis) and IH2 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):

Plz = exp i

4^a2 \b2

PZN = exp I

4^a2 \b2\ \ 1 + + (0.4a2 + 0.7) b-

where calculation of two parameters a2 and b2 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 DDI = -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 a2 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 b2 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 DDI = -103.3° with negligible dependence of DD2 angle. The effective coupling parameter a2 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 < a2 < 10.

3. 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)

20 30 40 50 60 70 80 DD2 (°)

Figure 3: The effective coupling parameter a2 along the seam line with DD1 = -103.3°.

potential energy surfaces at OBF-CI, and S0 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 S0 to Sj PES and then reaches the certain area of the present local cis (trans)-minimum on the 2D Sj 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 S0 to Sj 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 S0 to Sj 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 Sj state can have three ending ways; one is that it enters local cis-area on S0 (called cis-to-cis nonreactive trajectory), another is that it enters local-trans area on S0 (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 Tinit 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 Tinit, we have the following relations:

1 ( -2 -2 -2 -2\ rp c,

^mj (%! + ^ + x4 + ) = TinitS,

1m5 + >5 + + ¿e) = Tinit (1 - <5) ,

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).

OBF-CIa Initial from cis- area Initial from trans- area

Number of trajectories Cis-area (%) Trans-area (%) Trans-area (%) Cis-area (%)

100 46 (45) 54 (55) 46 (46) 54 (54)

500 35.4 (35.4) 64.6 (64.6) 47 (50) 53 (50)

1000 33 (36) 67 (64) 47 (49) 52.6 (50.7)

2000 33.4 (37) 66.6 (63) 47.8 (48.4) 52.2 (51.6)

Exp.b 50 50 50 50

aTotal 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. bReferences [59, 61].

where S is random number in the range of [0,1] and m1 (m5) is mass of carbon (hydrogen) atom. The x and y-components of velocities in (14) must satisfy

x1x1 + y1y1 = 0, x4x4 + y4y4 = 0,

X5X5 + 75% = 0 X6X6 + y&y& = 0

because of constraints of bond lengths and bond angles. Moreover, we choose the angular moment along z-axis for atoms 1 and 4 and for atoms 5 and 6 as follows:

(x1^1 - y1Xl) + (x4^4 - y4X4) = °

(X5^5 - y5X5) + (x6^6 - y6X6) =

This means that angular momentum along z-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 energybecause 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 knetic energy Tinit 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-to-cis isomerization, and the total energy is defined by the excitation energy 4.11 eV from S0 to S1 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 S0 to S1 state. Three types of ending trajectories are also collected as mentioned above (trans-to-trans nonreactive, trans-to-cis 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 Ak(t) = 0. Then, outcome of coordinates is inserted into the Jacobian determinant Jkki in (7) to have the new Xk(t) 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 a2 and effective collision energy b2 according to the method given in Appendix C. Then, we can evaluate one-passage nonadiabatic transition probability by (12) and (13). We obtain the seam line by calculating local minima of effective coupling a2 and found that they all are located at avoided crossing points. Comparing p with random number in between 0 and 1, we decide to hop trajectory from one adiabatic potential energy surface to another if p 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 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 trans-area. 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 trans-to-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 x 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 S0 to 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 x 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 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 b2 is well above the seam and this makestrajectoryhaveabigchancetohop down S0 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 b2 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-to-cis 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 small-ness 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 (a2 and b2) 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 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 state and increases on S0 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 S0 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-to-cis isomerization varies with time by starting at local trans-area with photo-excitation energy 3.98 eV and evolution time is 133 fs. Potential energies on S0 and 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 to S0 state happened at 50 fs, the second hopping from S0 to state at 71 fs, and the last hopping from to S0 state at 88 fs. After that, the trajectory is propagated on S0 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 cis-to-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 cis-area 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 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

ад 3

f 0 Рч

--- S,

40 60 80 100 120

Time (fs)

---DD2

60 Time (fs)

100 120

-a -100

-150 -200

0 20 40 60 80 100 120

Time (fs)

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.

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.

4. 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 a2 and b2 along trajectory moving on adiabatic potential energy surfaces. Effective coupling parameter a2 is not only used for determining transition probability but also for determining the seam line; in the present two-dimensional case this is

Time (fs) Time (fs)

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) DDI and DD2, and (c) D1 and D2 varying with time.

Time (fs) Time (fs)

Figure 6: The dihedral angels DDI and DD2 vary with time for 100 trajectories; (a) and (b) correspond to initially start from cis-area and trans-area, respectively.

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 a2 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-to-cis 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.

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 d^k2, bond angles ^k1k2k3, and dihedral angles ^k1k2k3ki in the following equations, respectively:

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

q with the derivations of the constraint equations, equation (5) can be rewritten as

= ?uc (t + At)-mkl(Atf

dgbk, (t)

V jk m ufk' (t)^ (0

!K> (t)

k'=1 N

a ^W' (t)

+ 1* i»dgi m

k'=1 N

dqkn (t)

in which

= sk'k,

dè (t) * dgk' «)

+ sk'k«

dqk[ (t) ^dqk, (t)

2 (<Vi- %,k2 ) (qk> -qk2 ),

^ = ^dq« (t) W'dqk, (t)

'k'kfl

dgQ (t) dq» (t)

= I Sk'kn - Sk'kf

)d_gkL

+ ISk'skls -Sk'2k^ dva'

0k(q) = ql1k2 -^klk2 = 0, k = \,...,Nb, (a.i) with

0k = cos0k1k2k3 \v\\v \ ■ v =0, k=1,

0\ | a a

kik2k, - V ,

d i I \^d\

0k = cos0k1k2k3k4 M lv I

d d ■ft ■ v = 0,

k = i,...,Nd

where Cartesian coordinates are qkik2 = (?k2-?kl) , H" = ?kl -

rk , = rk2-rkj ,ptd = axb, vd = bxc,fid ■ vd = (axb)-(bxc),

a = ?2-?i, b = ?3-?2 ,andc = r4-r3 and the subscript of each variable is atom number. In terms of the Cartesian coordinate

fl №.....

" cos0kik2k3 M V ,

S ÇgUt) S dgU^

dqh~ klk'dqk[ (t) k'2k?dqk, (t) , 3gd' (t) . dgd (t)

+ Sk3k- 3q~W) +Sk4k> sqTW)

q.n,k'3k'2

qm,k'k'

1njc[k'3 -qm,k'k'

dvm dvd

Inkc'kl

-qm,k'k'

Sk'kf Sk'kB

diabatic potential energy surfaces. The adiabatic parameter S that measures nonadiabatic transition strength is expressed in terms of the complex phase integrals:

5 = Im W [ ^2m(E-W~(Ry)

h JBo L >

2m (E-W+ (R))

qn,k[k'3 qn,k'2k'4 qn,k'3k'2 ~qm,k'4k'3 -qm,k'2k[ -qm,kl3k'2 j

Qlk'ak'p = qlk'« - ql,k'p '

(Imn) = (123), (231), (312) = (xyz), (yzx), (zxy).

For dihedral angles, we have

drf ~^Ykik2k,ki\^

Sk2kp Sk'2kp Sk'kR

= cos 0k1k2k3k4

m m d\^l - Vl>

¿9m h M m m

df = cos ^kik'k3k4 ]i ■

Finally, Jacobian determinants in (7) are given as

Jkk' =

dgpk (t + At)

t=1l=\

i=1l=1

dgk (t + At) dqUi (t + At) dqu (t + At) dXp,' (t)

dgk (t + At) dgk' (t) dqu (t + At) dqu (t)

where N = total atom number and I = x, y, z.

in which W+(R*) = W_(R*) at the complex crossing point R* (its real part R0 stands for transition point), m is reduced mass of system, and E 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 R0. Applying the exponential model (see (2.20) of [41]), we can convert (B.3) into

n^2m(AEo

h(d/dx)[W+(R) +W_(R)]\r=Ko

VE-Eo 1 + + AEo/(E - Eo)^d

W+ (Ro) - W- (Ro)

W+ (Ro) + W- (Ro)

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 S into the following two diabatic parameters:

B. Unified One-Passage Semiclassical Nonadiabatic Transition Probability

A unified one-passage semiclassical nonadiabatic transition probability was derived analytically as [73]

sinh [(d-1)S] P= sinh (dS) exp (-5)'

where d represents the ratio between adiabatic and diabatic potential energy gaps at transition point R0:

W+ (Ro)-W- (Ro) L ^2 (Ro)-V! (Ro)

where W+ and W- correspond to the excited state (upper state) and the ground state (lower state), while V2 and V1 are

[W+ (R) + W- (R)]\

b2 = E-Eo

11 - +

(E-Eo)

(B.6) (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 a2 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 b2 in (B.7) has an finite value at E = Eo. The parameter d 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 a2 and the Seam Line. The effective coupling parameter a2 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 a2, 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

2ID1 9D12 2ID2 dD22

~2I (Эх2 + ду2

where I = ^ID1ID2 is the reduced moment of inertia that corresponds to the reduced mass required in (B.6). The variables x and y are given by

q 0 0 c2

0 aI-D2

(D1 (D2

At conical configuration point (where D10 = -146.6° and D20 = -60.0°), we estimate from (10) and (11) that ID2 = 1342893.5 a.u. and ID1 = 3453.3 a.u. so that I = VWrn = 68098.6 a.u.

Let us start from conical configuration point (D10, D20) as reference point for the seam and they are transferred to variables x0 and y0 (see in Figure 8) by (C.2), and then we assume that x1 and y1 are a predefined point on the seam with s direction normal to distance direction between (x0, y0) and (x1,y1); thus the point (x,y) along s coordinate can be expressed as

x = s cos (в + + х1г y = s sin (в+^)+У1-

Derivative of potential energy surfaces defined in (B.6) with respect to s at (x1, y1) can be estimated as

dW±(x,y)

cos [в + — 2

sin [в + — 2

This equation has to be transformed to the original variables D1 and D2 in terms of which the potential energy surfaces are given by

9W±(D1, D2)

ds dW,

д M dW± 0 + — ) +

sin [в + — 2

where q and c2 are defined in (C.2), and в is given as bi - У0) C2 (D2 - D20)

(x1 -x0) q (D1- D10) '

Finally, we can have an expression for effective coupling parameter

[W+ + W_

161 (AE0)3 tg29 + 1 1 ( dW+ dW_

qV3D1 9D1

1 fdW, dW_\

c2 \ 9D2 9D2 ,

in which

АЕ0 =

W+ (D1, D2)-W_ (D1, D2) ~2 '

When 0 is rotating with respect to the conic point (x0,y0), the local maxima a2 can be obtained. In this way, we can determine one point on seam surface as well as a 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 a2 distribution on seam surface simultaneously.

C.2. The Effective Collision Energy b2 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

* s viV ■ + \ (xi,yi)

V 2 / e so Ni—^ e + w/2

O (xo,yo) \

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

hopping direction can be utilized. Let us write down classical kinetic energy in terms dihedral angles D1 and D2 as

1 ■ 9 1 -9

T=2Id1(D1)2 + -ID9 (D2)9

= +y2) = \l{Vl + V2),

where dot stands for derivative with respect to time t, and relation between (Dl, D2) and (x, y) is given in (C.2). The last equality in (C.9) is defined as

- sin 9 cos 0\ / cos 9 sin 0 )\y)'

(C.10)

where v1 is angular velocity along s direction and v2 is angular velocity along the direction normal to s as shown in Figure 8. If we assume that (v1y v2) change to (v\, v2 ) after hopping, then they should satisfy

v'1 = kv1,

(C.11)

where v2 is unchanged and k can be estimated from energy conservation as

4AE0 (iv2)

(C.12)

with AE0 = (W+ - W_)/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, k > 0 since the collision energy is quite large. We now need to convert (v\, v'2) to (x ,y ) that represent angular velocity change after hop along x and y direction shown in Figure 8. Equations (C.10) and (C.11) lead to the following relation:

y cos9 - x sin9 = k(y cos9 - xsin9), x cos 9 + y sin 9 = x cos 9 + y sin 9,

(C.13)

from which we can obtain

cos29 + k sin20 (1 - k) sin 9 cos 9\(x

(1-k) sin 9 cos 9 k cos20 + sin29 )\y

(C.14)

This can be converted to original dihedral angles D1 and D2 system as

^11 ^12 ^21 ^22

cos29 + k sin20 c92 (1-k) sin 9 cos 9 c2 (1-k) sin 9 cos 9 k cos2 9 + sin20

(C.15)

where (D1', D2') stand for angular velocity after hop in original dihedral angle system, and c1 and c2 are defined in (C.2). The effective collision energy for b2 in (B.7) is then given by collision energy

E=-I(v29) + W±

(C.16)

= 2I(C9D2 cos d-c1D1 sin 9) +W±.

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

cos (D1) - v = 0,

(C.17)

where fi = a xb and v = b xc (see Scheme 1 for definition of bond vectors a, b, and c). Taking derivative with respect to time t in (C.17) with consideration that all bond angles and bond distances are fixed, we have

D1 = -

PI sin (D1)]-1, (C.18)

[(a-b)(b-c)-(a-c){b-b)] 2 d(a-c)

dt (C.19)

= -\t>\ [(?2 -r1)-c + a-(r4-r3)]

= -p\ (¿2 -ri +rA - r3) -(c + a)

in which we utilized that all bond angles and bond distances are fixed again. Thus, we obtain

D1 = A - (r2 -r1 +r4 -r3),

(C.20)

(c + fl) jfcj

A = ' ^ '-. (C.21)

I a x fcI I fc x c I sin (D1)

In the same way as we derived (C.20), we can have

D2 = B-(?2 -?5 + ?6 -?3), (C.22)

where (see Scheme 1 for definition of bond vectors e and d)

(e + d) jfcj

je x fc j j d x fcj sin (D2)

(C.23)

Assume that velocities after hop in Cartesian coordinate

4' 4' 4' 4' 4' 4'

system are rp r2, r3, r4, r5, and r6, and then they are related to angular velocities after hop as

•>' •>' •>'

D1 = A • ( r2 - rj + r4 - r3),

, -+' 4' 4' 4'

D2 = ^ • ( r2 -r5+r6 -r3

(C.24)

in which A 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

A • ( r' - r,) + A • - r.

= rnA-[(?2 -?3) + (*4 -?i)]

+ Tl25- [(?2 -?3) + (?6 -?5)]

(?6 - ?5) + £• (r2 -r'

= T2iA^[(f2 -?3) + (r4 -?1)] + T225- [(?2 -?3) + (?6 -?5)]

where T^- 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

?4 - r'i = «j (?4 - ?i) x fö - rj + fcj (?4 - ?i) ,

?6 - ?5 = a2 (?6 - ?5) X (?6 - ) + fc2 (?6 - ?5) > (C.26)

?2 - ?3 = a3 (?2 - ?3) X (?2 - ?3) + fc3 (?2 - ?3) •

The first assumption is that all is equal to zero. Then, the relative motion (r3 - r2) does not change after hop (fc3 = 1)

(C.27)

as the rotation is around bond 2-3 as shown in Scheme 1. Therefore, we can derive

fci = (((Tn -1)A + Ti2 B)-(?2 -^3)

+ TnA-(?4 -?1 ) + Ti2B-(?6 -^5))

x(A-(?4 -?1)) 1,

^2 = ((^21 A+(T22 -l)5)-(?2 -^3)

+ T21A • (?4 - ?1 ) + T22B • (?6 - ?5) )

x(5-(?6 -?5))_1-

4' 4 4' 4 4'

Finally, we can assume after hop that r1 = fc1r1, r2 = r2, r3 =

^ = ^ = and = fc2?6.

The moment of inertia I = y/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

4' 44' 44' 4

atoms 1, 4, 5, and 6 as r1 = Afc1r1, r4 = Afc1r4, r5 = Afc2r5,

and r6 = Afc2 r6 to restore energy conservation in which A = Vl - Ae/T, where Ae is energy lose and T 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 A is close to unity from numerical calculation during surface hopping (A is exact unity if hopping is taking place right at conic point).

(C.25) 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. 100-2811-M-009-004. This work is supported by National Science Council of China under Grant no. 100-2113-M-009-005-MY3. 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 MOE-ATU Project of the National Chiao Tung University for support.

References

[1] J. C. Tully and R. K. Pkeston, "Trajectory surface hopping approach to nonadiabatic molecular collisions: the reaction of H+ with D2" The Journal of Chemical Physics, vol. 55, no. 2, pp. 562-572,1971.

[2] J. C. Tully, "Molecular dynamics with electronic transitions," The Journal of Chemical Physics, vol. 93, no. 2, pp. 1061-1071,1990.

[3] C. Zhu and H. Nakamura, "Theory of nonadiabatic transition for general two-state curve crossing problems. II. Landau-Zener case," The Journal of Chemical Physics, vol. 102, no. 19, pp. 74487461,1995.

[4] C. Zhu, H. Kamisaka, and H. Nakamura, "New implementation of the trajectory surface hopping method with use of the Zhu-Nakamura theory. II. application to the charge transfer processes in the 3D DH+ system," Journal of Chemical Physics, vol. 116, no. 8, pp. 3234-3247, 2002.

[5] E. Tapavicza, I. Tavernelli, U. Rothlisberger, C. Filippi, and M. E. Casida, "Mixed time-dependent density-functional theory/classical trajectory surface hopping study of oxirane photochemistry," The Journal of Chemical Physics, vol. 129, no. 12, Article ID 124108, 2008.

[6] M. Sanier, U. Manthe, and G. Stock, "Quantum-classical Liou-ville description of multidimensional nonadiabatic molecular dynamics," The Journal of Chemical Physics, vol. 114, no. 5, pp. 2001-2012, 2001.

[7] I. Burghardt and G. Parlant, "On the dynamics of coupled Bohmian and phase-space variables: a new hybrid quantum-classical approach," The Journal of Chemical Physics, vol. 120, no. 7, pp. 3055-3058, 2004.

[8] P. Pechukas, "Time-dependent semiclassical scattering theory. II. atomic collisions," Physical Review, vol. 181, no. 1, pp. 174-185, 1969.

[9] 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.

[10] M. Ben-Nun, J. Quenneville, and T. J. Martinez, "Ab initio multiple spawning: photochemistry from first principles quantum molecular dynamics," The Journal ofPhysical Chemistry A, vol. 104, no. 22, pp. 5172-5175, 2000.

[11] A. N. Alexandrova, J. C. Tully, and G. Granucci, "Photochemistry of DNA fragments via semiclassical nonadiabatic dynamics," The Journal ofPhysical Chemistry B, vol. 114, no. 37, pp. 12116-12128, 2010.

[12] B. F. E. Curchod, T. J. Penfold, U. Rothlisberger, and I. Tavernelli, "Local control theory in trajectory-based nonadiabatic dynamics," Physical Review A, vol. 84, no. 4, Article ID 042507, 2011.

[13] R. Crespo-Otero and M. Barbatti, "Cr(CO)6 photochemistry: semi-classical study of UV absorption spectral intensities and dynamics of photodissociation," The Journal of Chemical Physics, vol. 134, no. 16, Article ID 164305, 2011.

[14] 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.

[15] Z. Lan, E. Fabiano, and W. Thiel, "Photoinduced nonadia-batic dynamics of pyrimidine nucleobases: on-the-fly surface-hopping study with semiempirical methods," The Journal of Physical Chemistry B, vol. 113, no. 11, pp. 3548-3555, 2009.

[16] O. Weingart, Z. Lan, A. Koslowski, and W. Thiel, "Chiral pathways and periodic decay in cis-azobenzene photodynamics," The Journal ofPhysical ChemistryLetters, vol. 2, no. 13, pp. 15061509, 2011.

[17] 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. 1-2, pp. 145156, 2009.

[18] V. Leyva, I. Corral, F. Feixas et al., "A non-adiabatic quantum-classical dynamics study of the intramolecular excited state hydrogen transfer in ortho-nitrobenzaldehyde," Physical Chemistry Chemical Physics, vol. 13, no. 32, pp. 14685-14693, 2011.

[19] I. Tavernelli, B. F. E. Curchod, and U. Rothlisberger, "Mixed quantum-classical dynamics with time-dependent external fields: a time-dependent density-functional-theory approach," Physical Review A, vol. 81, no. 5, Article ID 052508, 2010.

[20] 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. Koppel, Eds., World Scientific, Singapore, 2010.

[21] 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.

[22] A. Gao, B. Li, P. Zhang, and K. Han, "Nonadiabatic ab initio molecular dynamics of photoisomerization in bridged azoben-zene," The Journal of Chemical Physics, vol. 137, no. 20, Article ID 204305, 2012.

[23] A. Gao, B. Li, P. Zhang, and J. Liu, "Photochemical dynamics simulations for trans-cis photoisomerizations of azobenzene and bridged azobenzen," Computational and Theoretical Chemistry, vol. 1031, pp. 13-21, 2014.

[24] 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 n-alkanes," Journal of Computational Physics, vol. 23, no. 3, pp. 327-341,1977.

[25] V. Krautler, W. F. van Gunsteren, and P. H. Hunenberger, "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.

[26] G. Ciccotti and J.-P. Ryckaert, "Molecular dynamics simulation of rigid molecules," Computer Physics Reports, vol. 4, no. 6, pp. 346-392, 1986.

[27] P. Gonnet, "P-SHAKE: a quadratically convergent SHAKE in O(n2)," Journal of Computational Physics, vol. 220, no. 2, pp. 740-750, 2007.

[28] S. Miyamoto and P. A. Kollman, "Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models," Journal ofComputational Chemistry, vol. 13, no. 8, pp. 952-962, 1992.

[29] 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.

[30] P. E. Gill, W. Murray, and M. H. Wright, Numerical Linear Algebra and Optimization, Addison-Wesley, Redwood City, Calif, USA, 1991.

[31] T. Chu, Y. Zhang, and K. Han, "The time-dependent quantum wave packet approach to the electronically nonadiabatic processes in chemical reactions," International Reviews in Physical Chemistry, vol. 25, no. 1-2, pp. 201-235, 2006.

[32] K. Han, J. Huang, S. Chai, S. Wen, and W. Deng, "Anisotropic mobilities in organic semiconductors," Nature Protocol Exchange, 2013.

[33] Y. Lei, C. Zhu, Z. Wen, and S. Lin, "New implementation of semi-classical dynamic simulation on the photoisomerization of cis- and trans- isomers of free stilbene," Acta Chimica Sinica, vol. 70, no. 17, pp. 1869-1876, 2012.

[34] A. J. Neukirch, L. C. Shamberger, E. Abad et al., "Nonadiabatic ensemble simulations of cis-stilbene and cis-azobenzene pho-toisomerization," Journal ofChemical Theory and Computation, vol. 10, no. 1, pp. 14-23, 2014.

[35] Y. Ootani, K. Satoh, A. Nakayama, T. Noro, and T. Taketsugu, "Ab initio molecular dynamics simulation of photoisomeriza-tion in azobenzene in the n n state," The Journal of Chemical Physics, vol. 131, no. 19, Article ID 194306, 2009.

[36] A. Toniolo, C. Ciminelli, M. Persico, and T. J. Martinez, "Simulation of the photodynamics of azobenzene on its first excited state: comparison of full multiple spawning and surface hopping treatments," The Journal ofChemical Physics, vol. 123, no. 23, Article ID 234308, 2005.

[37] P. Cattaneo and M. Persico, "An abinitio study of the photochemistry of azobenzene," Physical Chemistry Chemical Physics, vol. 1, pp. 4739-4743, 1999.

[38] M. Pederzoli, J. Pittner, M. Barbatti, and H. Lischka, "Nonadia-batic molecular dynamics study of the cis-trans photoisomer-ization of azobenzene excited to the Sj state," The Journal of Physical Chemistry A, vol. 115, no. 41, pp. 11136-11143, 2011.

[39] 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. 1218, 1968.

[40] 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.

[41] A. Gilbert and J. Baggott, Essentials of Molecular Photochemistry, Blackwell Science, Oxford, UK, 1991.

[42] R. S. H. Liu, "Photoisomerization by Hula-Twist: a funda-mentalsupramolecular photochemical reaction," Accounts of Chemical Research, vol. 34, no. 7, pp. 555-562, 2001.

[43] R. D. Sampedro, A. Cembran, M. Garavelli, M. Olivucci, and W. Fu^, "Structure of the conical intersections driving the cis-trans photoisomerization of conjugated molecules," Photochemistry and Photobiology, vol. 76, no. 6, pp. 622-633, 2002.

[44] Y. Amatatsu, "Ab initio study on the electronic structures of stilbene at the conical intersection," Chemical Physics Letters, vol. 314, no. 3-4, pp. 364-368,1999.

[45] J. Quenneville and T. J. Martinez, "Ab initio study of cis-trans photoisomerization in stilbene and ethylene," The Journal of Physical Chemistry A, vol. 107, no. 6, pp. 829-837, 2003.

[46] M. J. Bearpark, F. Bernardi, S. Clifford, M. Olivucci, M. A. Robb, andT. Vreven, "Cooperatingrings in cis-stilbenelead to anS0/S1 conical intersection," The Journal ofPhysical Chemistry A, vol. 101, no. 21, pp. 3841-3847,1997.

[47] J. H. Frederick, Y. Fujiwara, J. H. Penn, K. Yoshihara, and H. Petek, "Models for stilbene photoisomerization: experimental and theoretical studies of the excited-state dynamics of 1,2-diphenylcycloalkenes," The Journal of Physical Chemistry, vol. 95, no. 7, pp. 2845-2858,1991.

[48] 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. 306314, 1993.

[49] 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.

[50] 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. 2-3, pp. 383-391, 2002.

[51] W.-V. Chiang and J. Laane, "Fluorescence spectra and torsional potential functions for trans-stilbene in its S0 and electronic states," The Journal ofChemical Physics, vol. 100, no. 12, pp. 8755-8767,1994.

[52] C. W. Jiang, R. H. Xie, F. L. Li, and R. E. Allen, "Trans-to-cis isomerization of stilbene following an ultrafast laser pulse," Chemical Physics Letters, vol. 474, no. 4-6, pp. 263-267, 2009.

[53] Y. Dou and R. E. Allen, "Detailed dynamics of a complex photochemical reaction: cis-trans photoisomerization of stilbene," The Journal of Chemical Physics, vol. 119, no. 20, pp. 10658-10666, 2003.

[54] 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.

[55] W. Fuß, C. Kosmidis, W. E. Schmid, and S. A. Trushin, "The lifetime of the perpendicular minimum of cis-stibene observed by dissociative intense-laser field ionization," Chemical Physics Letters, vol. 385, pp. 423-430, 2004.

[56] W. Fuß, C. Kosmidis, W. E. Schmid, and S. A. Trushin, "The photochemical cis-trans isomerization of free stilbene molecules follows a hula-twist pathway," Angewandte Chemie International Edition, vol. 43, no. 32, pp. 4178-4182, 2004.

[57] D. C. Todd and G. R. Fleming, "Cis-stilbene isomerization: temperature dependence and the role of mechanical friction," The Journal of Chemical Physics,vol. 98, no. 1,pp. 269-279,1993.

[58] S. Abrash, S. Repinec, and R. M. Hochstrasser, "The viscosity dependence and reaction coordinate for isomerization of cis-stilbene," The Journal ofChemical Physics, vol. 93, no. 2, pp. 10411053, 1990.

[59] R. J. Sension, S. T. Repinec, A. Z. Szarka, and R. M. Hochstrasser, "Femtosecond laser studies of the cis-stilbene photoisomerization reactions," The Journal of Chemical Physics, vol. 98, no. 8, pp. 6291-6315,1993.

[60] 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 ofChemical Physics, vol. 93, no. 12, pp. 9185-9188, 1990.

[61] S. T. Repinec, R. J. Sension, A. Z. Szarka, and R. M. Hochstrasser, "Femtosecond laser studies of the cis-stilbene photoisomeriza-tion reactions: the cis-stilbene to dihydrophenanthrene reaction," The Journal of Physical Chemistry, vol. 95, no. 25, pp. 10380-10385, 1991.

[62] R. Improta and F. Santoro, "Excited-state behavior of trans and cis isomers of stilbene and stiff stilbene: a TD-DFT study," The Journal of Physical Chemistry A, vol. 109, no. 44, pp. 1005810067, 2005.

[63] C. D. Berweger, W. F. van Gunsteren, and F. Müller-Plathe, "Molecular dynamics simulation with an ab initio potential energy function and finite element interpolation: the photoiso-merization of cis-stilbene in solution," The Journal ofChemical Physics, vol. 108, no. 21, pp. 8773-8781,1998.

[64] J. Tatchen and E. Pollak, "Ab initio spectroscopy and photoin-duced cooling of the trans-stilbene molecule," The Journal of Chemical Physics, vol. 128, no. 16, Article ID 164303, 2008.

[65] 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.

[66] G. Karlstrom, R. Lindh, P.-A. Malmqvist et al., "MOLCAS: a program package for computational chemistry," Computational Materials Science, vol. 28, no. 2, pp. 222-239, 2003.

[67] Matlab, User Guide, The Mathworks, Torrance, Calif, USA, http://www.mathworks.com.

[68] J. K. Rice and A. P. Baronavski, "Ultrafast studies of solvent effects in the isomerization of cis-stilbene," The Journal of Physical Chemistry, vol. 96, no. 8, pp. 3359-3366,1992.

[69] C. Zhu, K. Nobusada, and H. Nakamura, "New implementation of the trajectory surface hopping method with use of the Zhu-Nakamura theory," The Journal of Chemical Physics, vol. 115, no. 7, pp. 3031-3044, 2001.

[70] J. J. Moré, B. S. Garbow, and K. E. Hillstrom, "User guide for MINPACK-1," Argonne National Laboratory Report ANL-80-74, Argonne, Ill, USA, 1980.

[71] M. J. D. Powell, Numerical Methods for Non-Linear Algebraic Equations, Gordon and Breach, New York, NY, USA, 1970, Edited by P. Rabinowitz.

[72] L. F. Shampine and H. A. Watts, Mathematical Software III, Academic Press, New York, NY, USA, 1977, Edited by J. R. Rice.

[73] C. Zhu, "Unified semiclassical theory for the two-state system: analytical solutions for scattering matrices," The Journal of Chemical Physics, vol. 105, no. 10, pp. 4159-4172,1996.

Copyright of International Journal of Photoenergy is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.