Intersystem crossing-branched excited-state intramolecular proton transfer for o-nitrophenol: An ab initio on-the-fly nonadiabatic molecular dynamic simulation

The 6SA-CASSCF(10, 10)/6-31G (d, p) quantum chemistry method has been applied to perform on-the-fly trajectory surface hopping simulation with global switching algorithm and to explore excited-state intramolecular proton transfer reactions for the o-nitrophenol molecule within low-lying electronic singlet states (S0 and S1) and triplet states (T1 and T2). The decisive photoisomerization mechanisms of o-nitrophenol upon S1 excitation are found by three intersystem crossings and one conical intersection between two triplet states, in which T1 state plays an essential role. The present simulation shows branch ratios and timescales of three key processes via T1 state, non-hydrogen transfer with ratio 48% and timescale 300 fs, the tunneling hydrogen transfer with ratios 36% and timescale 10 ps, and the direct hydrogen transfer with ratios 13% and timescale 40 fs. The present simulated timescales might be close to low limit of the recent experiment results.

Scientific RepoRts | 6:26768 | DOI: 10.1038/srep26768 argon matrices 23 and by laser-induced-fluorescence 24 . The rich intersystem crossing network makes o-nitrophenol from the first excited singlet state decaying to the triplet states rapidly, and the process occurs in the femtoseconds to a few picoseconds time scale dependent on the intermediate states [25][26][27] . The early experiment showed that time scale decaying to the first excited triplet state is lower than 50 ps by measuring the UV excitation of o-nitrophenol in benzene solvent 28 . However, the recent experiment showed the existence of unstable aci-nitrophenol isomers and ultrafast ISC to the triplet manifold on a subpicosecond time scale by both femtosecond transient absorption spectroscopy in solution and time-resolved photoelectron spectroscopy in the gas phase 29 .
Up to now, neither fluorescence nor phosphorescence has been reported for o-nitrophenol and this indicates the existence of pure ultrafast radiationless decay process 29 . Therefore, an ab initio on-the-fly nonadiabatic molecular dynamics simulation is necessary to be performed to quantitatively analyze ESIPT process of o-nitrophenol via internal conversion and intersystem crossing network. This is motivation of the present study. By using ab initio quantum chemistry method at 6SA-CASSCF(10, 10)/6-31G (d, p) and MRCI/cc-pVDZ level, we previously optimized geometries for all isomers and transition states on two low-lying singlet (S 0 and S 1 ) and two low-lying triplet (T 1 and T 2 ) electronic states and we computed all four-state potential energy profiles along ESIPT coordinates 30 . We confirmed the existence of unstable aci-nitrophenol isomers as observed in the early studies 23,24,29,30 . We found total five ISC zones with SOCs at ~10 and ~40 wavenumbers. By performing nonadiabatic molecular dynamics simulation in the present study, we actually found one new S 0 /T 1 ISC which governs decay process for non hydrogen transfer. The previous computational and experimental study has proposed two possible relaxation pathways 29 ; one is via the ISC following (O)H-O(NO) stretch and the other is via the S 0 /S 1 CI after overcoming potential energy barrier induced by hybrid torsion of HONO group and (O)H-O(NO) stretch. Actually, the new S 0 /T 1 ISC is located near Frank-Condon region while S 0 /S 1 CI is located in configuration after hydrogen transfer. The present simulation has shown that the new S 0 /T 1 ISC is very effective pathway and there is almost no trajectory going via S 0 /S 1 CI. We can expect that there must have strong intersystem crossing-branched ESIPT. Besides, the existence the certain potential barriers on S 1 and T 1 states along ESIPT coordinate can add tunneling effects as well. The present study should reveal entire picture of photoisomerization reaction and the photo-decay mechanisms including ESIPT for o-nitrophenol molecule.
Trajectory-based ab initio nonadiabatic molecular dynamics simulations have been successfully applied to photoisomerization and photoreaction processes involving intersystem crossings and conical intersections 31-48 . Recently, within Tully's fewest switching algorithm in diabatic representation dynamic simulations have been carried out for intersystem crossings [43][44][45][46][47][48] . We have developed trajectory-based ab initio nonadiabatic molecular dynamics simulations without involving calculation for nonadiabatic couplings 37 , and it has been applied for multiple-state azobenzene photoisomerization via conical intersections 38 . Now we extend our analytically global switching probability method to include intersystem crossing with spin-orbital couplings. Trajectories are calculated on on-the-fly ab initio singlet or triplet potential energy surfaces and then by applying global switching algorithm we treat trajectory surface hopping in a unified way for internal conversion and intersystem crossing processes. Intersystem-crossing dynamic simulation with global switching surface hopping algorithm is not really new with Landau-Zener switching probability; the spin-diabatic and spin-adiabatic dynamic simulations 31 were performed for the model system of spin-orbital couplings that are used as diabatic coupling parameter in Landau-Zener formula, and moreover good agreement with experimental observation for the triplet to singlet branching ratio has been achieved with Landau-Zener formula for performing real dynamic simulation of O( 3 P) + ethylene 32 . We utilize improved Landau-Zener formula, namely Zhu-Nakamura to treat intersystem-crossing dynamic simulation in the present report. Moreover, tunneling effects are treated by one-dimensional semilcassical method along ESIPT coordinate 49 .

