Generalized trajectory surface-hopping method for internal conversion and intersystem crossing

Trajectory-based fewest-switches surface-hopping (FSSH) dynamics simulations have become a popular and reliable theoretical tool to simulate nonadiabatic photophysical and photochemical processes. Most available FSSH methods model internal conversion. We present a generalized trajectory surface-hopping (GTSH) method for simulating both internal conversion and intersystem crossing processes on an equal footing. We consider hops between adiabatic eigenstates of the non-relativistic electronic Hamiltonian (pure spin states), which is appropriate for sufficiently small spin-orbit coupling. This choice allows us to make maximum use of existing electronic structure programs and to minimize the changes to available implementations of the traditional FSSH method. The GTSH method is formulated within the quantum mechanics (QM)/molecular mechanics framework, but can of course also be applied at the pure QM level. The algorithm implemented in the GTSH code is specified step by step. As an initial GTSH application, we report simulations of the nonadiabatic processes in the lowest four electronic states (S0, S1, T1, and T2) of acrolein both in vacuo and in acetonitrile solution, in which the acrolein molecule is treated at the ab initio complete-active-space self-consistent-field level. These dynamics simulations provide detailed mechanistic insight by identifying and characterizing two nonadiabatic routes to the lowest triplet state, namely, direct S1 → T1 hopping as major pathway and sequential S1 → T2 → T1 hopping as minor pathway, with the T2 state acting as a relay state. They illustrate the potential of the GTSH approach to explore photoinduced processes in complex systems, in which intersystem crossing plays an important role.


I. INTRODUCTION
9][40][41][42] To simulate these two types of nonadiabatic processes on an equal footing, it is desirable to reformulate the original FSSH method and to include spin-orbit interactions that may induce nonadiabatic intersystem crossings.a) Present address: Key Laboratory of Theoretical and Computational Photochemistry, Ministry of Education, College of Chemistry, Beijing Normal University, Beijing 100875, China.Electronic mail: ganglong.cui@bnu.edu.cn.b) Electronic mail: thiel@mpi-muelheim.mpg.de4][45][46][47][48][49][50][51][52][53] Direct "on-the-fly" trajectory-based FSSH simulations of photoinduced processes that account for SOC have however been reported only recently.Richter et al. 54,55 proposed a general FSSH method, in which the relevant physical quantities are subjected to unitary transformations between diabatic and adiabatic representations to enable the treatment of arbitrary couplings like SOC.Recently, Mai et al. 56 applied this method to study the nonadiabatic dynamics of SO 2 .Granucci et al. 57 formulated both spin-diabatic and spin-adiabatic FSSH methods and investigated their merits and limitations for model systems.Their subsequent applications at the semiempirical level addressed the photodynamics of acetone and 6thioguanine taking into account SOC in the spin-adiabatic framework. 58,59 abenicht and Prezhdo 60 presented a formulation within time-domain Kohn-Sham density functional theory and applied their approach to study the time-dependent evolution of singlet and triplet excited states of carbon nanotubes.
To our knowledge, corresponding generalized FSSH methods that can treat both internal conversion and intersystem crossing have not yet been reported so far within the quantum mechanics/molecular mechanics (QM/MM) framework. 61,62 uch a combined approach is expected to S1 S0 T1 ISC IC FIG. 1. Scope of the standard and the generalized FSSH method.The standard treatment covers nonadiabatic transitions (solid red line) leading to internal conversion (IC) between states with the same spin (black curves), for example, S 1 → S 0 .The generalized method also describes nonadiabatic transitions (dashed red line) giving rise to intersystem crossing (ISC) between states with different spin (black and blue curves), for example, S 1 → T 1 or T 1 → S 0 .be most useful for modeling nonadiabatic internal conversion and intersystem crossing processes in complex systems, e.g., molecules in solution and biological systems.
In this article we present a generalized FSSH method at the QM/MM level without invoking a diabatic representation of the electronic states.The exclusive use of adiabatic eigenstates of the non-relativistic electronic Hamiltonian offers the following benefits: (i) we can make maximum use of existing electronic structure programs; (ii) we minimize the need for code modifications in available implementations of the original FSSH method; and (iii) spin-orbit couplings can be computed "on-the-fly" (like the required potential energy surface information, i.e., energies, gradients, and derivative couplings).5][56][57] The scope of the original and the generalized FSSH methods is illustrated in Fig. 1.The treatment of intersystem crossing between states with different spin is added (red dashed line between blue and black curves), without modifying the treatment of internal conversion between states with the same spin (red solid line between black curves).
The article is structured as follows: In Sec.II, we present the formalism of the generalized FSSH method at the QM and QM/MM levels.Section III describes the actual algorithm and its implementation.In Sec.IV, we report an initial application to acrolein in vacuo and in acetonitrile solution.Sections V and VI offer some discussion and conclusions.

