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

Authors: Yibo Lei; Shaomei Wu; Chaoyuan Zhu; Zhenyi Wen; Sheng-Hsien Lin
Journal: International Journal of Photoenergy
 (2014)

---

## Abstract

Combining trajectory surface hopping (TSH) method with constraint molecular dynamics, we have extended TSH method from full to flexible dimensional potential energy surfaces. Classical trajectories are carried out in Cartesian coordinates with constraints in internal coordinates, while nonadiabatic switching probabilities are calculated separately in free internal coordinates by Landau-Zener and Zhu-Nakamura formulas along the seam. Two-dimensional potential energy surfaces of groundS
0 and excited S
1 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%.


---

## Body

## 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 theD
H
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 D
H
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–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 nonadiabatic 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=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 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 stateS
0 and the lowest excited state S
1. The dominant excitation from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) corresponds to the S
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, 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 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 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 S
0 and the first excited state S
1 are constructed in terms of the D1 and D2 angles, and the other internal coordinates are fixed at the OBF-CI configuration.Figure 1
Atom numbering of Stilbene at theS
1
/
S
0 conical intersection (OBF-CI) configuration in terms of two dihedral angles D1 and D2 at which D10 = −146.6° and D20 = −60.0°.The rest of this paper is organized as follows. Section2 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 withN
c 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 withN
c constraints in which the canonical Hamiltonian equations can be written as [24–30]
(1)
q
˙
i
=
∂
H
(
q
,
p
,
λ
)
∂
p
i
=
p
i
m
i
,
i
=
1,2
,
…
,
3
N
,
(2)
p
˙
i
=
-
∂
H
(
q
,
p
,
λ
)
∂
q
i
=
f
i
uc
+
f
i
c
=
-
∂
V
(
q
)
∂
q
i
-
∑
k
N
c
λ
k
∂
g
k
(
q
)
∂
q
i
,
where m
i, q
i, and p
i are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respectively. The g
k
(
q
)
=
(
k
=
1,2
,
…
,
N
c
) are constraint equations determined by N
c internal coordinates that are fixed at the certain structures of nuclear configuration all the time. The f
i
uc and f
i
c in (2) stand for unconstraint and constraint force acting on the corresponding atom, respectively. The Hamiltonian with N
c Lagrange multipliers λ
k
(
t
) is given as
(3)
H
(
q
,
p
,
λ
)
=
∑
i
3
N
p
i
2
2
m
i
+
V
(
q
)
+
∑
k
N
c
λ
k
g
k
(
q
)
,
where λ
k
(
t
) are to be determined. As g
k
(
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 solveq
i and p
i by the numerical integration of (1) and (2), and then the q
i 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]:
(4)
m
i
q
¨
i
(
t
)
=
-
∂
∂
q
i
[
V
(
q
)
+
∑
k
=
1
N
c
λ
k
(
t
)
g
k
(
q
)
]
.
Integrating both sides of (4) twice in time yields the constraint Cartesian coordinates at the time t
+
Δ
t:
(5)
q
i
(
t
+
Δ
t
)
=
q
i
uc
(
t
+
Δ
t
)
+
m
i
-
1
(
Δ
t
)
2
f
i
c
(
t
)
=
q
i
uc
(
t
+
Δ
t
)
-
m
i
-
1
(
Δ
t
)
2
∑
k
N
c
λ
k
∂
g
k
(
q
)
∂
q
i
,
where q
i
uc
(
t
+
Δ
t
) is solved by unconstrained force. The above q
i
(
t
+
Δ
t
) can be inserted to every type of the constraint equations to derive N
c nonlinear equations with respect to the corresponding N
c Lagrange multipliers. Although there are a number of algorithms to solve λ
k
(
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 N
c
=
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 g
k
=
0:
(6)
J
g
=
[
∂
g
1
∂
λ
1
∂
g
1
∂
λ
2
⋯
∂
g
1
∂
λ
N
c
∂
g
2
∂
λ
1
∂
g
2
∂
λ
2
⋯
∂
g
2
∂
λ
N
c
⋮
⋮
⋮
⋮
∂
g
N
c
∂
λ
1
∂
g
N
c
∂
λ
2
⋯
∂
g
N
c
∂
λ
N
c
]
,
where λ
k are then updated repeatedly by solving the system of linear equations iteratively:
(7)
(
J
g
|
λ
k
=
0
)
λ
=
-
g
|
λ
k
=
0
until the max
⁡
|
g
k
|
<
τ
(
k
=
1
,
…
,
N
c
) 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.


### 2.2. Two-Dimensional Potential Energy Surfaces for Stilbene

Quenneville and Martínez [45] have done extensive calculation for the critical S
1
/
S
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-3
1
G
*, and 6-31
G
*
*. Following their conclusion, we can safely utilize state-averaged CASSCF method abbreviated as SA-N-CAS(n/m), whereN refers to the number of states included in the average, whileN(=2),n(=2), andm(=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 onS
0 and S
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 DD
1
=
(
D
1
+
D
2
)
/
2 and DD
2
=
(
D
2
-
D
1
)
/
2 are better variables to describe the isomerization. If we set up interesting region for potential energy surface S
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 allab. 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:
(8)
W
(
x
,
y
)
=
c
0
exp
⁡
[
-
a
1
(
x
-
x
0
)
2
-
b
1
(
y
-
y
0
)
2
]
+
c
1
cos
⁡
(
x
2
)
+
c
2
cos
⁡
(
y
2
)
+
c
3
cos
⁡
(
x
)
+
c
4
cos
⁡
(
y
)
+
c
5
cos
⁡
(
2
x
)
+
c
6
cos
⁡
(
2
y
)
+
c
7
cos
⁡
(
3
x
)
+
c
8
cos
⁡
(
3
y
)
+
c
9
sin
[
x
+
y
2
]
+
c
10
sin
[
2
x
+
y
2
]
+
c
11
sin
[
x
+
2
y
2
]
+
c
12
sin
(
x
+
y
)
+
c
13
sin
(
2
x
+
y
)
+
c
14
sin
(
x
+
2
y
)
+
c
15
sin
(
2
x
+
2
y
)
+
c
16
sin
(
3
x
+
y
)
+
c
17
sin
(
x
+
3
y
)
+
c
18
sin
(
3
x
+
2
y
)
+
c
19
sin
(
2
x
+
3
y
)
+
c
20
sin
(
3
x
+
3
y
)
+
c
21
cos
⁡
[
x
+
y
2
]
+
c
22
cos
⁡
[
2
x
+
y
2
]
+
c
23
cos
⁡
[
x
+
2
y
2
]
+
c
24
cos
⁡
(
x
+
y
)
+
c
25
cos
⁡
(
2
x
+
y
)
+
c
26
cos
⁡
(
x
+
2
y
)
+
c
27
cos
⁡
(
2
x
+
2
y
)
+
c
28
cos
⁡
(
3
x
+
y
)
+
c
29
cos
⁡
(
x
+
3
y
)
+
c
30
cos
⁡
⁡
(
3
x
+
2
y
)
+
c
31
cos
⁡
(
2
x
+
3
y
)
+
c
32
cos
⁡
(
3
x
+
3
y
)
,
where x
=
DD
1 and y
=
DD
2, x
0
=
-
103.3° and y
0
=
43.3° are angles at OBF-CI. The fitting parameters a
1, b
1, and c
0 to c
32 are all given in Table 1. Potential energy surfaces of S
0 and S
1 calculated from (8) are plotted in Figure 2 with zero point of potential energy at W
(
x
0
,
y
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 DD
1
=
-
164.3° and DD
2
=
36.7° with energy 1.72 eV below OBF-CI, and local cis-isomer at DD
1
=
-
39.2° and DD
2
=
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 S
1 state is just right at conical intersection for the present two-dimensional potential energy surface.Table 1
The coefficients in (1) for two-dimensional S
0 and S
1 PESs around one-bond flip conical intersection (OBF-CI) (c
0
~
c
32 in units of eV) [33].
Index
a
1
b
1
c
0
c
1
c
2
S
0
35
3.5
0.495231
364.2436
166.2537
S
1
39
3.9
−0.52219
256.2847
−18.8715
Index
c
3
c
4
c
5
c
6
c
7
S
0
−63.2233
−228.305
−12.7773
25.79321
−3.78713
S
1
−74.255
−85.6398
−4.51131
16.61385
−1.65933
Index
c
8
c
9
c
10
c
11
c
12
S
0
−3.97288
9.467412
106.3829
−36.167
−72.9472
S
1
−2.44485
−219.158
171.3529
98.32809
−110.789
Index
c
13
c
14
c
15
c
16
c
17
S
0
5.205437
−0.08421
−4.58729
5.134825
3.857498
S
1
9.156298
−4.96308
−1.1299
0.059523
4.654567
Index
c
18
c
19
c
20
c
21
c
22
S
0
−6.20916
0.171133
2.090849
−538.156
−56.4412
S
1
−1.20414
−1.73863
0.658395
−277.529
87.56788
Index
c
23
c
24
c
25
c
26
c
27
S
0
254.7531
157.3086
17.19558
−58.19
−8.161
S
1
112.7934
47.53832
−7.42468
−41.5302
7.965267
Index
c
28
c
29
c
30
c
31
c
32
S
0
7.267131
4.132838
−3.08607
1.734321
0.088946
S
1
1.25399
3.375719
−0.42815
−0.71119
−0.09718Figure 2
Analytical two-dimensional PESs for the ground stateS
0 and the first excited state S
1 with respect to the combined internal angles DD1 and DD2 [33].


### 2.3. TSH Scheme for LZ and ZN Methods

Section2.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
(9)
T
=
-
ℏ
2
2
I
D
1
∂
2
∂
D
1
2
-
ℏ
2
2
I
D
2
∂
2
∂
D
2
2
,
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 Φ
1 and Φ
2 are rotation angles around z-axis for two phenyl rings, respectively. Thus, the moment of inertia I
D
1 is the reduced inertia from two phenyl rings
(10)
I
D
1
=
I
Φ
1
I
Φ
2
I
Φ
1
+
I
Φ
2
,
where I
Φ
1
=
∑
i
m
i
ρ
i
2 (m
i is atom mass in the one phenyl ring and ρ
i is distance from atom i to z-axis) and I
Φ
2 is estimated from the other phenyl ring. The same analysis can be applied for the moment of inertia I
D
2 with two hydrogen atoms:
(11)
I
D
2
=
I
H
1
I
H
2
I
H
1
+
I
H
2
,
where I
H
1
=
m
H
ρ
H
2 (m
H is hydrogen atom mass and ρ
H is distance from hydrogen atom H to z-axis) and I
H
2 is estimated from the other hydrogen atom.Scheme 1One-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):
(12)
p
LZ
=
exp
⁡
(
-
π
4
a
2
|
b
2
|
)
,
(13)
p
ZN
=
exp
⁡
(
-
π
4
a
2
|
b
2
|
2
1
+
1
+
(
0.4
a
2
+
0.7
)
b
-
4
)
,
where calculation of two parameters a
2 and b
2 is given in Appendix B.We can simply look at the present two-dimensional PESs in Figure2 to find that the seam line is appearing at the DD
1
=
-
10
3.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 a
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 b
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 atDD
1
=
-
103.3
° with negligible dependence of DD2 angle. The effective coupling parameter a
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
≤
a
2
≤
10.Figure 3
The effective coupling parametera
2 along the seam line with DD1 = −103.3°.


## 2.1. Constraint Molecular Dynamics

We review constraint molecular dynamics for system withN
c constraints in which the canonical Hamiltonian equations can be written as [24–30]
(1)
q
˙
i
=
∂
H
(
q
,
p
,
λ
)
∂
p
i
=
p
i
m
i
,
i
=
1,2
,
…
,
3
N
,
(2)
p
˙
i
=
-
∂
H
(
q
,
p
,
λ
)
∂
q
i
=
f
i
uc
+
f
i
c
=
-
∂
V
(
q
)
∂
q
i
-
∑
k
N
c
λ
k
∂
g
k
(
q
)
∂
q
i
,
where m
i, q
i, and p
i are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respectively. The g
k
(
q
)
=
(
k
=
1,2
,
…
,
N
c
) are constraint equations determined by N
c internal coordinates that are fixed at the certain structures of nuclear configuration all the time. The f
i
uc and f
i
c in (2) stand for unconstraint and constraint force acting on the corresponding atom, respectively. The Hamiltonian with N
c Lagrange multipliers λ
k
(
t
) is given as
(3)
H
(
q
,
p
,
λ
)
=
∑
i
3
N
p
i
2
2
m
i
+
V
(
q
)
+
∑
k
N
c
λ
k
g
k
(
q
)
,
where λ
k
(
t
) are to be determined. As g
k
(
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 solveq
i and p
i by the numerical integration of (1) and (2), and then the q
i 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]:
(4)
m
i
q
¨
i
(
t
)
=
-
∂
∂
q
i
[
V
(
q
)
+
∑
k
=
1
N
c
λ
k
(
t
)
g
k
(
q
)
]
.
Integrating both sides of (4) twice in time yields the constraint Cartesian coordinates at the time t
+
Δ
t:
(5)
q
i
(
t
+
Δ
t
)
=
q
i
uc
(
t
+
Δ
t
)
+
m
i
-
1
(
Δ
t
)
2
f
i
c
(
t
)
=
q
i
uc
(
t
+
Δ
t
)
-
m
i
-
1
(
Δ
t
)
2
∑
k
N
c
λ
k
∂
g
k
(
q
)
∂
q
i
,
where q
i
uc
(
t
+
Δ
t
) is solved by unconstrained force. The above q
i
(
t
+
Δ
t
) can be inserted to every type of the constraint equations to derive N
c nonlinear equations with respect to the corresponding N
c Lagrange multipliers. Although there are a number of algorithms to solve λ
k
(
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 N
c
=
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 g
k
=
0:
(6)
J
g
=
[
∂
g
1
∂
λ
1
∂
g
1
∂
λ
2
⋯
∂
g
1
∂
λ
N
c
∂
g
2
∂
λ
1
∂
g
2
∂
λ
2
⋯
∂
g
2
∂
λ
N
c
⋮
⋮
⋮
⋮
∂
g
N
c
∂
λ
1
∂
g
N
c
∂
λ
2
⋯
∂
g
N
c
∂
λ
N
c
]
,
where λ
k are then updated repeatedly by solving the system of linear equations iteratively:
(7)
(
J
g
|
λ
k
=
0
)
λ
=
-
g
|
λ
k
=
0
until the max
⁡
|
g
k
|
<
τ
(
k
=
1
,
…
,
N
c
) 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.


## 2.2. Two-Dimensional Potential Energy Surfaces for Stilbene

Quenneville and Martínez [45] have done extensive calculation for the critical S
1
/
S
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-3
1
G
*, and 6-31
G
*
*. Following their conclusion, we can safely utilize state-averaged CASSCF method abbreviated as SA-N-CAS(n/m), whereN refers to the number of states included in the average, whileN(=2),n(=2), andm(=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 onS
0 and S
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 DD
1
=
(
D
1
+
D
2
)
/
2 and DD
2
=
(
D
2
-
D
1
)
/
2 are better variables to describe the isomerization. If we set up interesting region for potential energy surface S
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 allab. 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:
(8)
W
(
x
,
y
)
=
c
0
exp
⁡
[
-
a
1
(
x
-
x
0
)
2
-
b
1
(
y
-
y
0
)
2
]
+
c
1
cos
⁡
(
x
2
)
+
c
2
cos
⁡
(
y
2
)
+
c
3
cos
⁡
(
x
)
+
c
4
cos
⁡
(
y
)
+
c
5
cos
⁡
(
2
x
)
+
c
6
cos
⁡
(
2
y
)
+
c
7
cos
⁡
(
3
x
)
+
c
8
cos
⁡
(
3
y
)
+
c
9
sin
[
x
+
y
2
]
+
c
10
sin
[
2
x
+
y
2
]
+
c
11
sin
[
x
+
2
y
2
]
+
c
12
sin
(
x
+
y
)
+
c
13
sin
(
2
x
+
y
)
+
c
14
sin
(
x
+
2
y
)
+
c
15
sin
(
2
x
+
2
y
)
+
c
16
sin
(
3
x
+
y
)
+
c
17
sin
(
x
+
3
y
)
+
c
18
sin
(
3
x
+
2
y
)
+
c
19
sin
(
2
x
+
3
y
)
+
c
20
sin
(
3
x
+
3
y
)
+
c
21
cos
⁡
[
x
+
y
2
]
+
c
22
cos
⁡
[
2
x
+
y
2
]
+
c
23
cos
⁡
[
x
+
2
y
2
]
+
c
24
cos
⁡
(
x
+
y
)
+
c
25
cos
⁡
(
2
x
+
y
)
+
c
26
cos
⁡
(
x
+
2
y
)
+
c
27
cos
⁡
(
2
x
+
2
y
)
+
c
28
cos
⁡
(
3
x
+
y
)
+
c
29
cos
⁡
(
x
+
3
y
)
+
c
30
cos
⁡
⁡
(
3
x
+
2
y
)
+
c
31
cos
⁡
(
2
x
+
3
y
)
+
c
32
cos
⁡
(
3
x
+
3
y
)
,
where x
=
DD
1 and y
=
DD
2, x
0
=
-
103.3° and y
0
=
43.3° are angles at OBF-CI. The fitting parameters a
1, b
1, and c
0 to c
32 are all given in Table 1. Potential energy surfaces of S
0 and S
1 calculated from (8) are plotted in Figure 2 with zero point of potential energy at W
(
x
0
,
y
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 DD
1
=
-
164.3° and DD
2
=
36.7° with energy 1.72 eV below OBF-CI, and local cis-isomer at DD
1
=
-
39.2° and DD
2
=
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 S
1 state is just right at conical intersection for the present two-dimensional potential energy surface.Table 1
The coefficients in (1) for two-dimensional S
0 and S
1 PESs around one-bond flip conical intersection (OBF-CI) (c
0
~
c
32 in units of eV) [33].
Index
a
1
b
1
c
0
c
1
c
2
S
0
35
3.5
0.495231
364.2436
166.2537
S
1
39
3.9
−0.52219
256.2847
−18.8715
Index
c
3
c
4
c
5
c
6
c
7
S
0
−63.2233
−228.305
−12.7773
25.79321
−3.78713
S
1
−74.255
−85.6398
−4.51131
16.61385
−1.65933
Index
c
8
c
9
c
10
c
11
c
12
S
0
−3.97288
9.467412
106.3829
−36.167
−72.9472
S
1
−2.44485
−219.158
171.3529
98.32809
−110.789
Index
c
13
c
14
c
15
c
16
c
17
S
0
5.205437
−0.08421
−4.58729
5.134825
3.857498
S
1
9.156298
−4.96308
−1.1299
0.059523
4.654567
Index
c
18
c
19
c
20
c
21
c
22
S
0
−6.20916
0.171133
2.090849
−538.156
−56.4412
S
1
−1.20414
−1.73863
0.658395
−277.529
87.56788
Index
c
23
c
24
c
25
c
26
c
27
S
0
254.7531
157.3086
17.19558
−58.19
−8.161
S
1
112.7934
47.53832
−7.42468
−41.5302
7.965267
Index
c
28
c
29
c
30
c
31
c
32
S
0
7.267131
4.132838
−3.08607
1.734321
0.088946
S
1
1.25399
3.375719
−0.42815
−0.71119
−0.09718Figure 2
Analytical two-dimensional PESs for the ground stateS
0 and the first excited state S
1 with respect to the combined internal angles DD1 and DD2 [33].


## 2.3. TSH Scheme for LZ and ZN Methods

Section2.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
(9)
T
=
-
ℏ
2
2
I
D
1
∂
2
∂
D
1
2
-
ℏ
2
2
I
D
2
∂
2
∂
D
2
2
,
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 Φ
1 and Φ
2 are rotation angles around z-axis for two phenyl rings, respectively. Thus, the moment of inertia I
D
1 is the reduced inertia from two phenyl rings
(10)
I
D
1
=
I
Φ
1
I
Φ
2
I
Φ
1
+
I
Φ
2
,
where I
Φ
1
=
∑
i
m
i
ρ
i
2 (m
i is atom mass in the one phenyl ring and ρ
i is distance from atom i to z-axis) and I
Φ
2 is estimated from the other phenyl ring. The same analysis can be applied for the moment of inertia I
D
2 with two hydrogen atoms:
(11)
I
D
2
=
I
H
1
I
H
2
I
H
1
+
I
H
2
,
where I
H
1
=
m
H
ρ
H
2 (m
H is hydrogen atom mass and ρ
H is distance from hydrogen atom H to z-axis) and I
H
2 is estimated from the other hydrogen atom.Scheme 1One-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):
(12)
p
LZ
=
exp
⁡
(
-
π
4
a
2
|
b
2
|
)
,
(13)
p
ZN
=
exp
⁡
(
-
π
4
a
2
|
b
2
|
2
1
+
1
+
(
0.4
a
2
+
0.7
)
b
-
4
)
,
where calculation of two parameters a
2 and b
2 is given in Appendix B.We can simply look at the present two-dimensional PESs in Figure2 to find that the seam line is appearing at the DD
1
=
-
10
3.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 a
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 b
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 atDD
1
=
-
103.3
° with negligible dependence of DD2 angle. The effective coupling parameter a
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
≤
a
2
≤
10.Figure 3
The effective coupling parametera
2 along the seam line with DD1 = −103.3°.


## 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) potential energy surfaces at OBF-CI, andS
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 S
0 to S
1 PES and then reaches the certain area of the present local cis (trans)-minimum on the 2D S
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 S
0 to S
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 S
0 to S
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 S
1 state can have three ending ways; one is that it enters local cis-area on S
0 (called cis-to-cis nonreactive trajectory), another is that it enters local-trans area on S
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 energyT
init to atoms 1, 4, 5, and 6 with their initial velocities analyzed in the following. As shown in Scheme 1, thez-component of velocity is assumed to be zero so that for a given kinetic energy T
init, we have the following relations:
(14)
1
2
m
1
(
x
˙
1
2
+
y
˙
1
2
+
x
˙
4
2
+
y
˙
4
2
)
=
T
init
δ
,
1
2
m
5
(
x
˙
5
2
+
y
˙
5
2
+
x
˙
6
2
+
y
˙
6
2
)
=
T
init
(
1
-
δ
)
,
where δ is random number in the range of [
0
,
1
] and m
1  (m
5) is mass of carbon (hydrogen) atom. The x and y-components of velocities in (14) must satisfy
(15)
x
1
x
˙
1
+
y
1
y
˙
1
=
0
,
x
4
x
˙
4
+
y
4
y
˙
4
=
0
,
x
5
x
˙
5
+
y
5
y
˙
5
=
0
,
x
6
x
˙
6
+
y
6
y
˙
6
=
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:
(16)
(
x
1
y
˙
1
-
y
1
x
˙
1
)
+
(
x
4
y
˙
4
-
y
4
x
˙
4
)
=
0
,
(
x
5
y
˙
5
-
y
5
x
˙
5
)
+
(
x
6
y
˙
6
-
y
6
x
˙
6
)
=
0
.
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 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 T
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-to-cis isomerization, and the total energy is defined by the excitation energy 4.11 eV fromS
0 to S
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 S
0 to S
1 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 λ
k
(
t
)
=
0. Then, outcome of coordinates is inserted into the Jacobian determinant J
k
k
′ in (7) to have the new λ
k
(
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 couplinga
2 and effective collision energy b
2 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 a
2 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 Table2 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 × 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 S
0 to S
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.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].Potential energy surface ofS
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 b
2 is well above the seam and this makes trajectory have a big chance to hop down S
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 b
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-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 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 Table2 still stand. This is because effective coupling and collision energy parameters (a
2 and b
2) 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. Figure4 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 S
1 state and increases on S
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 S
1 to S
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.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.
(a)
(b)
(c)Figure5 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 S
0 and S
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 S
1 to S
0 state happened at 50 fs, the second hopping from S
0 to S
1 state at 71 fs, and the last hopping from S
1 to S
0 state at 88 fs. After that, the trajectory is propagated on S
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 cis-to-trans trajectory in Figure 4.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.
(a)
(b)
(c)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 Figure6; 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 S
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.The dihedral angels DD1 and DD2 vary with time for 100 trajectories; (a) and (b) correspond to initially start from cis-area and trans-area, respectively.
(a)
(b)


## 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 Table2 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 × 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 S
0 to S
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.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].Potential energy surface ofS
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 b
2 is well above the seam and this makes trajectory have a big chance to hop down S
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 b
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-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 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 Table2 still stand. This is because effective coupling and collision energy parameters (a
2 and b
2) 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. Figure4 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 S
1 state and increases on S
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 S
1 to S
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.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.
(a)
(b)
(c)Figure5 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 S
0 and S
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 S
1 to S
0 state happened at 50 fs, the second hopping from S
0 to S
1 state at 71 fs, and the last hopping from S
1 to S
0 state at 88 fs. After that, the trajectory is propagated on S
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 cis-to-trans trajectory in Figure 4.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.
(a)
(b)
(c)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 Figure6; 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 S
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.The dihedral angels DD1 and DD2 vary with time for 100 trajectories; (a) and (b) correspond to initially start from cis-area and trans-area, respectively.
(a)
(b)


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



---

*Source: 132149-2014-07-07.xml*