Results
Global switching algorithm. Global switching algorithm makes trajectory hopping at the time t where reaches local maximum (U 2 and U 1 are adjacent two adiabatic potential energy surfaces, while V 2 and V 1 are corresponding two diabatic potential energy surfaces). Analytical switching probability can be generally expressed in terms of d(t) as 50 in which δ is estimated from two potential energy surfaces and its gradients at hopping spot along on-the-fly running trajectory. In the case of d  1, it basically goes to Landau-Zener or improved version Zhu-Nakamura formula 51,52 , bewteen two adiabatic potentials 1 exp( 2 ) bewteen two diabatic potentials , (2) where (3) 2 2 4 in which the effective coupling parameter a 2 and effective collision energy b 2 are given by  where F 1 and F 2 are forces on two diabatic potential energy surfaces, V 12 is diabatic coupling, μ is reduced mass and E X is energy at crossing point and E t is potential energy plus kinetic energy component in direction of hopping direction. All those quantities in equation 4 are calculated along on-the-fly running trajectory at hopping spot 37 and its details are also given in the Supplementary Note 1. The hopping direction in the present theory is defined based on sort of the local modes and this agrees with the normal modes constructed from the regularized diabatic states 53 . Global switching algorithm has been extensively compared with the Tully's fewest switching algorithm 37,40 and two algorithms are basically similar for highly averaged quantities like quantum yields and lifetimes. In the present on-the-fly simulation, along three consecutive time steps we detect minimum energy separation between two singlet-states (S 1 and S 0 ) or two triplet-states (T 2 and T 1 ) and this gap is considered as diabatic coupling 2V 12 in equation 4. In this case, the switching is considered as in between two adiabatic potential energy surfaces. On the other hand, we detect the local maximum d(t) between the singlet and triplet states. The diabatic coupling is considered as spin-orbital coupling as V 12 = SOC in equation 4 as proposed in ref. 43., and in fact, singlet state V 1 and triplet state V 2 can be transformed into adiabatic representation according to However, in this case, we can do the switching directly between the singlet and triplet states as the spin-orbital coupling is known.
Intersystem crossings with ab initio dynamics. We have previously performed ab initio quantum chemistry calculations at 6SA-CASSCF(10, 10)/6-31G (d, p) level and have compared with MRCI/cc-pVDZ energy corrections for two singlet states (S 0 and S 1 ) and two triplet states (T 1 and T 2 ). It is well-known that the CASSCF method is lack of dynamic correlation. Therefore, we carried out energy correction by multi-reference configuration interaction (MRCI) calculation at CASSCF optimized geometries. Both vertical excitation energies (see Fig. S1(a) at CASSCF and Fig. 1S(b) at MRCI given in Supplementary Note 2) and adiabatic energies (see Fig. S2(a) at CASSCF and Fig. 2S(b) at MRCI given in Supplementary Note 2) at all key geometries show the same energy sequences. Relative energy differences between CASSCF and MRCI are mostly smaller than 0.1 eV, especially for energy gaps at six intersystem crossings and one conical intersection being even smaller (see Table S1). Therefore, we think that the present on-the-fly dynamical simulation based on CASSCF level should present the almost same results as it is based on MRCI level. We have confirmed 6SA-CASSCF(10, 10)/6-31G (d, p) is suitable choice for the present dynamic simulation.
We found there are two intersystem crossings between S 1 and T 2 (S 1 T 2 -IC1 and S 1 T 2 -IC2) around which SOC can be almost considered as constant 10 cm −1 , one between S 1 and T 1 (S 1 T 1 -IC) around which SOC is 40 cm −1 , three between S 0 and T 1 (S 0 T 1 -IC1, S 0 T 1 -IC2 and S 0 T 1 -ICX) around which SOCs are close to 40 cm −1 . The S 0 T 1 -ICX is newly found in the present study by trajectory surface hopping dynamic simulation, and the other five were optimized previously 30 . Furthermore, the present dynamic simulation shows that only three of six are active for ESIPT reaction, namely S 1 T 2 -IC1, S 0 T 1 -ICX, and S 0 T 1 -IC1 as shown in Fig. 1. The S 1 T 1 -IC is less active, and S 1 T 2 -IC2 and S 0 T 1 -IC2 are not active at all (key geometries and Cartesian coordinates for all six intersystem crossings are given in Supplementary Note 2 and Note 4). Figure 1 shows that the dihedral angels C 4 C 5 O 11 H 12 are zero for S 0 T 1 -ICX and S 1 T 2 -IC1, and the O 11 H 12 bond lengths are 0.945 Å and 0.955 Å, respectively. However, unlike to the planar geometry of S 1 T 2 -IC1, the nitro group of S 0 T 1 -ICX and S 0 T 1 -IC1 show a similarity although the H 12 bonded to O 11 in S 0 T 1 -ICX instead of O 15 in S 0 T 1 -IC1. The nitro group of them is out of the aromatic skeleton. Those geometry differences essentially govern the three intersystem crossings undergoing distinct ESIPT reaction pathways. Interestingly, the structure of S 0 T 1 -ICX is close Franck-Condon geometry and the structure of S 0 S 1 -CI is close to S 0 T 1 -IC1 geometry (see Fig. S3 given in Supplementary Note 2). The S 0 T 1 -ICX is responsible for relaxing pathway of non-hydrogen transfer, while the S 0 T 1 -IC1 responsible for relaxing pathway of after-hydrogen transfer.
The initial condition of trajectories is started from Franck-Condon region of o-nitrophenol. We perform frequency calculation at the ground state equilibrium geometry of S 0 to obtain the normal mode coordinates. Initial normal-mode coordinates and velocities of trajectory are selected according to the Wigner distribution on S 0 state. Finally, these initial normal-mode coordinates and velocities are converted into Cartesian coordinates and velocities plus vertical excitation energy to excited S 1 state. The thermal kinetic energy with T = 300K is added to all sampling trajectories with randomly distributing into initial Wigner velocities. However, such equally distributed initial conditions are not suitable for stimulating ESIPT reaction, there are the certain vibronic modes enhanced more than the others as molecule absorbed light to be excited to S 1 state. We found that four normal modes (O 11 H 12 stretch (4172 cm −1 ), C 5 O 11 H 12 bend (1473 cm −1 ), and two C 5 O 11 H 12 -C 4 N 13 O 15 scissor (408cm −1 and 310 cm −1 )) involved in the ESIPT reaction path have stronger vibronic couplings than the others. Therefore, we added extra kinetic energy (equivalent to T = 500 K) to these four normal modes. Along an on-the-fly running trajectory, the nuclear coordinates and velocities in Cartesian coordinates are propagated by numerically integrating the Newtonian equation of equation of motion with the velocity-Verlet method 54 . We determine the minimum potential energy gap between S 1 and S 0 states, and between T 1 and T 2 states at the conical intersection zones, and determine maximum d (t) between singlet and triplet states at intersystem crossing zones, at which we compute the effective coupling parameter a 2 and the effective collision energy b 2 . We found SOC varies very slowly at intersystem crossing zones, so that we choose constants accordingly 10 and ~40 cm −1 in simulation (when the energy gap is smaller than 0.25 eV, we start to check attempted trajectory surface hopping). We have run several expensive on-the-fly SOC trajectories in comparison with fixed SOC trajectories and results show small difference. We have checked for four active intersystem crossings, relative spin-orbital coupling variations (δ SOC/SOC) are about less than 2% around crossing zones in which trajectory hops take place. This fixed SOC technique was also utilized for performing real dynamic simulation of O( 3 P) + ethylene 32 .
Due to the present dynamics involving OH stretch vibration (over 3500 cm −1 ) along ESIPT reaction path, we made test runs with time steps of 0.05 fs, 0.1 fs and 0.15 fs, and finally we set up the 0.1 fs time step from the beginning to the end for the entire dynamic simulation. The dynamics simulation time is set up as 500 fs and as we increase up to 1000 fs, there is no notable difference. All the quantum chemical calculations of on-the-fly potential energy surfaces and its gradients are carried out at 6SA-CASSCF(10, 10)/6-31G (d, p) level by using the quantum chemistry package MOLPRO 2009.1. 55 and the dynamic simulation is carried by our own code.
On-the-fly trajectory analysis. We have run total 280 trajectories and the outcomes show that 6 stay on S 1 state regarded as resonance, 6 go via S 1 T 1 -IC as direct hydrogen transfer via T 1 state, and 268 go via S 1 T 2 -IC1 initially as shown in Fig. 2 (this is basically Fig. 7 in ref. 30). After going via S 1 T 2 -IC1, 265 trajectories hop to T 1 states via conical intersection between T 1 and T 2 states, and 3 go direct hydrogen transfer via T 2 state. Further three bifurcations take place for the 265 trajectories on T 1 state, 101 regarded as resonance on T 1 state, 133 going via S 0 T 1 -ICX as non-hydrogen transfer, and 31 as direct hydrogen transfer via T 1 state. This is overview how trajectories decay through competing nonadiabatic pathways in the singlet and triplet excited-state manifold. The percentage distribution of various nonadiabatic transition pathways is shown in Fig. 3a. The present ab initio dynamic simulation shows that the S 1 → T 2 → T 1 → S 0 (61.29%) and S 1 → T 2 → T 1 (34.77%) are the dominant processes, while the S 1 → T 1 → S 0 (0.72%) and S 1 → T 1 (1.43%) are kind of rare cases. Initially starting from Frank-Condon region on S 1 state, the majority sampling trajectories run in the electronic configuration where three excited states (S 1 , T 2 and T 1 ) are energetically close together, T 2 state always keeps in between the S 1 and T 1 states for a while. On the other hand, the present ab initio dynamic simulation confirms that intersystem crossing zone between the S 1 and T 2 states with small SOC (10 cm −1 ) are much wider than intersystem crossing zone between the S 1 and T 1 states with large SOC (40 cm −1 ). Therefore, without performing ab initio dynamic simulation, we could not realize that S 1 → T 2 intersystem crossing (ISC-S 1 T 2 ) channel is paramount relaxing pathway. Actually, similar situation is found that the efficient ultrafast ISC dynamic can still be taken place with a rather small spin-orbit coupling due to the counterbalance mechanism 56,57 . Further detailed distribution is shown in Fig. 3b for the sampling trajectories passing through the isomerization of the o-nitrophenol and aci-nitrophenol. In the present simulations, we defined the aci-nitrophenol isomers as aci-isomer1 and aci-isomer2, respectively. For aci-isomer1, it still exhibits intramolecular hydrogen bond with ortho oxygen atom of benzene and while for aci-isomer2, the hydrogen rotates with nitro group and far away from the ortho oxygen atom. Diverse distributions are shown, for instance, that o-nitrophenol overcomes barrier of hydrogen transfer and reaches the aci-isomer1 (13.33%), after getting over the hydrogen transfer barrier it passes to aci-isomer1 and then back to o-nitrophenol (42.22%), the o-nitrophenol converts to aci-isomer1 and then continue to produce aci-isomer2 with the inversion of O 15 H 12 (13.33%), the o-nitrophenol changes to aci-isomer1and then regenerated again after going via aci-isomer2 (13.33%), and so on. Starting from Frank-Condon region on S 1 state, among 280 sampling trajectories there are 3, 37, and 5 trajectories undergoing direct hydrogen transfer on T 2 , T 1 and S 0 states, respectively. However there are 6 and 101 trajectories undergoing tunneling hydrogen transfer on S 1 and T 1 states, respectively. Although direct hydrogen transfer on T 2 is rare case, it presents a new perspective for general ESIPT reactions. In most of situation, the T 2 state looks more like a "hub", trajectories stop on it for a moment and then quickly decay to the reactive T 1 state where is major channel for direct hydrogen transfer. The barrier for hydrogen transfer on T 1 state is lower than that of on S 1 state as shown in Fig. 2. On the other hand, the present dynamic simulation shows that kinetic energy from the four vibartional normal modes responding for hydrogen transfer easily dissipates to the other vibrational modes on S 1 state much faster than on T 1 state. This is part of reason that direct hydrogen transfer does not occur on S 1 state. Hydrogen transfer has a peculiar property on S 0 state which differs from that of on T 1 state, trajectory reaches to S 0 state via multi-steps continuous hops in a relatively short period of time and then it runs on the ground state, eventually the hydrogen transfer reaction occured and aci-isomer formed. It should be noted that the aci-isomer can undergo backward hydrogen transfer reaction to o-nitrophenol due to very low barrier (nearly barrier less). This is reason that hydrogen transfer occurs on S 0 state not very often. The newly found intersystem crossing S 0 T 1 -ICX has its geometry close to S 0 -NP, so that 133 out of 280 go via it and those trajectories relax to ground state by non-hydrogen transfer pathway. In brief summary, we conclude that S 0 T 1 -ICX serves as the center player of intersystem crossing for three competing processes, 45 trajectories results in direct hydrogen transfer, 107 in tunneling hydrogen transfer, and 128 in non-hydrogen transfer.
Although every trajectory has its own lifetime when surface hopping occurs after photoexcitation, and the final result is estimated by the average assemble of all sampling trajectories. Figure 4 shows average population decay with respect of time from excited states S 1 , T 2 and T 1 , respectively. A single exponential function curve fitting is applied to calculate average lifetimes on those excited states. The lifetime of S 1 state was estimated to be about 8 fs as an ultrafast decay process and population is almost diminished after 25 fs. This is understandable in common sense that the first excited state decay more likely undergoes via intersystem crossing rather than vibrational cooling 58,59 . The decaying process is also an ultrafast on T 2 state and the lifetime of T 2 state is estimated to be 14 fs. The population on T 2 state is quickly raised first and it reaches maximum at 10 fs and than it almost diminishes after 40 fs. The most of population via T 2 state transfers very quickly to T 1 state via conical intersection between T 1 and T 2 states. The population on T 1 state reaches its maximum (∼ 93%) after 40 fs and the lifetime of T 1 state is estimated to be around 1000 fs. This is the most likely to contribute to the long-lived excited species that should be observed in both gas and liquid phase because the population on T 1 state decays very slowly. However, the real lifetimes on excited states from experimental observation might be longer than the present estimation due to some extra kinetic energy being initially put in four vibrational modes which are enhancing hydrogen transfer dynamics.
The triplet T 1 state plays an essential role for intersystem crossing-branched ESIPT photoisomerization dynamics in which 97% sampling trajectories bifurcate into three major distinct decay channels. The first branch is the non-hydrogen transfer that counts for 47% sampling trajectories, and timescale is about 100 fs to 500 fs (average is about 300 fs). For trajectories running on T 1 state, once they form out-of plane synchronous motion by O 11 H 12 vibrating around C 5 O 11 bond and the nitro group vibrating around C 4 N 13 bond simultaneously, trajectories can reach S 0 T 1 -ICX to decay to ground S 0 state. The present simulation shows if this out-of plane synchronous motion form coherent motion before trajectories reach S 0 T 1 -ICX (we checked those trajectories from 500 fs to 1000 fs), they stay on T 1 state as resonance trajectories. This is the second branch as the tunneling hydrogen transfer that counts for 36% sampling trajectories. We have roughly estimated thermal and microcanonical rate constants of nonadiabatic chemical reaction along ESIPT coordinate as shown in Fig. 2 r in which β is thermal energy (temperature is set up T = 300 K), Z γ is reactant partition function and P(E) is nonadiabatic transmission probability for the tunneling trajectories. The average tunneling timescale estimated from equation 6 is about 10 ps. This tunneling timescale is about the same as tunneling proton transfer reaction for tropolone molecule 60 . The third branch is the direct hydrogen transfer that counts for 13% sampling trajectories, and timescale is from 20 fs to 100 fs (average is about 40 fs). In this case, trajectories can go the direct hydrogen transfer before the out-of plane synchronous motion occurs. Then, the trajectories continue to run on the T 1 state until they reach S 0 T 1 -IC1 accompanying with O 14 N 13 O 15 H 12 group out-of-plane rotating and finally hop to S 0 state.
Typical trajectories. Figure 5 shows the non-hydrogen transfer for a trajectory that hops to the S 0 state via  Figure 6 shows the direct hydrogen transfer on T 1 state for a trajectory that takes the aci-nitrophenol isomerization in the ground state finally. This trajectory hops from the S 1 to T 2 state via S 1 T 2 -IC1 at the 16.6 fs and from the T 2 to T 1 state via conical intersection at the 18.2 fs. Sequentially, the hydrogen migration is occurred around   rotates slightly around the C 4 N 13 bond and the transferred H 12 turns to away from the original O 11 , finally the aci-isomer2 is observed in this trajectory. Figure 7 shows that the trajectory makes the hydrogen transfer on the T 1 and back hydrogen transfer on S 0 state. Staring from Frank-Condon region on S 1 state, the trajectory hops from the S 1 to T 2 state via S 1 T 2 -IC1 at the 7.5 fs and from the T 2 to T 1 state via conical intersection at the 9.1 fs. Then, the hydrogen transfer is complete and the aci-nitrophenol species is formed at 35 fs on T 1 state. During trajectory evolution, the distance between O 15 and H 12 is shortened gradually until bond forming, and bond angles of C 5 O 11 H 12 and C 5 C 4 N 13 change drastically while dihedral angles of C 4 C 5 O 11 H 12 and C 5 C 4 N 13 O 15 fluctuate smoothly. At 86.1 fs, trajectory decays to S 0 state via S 0 T 1 -IC1. Soon after 130.0 fs, the back hydrogen transfer takes place on S 0 state and regenerates the o-nitrophenol. In the rest of evolution, the trajectory shows a periodic fluctuation with H 12 vibrates around the O 11 H 12 bond. We have plotted more cases in Supplementary Note 3 for various more complicated fast and slow hydrogen transfer processes.
The present trajectory simulation indicates that the aci-nitro isomers are mostly generated on the triplet states, especially on T 1 state and this confirms the previous assumption in which the HONO split-off motion is taken place in the triplet manifold 29,61,62 . The present simulation also indicates that the o-NP dynamic decay via ISC ESIPT is mostly in time scale of subpicosecond (expect tunneling ESIPT) and this is in close agreement with the observations for a large number of nitrated polycyclic aromatic compounds [63][64][65] . The present simulation shows that there excites the wide Franck-Condon region with very small energy gap between the first excited singlet S 1 and the second triplet T 2 states, and this can facilitate the fast ISC radiationless process. Therefore, trajectories hop from S 1 to T 2 state along the (O)H-O(NO) barrierless stretching pathway and in the tunneling zone the ESIPT pathways can generally have lower barrier via triplet states than via S 1 state. The present simulation shows that the (O)H-O(NO) stretch related vibrational modes can enhance ESIPT and the torsion motion of nitro group would hinder the hydrogen transfer reaction.

Conclusion
By applying ab initio on-the-fly nonadiabatic molecular dynamic simulation for intersystem crossing-branched ESIPT for o-nitrophenol with use of global switching trajectory surface hopping algorithm, we have estimated reaction branch ratios and timescales of various ESIPT processes. As is summarized in Fig. 8, there are three decay branches, the first one is non-hydrogen transfer via newly found S 0 T 1 -ICX with ratio 47.5% and average timescale 300 fs, the second branch is tunneling hydrogen transfer with ratios 36.08% and 2.14% taken place on T 1 and S 1 states, respectively and timescale 10 ps, and the third branch is direct hydrogen transfer with ratios 13.21% and 1.07% taken place on T 1 and T 2 states, respectively and timescale 40 fs. However, the real lifetimes on excited states from experimental observation might be longer than the present simulated timescales due to some initial extra kinetic energy maybe enhance hydrogen transfer dynamics 29 . Finally, we believe that the present trajectory-based on-the-fly nonadiabatic molecular dynamic simulation method can be generally applied to ultrafast photophysical and photochemical processes in a unified way involved in conical intersections and intersystem crossings for large-scale simulation.