II. GENERALIZED TRAJECTORY SURFACE-HOPPING (GTSH) METHOD
The original FSSH method 1, 2 was devised for describing radiationless transitions between same-spin states.It employs the non-relativistic electronic Hamiltonian operator which depends on the electronic coordinates, r, and parametrically on the time-dependent nuclear coordinates, R(t).T (r) denotes the kinetic energy operator of the electrons, and the potential energy operator V (r, R(t)) includes inter-electronic and inter-nuclear repulsions as well as electronnucleus attractions.The basic equations of the original FSSH method are documented in the supplementary material for easy reference. 63

A. The generalized FSSH method
In order to treat intersystem crossing, the electronic Hamiltonian must include the spin-orbit interaction, which drives radiationless transitions between different spin states.We thus add the spin-orbit operator to the zero-order electronic Hamiltonian and consider it as a perturbation, where r and s correspond to spatial and spin coordinates of electrons.The time-dependent Schrödinger equation now reads The time-dependent electronic wavefunction can be expressed in terms of adiabatic zero-order electronic wavefunctions, These zero-order wavefunctions are eigenfunctions of Ĥ0 and of the spin operators S 2 and S z .They can be written as a sum of products of space and spin parts, which as a whole must transform like the totally antisymmetric irreducible representation (irrep) under simultaneous permutation of the electronic space-spin coordinates, with the spin parts forming a basis for an irrep in spin space. 64,65 herefore, matrix elements of spin-free operators contain a spin factor ji that is one (zero) between states of the same (different) spin.The space part 0 i (r, R(t)) is an eigenfunction of the zero-order Hamiltonian Ĥ0 (r, R(t)) (not Ĥ (r, R(t), s)) at nuclear coordinates R(t).
After inserting Eq. ( 4) into Eq.( 3), multiplying by 0 j (r, R(t), s)| from the left-hand side, and integrating over electronic spatial and spin coordinates, we obtain where is the spin-orbit matrix element between electronic states j and i at nuclear coordinates R(t), and ji (t is the nonadiabatic coupling term, i.e., the dot product of the nuclear velocities v(t) and the derivative coupling vector d ji (R(t)) (see the supplementary material 63 for details).We thus arrive at our final generalized working equation This is the central equation of the generalized FSSH method.It describes radiationless transitions between states with the same spin and with different spins, i.e., internal conversion and intersystem crossing.
For radiationless transitions between same-spin states, H so ji (t) is assumed to vanish. 66Hence, Eq. ( 6) simplifies to the working expression of the traditional FSSH approach (see the supplementary material 63 ). 1,2 or radiationless transitions between different spin states, the derivative coupling contribution vanishes because of ji = 0, so Eq. ( 6) becomes

B. The generalized fewest-switches criterion
In the generalized FSSH method, the fewest-switches criterion must also account for hops caused by spin-orbit coupling.The derivation is analogous to the one given by Tully (see the supplementary material 63 ). 1,2 he main difference is that the use of Eq. ( 6) introduces an additional spin-orbit term into the transition probability from state i to state j:

C. Spin-orbit couplings
Our generalized FSSH method is a spin-diabatic approach in the notation of Granucci et al. 57 They stress in their analysis that the total intermultiplet transition probability between two states of different spin should be rotationally invariant in such a scheme. 57This can be achieved by assuming that each spin multiplet can be considered as a single electronic state -an approximation that has commonly been used in previous studies of SOC-mediated processes. 43,45,48,50,52,53,57 We adpt this approximation and thus need to determine an effective SOC value between the two interacting multiplets.
In the literature, several approximate schemes are available for this purpose.9]57 Another option is to evaluate the SOC value from half the energy difference between full calculations with and without the spinorbit coupling term in the Hamiltonian. 60,67 n a semiempirical configuration interaction (CI) framework, a one-electron effective spin-orbit Hamiltonian 68,69 can be used to compute approximate SOC values. 58,59 n our approach, we calculate the SOC matrix elements at each step of the simula-tion ("on-the-fly") at the ab initio complete active space selfconsistent field (CASSCF) level using the atomic mean-field approximation 70,71 and then condense them into an effective SOC value for a given pair of states.For the latter step, several procedures are possible.In the present pilot application (see below), we need effective SOC values for transitions between singlet and triplet states, which we approximate by averaging over the interactions between the singlet and the three triplet substates (i.e., as one-third of the corresponding norm).In schematic notation, With this choice, the effective SOC values for singlet-triplet transitions are the same in either direction and of similar magnitude as those obtained from previous formulations. 57urther work is needed in the future to explore better representations of the effective spin-orbit coupling in GTSH simulations.

