Research Article
Constraint Trajectory Surface-Hopping Molecular Dynamics
Simulation of the Photoisomerization of Stilbene
Yibo Lei,
1,2Shaomei Wu,
1Chaoyuan Zhu,
2Zhenyi Wen,
1and Sheng-Hsien Lin
2 1Key Laboratory of Synthetic and Natural Functional Molecule Chemistry of Ministry of Education,College of Chemistry & Materials Science, Shaanxi Key Laboratory of Physico-Inorganic Chemistry, Northwest University, Xi’an 710069, China
2Department of Applied Chemistry, Institute of Molecular Science and Center for Interdisciplinary Molecular Science, National Chiao-Tung University, Hsinchu 300, Taiwan
Correspondence should be addressed to Yibo Lei; [email protected] and Chaoyuan Zhu; [email protected] 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𝑆0and excited𝑆1states 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. Trajec-tories 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 equa-tions 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 cou-pling vectors while the transition probabilities calculated from LZ and ZN formulas require information only from two adiabatic potential energy surfaces around the crossing
Volume 2014, Article ID 132149, 18 pages http://dx.doi.org/10.1155/2014/132149
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 ener-gies, energy gradients, and nonadiabatic coupling vectors are calculated on-the-fly for trajectories with full degree of free-doms 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 tra-jectory 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 coor-dinates 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 quan-tum yields of the photoisomerization were measured exper-imentally [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 accompa-nying electronically nonadiabatic transition among various electronic states. Traditionally, the isomerization of cis- and trans-conformations has been investigated by computing the ground state 𝑆0 and the lowest excited state 𝑆1. The dominant excitation from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) corresponds to the 𝑆1 𝜋𝜋∗ transition which was confirmed to be a bright state [45]. In the recent years, the photoisomerization involving the conical intersections plus fluorescence spectra has been extensively investigated by both experimentalists and theoreticians [34,35,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 mecha-nism [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 investi-gated 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
1 11 10 8 9 2 12 6 5 3 17 7 21 22 23 24 26 25 20 4 18 16 14 13 15 19 D1: C11-C8-C1-C2 D2: H10-C8-C1-H9
Figure 1: Atom numbering of Stilbene at the𝑆1/𝑆0conical 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 dynam-ical 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 inFigure 1). The potential energy surfaces of ground-state 𝑆0 and the first excited state𝑆1 are constructed in terms of the D1 and D2 angles, and the other internal coordinates are fixed at the OBF-CI configuration.
The rest of this paper is organized as follows.Section 2
first summarizes constrained molecular dynamics method and constructs two-dimensional potential energy surfaces analytically for Stilbene. Trajectory surface hopping with one-passage nonadiabatic transition probability is carried out by LZ and ZN formula.Section 3presents the implementations of the methods mentioned inSection 2, and results and dis-cussions about branching ratio of isomerization and picture of evolution trajectory are also given therein. Concluding remarks are mentioned inSection 4followed by three Appen-dicesA,B, andCin which the detailed constraint relations, nonadiabatic transition probability, and hopping scheme are presented.
2. Constraint TSH Molecular
Dynamics with LZ and ZN Formulas for
Nonadiabatic Transitions
We first introduce constraint molecular dynamics for system with𝑁𝑐constraints, and then we construct 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𝑁𝑐constraints in which the canonical Hamiltonian equations can be written as [24–
30] ̇𝑞 𝑖= 𝜕𝐻 (𝑞, 𝑝, 𝜆)𝜕𝑝 𝑖 = 𝑝𝑖 𝑚𝑖, 𝑖 = 1, 2, . . . , 3𝑁, (1) ̇𝑝 𝑖= −𝜕𝐻 (𝑞, 𝑝, 𝜆)𝜕𝑞 𝑖 = 𝑓 uc 𝑖 + 𝑓𝑖𝑐= − 𝜕𝑉 (𝑞) 𝜕𝑞𝑖 − 𝑁𝑐 ∑ 𝑘 𝜆𝑘𝜕𝑔𝑘(𝑞) 𝜕𝑞𝑖 , (2) where𝑚𝑖,𝑞𝑖, and𝑝𝑖are the mass, the Cartesian coordinate, and conjugate momentum for each atom in molecule, respec-tively. The𝑔𝑘(𝑞) = (𝑘 = 1, 2, . . . , 𝑁𝑐) are constraint equations determined by𝑁𝑐internal coordinates that are fixed at the certain structures of nuclear configuration all the time. The 𝑓uc
𝑖 and 𝑓𝑖𝑐 in (2) stand for unconstraint and constraint
force acting on the corresponding atom, respectively. The Hamiltonian with𝑁𝑐Lagrange multipliers𝜆𝑘(𝑡) is given as
𝐻 (𝑞, 𝑝, 𝜆) =3𝑁∑ 𝑖 𝑝2𝑖 2𝑚𝑖 + 𝑉 (𝑞) + 𝑁𝑐 ∑ 𝑘 𝜆𝑘𝑔𝑘(𝑞) , (3) where𝜆𝑘(𝑡) are to be determined. As 𝑔𝑘(𝑞) in (3) are zero all the time, the energy conversation is satisfied.
The Lagrange multipliers are solved iteratively. With the given initial Lagrange multipliers, we solve 𝑞𝑖 and 𝑝𝑖 by the numerical integration of (1) and (2), and then the 𝑞𝑖 is inserted into the constraint equations (including bond lengths, bond angles, and dihedral angles) to get the Lagrange multipliers. This iterative process is eventually convergent. The technique to implement this iterative procedure is start-ing from the classical (Newtonian) equations [27]:
𝑚𝑖 ̈𝑞𝑖(𝑡) = −𝜕𝑞𝜕 𝑖[𝑉 (𝑞) + 𝑁𝑐 ∑ 𝑘=1 𝜆𝑘(𝑡) 𝑔𝑘(𝑞)] . (4)
Integrating both sides of (4) twice in time yields the con-straint Cartesian coordinates at the time𝑡 + Δ𝑡:
𝑞𝑖(𝑡 + Δ𝑡) = 𝑞𝑖uc(𝑡 + Δ𝑡) + 𝑚−1𝑖 (Δ𝑡)2𝑓𝑖𝑐(𝑡) = 𝑞uc 𝑖 (𝑡 + Δ𝑡) − 𝑚−1𝑖 (Δ𝑡)2 𝑁𝑐 ∑ 𝑘 𝜆𝑘𝜕𝑔𝜕𝑞𝑘(𝑞) 𝑖 , (5)
where𝑞𝑖uc(𝑡 + Δ𝑡) is solved by unconstrained force. The above 𝑞𝑖(𝑡 + Δ𝑡) can be inserted to every type of the constraint equations to derive𝑁𝑐nonlinear equations with respect to the corresponding𝑁𝑐Lagrange multipliers. Although there are a number of algorithms to solve𝜆𝑘(𝑡), they only differ on how to solve these nonlinear equations. The SETTLE algorithm [28,29] obtained the solutions of these nonlinear equations analytically for𝑁𝑐= 3 constraints, but it cannot be extended to the large system. The originally simple SHAKE algorithm [24,25] was developed to solve bond length constraints which are converged linearly. Later advanced M-SHAKE Newton method [24, 25, 30] and quasi-Newton method [30] were straight-forward and fast to solve the nonlinear equations. All of these methods involved the calculations of the Jacobian determinant about the constraint equations𝑔𝑘= 0:
𝐽𝑔= [ [ [ [ [ [ [ [ [ [ [ [ 𝜕𝑔1 𝜕𝜆1 𝜕𝑔1 𝜕𝜆2 ⋅ ⋅ ⋅ 𝜕𝑔1 𝜕𝜆𝑁𝑐 𝜕𝑔2 𝜕𝜆1 𝜕𝑔2 𝜕𝜆2 ⋅ ⋅ ⋅ 𝜕𝑔2 𝜕𝜆𝑁𝑐 .. . ... ... ... 𝜕𝑔𝑁𝑐 𝜕𝜆1 𝜕𝑔𝑁𝑐 𝜕𝜆2 ⋅ ⋅ ⋅ 𝜕𝑔𝑁𝑐 𝜕𝜆𝑁𝑐 ] ] ] ] ] ] ] ] ] ] ] ] , (6)
where𝜆𝑘are then updated repeatedly by solving the system of linear equations iteratively:
(𝐽𝑔𝜆
𝑘=0) 𝜆 = −𝑔𝜆𝑘=0 (7)
until the max|𝑔𝑘| < 𝜏 (𝑘 = 1, . . . , 𝑁𝑐) in which 𝜏 is small number with the prescribed tolerance of the constraints. The Jacobian determinant for the quasi-Newton method is only computed once in the first iteration with the linear con-vergence, 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 inAppendix A.
2.2. Two-Dimensional Potential Energy Surfaces for Stilbene.
Quenneville and Mart´ınez [45] have done extensive calcula-tion for the critical𝑆1/𝑆0 conical intersection OBF-CI and potential energy surfaces in terms of two internal angles: phenyl rotation and the pyramidalization coordinates. They have done detailed comparison for the results calculated by complete active space self-consistent field (CASSCF) and multireference perturbation theory (CASPT2) methods with
basis sets of 6-31G., 6-31G∗, and 6-31G∗∗. Following their
conclusion, we can safely utilize state-averaged CASSCF method abbreviated as SA-N-CAS(n/m), where N refers to the number of states included in the average, while
N(=2), n(=2), and m(=2) are the number of active electrons
and orbitals, respectively. This means that just 𝜋 and 𝜋∗ orbitals corresponding to HOMO and LUMO are involved in CASSCF calculation. Molcas 7.5 [66] program packages were employed for all calculations.
We have first used SA-2-CAS(2/2)/6-31G method to reproduce geometries of conical intersection OBF-CI, the local minima on𝑆0 and𝑆1 potential energy surfaces given by [45]. Then, we concentrate on configuration of conical intersection at which the torsion angle D1 (C11-C8-C1-C2) and dihedral angle D2 (H10-C8-C1-H9) are −146.6∘ and −60.0∘, respectively, as shown inFigure 1. We have carried out
a preliminary test calculation for two-dimensional potential energy surfaces with variables D1 and D2 (all bond lengths, all bond angles, and the other dihedral angles are fixed at OBF-CI configuration). We found that DD1 = (D1 + D2)/2 and DD2 = (D2 − D1)/2 are better variables to describe the isomerization. If we set up interesting region for potential energy surface𝑆1 with energy not higher than 5 eV above OBF-CI (where it is considered as zero point of potential energy), DD2 can be chosen as in [−20∘, 90∘] and DD1 is
chosen as in [0∘, −180∘] that is sufficient to describe the present isomerization process. We use equal spacing 5∘ for grid points: DD1 is0, −5, . . . , −180 (37 grid points) and DD2 is −20, −15, . . . , 90 (23 grid points). Total configurations are 851 plus OBF-CI which all are computed by SA-2-CAS(2/2)/6-31G method, and then all ab. initio data are fitted into analytical function of potential energy surfaces by the least square method from Matlab 7.5 package [67]. Finally, we have adiabatic potential energy surfaces in the analytical form:
𝑊 (𝑥, 𝑦) = 𝑐0exp[−𝑎1(𝑥 − 𝑥0)2− 𝑏1(𝑦 − 𝑦0)2] + 𝑐1cos(𝑥
2) + 𝑐2cos( 𝑦
2) + 𝑐3cos(𝑥) + 𝑐4cos(𝑦) + 𝑐5cos(2𝑥) + 𝑐6cos(2𝑦) + 𝑐7cos(3𝑥) + 𝑐8cos(3𝑦) + 𝑐9sin[𝑥 + 𝑦2 ]
+ 𝑐10sin[2𝑥 + 𝑦 2 ] + 𝑐11sin[𝑥 + 2𝑦2 ] + 𝑐12sin(𝑥 + 𝑦) + 𝑐13sin(2𝑥 + 𝑦) + 𝑐14sin(𝑥 + 2𝑦) + 𝑐15sin(2𝑥 + 2𝑦) + 𝑐16sin(3𝑥 + 𝑦) + 𝑐17sin(𝑥 + 3𝑦) + 𝑐18sin(3𝑥 + 2𝑦) + 𝑐19sin(2𝑥 + 3𝑦) + 𝑐20sin(3𝑥 + 3𝑦) + 𝑐21cos[𝑥 + 𝑦 2 ] + 𝑐22cos[2𝑥 + 𝑦2 ] + 𝑐23cos[𝑥 + 2𝑦2 ]
0 0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 P o te n tial ener g y (eV) 10 20 −180 −150 −120 −90 −60 −30 −2 −1 DD2 (∘) DD 1 (∘)
Figure 2: Analytical two-dimensional PESs for the ground state𝑆0 and the first excited state𝑆1with respect to the combined internal angles DD1 and DD2 [33]. + 𝑐24cos(𝑥 + 𝑦) + 𝑐25cos(2𝑥 + 𝑦) + 𝑐26cos(𝑥 + 2𝑦) + 𝑐27cos(2𝑥 + 2𝑦) + 𝑐28cos(3𝑥 + 𝑦) + 𝑐29cos(𝑥 + 3𝑦) + 𝑐30cos(3𝑥 + 2𝑦) + 𝑐31cos(2𝑥 + 3𝑦) + 𝑐32cos(3𝑥 + 3𝑦) , (8) where𝑥 = DD1 and 𝑦 = DD2, 𝑥0 = −103.3∘and𝑦0 = 43.3∘ are angles at OBF-CI. The fitting parameters 𝑎1, 𝑏1, and𝑐0 to 𝑐32 are all given in Table 1. Potential energy surfaces of 𝑆0 and 𝑆1 calculated from (8) are plotted in Figure 2with
zero point of potential energy at𝑊(𝑥0, 𝑦0) = 0 in a unit of electron volt. We computed the mean absolute error between calculated results from (8) and numerical data computed from method SA-2-CAS(2/2)/6-31G for about 100 randomly picked up configuration points and we found that this error is about 2.4 kcal/mol for the both potential energy surfaces, and less than 1.0 kcal/mol around OBF-CI region. From analytical PES in (8), we have estimated local trans-isomer at DD1 =
−164.3∘ and DD2 = 36.7∘ with energy 1.72 eV below
OBF-CI, and local cis-isomer at DD1 = −39.2∘ and DD2 =
50.9∘ with energy 1.28 eV below OBF-CI. Furthermore, the
vertical excitation energies have been calculated to be 4.11 eV and 4.45 eV at local trans- and cis-minima, respectively, which are accidently in consistence with the experimental measurement 4.13 eV and 4.59 eV [68]. We confirmed that a global minimum of𝑆1state 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
imple-ments classical trajectory simulation on a single potential energy within the Cartesian coordinate system plus con-straints in internal coordinate systems. This not only sim-plified numerical calculation for integration of classical tra-jectory in molecular Hamiltonian but also avoided searching complicated form of kinetic operators in internal coordinates. However, nonadiabatic transition between two PESs is com-pletely 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𝑆0and𝑆1PESs around one-bond flip conical intersection (OBF-CI) (𝑐0 ∼ 𝑐32in units of eV) [33]. Index 𝑎1 𝑏1 𝑐0 𝑐1 𝑐2 𝑆0 35 3.5 0.495231 364.2436 166.2537 𝑆1 39 3.9 −0.52219 256.2847 −18.8715 Index 𝑐3 𝑐4 𝑐5 𝑐6 𝑐7 𝑆0 −63.2233 −228.305 −12.7773 25.79321 −3.78713 𝑆1 −74.255 −85.6398 −4.51131 16.61385 −1.65933 Index 𝑐8 𝑐9 𝑐10 𝑐11 𝑐12 𝑆0 −3.97288 9.467412 106.3829 −36.167 −72.9472 𝑆1 −2.44485 −219.158 171.3529 98.32809 −110.789 Index 𝑐13 𝑐14 𝑐15 𝑐16 𝑐17 𝑆0 5.205437 −0.08421 −4.58729 5.134825 3.857498 𝑆1 9.156298 −4.96308 −1.1299 0.059523 4.654567 Index 𝑐18 𝑐19 𝑐20 𝑐21 𝑐22 𝑆0 −6.20916 0.171133 2.090849 −538.156 −56.4412 𝑆1 −1.20414 −1.73863 0.658395 −277.529 87.56788 Index 𝑐23 𝑐24 𝑐25 𝑐26 𝑐27 𝑆0 254.7531 157.3086 17.19558 −58.19 −8.161 𝑆1 112.7934 47.53832 −7.42468 −41.5302 7.965267 Index 𝑐28 𝑐29 𝑐30 𝑐31 𝑐32 𝑆0 7.267131 4.132838 −3.08607 1.734321 0.088946 𝑆1 1.25399 3.375719 −0.42815 −0.71119 −0.09718 D1 D2 H H 1 1 5 6 4 z-axis 2 3 → a =→r2−→r1 → b =→r3− → r2 → c =→r4− → r3 → d=→r2−→r5 → e =→r6− → r3 Φ1 Φ2 Scheme 1
dihedral angles D2 (transferred into DD1 and DD2 for PESs depicted inFigure 2) and also present the following regular form of kinetic operators explained inScheme 1, with
𝑇 = −2𝐼ℎ2 D1 𝜕2 𝜕D12 − ℎ2 2𝐼D2 𝜕2 𝜕D22, (9)
where the dihedral angle D1 represents the relative torsion motion between two phenyl rings and the D2 represents the relative torsion motion between two hydrogen atoms. We assume that atoms 2 and 3 are fixed and two phenyl rings and two hydrogen atoms rotate on the circles around𝑧-axis as shown inScheme 1. For instance, we can think thatΦ1 and Φ2 are rotation angles around 𝑧-axis for two phenyl rings, respectively. Thus, the moment of inertia𝐼D1is the reduced inertia from two phenyl rings
𝐼D1= 𝐼𝐼Φ1𝐼Φ2
Φ1+ 𝐼Φ2, (10)
where𝐼Φ1= ∑𝑖𝑚𝑖𝜌𝑖2(𝑚𝑖is atom mass in the one phenyl ring and𝜌𝑖is distance from atom𝑖 to 𝑧-axis) and 𝐼Φ2is estimated
from the other phenyl ring. The same analysis can be applied for the moment of inertia𝐼D2with two hydrogen atoms:
𝐼D2= 𝐼𝐼𝐻1𝐼𝐻2
𝐻1+ 𝐼𝐻2, (11)
where 𝐼𝐻1 = 𝑚𝐻𝜌𝐻2 (𝑚𝐻 is hydrogen atom mass and 𝜌𝐻 is distance from hydrogen atom 𝐻 to 𝑧-axis) and 𝐼𝐻2 is estimated from the other hydrogen atom.
One-passage semiclassical nonadiabatic transition prob-ability (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):
𝑝LZ= exp (− 𝜋 4√𝑎2𝑏2) , (12) 𝑝ZN= exp (− 𝜋 4√𝑎2𝑏2√ 2 1 + √1 + (0.4𝑎2+ 0.7) 𝑏−4) , (13) where calculation of two parameters𝑎2 and 𝑏2 is given in
Appendix B.
We can simply look at the present two-dimensional PESs
in Figure 2 to find that the seam line is appearing at the
DD1 = −103.3∘. Seam line is defined as the trace of local
minimum gap between two adiabatic PESs [69]. However, in the present work we propose to calculate the local maxima of the effective coupling parameter 𝑎2 as introduced in
Appendix C. Once the seam line is found, the nonadiabatic
coupling vector is assumed to be perpendicular to the seam line and it appears in DD1 direction in the present case. Therefore, the effective collision energy parameter𝑏2is calcu-lated 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 inAppendix C.
The calculation shows that the seam line appears at DD1 = −103.3∘ with negligible dependence of DD2 angle. The effective coupling parameter𝑎2 along the seam line is calculated and depicted inFigure 3.Figure 3shows 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 with0.1 ≤ 𝑎2≤ 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 sur-faces. We construct the reduced two-dimensional (2D)
0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 7 8 −20 −10 −2 −1 log a 2 DD2 (∘)
Figure 3: The effective coupling parameter𝑎2along the seam line with DD1 =−103.3∘.
potential energy surfaces at OBF-CI, and𝑆0potential 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 𝑆0to𝑆1PES and then reaches the certain area of the present local cis (trans)-minimum on the 2D 𝑆1 PES, from where trajectories start isomerization. Then, question is how to determine such area for both local cis- and trans-isomers. Let us start with the test local cis-area and local trans-area in rectangular region of [−19.2∘ < DD1 < −59.2∘, 30.9∘ < DD2< 70.9∘] and [−180.0∘ < DD1 < −144.3∘, 16.7∘ < DD2 < 56.7∘], respectively, for both initial and final conditions
of trajectories. Total energy for cis-to-trans isomerization is defined by the excitation energy 4.45 eV from𝑆0to𝑆1state estimated at local cis-minimum (DD1 =−39.2∘, DD2 = 50.9∘). The initial conditions of DD1 and DD2 for trajectories are randomly picked up in the local cis-area vertically excited from𝑆0 to𝑆1 state; if this trajectory has vertical excitation energy which is larger than 4.45 eV, it is abandoned, and if this trajectory has vertical excitation energy which is smaller than 4.45 eV, it can be proceeded. Then, we specify initial kinetic energy for proceeding trajectory withScheme 1by choosing atoms 1, 4, 5, and 6 with nonzero velocities and all atoms in molecule with zero velocities. The proceeding trajectory which starts from local cis-area on𝑆1 state can have three ending ways; one is that it enters local cis-area on𝑆0(called cis-to-cis nonreactive trajectory), another is that it enters local-trans area on𝑆0(called cis-to-trans reactive trajectory), and the third is that it enters outside of 2D PES boundary region depicted in Figure 2; for instance, DD1 is smaller than −180∘ or larger than 0∘ (called unreactive trajectory). An unreactive trajectory is due to that 2D PES has periodic properties in terms of DD1 and DD2 angles and we do count it as nonreactive trajectory.
Now we present how to distribute a given kinetic energy 𝑇initto atoms 1, 4, 5, and 6 with their initial velocities analyzed
in the following. As shown in Scheme 1, the z-component of velocity is assumed to be zero so that for a given kinetic energy𝑇init, we have the following relations:
1
2𝑚1( ̇𝑥21+ ̇𝑦21+ ̇𝑥24+ ̇𝑦42) = 𝑇init𝛿,
1
2𝑚5( ̇𝑥25+ ̇𝑦25+ ̇𝑥26+ ̇𝑦62) = 𝑇init(1 − 𝛿) ,
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𝛿 is random number in the range of [0, 1] and 𝑚1 (𝑚5) is mass of carbon (hydrogen) atom. The𝑥 and 𝑦-components of velocities in (14) must satisfy
𝑥1 ̇𝑥1+ 𝑦1 1̇𝑦 = 0, 𝑥4 ̇𝑥4+ 𝑦4 4̇𝑦 = 0,
𝑥5 ̇𝑥5+ 𝑦5 5̇𝑦 = 0, 𝑥6 ̇𝑥6+ 𝑦6 6̇𝑦 = 0,
(15) because of constraints of bond lengths and bond angles. Moreover, we choose the angular moment along𝑧-axis for atoms 1 and 4 and for atoms 5 and 6 as follows:
(𝑥1 ̇𝑦1− 𝑦1 ̇𝑥1) + (𝑥4 4̇𝑦 − 𝑦4 ̇𝑥4) = 0,
(𝑥5 ̇𝑦5− 𝑦5 ̇𝑥5) + (𝑥6 6̇𝑦 − 𝑦6 ̇𝑥6) = 0.
(16) This means that angular momentum along𝑧-axis for the twist of atoms 1 and 4 is cancelled out from each other, so does for atoms 5 and 6. Equations (14), (15), and (16) all together have the eight related equations for the eight velocity components to be determined. However, the above procedure slightly violates conservation of total energy because two phenyl rings associated with atoms 1 and 4 inScheme 1can have induced velocities through the constraints for entire molecule. We can resolve this problem by choosing smaller kinetic energy 𝑇init in (14) than the previously specified one, and this can
be obtained by subtracting induced kinetic energy from two phenyl rings. Of course, this has to be done iteratively with the constraint Hamiltonian equations (1) and (2) in the first step of classical trajectory integration.
We can do the same as mentioned above for trans-to-cis isomerization, and the total energy is defined by the excitation energy 4.11 eV from𝑆0to𝑆1state 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𝑆0to 𝑆1state. Three types of ending trajectories are also collected as mentioned above (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 con-jugate moments, we integrate (1) and (2) with initial
the Lagrange multipliers𝜆𝑘(𝑡) = 0. Then, outcome of coor-dinates is inserted into the Jacobian determinant𝐽𝑘𝑘 in (7)
to have the new𝜆𝑘(𝑡) that are inserted into (1) and (2) again. We do this iteratively until all Lagrange multipliers converge together with the update coordinates and momentums. For the second step of integration, the Lagrange multipliers are set up as the previous values, and in this way, we propagate trajectory until it ends with use of the fourth-order Runge-Kutta method (RK4) [72] by equal time step size 2.07 a.u. (about 0.05 fs).
Classical trajectory is propagated on a single adiabatic potential energy surface all the time. However, once tra-jectory is across the seam line, we compute effective cou-pling𝑎2 and effective collision energy 𝑏2 according to the method given inAppendix 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 𝑎2 and found that they all are located at avoided crossing points. Comparing𝑝 with random number in between 0 and 1, we decide to hop trajectory from one adiabatic potential energy surface to another if𝑝 is larger than the random number; otherwise, it stays in the original potential energy surface. As trajectory hops to the other PES vertically, and it means that all coordinates are unchanged, but the moments are changed with method introduced in
Appendix C.
3.1. Branching Ratios of Cis-to-Trans and Trans-to-Cis Iso-merizations. 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× 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 simula-tion. There are two pathways when Stilbene is photo-excited from the global cis-isomer 𝑆0 to 𝑆1 state (see [59, 61] and references therein); one pathway leads to the OBF-CI with branching ratio 70% and another goes to the other conical intersection (branching ratio 30%) which is for side reaction of a ring closing to form DHP. If we take this branching ratio into consideration, we can estimate the quantum yield for cis-to-trans isomerization as 0.7 × 67% = 47% which is larger than experimental measurement with 35%. The present simulation shows isomerization in favor of reactive trajectories regardless of initial starting from local trans-area or local cis-area and this means that branching ratio is not 50 : 50 at OBF-CI.
Potential energy surface of𝑆1state as shown inFigure 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𝑏2is well above the seam and this makes trajectory have a big chance to hop down𝑆0state 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 inTable 2. Meanwhile, initial vertical excitation energy is 3.17 eV above OBF-CI at cis-area and is 2.39 eV above OBF-CI at trans-area, so that effective collision energy𝑏2 is higher for cis-to-trans than for trans-to-cis. This explains why branching ratio is more in favor of reactive cis-to-trans than trans-to-cis. We know that the conical intersection Hula-Twist (HT-CI) is higher in energy than in OBF-CI, and these two CIs might be in close configuration. Therefore, we suspect that there is a big chance for cis-to-trans trajectory to access HT-CI with high vertical excitation energy, but this is not the case for trans-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
inTable 2still stand. This is because effective coupling and
collision energy parameters (𝑎2 and𝑏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.Figure 4shows 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 inFigure 4(a); it can be seen that potential energy steeply decreases on𝑆1state and increases on𝑆0 state with crossing seam line around 52 fs first and then crossing again shortly, but no hopping occurs due to the fact that nonadiabatic transition probability is smaller than
random number there. Trajectory hopping from𝑆1to𝑆0state 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 inFigure 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
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 𝑆0 and 𝑆1 states vary with time with four selected snapshots of structure as shown
in Figure 5(a) for a reactive trans-to-cis trajectory. There
are three-time surface hopping from one potential energy surface to another during evolution of trajectory; the first hopping from𝑆1 to𝑆0 state happened at 50 fs, the second hopping from 𝑆0 to 𝑆1 state at 71 fs, and the last hopping from 𝑆1 to 𝑆0 state at 88 fs. After that, the trajectory is propagated on 𝑆0 state until its end at 133 fs. Figure 5(b)
shows that both DD1 and DD2 angles change with oscillation against time; DD1 rises considerably while DD2 keeps almost constant oscillation.Figure 5(c)shows that the torsion angle D1 increases monotonically from−188∘(at local trans-area) to −97∘(at local cis-area), and this shows one-bond flip picture
around the central C=C bond again. Dihedral angle D2 increases from the beginning at−128∘to the end at−5∘with four-circle oscillations. The reactive trans-to-cis trajectory in
Figure 5evolutes longer time period than that of reactive
cis-to-trans trajectory inFigure 4.
With randomly choosing 100 initial conditions as men-tioned before in both trans-area and cis-area, we have sam-pled 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 inFigure 6;
Figure 6(a)shows the most of trajectories started from
cis-area end around 120 fs in average, whileFigure 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𝑆1state 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
0 20 40 60 80 100 120 0 1 2 3 4 5 6 Time (fs) P o te n tial ener g y (eV) −2 −1 S0 S1 (a) 0 20 40 60 80 100 120 0 50 100 Time (fs) DD1 DD2 −200 −150 −100 −50 Dihedral a n gle ( ∘ ) (b) 0 20 40 60 80 100 120 −100 −50 0 50 Time (fs) D1 D2 −200 −150 Dihedral a n gle ( ∘) (c)
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 trajec-tory, in which D1 changes monotonically from reactant to product for both cis-to-trans and trans-to-cis isomerizations. Trajectory surface hopping simulation based on the present two-dimensional potential energy surfaces can present quite quantitative information for isomerization of Stilbene. We have proposed to calculate two effective parameters𝑎2 and 𝑏2 along trajectory moving on adiabatic potential energy surfaces. Effective coupling parameter𝑎2is not only used for determining transition probability but also for determining the seam line; in the present two-dimensional case this is
0 20 40 60 80 100 120 140 0 1 2 3 4 5 6 Time (fs) P o te n tial ener g y (eV) S0 S1 −2 −1 (a) 0 20 40 60 80 100 120 140 0 50 100 Time (fs) −200 −150 −100 −50 Dihedral a n gle ( ∘) DD1 DD2 (b) 0 20 40 60 80 100 120 140 0 50 Time (fs) −200 −150 −100 −50 Dihedral a n gle ( ∘) D1 D2 (c)
Figure 5: A typical reactive trajectory for trans-to-cis isomerization starts from local trans-area with photo-excitation energy 3.98 eV. (a) Snapshots of potential energy surfaces and structures varying with time, (b) DD1 and DD2, and (c) D1 and D2 varying with time.
0 50 100 150 200 250 300 350 400 Time (fs) −100 −50 0 100 50 −200 −150 Dihedral a n gle ( ∘ ) DD1 DD2 (a) 0 100 200 300 400 500 600 Time (fs) −100 −50 0 100 50 −200 −150 Dihedral a n gle ( ∘) DD1 DD2 (b)
Figure 6: 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.
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𝑎2along 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 exci-tation 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𝑑𝑘1𝑘2, bond angles 𝜃𝑘1𝑘2𝑘3, and dihedral angles𝜙𝑘1𝑘2𝑘3𝑘4 in the following equa-tions, respectively: 𝑔𝑏𝑘(𝑞) = 𝑞2𝑘1𝑘2− 𝑑2𝑘1𝑘2 = 0, 𝑘 = 1, . . . , 𝑁𝑏𝑐, (A.1) 𝑔𝑘𝑎= cos 𝜃𝑘1𝑘2𝑘3 ⃗𝜇𝑎 ⃗]𝑎 − ⃗𝜇𝑎⋅ ⃗]𝑎= 0, 𝑘 = 1, . . . , 𝑁𝑐𝑎, (A.2) 𝑔𝑘𝑑= cos 𝜙𝑘1𝑘2𝑘3𝑘4 ⃗𝜇 𝑑 ⃗]𝑑 − ⃗𝜇𝑑⋅ ⃗]𝑑= 0, 𝑘 = 1, . . . , 𝑁𝑐𝑑, (A.3) where Cartesian coordinates are𝑞2𝑘
1𝑘2 = ( ⃗𝑟𝑘2− ⃗𝑟𝑘1) 2, ⃗𝜇𝑎= ⃗𝑟 𝑘1− ⃗𝑟 𝑘2, ⃗] 𝑎= ⃗𝑟 𝑘2− ⃗𝑟𝑘3, ⃗𝜇 𝑑= ⃗𝑎 × ⃗𝑏, ⃗V𝑑= ⃗𝑏× ⃗𝑐, ⃗𝜇𝑑⋅ ⃗]𝑑= ( ⃗𝑎× ⃗𝑏)⋅( ⃗𝑏× ⃗𝑐),
⃗𝑎 = ⃗𝑟2− ⃗𝑟1, ⃗𝑏 = ⃗𝑟3− ⃗𝑟2, and ⃗𝑐 = ⃗𝑟4− ⃗𝑟3and the subscript of each
variable is atom number. In terms of the Cartesian coordinate
1
2 3
4 a
b c
Figure 7: Schematic diagram for relations between Cartesian and internal coordinates for four atoms.
𝑞 with the derivations of the constraint equations, equation (5) can be rewritten as 𝑞𝑘𝛽(𝑡 + Δ𝑡) = 𝑞uc 𝑘𝛽(𝑡 + Δ𝑡) − 𝑚 −1 𝑘𝛽(Δ𝑡) 2[ [ 𝑁𝑏 𝑐 ∑ 𝑘=1 𝜆𝑏𝑘(𝑡) 𝜕𝑔𝑏 𝑘(𝑡) 𝜕𝑞𝑘𝛽(𝑡) + 𝑁𝑎 𝑐 ∑ 𝑘=1 𝜆𝑎𝑘(𝑡)𝜕𝑔 𝑎 𝑘(𝑡) 𝜕𝑞𝑘𝛽(𝑡) + 𝑁 𝑑 𝑐 ∑ 𝑘=1 𝜆𝑑𝑘(𝑡) 𝜕𝑔 𝑑 𝑘(𝑡) 𝜕𝑞𝑘𝛽(𝑡)] ] (A.4) in which 𝜕𝑔𝑏 𝑘 𝜕𝑞𝑘𝛽 = 𝛿𝑘 1𝑘𝛽 𝜕𝑔𝑏 𝑘(𝑡) 𝜕𝑞𝑘 1(𝑡) + 𝛿𝑘 2𝑘𝛽 𝜕𝑔𝑏 𝑘(𝑡) 𝜕𝑞𝑘 2(𝑡) = 2 (𝛿𝑘𝛽,𝑘1− 𝛿𝑘𝛽,𝑘2) (𝑞𝑘1− 𝑞𝑘2) , 𝜕𝑔𝑎𝑘 𝜕𝑞𝑘𝛽 = 𝛿𝑘1𝑘𝛽 𝜕𝑔𝑘𝑎(𝑡) 𝜕𝑞𝑘 1(𝑡) + 𝛿𝑘 2𝑘𝛽 𝜕𝑔𝑘𝑎(𝑡) 𝜕𝑞𝑘 2(𝑡) + 𝛿𝑘 3𝑘𝛽 𝜕𝑔𝑎 𝑘(𝑡) 𝜕𝑞𝑘 3(𝑡) = (𝛿𝑘 1𝑘𝛽− 𝛿𝑘2𝑘𝛽) 𝜕𝑔𝑎 𝑘 𝜕𝜇𝑎 + (𝛿𝑘3𝑘𝛽− 𝛿𝑘2𝑘𝛽) 𝜕𝑔𝑎 𝑘 𝜕]𝑎, (A.5) with 𝜕𝑔𝑎 𝑘 𝜕𝜇𝑎 𝑙 = cos 𝜃𝑘1𝑘2𝑘3 ⃗𝜇 ⃗]𝑎𝑎𝜇𝑎𝑙 − ]𝑎𝑙, 𝜕𝑔𝑎𝑘 𝜕]𝑎 𝑙 = cos 𝜃𝑘1𝑘2𝑘3 ⃗𝜇 ⃗]𝑎𝑎]𝑎𝑙 − 𝜇𝑙𝑎, 𝜕𝑔𝑑 𝑘 𝜕𝑞𝑘𝛽 = 𝛿𝑘 1𝑘𝛽 𝜕𝑔𝑑 𝑘(𝑡) 𝜕𝑞𝑘 1(𝑡) + 𝛿𝑘 2𝑘𝛽 𝜕𝑔𝑑 𝑘(𝑡) 𝜕𝑞𝑘 2(𝑡) + 𝛿𝑘 3𝑘𝛽 𝜕𝑔𝑑𝑘(𝑡) 𝜕𝑞𝑘 3(𝑡) + 𝛿𝑘 4𝑘𝛽 𝜕𝑔𝑘𝑑(𝑡) 𝜕𝑞𝑘 4(𝑡)
= [𝜕𝑔𝑑𝑘 𝜕𝜇𝑑 𝑚 𝜕𝑔𝑑 𝑘 𝜕𝜇𝑑 𝑛 ] ⋅ [ 𝑞𝑛,𝑘 3𝑘2 𝑞𝑛,𝑘1𝑘3 𝑞𝑛,𝑘2𝑘1 −𝑞𝑚,𝑘 3𝑘2 −𝑞𝑚,𝑘1𝑘3 −𝑞𝑚,𝑘2𝑘1 ] ⋅ [[ [ 𝛿𝑘 1𝑘𝛽 𝛿𝑘 2𝑘𝛽 𝛿𝑘 3𝑘𝛽 ] ] ] + [𝜕𝑔𝑘𝑑 𝜕]𝑑 𝑚 𝜕𝑔𝑑 𝑘 𝜕]𝑑 𝑛 ] ⋅ [ 𝑞𝑛,𝑘4𝑘3 𝑞𝑛,𝑘2𝑘4 𝑞𝑛,𝑘3𝑘2 −𝑞𝑚,𝑘 4𝑘3 −𝑞𝑚,𝑘2𝑘4 −𝑞𝑚,𝑘3𝑘2 ] ⋅ [[ [ 𝛿𝑘 2𝑘𝛽 𝛿𝑘 2𝑘𝛽 𝛿𝑘 2𝑘𝛽 ] ] ] , 𝑞𝑙,𝑘 𝛼𝑘𝛽 = 𝑞𝑙,𝑘𝛼− 𝑞𝑙,𝑘𝛽, (𝑙𝑚𝑛) = (123) , (231) , (312) = (𝑥𝑦𝑧) , (𝑦𝑧𝑥) , (𝑧𝑥𝑦) . (A.6) For dihedral angles, we have
𝜕𝑔𝑑 𝑘 𝜕𝜇𝑑 𝑙 = cos 𝜙𝑘1𝑘2𝑘3𝑘4] 𝑑 𝜇𝑑𝜇𝑑𝑙 − ]𝑑𝑙, 𝜕𝑔𝑑𝑘 𝜕]𝑑 𝑙 = cos 𝜙𝑘1𝑘2𝑘3𝑘4𝜇 𝑑 ]𝑑]𝑑𝑙 − 𝜇𝑑𝑙. (A.7)
Finally, Jacobian determinants in (7) are given as 𝐽𝑘𝑘= 𝜕𝑔 𝑝 𝑘(𝑡 + Δ𝑡) 𝜕𝜆𝑝𝑘(𝑡) =∑𝑁 𝑖=1 3 ∑ 𝑙=1 [𝜕𝑔 𝑝 𝑘(𝑡 + Δ𝑡) 𝜕𝑞𝑙,𝑖(𝑡 + Δ𝑡) 𝜕𝑞𝑙,𝑖(𝑡 + Δ𝑡) 𝜕𝜆𝑝𝑘(𝑡) ] =∑𝑁 𝑖=1 3 ∑ 𝑙=1 [𝜕𝑔 𝑝 𝑘(𝑡 + Δ𝑡) 𝜕𝑞𝑙,𝑖(𝑡 + Δ𝑡) 𝜕𝑔𝑝𝑘(𝑡) 𝜕𝑞𝑙,𝑖(𝑡)] , (A.8) where𝑁 = total atom number and 𝑙 = 𝑥, 𝑦, 𝑧.
B. Unified One-Passage Semiclassical
Nonadiabatic Transition Probability
A unified one-passage semiclassical nonadiabatic transition probability was derived analytically as [73]
𝑝 = sinh[(𝑑 − 1) 𝛿]
sinh(𝑑𝛿) exp(−𝛿) , (B.1)
where𝑑 represents the ratio between adiabatic and diabatic potential energy gaps at transition point𝑅0:
𝑑 = [𝑊+(𝑅0) − 𝑊−(𝑅0) 𝑉2(𝑅0) − 𝑉1(𝑅0) ]
2
, (B.2)
where 𝑊+ and 𝑊− correspond to the excited state (upper state) and the ground state (lower state), while𝑉2and𝑉1are
diabatic potential energy surfaces. The adiabatic parameter𝛿 that measures nonadiabatic transition strength is expressed in terms of the complex phase integrals:
𝛿 = Im1 ℎ∫ 𝑅∗ 𝑅0 [√2𝑚 (𝐸 − 𝑊−(𝑅)) −√2𝑚 (𝐸 − 𝑊+(𝑅))] 𝑑𝑅, (B.3)
in which𝑊+(𝑅∗) = 𝑊−(𝑅∗) at the complex crossing point 𝑅∗(its real part𝑅0stands for transition point),𝑚 is reduced mass of system, and𝐸 is collision energy. In order to apply the above formula to multidimensional case, we generalize (B.3) in terms of adiabatic potential energy surfaces at the transition point 𝑅0. Applying the exponential model (see (2.20) of [41]), we can convert (B.3) into
𝛿 = 𝜋√2𝑚(Δ𝐸0) 2 ℎ(𝜕/𝜕𝑥)[𝑊+(𝑅) + 𝑊−(𝑅)]𝑅=𝑅 0 × 1 √𝐸 − 𝐸0 1 1 + √1 + Δ𝐸0/ (𝐸 − 𝐸0) √𝑑, (B.4) where Δ𝐸0= 𝑊+(𝑅0) − 𝑊−(𝑅0) 2 , 𝐸0= 𝑊+(𝑅0) + 𝑊−(𝑅0) 2 . (B.5)
It should be noted that (B.4) is exact for the exponential model and it becomes an approximation when it is general-ized to the general two-state system. In comparison with the Landau-Zener transition formula in (12), we can decompose the adiabatic parameter 𝛿 into the following two diabatic parameters: 𝑎2= ℎ2 16𝑚 1 (Δ𝐸0)3[ 𝜕 𝜕𝑅[𝑊+(𝑅) + 𝑊−(𝑅)]𝑅=𝑅0] 2 , (B.6) 𝑏2= 𝐸 − 𝐸0 2Δ𝐸0 [ 1 2+ 1 2√1 + Δ𝐸0 (𝐸 − 𝐸0) √𝑑] 2 . (B.7)
Equations (B.6) and (B.7) are similar to the expressions derived by Zhu and Nakamura [3], but adiabatization form is much simpler. An effective coupling parameter𝑎2in (B.6) is very useful to determine nonadiabatic transition zone and the seam surface for multidimensional potential energy surfaces, and an effective collision energy𝑏2in (B.7) has an finite value at𝐸 = 𝐸0. The parameter𝑑 defined in (B.2) which represents general type of nonadiabatic transition is only included in (B.7) for the effective collision parameter.
C. The Two Diabatic Parameters and
Trajectory Surface Hopping along
the Seam Surface
C.1. The Effective Coupling𝑎2and the Seam Line. The
effec-tive coupling parameter 𝑎2 in (B.6) is defined along the certain one-dimensional direction, and the best direction is represented by the nonadiabatic coupling vector. However, we know that the direction perpendicular to the seam surface can be approximately considered as the direction of the nonadiabatic coupling vector. The seam surface can be calculated from two adiabatic potential energy surfaces without knowing any information of nonadiabatic coupling [69], in which the method requires to solve the certain partial differential equation. This is not easy to be generalized to the high-multidimensional case. We propose in the present work to calculate the local maxima of the effective coupling parameter𝑎2, and then the collection of all the local maxima constructs the seam surface. The kinetic part of Hamiltonian in terms of two dihedral angles D1 and D2 as shown in (9) can be rewritten as 𝑇 = −2𝐼ℎ2 D1 𝜕2 𝜕D12 − ℎ2 2𝐼D2 𝜕2 𝜕D22 = −ℎ2 2𝐼( 𝜕2 𝜕𝑥2 + 𝜕2 𝜕𝑦2) , (C.1)
where𝐼 = √𝐼D1𝐼D2 is the reduced moment of inertia that
corresponds to the reduced mass required in (B.6). The variables𝑥 and 𝑦 are given by
(𝑥𝑦) = (𝑐1 0 0 𝑐2) (D1D2) = ( √ 𝐼D1 𝐼 0 0 √ 𝐼D2 𝐼 ) (D1D2) . (C.2)
At conical configuration point (where D10 = −146.6∘ and D20 =−60.0∘), we estimate from (10) and (11) that 𝐼D2 =
1342893.5 a.u. and 𝐼D1 = 3453.3 a.u. so that 𝐼 = √𝐼D1𝐼D2 =
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𝑥0and𝑦0 (see inFigure 8) by (C.2), and then we assume that𝑥1 and𝑦1 are a predefined point on the seam with𝑠 direction normal to distance direction between (𝑥0, 𝑦0) and(𝑥1, 𝑦1); thus the point (𝑥, 𝑦) along 𝑠 coordinate can be expressed as
𝑥 = 𝑠 cos (𝜃 +𝜋2) + 𝑥1,
𝑦 = 𝑠 sin (𝜃 +𝜋2) + 𝑦1.
(C.3)
Derivative of potential energy surfaces defined in (B.6) with respect to𝑠 at (𝑥1, 𝑦1) can be estimated as
𝜕𝑊±(𝑥, 𝑦) 𝜕𝑠 𝑠=0 = 𝜕𝑊± 𝜕𝑥 𝑠=0cos(𝜃 + 𝜋 2) + 𝜕𝑊𝜕𝑦± 𝑠=0 sin(𝜃 +𝜋 2) . (C.4)
This equation has to be transformed to the original variables D1and D2in terms of which the potential energy surfaces are given by 𝜕𝑊±(D1, D2) 𝜕𝑠 𝑠=0 = 𝜕𝑊± 𝑐1𝜕D1𝑠=0cos(𝜃 + 𝜋 2) + 𝜕𝑊± 𝑐2𝜕D2𝑠=0sin(𝜃 + 𝜋 2) , (C.5) where𝑐1and𝑐2are defined in (C.2), and𝜃 is given as
𝑡𝑔𝜃 = (𝑦1− 𝑦0) (𝑥1− 𝑥0) =
𝑐2(D2 − D20)
𝑐1(D1 − D10). (C.6)
Finally, we can have an expression for effective coupling parameter 𝑎2=16𝐼ℎ2 1 (Δ𝐸0)3( 𝜕 𝜕𝑠[𝑊++ 𝑊−]) 2 𝑠=0 =16𝐼ℎ2 1 (Δ𝐸0)3 1 𝑡𝑔2𝜃 + 1 × (−1 𝑐1( 𝜕𝑊+ 𝜕D1+ 𝜕𝑊− 𝜕D1) 𝑡𝑔𝜃 + 1 𝑐2( 𝜕𝑊+ 𝜕D2+ 𝜕𝑊− 𝜕D2)) 2 (C.7) in which Δ𝐸0= 𝑊+(D1, D2) − 𝑊−(D1, D2) 2 . (C.8)
When𝜃 is rotating with respect to the conic point (𝑥0, 𝑦0), the local maxima 𝑎2 can be obtained. In this way, we can determine one point on seam surface as well as𝑎2 at that point. Then we can use this point as reference point to search next point on seam surface. Finally, we can determine whole seam surface and 𝑎2 distribution on seam surface simultaneously.
C.2. The Effective Collision Energy 𝑏2 and the Momentum Changes during Hopping. Classical trajectory is running in
the Cartesian coordinate system, while the surface hopping is taking place in the system of internal dihedral angles D1 and D2. We need project internal coordinates to Cartesian coordinates after hopping; besides we cannot use all kinetic energy for hopping and only part of kinetic energy along
y x O (x0, y0) s0 1 2 (x1, y1) x y s 𝜃 𝜃 + 𝜋/2 𝜃 + 𝜋/2
Figure 8: Schematic diagram in which𝑠 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
𝑇 = 12𝐼D1( ̇D1)2+12𝐼D2( ̇D2)2 =1 2𝐼 ( ̇𝑥2+ ̇𝑦2) = 1 2𝐼 (V21+ V22) , (C.9)
where dot stands for derivative with respect to time𝑡, and relation between( ̇D1, ̇D2) and ( ̇𝑥, ̇𝑦) is given in (C.2). The last equality in (C.9) is defined as
(V1
V2) = (−cossin𝜃 sin 𝜃) (𝜃 cos 𝜃 ̇𝑥̇𝑦) , (C.10) whereV1is angular velocity along𝑠 direction and V2is angular velocity along the direction normal to𝑠 as shown inFigure 8. If we assume that(V1, V2) change to (V1, V2) after hopping, then they should satisfy
V1= 𝑘V1, V2= V2, (C.11) whereV2is unchanged and𝑘 can be estimated from energy conservation as
𝑘 = √1 ±4Δ𝐸0 (𝐼V2
1)
(C.12) withΔ𝐸0 = (𝑊+ − 𝑊−)/2, in which “+” (“−”) represents hopping from upper (lower) potential to lower (upper) potential. We discuss the classical allowed vertical hop in the present system, namely,𝑘 ≥ 0 since the collision energy is quite large. We now need to convert(V1, V2) to ( ̇𝑥, ̇𝑦) that represent angular velocity change after hop along 𝑥 and 𝑦 direction shown inFigure 8. Equations (C.10) and (C.11) lead to the following relation:
̇𝑦
cos𝜃 − ̇𝑥sin𝜃 = 𝑘 ( ̇𝑦 cos 𝜃 − ̇𝑥 sin 𝜃) ,
̇𝑥
cos𝜃 + ̇𝑦sin𝜃 = ̇𝑥 cos 𝜃 + ̇𝑦 sin 𝜃, (C.13)
from which we can obtain ( ̇𝑥̇𝑦) = ( cos
2𝜃 + 𝑘 sin2𝜃 (1 − 𝑘) sin 𝜃 cos 𝜃
(1 − 𝑘) sin 𝜃 cos 𝜃 𝑘 cos2𝜃 + sin2𝜃 ) ( ̇𝑥̇𝑦) .
(C.14) This can be converted to original dihedral angles D1 and D2 system as
( ̇D1̇D2) = (𝑇𝑇11 𝑇12 21 𝑇22) (
̇D1 ̇D2)
= ( cos2𝜃 + 𝑘 sin2𝜃 𝑐22(1 − 𝑘) sin 𝜃 cos 𝜃
𝑐12(1 − 𝑘) sin 𝜃 cos 𝜃 𝑘 cos2𝜃 + sin2𝜃 ) × ( ̇D1̇D2) ,
(C.15) where ( ̇D1, ̇D2) stand for angular velocity after hop in original dihedral angle system, and𝑐1 and𝑐2 are defined in (C.2). The effective collision energy for 𝑏2 in (B.7) is then given by collision energy
𝐸 = 1 2𝐼 (V21) + 𝑊± = 1 2𝐼(𝑐2 ̇D2 cos 𝜃 − 𝑐1 ̇D1 sin 𝜃) 2+ 𝑊 ±. (C.16)
Next, we derive relation between internal angular veloc-ities and their corresponding velocveloc-ities in Cartesian coordi-nate 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)⇀𝜇⇀] − ⃗𝜇 ⋅ ⃗] = 0, (C.17) where ⃗𝜇 = ⃗𝑎 × ⃗𝑏 and ⃗V = ⃗𝑏 × ⃗𝑐 (seeScheme 1for definition of bond vectors ⃗𝑎, ⃗𝑏, and ⃗𝑐). Taking derivative with respect to time𝑡 in (C.17) with consideration that all bond angles and bond distances are fixed, we have
̇D1 = −𝑑 ( ⃗𝜇 ⋅𝑑𝑡 ⃗])[ ⃗𝜇⇀]sin(D1)]−1, (C.18) where 𝑑 ( ⃗𝜇 ⋅ ⃗]) 𝑑𝑡 = 𝑑 𝑑𝑡[( ⃗𝑎 ⋅ ⃗𝑏) ( ⃗𝑏 ⋅ ⃗𝑐) − ( ⃗𝑎 ⋅ ⃗𝑐) ( ⃗𝑏 ⋅ ⃗𝑏)] = − ⃗𝑏2𝑑 ( ⃗𝑎 ⋅ ⃗𝑐)𝑑𝑡 = − ⃗𝑏2[( ̇⃗𝑟2− ̇⃗𝑟1) ⋅ ⃗𝑐 + ⃗𝑎 ⋅ ( ̇⃗𝑟4− ̇⃗𝑟3)] = − ⃗𝑏2( ̇⃗𝑟2− ̇⃗𝑟1+ ̇⃗𝑟4− ̇⃗𝑟3) ⋅ ( ⃗𝑐 + ⃗𝑎) (C.19)
in which we utilized that all bond angles and bond distances are fixed again. Thus, we obtain