D. QM/MM implementation
In the following, we present our implementation of the GTSH approach in the QM/MM framework 61,62 (for details see the supplementary material 63 ).In an additive QM/MM scheme, the energy of electronic state I is given by (10)   where R qm and R mm are the atomic coordinates of the QM and MM subsystems; E I,QM (R qm ), E I,MM (R mm ), and E I,QM−MM (R qm , R mm ) denote the QM energy of the QM subsystem, the MM energy of the MM subsystem, and the QM-MM interaction energy, respectively.In our approach, we use a standard force field for calculating E I,MM (R mm ), without distinguishing between the different electronic states I.For sufficiently large QM subsystems, the approximation E I,MM (R mm ) E MM (R mm ) is expected to be valid because the electronic excitation and the relevant photoinduced processes will be well localized in the QM region.For analogous reasons, standard MM force fields are considered adequate for evaluating the bonded and van der Waals QM-MM interaction terms in all states.Adopting the usual electronic embedding, 72 the electrostatic QM-MM interaction terms are calculated explicitly at the QM level by incorporating the MM charges into the one-electron Hamiltonian.If covalent bonds are broken in QM/MM setups, the QM-MM boundary can be treated by a variety of approaches including the link atom scheme, 62 the generalized hybrid orbital scheme, 73 and the pseudobond scheme. 74In our QM/MM implementation, we use the link atom approach in combination with a charge-shift scheme as implemented in the ChemShell3.5 package. 75he GTSH method requires derivative coupling vectors (as in the original FSSH method 1,2 ) and spin-orbit couplings.In the QM/MM framework, derivative coupling vectors between electronic states I and J with the same spin can be decomposed into QM and MM components resulting from differentiation with respect to the nuclear coordinates in the QM and MM regions, respectively.In our implementation, only the QM components of the derivative coupling vectors are used in QM/MM GTSH dynamics simulations because the MM components are expected to be very small; in addition, they are available only numerically so that their evaluation is expensive.Similar approximations have been used in previous QM/MM FSSH applications 7,12,14,24,[76][77][78] and in QM/MM optimizations of conical intersections. 79,80 n the same spirit, the spin-orbit couplings between electronic states I and J with different spin are determined from the QM wavefunctions obtained in the QM/MM calculations with electronic embedding,

III. ALGORITHM AND IMPLEMENTATION
The algorithm for running GTSH simulations at the QM level can be summarized as follows: In GTSH simulations at the QM/MM level, the MM contributions to the energy and forces are computed and added to the corresponding QM terms after steps 3 and 5.
We have implemented this algorithm in the GTSH package that is interfaced with our local version of the ChemShell3.5 package. 75Codes used for data input, output, and exchange are written in the C language, while computationally more demanding parts such as the electronic propagation with complex numbers are written in Fortran77.Electronic energies, gradients, derivative couplings, and spinorbit couplings are computed using external electronic structure programs, e.g., MOLPRO2010 82 and MNDO99. 83In QM/MM applications, the MM contributions to the energy and forces are computed using the DL-POLY code incorporated into ChemShell.The data exchange between the GTSH package and external programs is realized by flexible Linux scripts.The overall structure of the GTSH package is shown in Fig. 2.

IV. GTSH SIMULATION OF ACROLEIN
As an initial illustrative application, we present GTSH simulations of acrolein in vacuo (QM) and in acetonitrile solution (QM/MM) to demonstrate the feasibility of our approach.In these simulations, we consider only the lowest four electronic states (S 0 , S 1 , T 1 , and T 2 ).

A. Computational details
We applied the CASSCF method for the QM computations, which is capable of providing the derivative coupling vectors and spin-orbit couplings efficiently and with reasonable accuracy. 84Computational efficiency is very important in GTSH simulations because hundreds of "on-the-fly" trajectories are needed to achieve statistically meaningful results.In all CASSCF calculations of acrolein in vacuo and in solution, we adopted an active space of 8 electrons in 7 orbitals, which has also been employed in previous theoretical work. 85The 6-31G basis set was used throughout. 86

B. QM simulations
Initial conditions (velocities and coordinates) of 250 trajectories were generated using Wigner sampling. 81Whenever the energy gap was less than 15 kcal/mol, the generalized FSSH criterion, Eq. ( 8), was applied to decide whether to hop.We chose time steps of 1.0 fs for nuclear propagation and 0.002 fs for electronic propagation.The CASSCF GTSH dynamics simulations for acrolein in vacuo were carried out using the GTSH package. 87The energies, gradients, derivative coupling vectors, and spin-orbit couplings in atomic meanfield approximation (AMFI) 70,71 were computed "on-the-fly" using MOLPRO2010. 82

C. QM/MM simulations
Acrolein was first solvated in a spherical acetonitrile box of 25 Å radius (599 solvent molecules; 3602 atoms in total) using Packmol scripts. 88The solvent molecules were then FIG. 4. Acrolein in acetonitrile solution.Color code: QM acrolein in yellow, MM atoms allowed to move in blue, and frozen MM atoms in red.The atomic numbering of acrolein is same as in Fig. 3. See text for further details.minimized (2000 steps) and equilibrated (5 ps MD simulations at 300 K) at the MM level.Missing parameters were fitted through the SwissParam website service. 89The final classical MD snapshot was the starting point for a subsequent 2 ps QM(OM2)/MM MD simulation. 90Snapshots from this QM/MM trajectory were used to prepare the initial conditions for 200 QM/MM GTSH nonadiabatic dynamics runs.In all QM/MM computations, only acrolein and solvent molecules within 15 Å from its center of mass were allowed to move (1358 atoms), whereas all other atoms were frozen at the positions adopted at the end of the initial classical 5 ps MD simulations.Fig. 4 illustrates the QM/MM setup.
All MM computations were carried out using CHARMM. 91All QM (OM2)/MM MD computations were performed using ChemShell3.5. 75All QM (CASSCF)/MM FSSH dynamics simulations for acrolein in acetonitrile solution were done with the GTSH package; 87 the QMrelated energies, gradients, derivative coupling vectors, and spin-orbit couplings were computed by MOLPRO2010, 82 and the MM-related energies and gradients were obtained using a locally modified ChemShell3.5 version. 75he GTSH nonadiabatic dynamics were run in an analogous manner at the QM and QM/MM levels, except that the initial conditions were sampled from the Wigner distribution of the ground vibrational state in vacuo (QM) and from a thermal MD trajectory in acrolein (QM/MM), since the use of the preferred Wigner distribution is not practical in the QM/MM case.Hence, the initial excess of vibrational energy is expected to be larger in vacuo than in acrolein solution.

D. Time-dependent state populations
Fig. 5 shows the in vacuo time-dependent state populations of the lowest four electronic states (S 0 , S 1 , T 1 , and T 2 ) as determined from the 162 successful GTSH dynamics trajectories at the CASSCF level.The S 1 population monotonously decreases from 1 to 0.10 at 500 fs, whereas the S 0 and T 1 populations increase from zero to 0.38 and 0.52, respectively.The T 1 population rises notably faster than the S 1 population, and the T 2 population remains negligible in the whole simulation time.
The acetonitrile solvent influences the time-dependent state populations of acrolein considerably (Fig. 6).In this case, 165 of the 200 QM/MM GTSH trajectories finished successfully.Within 500 fs, the S 0 population decreases to 0.25, less so than in vacuo (0.10), indicating that the solvent slows down the excited-state deactivation.Correspondingly, the S 1 and T 1 populations increase more slowly, to 0.28 and 0.47, respectively.The T 2 population again remains small during the whole simulation time of 500 fs; nevertheless, the T 2 state plays an important role for S 1 deactivation (see below).

E. Distribution of hopping points
An analysis of all hops between the S 0 , S 1 , T 1 , and T 2 states in all GTSH trajectories shows that the T 2 state is often  involved in the deactivation of the S 1 state.We first discuss the results in vacuo (see Table I and Fig. 7).
The S 1 → T 2 and T 2 → S 1 intersystem crossings occur in 73 (23%) and 51 (16%) trajectories, respectively; the T 2 → T 1 internal conversion happens 19 times (6%).Hence, the T 2 state can act as a relay electronic state in the S 1 → T 2 → T 1 deactivation channel.This dynamic feature originates from the special topology of the S 1 , T 1 , and T 2 potential energy surfaces.Previous high-level electronic structure calculations 85 established that the S 1 and T 1 states are of nπ * character, while the T 2 state is of ππ * character in the S 1 /T 1 /T 2 three-state intersection region.According to El-Sayed's rule, 92 the S 1 → T 2 intersystem crossing will be strongly favored over the S 1 → T 1 intersystem crossing in this region.Therefore, the T 2 state can act as a relay state, which then undergoes a fast T 2 → T 1 internal conversion to the T 1 state.
The frequent S 1 → T 2 and T 2 → S 1 hops are facilitated by the energetic proximity of these two states and the fact that acrolein remains almost planar in these states.By contrast, we did not see any T 1 → S 1 hop in any of the 500 fs trajectories.This may be due to the fast relaxation of acrolein in the T 1 state towards its minimum geometry with a twisted C=C double bond.Once this happens, the T 1 → S 1 hop is energetically forbidden, because the relaxed T 1 state is much lower in energy than the S 1 and T 2 states.
The direct S 1 → T 1 nonadiabatic hops occur most often, 116 times (37%).Closer examination reveals that acrolein initially stays in a nearly planar conformation where the S 1 , T 1 , and T 2 states are energetically close to each other (see the trajectories below).In certain of these initially visited regions, the T 1 state acquires some ππ * character that makes the S 1 → T 1 nonadiabatic transition become efficient.
There are three types of hops that can populate the S 0 state.Among them, the T 1 → S 0 hops are predominant (51 times, 90%), whereas the S 1 → S 0 and T 2 → S 0 hops are very rare (together only 7 times in all trajectories).
In acetonitrile solution, the total number of hops is markedly reduced, from 317 in vacuo to 231 (see Table I and Fig. 7).This is probably due to the steric interactions between acrolein and surrounding solvent molecules that make it harder for acrolein to reach the distorted geometries in the conical intersection region with small energy gaps, which will decrease the chances for hopping events.The back-and-forth hops between T 2 and S 1 are partly suppressed in solution.However, the relative frequency of the different nonadiabatic hops does not change much when going from the gas phase to acetonitrile solution (see Fig. 7).

F. Four typical trajectories
We first address a typical example of a direct S 1 → T 1 intersystem crossing.Fig. 8 shows the time evolution of the nonadiabatic couplings, spin-orbit couplings, potential energies, and several key geometric parameters in a representative trajectory.The S 1 → T 1 hop takes place at 308 fs when there is a large S 1 -T 1 spin-orbit coupling of 0.1 kcal/mol (35 cm −1 , blue line in the middle-left panel) and a small S 1 -T 1 energy gap of 1.5 kcal/mol.Acrolein does not access the T 2 relay state in this trajectory, and, therefore, the large T 2 /T 1 nonadiabatic coupling at this hopping point (500 ps −1 , green line in the middle-right panel) has no relevance.As shown in the bottom panel, there are no pronounced structural changes in this trajectory.
The second example shows the importance of the S 1 /T 2 /T 1 three-state intersection region for S 1 excited-state deactivation.4][95][96][97] Here, we directly see their role in the excited-state nonadiabatic dynamics.In the trajectory depicted in Fig. 9, a S 1 → T 2 intersystem crossing takes place after 150 fs, when the spin-orbit coupling is computed to be 0.1 kcal/mol (35 cm −1 ).After another 9 fs, there is a T 2 → T 1 internal conversion triggered by a large nonadiabatic coupling of 765 ps −1 .During this small time interval, the three involved states are almost degenerate in energy (see the top panel of Fig. 9), which facilitates the observed sequential nonadiabatic transitions.Evidently, the T 2 state acts overall as a relay electronic state and mediates the S 1 → T 1 nonadiabatic transition.After another 190 fs (i.e., at 350 fs), acrolein in the T 1 state starts to twist around the H3-C1-C2-C6 dihedral angle, accompanied by a shortening of the C6-O7 bond and an elongation of the C1-C2 bond (see the bottom panel).These structural changes reflect the change of the electronic character of the T 1 state from nπ * * .
The third example (see Fig. 10) concerns an indirect S 1 → S 0 nonadiabatic transition mediated by the T 1 state.In this trajectory, the first S 1 → T 1 transition occurs at 29 fs due to a large spin-orbit coupling of 0.1 kcal/mol.After 70 fs the T 1 state starts to relax from its initial planar conformation towards its twisted minimum.At 130 fs, this T 1 minimum is reached.It has a broken C1=C2 double bond and a nearly perpendicular H3-C1-C2-C6 arrangement (see the bottom panel), indicating ππ * character of the T 1 state at this point.After another 80 fs (i.e., at 210 fs), the second T 1 → S 0 intersystem crossing happens at a geometry where the T 1 -S 0 energy gap is small (see the top panel) and the spin-orbit coupling is also surprisingly small (3 cm −1 ).Once the S 0 state is reached, both the C1-C2 and C6-O7 bond lengths decrease quickly and then oscillate around their ground-state equilibrium values; the H3-C1-C2-C6 dihedral angle twists back gradually, and internal rotation around the C2-C6 single bond starts.
The final example is a direct S 1 → S 0 internal conversion (see Fig. 11), which is only found rarely (in 6 trajectories).Here, the only nonadiabatic hopping takes place at 157 fs when there is a large nonadiabatic coupling value of 739 ps −1 and a small energy gap of 1.8 kcal/mol.After returning to the S 0 state, the C6=O7 bond length quickly relaxes to typical ground-state values around 1.25 Å (see the bottomleft panel), and acrolein evolves towards the most stable trans conformation (see C1-C2-C6-O7 dihedral angle).
The preceding four examples were all taken from QM GTSH simulations of acrolein in vacuo.The QM/MM GTSH trajectories of acrolein in acetonitrile solution show qualitatively similar behavior and will therefore not be discussed here.

G. Photophysics of acrolein: Summary
Previous electronic structure calculations 85,98 and the present nonadiabatic dynamics simulations indicate that once the S 1 state is populated, there are two nonadiabatic pathways to populate the lowest T 1 triplet state.The major one is direct S 1 → T 1 hopping, which accounts for 37% and 44% of all hopping events in vacuo and solution, respectively; the minor one is the sequential S 1 → T 2 → T 1 pathway, in which the T 2 state acts as a relay electronic state.In addition, we have observed a few trajectories that decay directly from the S 1 to the S 0 state.Within 500 fs, the S 0 state is reached by 38% (25%) of the trajectories in vacuo (in acetonitrile solution).In both cases, the pathway via the T 1 state is favored.
The relevance of the S 1 /T 2 /T 1 three-state intersection region for the formation of the lowest triplet state T 1 has been  shown here for the first time by direct nonadiabatic dynamics simulations.However, it should be stressed that the sequential S 1 → T 2 → T 1 hopping mediated by the relay T 2 state is the minor channel for populating the T 1 state compared with the direct S 1 → T 1 intersystem crossing.Finally, we note that the 193 nm photodynamics of acrolein in acetonitrile solution has been studied experimentally using time-resolved Fourier transform infrared spectroscopy. 99The dominant reaction channel was found to be 1,3-hydrogen migration in the T 1 state (78%), accompanied by some α C-H cleavage in the S 1 state (12%).This is compatible with the results from our GTSH dynamics simulations of acrolein in acetonitrile solution, which give preferential population of the T 1 state after 500 ps (Fig. 6).

V. DISCUSSION
An alternative to spin-diabatic treatments (such as GTSH) is the spin-adiabatic approach, 57 which considers hops between the eigenstates of the full Hamiltonian Ĥ (r, R(t), s) with spin-orbit terms, Eq. ( 2).This approach is formally attractive for surface hopping simulations since FIG.9. A typical trajectory with a S 1 → T 2 intersystem crossing (first vertical line), which is immediately followed by a T 2 → T 1 internal conversion (second vertical line).See text for discussion.
spin-orbit effects are properly included at the level of the spin-coupled eigenstates, and transitions between these eigenstates are thus mediated solely by the nonadiabatic coupling (as in the original FSSH method). 57Hence, there is no formal need to introduce effective spin-orbit interactions.It has also been shown for one-dimensional model systems that the spin-adiabatic approach is able to reproduce the results from full quantum calculations quite well even in difficult situations where spin-diabatic treatments fail, and it has been recommended as the method of choice when more than one intermultiplet crossing is important and/or when the spinorbit coupling is large. 57However, surface hopping simulations in the spin-adiabatic framework require the derivatives of the eigenenergies of the full Hamiltonian Ĥ (r, R(t), s), which have up to now only been implemented at the semiempirical CI level using a one-electron effective spin-orbit Hamiltonian, 68,69 and corresponding applications have thus been done only at this level. 58,59 ence, in the absence of a viable ab initio implementation of the spin-adiabatic approach, spin-diabatic treatments are currently the only practical option for ab initio trajectory surface hopping simulations that aim at covering both internal conversion and intersystem crossing.One crucial approximation in spin-diabatic treatments is the definition of the effective spin-orbit coupling (see Sec. II C).Our choice is based on the norm of the CASSCFcalculated SOC matrix elements, Eq. ( 9); the effective SOC value is real and contains no phase factor.The relevance of including a phase factor in spin-diabatic surface hopping simulations has been studied for singlet-triplet ISC processes in a one-dimensional model system through comparison with full quantum calculations; inclusion of the phase factor was found to be important when the SOC changes sign at the position of the singlet-triplet crossing, while it could safely be neglected when the SOC changes sign far away from this crossing. 57We also note in this context that our choice of the effective spinorbit interaction, Eq. ( 9), is not unique and that the adopted AMFI procedure 70,71   In our current simulations, we allow ISC hops between states only if the corresponding energy gap is less than 15 kcal/mol.There is no rigorous justification for this convention, which however seems physically reasonable in view of CASSCF-computed effective SOC values of the order of 0.1 kcal/mol or less (see Figs. 8-11).An obvious task for future work is to check and possibly adjust this criterion.
From a more general point of view, it is commonly accepted that FSSH simulations are well justified if the coupling between the relevant states is well localized in space.This condition is ideally satisfied by the nonadiabatic vibronic cou-pling, whereas spin-orbit coupling is usually considered to be quite nonlocal.At first sight, the fewest-switches condition in surface hopping and the nonlocal character of spin-orbit coupling would seem hardly compatible.Our GTSH simulations for acrolein provide some evidence to the contrary.First, the CASSCF-computed SOC values are not constant along the actual trajectories but often fluctuate quite strongly, albeit in a less spiky manner than the nonadiabatic vibronic couplings (see Figs. 8-11) Second, we indeed find only very few hops in any single trajectory (see Table I and Figs geometrically similar T 2 and S 1 states, i.e., once in 31% (18%) of the QM (QM/MM) trajectories.Apparently, the acrolein molecule relaxes quite fast after each hop such that the back hop becomes unfavorable because of the quickly increasing energy gap; this actually leads to a "fewest-switches scenario" even for SOC-mediated intersystem crossing.Fast ISC processes in real molecules may thus be more amenable to surface hopping simulations in a spin-diabatic framework than one might expect from formal considerations or model studies.We anticipate that this may hold in particular for systems in which the spin-orbit coupling is rather small (e.g., organic molecules) so that hops will occur mostly between nearly degenerate states (in regions that tend to be well localized in space).

VI. CONCLUSIONS
We present a GTSH method that can simulate internal conversion and intersystem crossing processes on an equal footing, both at the QM and QM/MM levels.The chosen approach allows us to make maximum use of electronic structure programs and to minimize the required changes in existing implementations of the traditional FSSH method.An initial application with ab initio GTSH simulations of acrolein in vacuo (CASSCF) and in acetonitrile solution (CASSCF/MM) demonstrates the potential of the GTSH approach for studying the full set of photoinduced processes including both internal conversions and intersystem crossings.The current developments pave the way for future photodynamics simulations of complex systems.

FIG. 3 .
FIG. 3. Acrolein (CH 2 CHCHO).Color code: Carbon atoms in yellow, oxygen atom in red, and hydrogen atoms in gray.Also shown is the atomic numbering.

FIG. 6 .
FIG.6.Time-dependent state populations of the lowest four electronic states (S 0 , S 1 , T 1 , and T 2 ) of acrolein in acetonitrile solution obtained from the present GTSH dynamics simulations.

FIG. 11 .
FIG.11.A typical trajectory with a single S 1 → S 0 internal conversion (first vertical line).See text for discussion.
. 5-7), and the only back hops occur between the energetically and Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions.Downloaded to IP: 192.108.70.170On: Thu, 09 Jun 2016 10:05:23 + dt) using an external electronic structure program; 6. interpolate relevant physical quantities (d ij and H so ji ) linearly between t and t + dt; 7. propagate coefficients of electronic states between t and t + dt based on the electronic equation of motion Eq.
1. generate initial atomic velocities v(t) and coordinates x(t) by means of Wigner sampling 81 or ground-state molecular dynamics (MD) simulations; 2. initialize propagation coefficients of the electronic states involved and random-number seeds; 3. compute energies, forces, derivative coupling vectors d ij (t), and spin-orbit couplings H so ji (t) at x(t) using an external electronic structure program; 4. propagate atomic coordinates from t to t + dt according to the velocity Verlet method; 5. compute energies, forces, derivative coupling vectors d ij (t + dt), and spin-orbit couplings H so ji (t + dt) at x(t

TABLE I .
Number of nonadiabatic hops between the lowest four electronic states (S 0 , S 1 , T 1 , and T 2 ) of acrolein: Results from 162 QM (1st row) and 165 QM/MM (2nd row) trajectoriesTotal S 1 → T 1 S 1 → T 2 S 1 → S 0 T 2 → T 1 T 2 → S 1 T 2 → S 0 T 1 → S 0 is just one of several possible quantumchemical methods for computing SOC values.It is obviously desirable in future work to explore other options in this regard.