We have performed trajectory surface hopping dynamics for cis,cis-1,3-cyclooctadiene to investigate the photochemical pathways involved after being excited to the S1 state. Our calculations reveal ultrafast decay to the ground state, facilitated by conical intersections involving distortions around the double bonds. The main distortions are localized on one double bond, involving twisting and pyramidalization of one of the carbons of that double bond (similar to ethylene), while a limited number of trajectories decay via delocalized (non-local) twisting of both double bonds. The interplay between local and non-local distortions is important in our understanding of photoisomerization in conjugated systems. The calculations show that a broad range of the conical intersection seam space is accessed during the non-adiabatic events. Several products formed on the ground state have also been observed.

Light-induced isomerization plays a key role in a myriad of fundamental natural processes.1 Such isomerization can also be used in molecular switches and to create molecular machines, where light can aid in conversion between two different stable forms or configurations of a molecule.2,3 The simplest and the most basic unit to gain insight into photoisomerization is ethylene, an unsaturated hydrocarbon having a single double bond. Hence, ethylene’s dynamics after photoexcitation are well understood as a result of both experimental and theoretical studies.4–9 In order to gain further insight into the excited state dynamics of extended conjugated systems, a multitude of studies have been performed on linear and cyclic dienes, such as butadiene (BD)10–17 and 1,3-cyclohexadiene (CHD).18–28 Although the vertical excitation character for both molecules is similar, the nuclear dynamics via conical intersections (CIs) and the photochemistry are different. UV excitation of BD leads to a competition between cistrans isomerization, electrocyclic ring closure to cyclobutene, and other products. On the other hand, photoexcitation of CHD leads to opening of the ring to form 1,3,5-hexatriene (HT). The excited state dynamics are also driven by whether the excitation is localized on one double bond resembling ethylene, which is the case in BD, vs delocalized on both double bonds, as seen in CHD.

One of the broader questions that we aim to address with this work is how the structure of conjugated systems affects their excited state dynamics. How does the number and position of double bonds and the remaining system affect the isomerization products and the lifetimes of excited states? How does the structure affect the excess kinetic energy released after returning to the electronic ground state? Are the excited state dynamics localized on individual double bonds or delocalized on both bonds? This last question has been discussed for BD and CHD in previous studies. The effect of localized vs delocalized dynamics has been examined in more detail as well by Schalk et al. by comparing CHD to 1,4-cyclohexadiene (1,4-CHD) and cyclohexene, leading them to coin the term “dynamophore.”29 1,4-CHD shows very similar dynamics and CIs to ethylene, cyclohexene, and BD.29Cis,cis-1,3-cyclooctadiene (cc-COD) is a cyclic π-conjugated diene, similar to CHD, but larger and more flexible. We wanted to compare the dynamics of cc-COD to those of CHD and BD. The structure of the three molecules is shown in Fig. 1. There are two alternating double bonds in all systems. Nevertheless, the larger size of the cc-COD cyclic system compared to CHD may have important implications for its dynamics. In BD, it was found that a large amplitude motion allows the localization of dynamics on one double bond rather than two (opposite of CHD). We investigate whether localization is possible in cc-COD, which is more flexible than CHD, but more constrained than BD. The general similarities with some distinctions in the electronic structure of the three systems allow for an important comparison of their dynamics.

FIG. 1.

(a) Structures of dienes: BD, CHD, and cc-COD (from the left to the right) (b) Main configurations describing the ground and excited states in cc-COD. For dynamics simulations, the HOMO-1, HOMO, and LUMO orbitals have been incorporated in the (4,3) active space. The active space at the CASSCF(4,3)/cc-pVDZ level is shown on the right side of panel (b).

FIG. 1.

(a) Structures of dienes: BD, CHD, and cc-COD (from the left to the right) (b) Main configurations describing the ground and excited states in cc-COD. For dynamics simulations, the HOMO-1, HOMO, and LUMO orbitals have been incorporated in the (4,3) active space. The active space at the CASSCF(4,3)/cc-pVDZ level is shown on the right side of panel (b).

Close modal

Unlike CHD and BD, there has not been much work on the excited state dynamics of cc-COD. Fuß et al. performed UV–IR femtosecond pump–probe measurements on cc-COD and predicted yields for possible photochemical products.30 These products are cis,trans-1,3-cyclooctadiene (ct-COD) that can form via a cistrans photoisomerization and cis-bicyclo[4.2.0]oct-7-ene (BCO) that can form via a photochemical electrocyclic ring closure. They also predicted the quantum yields to be 0.28 and 0.01, respectively. Nanosecond UV laser irradiation of cc-COD has also been performed to drive selective cyclization to BCO.31 Although these studies suggested certain photochemical mechanisms, there has been no theoretical work to explain the findings. A simulation of the non-adiabatic excited state dynamics can help evaluate these suggested mechanisms. To the best of our knowledge, there has only been one theoretical study related to cc-COD, which predated the experimental work and considered conical intersections of a methylated cc-COD.32 

In earlier work, we performed UV pump (4.77 eV) and VUV probe (7.94 eV) measurements of internal conversion of CHD and cc-COD and demonstrated that CHD shows substantial “hot” ground state ionization after internal conversion to the ground state, whereas cc-COD does not show below-threshold ionization.24 In that work, we located CIs for cc-COD and CHD and connected them to the Franck–Condon (FC) geometry (S0 minimum) and to the products using linear interpolations in order to demonstrate the presence of ultrafast excited state dynamics in both cc-COD and CHD. The detailed dynamics of cc-COD were neither examined nor compared to the dynamics in CHD. In this work, we report excited state dynamics simulations of cc-COD using the trajectory surface hopping methodology, in order to understand the non-adiabatic dynamics facilitated by CIs and compare and contrast its dynamics with similar linear and cyclic dienes.

cc-COD is substantially more flexible compared to CHD because of the presence of a —C—C—C—C— skeleton, which guarantees the presence of several conformers. Initial conformers were obtained using a conformational search using molecular mechanics, as implemented in Spartan.33 These conformers were then optimized at the Density Functional Theory (DFT)34,35 level using the B3LYP36–39 functional and the 6-31G(d)40 basis set. Single point (SP) energy calculations at the coupled-cluster with single, double, and perturbative triple excitations [CCSD(T)]41,42/cc-pVDZ43 level of theory were performed to achieve more accurate energies of the conformers to calculate their Boltzmann distribution at room temperature (298.15 K). Two conformers were found to populate the ground state with probabilities of 66% and 34%. The optimizations and SP calculations were performed using the Gaussian09 package.44 

The FC geometry of the lowest energy conformer was selected for the investigation of the excited states in cc-COD. Vertical excitation energies, oscillator strengths, gradients, and non-adiabatic couplings were calculated at the complete active space self-consistent field (CASSCF)45 and extended multi-state complete active space with second order perturbation (XMS-CASPT2) theory46–48 levels with the cc-pVDZ basis set. A Multi-configurational self-consistent field (MCSCF), of which the CASSCF is a specific type, combines Hartree–Fock or SCF theory with the configuration interaction method.49 In the CASSCF, a subset of orbitals important to the problem at hand are selected. These orbitals are called the active space, and all possible combinations of them are considered for optimizing the orbitals and coefficients. In doing so, the CASSCF (and MCSCF) ensures a correct description of the static correlation of the electrons. However, this method cannot treat the dynamic correlation which involves excitation of electrons among these configurations. We considered the multi-state version (specifically the “extended” version) of the CASPT2 method which takes into account the dynamic correlation using second order perturbation between the involved electronic states. The extended version (XMS-CASPT2) improves on the drawbacks of the MS-CASPT2 version when strongly mixed states are present.48 We have performed CASSCF and XMS-CASPT2 calculations, using the same cc-pVDZ basis set, but changing the active space in order to evaluate the proper state ordering in the FC region. Three states were averaged in each case. The XMS-CASPT2 calculations were performed using the corresponding CASSCF reference wavefunction with an imaginary shift of 0.2 au. The single-state single-reference (SS-SR)50 contraction scheme was employed for these calculations. Several different active spaces have been used. These are shown in Table II in a format (e,o) where e and o represent the number of electrons and orbitals in the active space, respectively. The CASSCF and XMS-CASPT2 calculations were performed with the Columbus 7.051–53 and Bagel54,55 packages, respectively.

The dynamics were calculated at the SA3-CASSCF(4,3)/cc-pVDZ level, and details about the dynamics will be described in Sec. II B. However, once the dynamics were over, several CI optimizations were started from S1/S0 hop geometries collected from different trajectories. These calculations were all performed at SA3-CASSCF(4,3)/cc-pVDZ level of theory using Columbus 7.0. CI searches between S1 and S2, as well as minimizations on S1, were done using SA3-CASSCF(4,3)/cc-pVDZ as well. We also performed linear interpolations in internal coordinates (LIICs) to some of the important CIs at CASSCF(4,3)/cc-pVDZ and XMS-CASPT2/CAS(6,6)/cc-pVDZ levels of theory to show the presence of pathways from the FC geometry to the CIs, after excitation to the S1 state. LIICs were also performed from the CIs to most of the products using the SA3-CASSCF(4,3)/cc-pVDZ level, which exhibits the presence of monotonic ground state pathways to the products. For some products, the CASSCF may not be the most appropriate method to describe them, and hence, we used DFT/B3LYP/6-31G(d) level of theory for the LIIC paths to confirm the dynamics. All products were optimized at the B3LYP/6-31G(d) level.

The ground state of all conformers of cc-COD was optimized at the DFT level using the B3LYP functional and 6-31G(d) basis set using the Gaussian09 package. The frequencies and normal modes were calculated at the same level of theory. Sampling was performed using a harmonic oscillator Wigner distribution in Newton-X56,57 to generate initial conditions (nuclear coordinates and velocities) based on the optimized geometry and the normal modes from the previous calculation. 200 initial conditions were generated for each of the conformers of cc-COD at 298 K. Vertical excitation energies and oscillator strengths were calculated for all initial conditions at the CASSCF level using the cc-pVDZ basis set. An active space of 4 electrons in 3 orbitals was employed and three states were averaged. Excitation energies and oscillator strengths of the initial conditions were used to calculate the absorption cross section and simulate the first absorption band of cc-COD. A Lorentzian line shape, temperature of 298 K, and a phenomenological broadening (δ) value of 0.2 eV were employed.

We performed non-adiabatic excited state dynamics simulations using trajectory surface hopping (TSH) in Newton-X on CASSCF(4,3)/cc-pVDZ full-dimensional potential energy surfaces58–62 (PESs) calculated on-the-fly using the Columbus 7.0 package. The trajectories were propagated starting from the S1 state which is the brighter of the two states (S1 and S2) in the FC region with the active space chosen for this study. The fewest switches surface hopping (FSSH)63 algorithm was employed to take into account non-adiabatic coupling (NAC) between the S2, S1, and S0 states. The velocity verlet algorithm was used to integrate Newton’s equations of motion with a time step of 0.5 fs. The semiclassical Schrödinger equation was integrated using fifth-order Butcher’ s algorithm with a time step of 0.005 fs. The simulations for both conformers were performed for 500 fs using XSEDE’s computational resources.64 

There are some issues in TSH that have to be addressed in order to properly run and interpret the trajectories. Some of the basic issues and how we dealt with them are outlined here. When a wavepacket splits between two different potential energy surfaces, the different branches move away in space leading to loss of coherence. This is not correctly described in TSH, however, since the propagation of every trajectory occurs on one PES. Hence, various decoherence corrections have been developed.65–68 Here, the approach of non-linear decay of mixing by Granucci and Persico66 has been used. In this approach, a decoherence time is defined as

τKL=EKEL1+αEkin,
(1)

where EL and EK are the energies of the current state, L, and any other (inactive) state, K; Ekin is the nuclear kinetic energy; and α is an empirical parameter with a recommended value of 0.1 hartree (which is the value used here).65 τKL is used in Eq. (2) to damp all the electronic coefficients of the inactive states (CK) at every time step, Δt. This leads to the new inactive state coefficient CK, which is then used to rescale the current state electronic coefficient (CL), as shown in Eq. (3),

CK=CKexpΔt/τKLKL,
(2)
CL=CL1KLCK2CL21/2.
(3)

Another obvious issue is that after a hop, electronic energy has to convert to nuclear kinetic energy conserving the total energy and/or momentum. The approaches that have been developed to address this involve rescaling of the momentum vector of all nuclei in order to conserve energy. An approach that was shown to perform better involves rescaling the momentum only along the derivative coupling vector,69–72 and it has been chosen here. Another problem occurs when a hop is predicted, but there is not enough kinetic energy in the system for it to occur. This is called a frustrated hop, since it is classically not allowed. To deal with this, the momentum direction could be left unaltered or be reflected at the frustrated hop. Here, we kept the momentum unaltered in its direction. The dynamics are affected by how each one of the above issues is addressed, and it is not always clear which is the best approach. A comprehensive study on these issues was published recently.73 The Wigner distribution used to obtain the initial conditions also has some issues, since it cannot treat anharmonicities and a distribution over many local minima with similar energies. The latter problem has been treated here by doing separate distributions and dynamics for the two conformers (local minima) of cc-COD.

cc-COD has 54 vibrational degrees of freedom. Due to the flexible nature of the —C—C—C—C— subunit in cc-COD, rotations around those sp3 carbons are possible, which lead to more than one conformer in the ground state. Using Spartan, we were able to find three conformers, the most stable of which has C2 symmetry. These conformers were optimized at the DFT/B3LYP/6-31G(d) level, and those optimized geometries were used to calculate better and more accurate energies at the CCSD(T)/cc-pVDZ level and calculate their Boltzmann distributions in the ground state at room temperature.

Table I shows that conformers 1 and 2 are both populated in the ground state of cc-COD at room temperature. Figure 2 shows the structure of the two conformers. As both are populated in the ground state, we decided to use both conformers for generating the initial conditions and, subsequently, running dynamics simulations.

TABLE I.

Energies at the CCSD(T)/cc-pVDZ level and Boltzmann distribution of conformers of cc-COD in the ground state at 298.15 K

ConformersRelative Energy (eV)Boltzmann Distribution (%)
0.00 66.056 
0.02 33.942 
0.27 0.002 
ConformersRelative Energy (eV)Boltzmann Distribution (%)
0.00 66.056 
0.02 33.942 
0.27 0.002 
FIG. 2.

Absorption spectrum calculated at the SA3-CASSCF(4,3)/cc-pVDZ level using 200 geometries for each conformer and a Lorentzian line shape with δ = 0.2 eV. The spectra from each conformer were combined using their relative populations. The structures of the two conformers are also shown here. Conformer one is also shown with atom numbering which will be used throughout this paper (also provided in the supplementary material).

FIG. 2.

Absorption spectrum calculated at the SA3-CASSCF(4,3)/cc-pVDZ level using 200 geometries for each conformer and a Lorentzian line shape with δ = 0.2 eV. The spectra from each conformer were combined using their relative populations. The structures of the two conformers are also shown here. Conformer one is also shown with atom numbering which will be used throughout this paper (also provided in the supplementary material).

Close modal

The energetic ordering of excited states in the FC region for linear and cyclic polyenes has been a contentious topic for a long time. Figure 1 shows an orbital energy diagram describing the basic characteristics of the first two excited states in dienes (such as cc-COD, BD, and CHD). The main excited state is a singly excited (1B) bright state corresponding to a ππ* excitation with B symmetry (using the C2 point group). A second excited state exists as well with similar energy, which is multi-configurational and has a significant doubly excited character (A symmetry), and is expected to be dark. The energetic ordering of these states, as well as their gap, is difficult to obtain both theoretically and experimentally. The dark state cannot be detected in absorption, making its experimental determination difficult. Consequently, the gap between these two states is uncertain.

We have used the FC geometry of the lowest energy conformer to decide the correct energetic ordering in cc-COD. Determining the correct energetic ordering is very important for the dynamics. In BD and CHD, it has been determined that the ordering is as shown in Fig. 1 [S1(B) and S2(A)], and high level calculations for cc-COD show the same ordering. Table II shows the results for vertical excitation energies and oscillator strengths at all the different levels of theory that have been employed for the FC geometry of the lowest energy conformer of cc-COD. XMS-CASPT2 calculations show that S1 is the bright state, while S2 has a smaller oscillator strength and has mostly a doubly excited character. All XMS-CASPT2 expansions yield similar results, with the best energy for the bright state (5.82 eV) being at the XMS-CASPT2 level using CAS(4,4) which is still 0.38 eV blue shifted from the experimental maxima at 5.44 eV (228 nm).31 Nonetheless, the experimental spectrum is in the solution phase, and it is not clear what the solvation effects are. CASSCF/cc-pVDZ calculations with all the active spaces other than (4,3) predict the wrong ordering of states. CASSCF/cc-pVDZ with CAS(4,3) gives the correct ordering and an energy gap between S1 and S2 similar to XMS-CASPT2/CAS(4,4)/cc-pVDZ. It predicts different ratios of oscillator strengths, but we expect this to have a smaller effect on the dynamics. In addition to switching state characters, CASSCF with CAS(4,4) and CAS(6,6) gives an energy difference of over 1 eV between the two states. Increasing the active space reduces the energy separation of the two states, but does not give the correct state character. In addition, the extra σ orbitals to be used in these larger active spaces can be very ambiguous. More information about the orbitals in each case is given in the supplementary material. Similar problems with switching of the state ordering have also been seen in CHD and BD.15,16,21,74

TABLE II.

Vertical excitation energies in electron volt, oscillator strengths (f), magnitude of gradient vectors for S1 and S2 states in hartree bohr−1, and non-adiabatic coupling (NAC) between S1 and S2 in bohr−1, for FC geometry for the lowest energy conformer of cc-COD. The cc-pVDZ basis set was used in all cases. The experimental value is taken from Ref. 31. The methods that predict the correct ordering (S1 state having B symmetry and S2 having A symmetry) are highlighted with boldface.

MethodE(S1) (f)E(S2) (f)g(S1)g(S2)NAC(S1S2)
SA3-CASSCF(4,3) 7.51 (0.69) 8.19 (0.29) 0.18 0.29 3.47 
SA3-CASSCF(4,4) 7.11 (1.4 × 10−28.39 (0.56) 0.37 0.19 0.47 
SA3-CASSCF(6,6) 6.89 (3.7 × 10−38.09 (0.51) 0.35 0.16 0.50 
SA3-CASSCF(8,8) 7.13 (6.8 × 10−38.00 (0.46)    
SA3-CASSCF(10,10) 7.11 (8.4 × 10−37.94 (0.44)    
XMS-CASPT2/CAS(4,3) 5.96 (0.35) 6.67 (0.16) 0.21 0.24 1.33 
XMS-CASPT2/CAS(4,4) 5.82 (0.27) 6.72 (2.4 × 10−30.22 0.30 0.79 
XMS-CASPT2/CAS(6,6) 5.91 (0.27) 6.72 (2.8 × 10−30.22 0.29 0.79 
XMS-CASPT2/CAS(8,8) 6.02 (0.28) 6.66 (5.5 × 10−3   
Exp 5.44     
MethodE(S1) (f)E(S2) (f)g(S1)g(S2)NAC(S1S2)
SA3-CASSCF(4,3) 7.51 (0.69) 8.19 (0.29) 0.18 0.29 3.47 
SA3-CASSCF(4,4) 7.11 (1.4 × 10−28.39 (0.56) 0.37 0.19 0.47 
SA3-CASSCF(6,6) 6.89 (3.7 × 10−38.09 (0.51) 0.35 0.16 0.50 
SA3-CASSCF(8,8) 7.13 (6.8 × 10−38.00 (0.46)    
SA3-CASSCF(10,10) 7.11 (8.4 × 10−37.94 (0.44)    
XMS-CASPT2/CAS(4,3) 5.96 (0.35) 6.67 (0.16) 0.21 0.24 1.33 
XMS-CASPT2/CAS(4,4) 5.82 (0.27) 6.72 (2.4 × 10−30.22 0.30 0.79 
XMS-CASPT2/CAS(6,6) 5.91 (0.27) 6.72 (2.8 × 10−30.22 0.29 0.79 
XMS-CASPT2/CAS(8,8) 6.02 (0.28) 6.66 (5.5 × 10−3   
Exp 5.44     

Table II also compares the magnitude of the gradients for the S1 and S2 surfaces as well as the non-adiabatic coupling vectors. The components of these vectors are given by gA(SI) = ∇A⟨ΨI|HI⟩ for the gradients and NACA(SISJ) = ⟨ΨI|∇AΨJ⟩ for the non-adiabatic (derivative) vector, where ΨI is the wavefunction of the I state and ∇A is the derivative with respect to nuclear coordinate A. CASSCF predicts a steeper S2 surface compared to XMS-CASPT2 which predicts more parallel S1 and S2 surfaces. The CASSCF method also gives the highest magnitude for the derivative coupling, indicating a very fast S1S2 transition. These values provide some evidence for differences between CASSCF and XMS-CASPT2, although we do not know exactly how these would manifest themselves in the dynamics.

Considering all of these factors, it became quite clear that a study at the XMS-CASPT2 level with CAS(4,4) or CAS(6,6) can provide the most accurate description of the dynamics. However, a TSH dynamics simulation for a molecule with 54 nuclear degrees of freedom at the XMS-CASPT2 level of theory with a significant number of trajectories is computationally extremely expensive. Hence, we selected the SA3-CASSCF(4,3)/cc-pVDZ level of theory for our TSH simulation, since it can reproduce the correct state ordering in the FC region and is computationally tractable with TSH. A comparison has been made previously for BD dynamics where both CASSCF(4,3) and MS-CASPT2 were used, and it was found that the dynamics were qualitatively correct at the CASSCF(4,3) level.15,16

Excitation energies and oscillator strengths were calculated for all initial conditions for both conformers as mentioned before in order to calculate the first absorption band. The results from both conformers were combined keeping in mind their Boltzmann distributions in the ground state. Figure 2 shows the normalized absorption spectrum of cc-COD. The theoretical absorption peak is at 7.02 eV, 1.58 eV higher than the experimental absorption peak at 5.44 eV.31 This is very typical of CASSCF, as it always overestimates excited state energies in the FC region since it lacks any treatment of dynamic correlation. It can also be seen from Table II that the brighter S1 state has an energy of 7.51 eV at the FC geometry at the same level of theory, which is also 0.49 eV higher than the theoretical maximum. This suggests that the experimental absorption maximum may not correspond to the vertical excitation energy because of FC factors.

Figure 3 shows the population of the three states for cc-COD for 500 fs. The population plot is a weighted average of the results from the two conformers for all valid trajectories. An important problem in TSH using multi-reference methods is that sometimes (especially after hopping to the S0 state) trajectories fail because the active space does not converge. This can happen if the proper antibonding orbitals are not present in the active space or if the molecule distorts significantly. This can engender a change in the orbitals in the active space leading to an energy conservation failure for a particular trajectory. Trajectories typically fail more often when an unbalanced active space is employed. Since we are working with a very flexible molecule and have employed a (4,3) active space which does not have an antibonding π* orbital corresponding to its bonding π orbital, a portion of the trajectories failed before the end of the calculation (500 fs), mostly after reaching S0. It is crucial to deal with this issue since the fraction of trajectories in each state depends sensitively on the fraction of trajectories that crash in each state. The default choice is that trajectories are excluded from the counting once they fail. However, in this case, when a significant portion of trajectories on the ground state fail, the excited state decay will appear to be slower than it is. Furthermore, after hopping to the S0 state, the potential energy surfaces S1 and S0 typically separate in energy by a substantial amount, rendering the probability of a back-hop low. So, it is reasonable to expect that a failed trajectory will remain on S0 and include it in the population count as such. This is the approach we used in Fig. 3. A small percentage of trajectories also failed while propagating on S1 (and S2) due to failure of energy conservation. Since it is not possible to establish with confidence the fate of these trajectories, we exclude them from the population count after their failure.

FIG. 3.

Calculated population dynamics of S2, S1, and S0 states of cc-COD at the SA3-CASSCF(4,3)/cc-pVDZ level for a simulation window of 500 fs. The populations from the two conformers were combined using the weighted average of their Boltzmann probabilities.

FIG. 3.

Calculated population dynamics of S2, S1, and S0 states of cc-COD at the SA3-CASSCF(4,3)/cc-pVDZ level for a simulation window of 500 fs. The populations from the two conformers were combined using the weighted average of their Boltzmann probabilities.

Close modal

In the TSH simulations, all trajectories were initiated on the S1 state. For all initial conditions, the two excited states are in close proximity to each other, so there is a rapid decay of the population on S1 and an increase in the population of S2 within the first ∼14 fs because of non-adiabatic transitions to S2. The population transfers back to S1 rather quickly leading to some increase in its population. Eventually, decay to S0 dominates leading to a rough exponential decay of S1 to S0 after about 50 fs. The half-life of the excited state is ∼96 fs, and the populations do not change much after 300 fs—with only 6% left in the excited states by this time.

Figure 4 gives more detailed information about the initial vibrational motion driving the dynamics. Figure 4(a) shows the gradient vectors on the S1 surface at the FC geometry. The gradient involves mainly a symmetric stretching of the two C=C double bonds. Figure 4(b) shows the bond alternation coordinate averaged over all trajectories as a function of time. There is a large amplitude motion for the first 50 fs, confirming that this motion plays a significant role in the dynamics. Similar behavior has also been observed in BD in both resonance Raman spectra75 and calculated dynamics.15 The figure also shows the derivative coupling vector between S1 and S2 in the FC region. This vector involves an asymmetric vibration of the two double bonds and is perpendicular to the gradient. Since the initial geometry has C2 symmetry and the two states have A and B symmetry, it is expected that the coupling vibration should have B symmetry and will be perpendicular to the totally symmetric gradient. This means that the initial motion drives the two states together along the bond alternation coordinate, but an asymmetric stretch vibration is needed to couple them. It is obvious from the ultrafast transition in the dynamics that this occurs very rapidly.

FIG. 4.

(a) (Left) Gradient of the S1 state and (right) derivative coupling vector between S2 and S1 at the FC geometry. The magnitude of the vectors is scaled in order to clearly show the direction. (b) Bond alternation coordinate [defined as (RC1=C8 − RC8C7 + RC7=C6)] averaged over all trajectories as a function of time.

FIG. 4.

(a) (Left) Gradient of the S1 state and (right) derivative coupling vector between S2 and S1 at the FC geometry. The magnitude of the vectors is scaled in order to clearly show the direction. (b) Bond alternation coordinate [defined as (RC1=C8 − RC8C7 + RC7=C6)] averaged over all trajectories as a function of time.

Close modal

We now focus on the non-adiabatic transitions. As expected, CIs play a crucial role in the fast dynamics we observe. A CI between S2 and S1 was found which preserves C2 symmetry. Its energy is at 6.51 eV, about 1 eV below the vertical excitation energy. The structure is given in the supplementary material. Because of the flexibility of the molecule, several local minima on the S1S0 seam have been found. The details of these structures can be found in the supplementary material, and we discuss them further in relation to the hopping geometries below. The lowest energy CI has an energy of 4.68 eV above the S0 minimum, while the others are within 1 eV of that energy.

Since there are two double bonds, the S1/S0 CIs involve either only one double bond (local) twisting or both double bonds twisting at the same time. Most of the CIs located involve a local twisting of one double bond, while only one was found with both double bonds twisted. This is also reflected in the hopping geometries, as shown in Fig. 5. This figure shows a correlation plot between the two dihedral angles around the double bonds. The cross pattern shows that most hops occur via a localized twist where one bond twists while the other double bond has a dihedral value close to zero, rather than a concerted motion with both double bonds twisting at the same time. Using the absolute values of the dihedral angles with a threshold of 30°, we can estimate that 16% of the non-adiabatic transitions occur through a delocalized distortion. In BD, it was found that this distortion was responsible for 23% of the non-adiabatic transitions, higher than the value we see here.16 It was also shown that the participation of this delocalized distortion is directly related to the S1S2 gap. The smaller the gap is the larger the participation. Since the gap here is larger than that in BD, we expect to see a smaller contribution.

FIG. 5.

Correlation between the two dihedral angles defining twisting around the double bonds at all the S1 − >S0 hopping geometries. Twist 1 and twist 2 are defined as the dihedral angles H12—C1=C8—H11 and H10—C7=C6—H9, respectively. The pattern indicates that in most cases, only one double bond twists. Only 16% of the hops have both dihedrals twisted, using a threshold of ±30°.

FIG. 5.

Correlation between the two dihedral angles defining twisting around the double bonds at all the S1 − >S0 hopping geometries. Twist 1 and twist 2 are defined as the dihedral angles H12—C1=C8—H11 and H10—C7=C6—H9, respectively. The pattern indicates that in most cases, only one double bond twists. Only 16% of the hops have both dihedrals twisted, using a threshold of ±30°.

Close modal

Figure 5 establishes that in most of the hops, the twisting is local, but it also indicates that there is a wide spread of twisting values. The variation is highlighted even more in Fig. 6 which shows a bar graph of the values of the twisting angle at the hop geometries, for either double bond. Using the thresholds of <30° as no twisting, 30°–90° as small twisting, and >90° as large twisting, we see that roughly half of the hops twist a particular double bond with 30% having small twist values and 20% large values. Thus, a smaller twist is more likely than a larger one. Both double bonds give the same picture, as one would expect given the symmetry of the system.

FIG. 6.

Bar graphs of the twist of two double bonds, (a) C1=C8 and (b) C7=C6. Twist 1 is defined as the dihedral angle H12—C1=C8—H11, whereas twist 2 is defined as the dihedral angle H10—C7=C6—H9. The two twists are also shown in panel (a). The fraction of S1 − >S0 hops that took place with different twist angles are: 47.5% for no twist 1(<30°), 32.5% for small twist 1 (30–90°), and 20% for large twist 1 (>90°). Similarly for twist 2, the fraction of S1 − >S0 hops are: 46.2% for no twist, 34.4% for small twist, and 19.4% for large twist. Absolute values of the dihedral angles have been used to simplify the plots.

FIG. 6.

Bar graphs of the twist of two double bonds, (a) C1=C8 and (b) C7=C6. Twist 1 is defined as the dihedral angle H12—C1=C8—H11, whereas twist 2 is defined as the dihedral angle H10—C7=C6—H9. The two twists are also shown in panel (a). The fraction of S1 − >S0 hops that took place with different twist angles are: 47.5% for no twist 1(<30°), 32.5% for small twist 1 (30–90°), and 20% for large twist 1 (>90°). Similarly for twist 2, the fraction of S1 − >S0 hops are: 46.2% for no twist, 34.4% for small twist, and 19.4% for large twist. Absolute values of the dihedral angles have been used to simplify the plots.

Close modal

As we have learned from ethylene, twisting alone is not enough to make the ground state degenerate with S1. Pyramidalization of one of the carbons is also involved as a result of charge transfer along the double bond. In this case, either the inner carbon (C8 or C7) or the outer one (C1 or C6) pyramidalizes. This has also been seen in BD with inner defined as Me+ and outer as Me.15,16 Both these types of CIs are found here. CIs with outer C pyramidalization have higher energy compared to the ones with inner. Figure 7 shows bar graphs of the dihedral angles describing pyramidalization around each of the four carbons. The top panels (a) and (b) show pyramidalization of the outer carbons while the bottom panels (c) and (d) show pyramidalization of the inner carbons. In general, there are more hops where the inner carbons are pyramidalized (Me+ in BD). About 40%–50% of the hopping geometries have an inner carbon pyramidalized, while 27%–29% have an outer carbon. This is consistent with the fact that the lowest energy CI has the inner carbon pyramidalized. It is, however, different from what was seen in BD, where many more trajectories showed the outer carbon (Me) pyramidalized. This can be rationalized if we consider the very different constraints of the outer carbons in cc-COD, which limit the flexibility of these carbons.

FIG. 7.

Bar graphs of pyramidalization angles (a) pyr out 1, (b) pyr out 2, (c) pyr in 1, and (d) pyr in 2. Outer pyramidalization angles pyr out 1 (pyramidalization at C1) and pyr out 2 (pyramidalization at C6) have been defined in terms of dihedral angles C2—C1—H12—C8 and C7—C6—H9—C5, respectively. Inner pyramidalization angles pyr in 1 (pyramidalization at C8) and pyr in 2 (pyramidalization at C7) have been defined in terms of dihedral angles C1—C8—H11—C7 and C8—C7—H10—C6, respectively. Panel (a) also shows the position of all the pyramidalization. For the geometries at S1 − >S0 hops for outer C pyramidalization, we found 71% of the geometries non-pyramidalized at C1 [panel (a)] and 73% non-pyramidalized at C6 [panel(b)]. For the geometries at S1 − >S0 hops for inner C pyramidalization, we found 60% non-pyramidalized at C8 [panel (c)] and 49% non-pyramidalized at C7 [panel (d)]. An angular range of 150°–180° was used to classify hops as non-pyramidalized.

FIG. 7.

Bar graphs of pyramidalization angles (a) pyr out 1, (b) pyr out 2, (c) pyr in 1, and (d) pyr in 2. Outer pyramidalization angles pyr out 1 (pyramidalization at C1) and pyr out 2 (pyramidalization at C6) have been defined in terms of dihedral angles C2—C1—H12—C8 and C7—C6—H9—C5, respectively. Inner pyramidalization angles pyr in 1 (pyramidalization at C8) and pyr in 2 (pyramidalization at C7) have been defined in terms of dihedral angles C1—C8—H11—C7 and C8—C7—H10—C6, respectively. Panel (a) also shows the position of all the pyramidalization. For the geometries at S1 − >S0 hops for outer C pyramidalization, we found 71% of the geometries non-pyramidalized at C1 [panel (a)] and 73% non-pyramidalized at C6 [panel(b)]. For the geometries at S1 − >S0 hops for inner C pyramidalization, we found 60% non-pyramidalized at C8 [panel (c)] and 49% non-pyramidalized at C7 [panel (d)]. An angular range of 150°–180° was used to classify hops as non-pyramidalized.

Close modal

The above discussion shows that the non-adiabatic transitions occur over a wide distribution of internal distortions, and the molecule accesses a wide part of the S1S0 seam. Figure 8 shows the energetic spread of the hops as well. The hops occur with relaxation energies that spread over 5 eV, as shown in Fig. 8(a), while they are still near a seam as shown in the values of the S1S0 energy gap shown in Fig. 8(b).

FIG. 8.

(a) Initial excitation energy of trajectories correlated with the relaxation energy at the hopping geometry. (b) Histogram of the energy gap S1S0 at the hopping geometries.

FIG. 8.

(a) Initial excitation energy of trajectories correlated with the relaxation energy at the hopping geometry. (b) Histogram of the energy gap S1S0 at the hopping geometries.

Close modal

Since we are not able to do the dynamics at the XMS-CASPT2 level, we compare LIICs between the FC region and CIs to see if we can learn anything more about the accuracy of the dynamics. Figure 9 shows two LIICs at both the CASSCF and XMS-CASPT2 levels: one between the FC region and a local CI, and one between FC and the non-local CI. LIICs to other CIs have also been constructed (shown in the supplementary material), and they show similar behavior. The curves show that initially the two excited states stay close in energy, and this facilitates transitions back and forth between them, as is seen in the dynamics. After that initial behavior, the S1 state stabilizes rapidly, leading to the CIs with the ground state and rapid radiationless decay. In the local CI, the overall pathway is barrierless, while the pathway to the non-local CI has a barrier on the S1 surface. The LIICs only give an upper bound to any barrier, and this upper bound is 0.38 eV at the CASSCF level and 0.62 eV at the XMS-CASPT2 level. The LIICs provide a simple explanation for the ultrafast decay to the ground state observed in the dynamics, as well as the fact that the local CIs are dominant. Furthermore, the comparison between the two methods in the figure highlights that the overall shape of the curves is qualitatively similar for the two methods.

FIG. 9.

Energies of the first three excited states along a linear interpolation path connecting the FC geometry to [(a) and (b)] a local CI with inner C pyramidalization (CI10 3) at (a) CASSCF(4,3)/cc-pVDZ and (b) XMS-CASPT2/CAS(6,6)/cc-pVDZ level and to [(c) and (d)] the non-local CI (CI10 10) at (c) CASSCF(4,3)/cc-pVDZ and (d) XMS-CASPT2/CAS(6,6)/cc-pVDZ levels. The x-axis is dimensionless in all the plots. Refer to the supplementary material for the structures of the CIs.

FIG. 9.

Energies of the first three excited states along a linear interpolation path connecting the FC geometry to [(a) and (b)] a local CI with inner C pyramidalization (CI10 3) at (a) CASSCF(4,3)/cc-pVDZ and (b) XMS-CASPT2/CAS(6,6)/cc-pVDZ level and to [(c) and (d)] the non-local CI (CI10 10) at (c) CASSCF(4,3)/cc-pVDZ and (d) XMS-CASPT2/CAS(6,6)/cc-pVDZ levels. The x-axis is dimensionless in all the plots. Refer to the supplementary material for the structures of the CIs.

Close modal

Figure 10 shows all the products observed after internal conversion to the ground state. The yields are derived using dynamics from both conformers taking into account their Boltzmann distribution. Our simulation shows that after internal conversion to the ground state, 38% of the population reverts back to the original cc-COD, whereas another 38% undergoes cistrans isomerization to form ct-COD and 1% undergoes a disrotatory Woodword–Hoffmann type electrocyclic ring closure to form BCO. The yield of BCO is the same as that predicted by Fuß et al.,30 but the yield of ct-COD is slightly higher than their prediction. However, we found more products where biradicals form due to the [1,2] H-shift29 and [1,3] C-shift.76 H-migration, i.e., [1,2] H-shift, can take place from one C at one end of a double bond to the other end forming an ethylidene type of biradical. On the other hand, the [1,3] C-shift can break the C—C single bond between the double bonds to form a biradical with a 7-member ring, which can serve as an intermediate to form 3-methylenecycloheptene (3-MCH). All products are shown in Fig. 10. Linear interpolations connecting the CIs to the products are given in the supplementary material.

FIG. 10.

Photochemical products: cc-COD and ct-COD via cistrans isomerization, biradicals via [1,2] H-shift, biradical via [1,3] C-shift and 3-MCH, and BCO via disrotatory Woodward-Hoffmann ring closure. Yields are given below the arrows. The energies of the products [calculated at the DFT/B3LYP/6-31G(d) level] relative to the lowest energy FC geometry (reactant) of conformer 1 are shown in the supplementary material.

FIG. 10.

Photochemical products: cc-COD and ct-COD via cistrans isomerization, biradicals via [1,2] H-shift, biradical via [1,3] C-shift and 3-MCH, and BCO via disrotatory Woodward-Hoffmann ring closure. Yields are given below the arrows. The energies of the products [calculated at the DFT/B3LYP/6-31G(d) level] relative to the lowest energy FC geometry (reactant) of conformer 1 are shown in the supplementary material.

Close modal

Figure 11 presents the internal coordinates which serve to determine the branching into all of the different products discussed in the previous paragraph. These plots are provided only for trajectories starting from conformer 1 for simplicity. The biradical that forms after the [1,3] C-shift can be recognized by plotting the time-evolution of the C7—C8 single bond separating the two double bonds, as this bond breaks to form a C6—C8 bond (or C1—C7 bond) which creates a 7-membered ring. Figure 11(a) shows that all of the trajectories start very close to 1.47 Å, close to the C7—C8 single bond distance at the S0 minimum. However, after S1 − >S0 hops, the C7—C8 bond breaks in a few trajectories and the bond distance increases to around 2.5 Å. The bond distance vibrates around 1.5 Å for the rest of the trajectories for the whole simulation window. The ensemble averages of the two separate branches of trajectories show that the averages separate after ∼50 fs.

FIG. 11.

(a) Formation of a biradical and, subsequently, 3-MCH via [1,3] C-shift is demonstrated using C7—C8 bond breaking. (b) Formation of a biradical following the [1,2] H-shift is shown for C1—H12 bond breaking where H12 migrates to C8 from C1. (c) Formation of ct-COD via cistrans isomerization is shown using H10—C7=C6—H9 dihedral angle.

FIG. 11.

(a) Formation of a biradical and, subsequently, 3-MCH via [1,3] C-shift is demonstrated using C7—C8 bond breaking. (b) Formation of a biradical following the [1,2] H-shift is shown for C1—H12 bond breaking where H12 migrates to C8 from C1. (c) Formation of ct-COD via cistrans isomerization is shown using H10—C7=C6—H9 dihedral angle.

Close modal

Figure 11(b) shows the time-evolution of the [1,2] H-shift for one of the C—H bonds (C1—H12) corresponding to a double bond. Here, the H12 migrates from C1 to C8, forming a biradical. Similarly, H11 can also migrate to C1 from C8 (not shown here). The same phenomenon can happen at the other double bond too. Very interestingly, the formation of these products starts at later times, after approximately 100 fs.

As the yield was very low for BCO, no internal coordinates were plotted in this case. Once we separate the trajectories which go through the [1,2] H-shift, [1,3] C-shift, and electrocyclic ring closure, we are left with the trajectories that either revert back to cc-COD or form ct-COD via cistrans isomerization. Figure 11(c) shows the time-evolution of the cis or trans position of H’s with respect to each other at the C7=C6 double bond. If they are cis to each other, then the dihedral angle H10—C7=C6—H9 will be close to 0°, whereas when they are trans to each other, the same angle will be close to 180° or −180°. It can be seen that all the trajectories start initially close to 0°, but they start to move toward higher values of the angle ∼ between 20 fs and 150 fs. However, after S1 − >S0 hops, this angle decreases to vibrate around 0° for a set of trajectories (cis), whilst for the other set of trajectories (trans), the same angle increases to vibrate around 180° or −180°. The same is true for the dihedral angle H12—C1=C8—H11 (which is not shown here) too. By adding the number of trajectories that ended up in the trans branch for both dihedral angles, we can calculate the number of trajectories that formed ct-COD after the cistrans isomerization, while the rest constitutes the trajectories that reverted back to cc-COD.

There were also some trajectories which failed right after the S1 − >S0 hop (failed while still in the S1/S0 seam), and hence, no decision could be taken regarding the product for these trajectories. Similarly, no decision could be taken about the trajectories that failed in S1 and S2. Hence, they are excluded from the calculation of yields. This could be a reason why our calculated yield for ct-COD is higher than that predicted by Fuß et al.30 

In this study, we have performed trajectory surface hopping dynamics for cc-COD to investigate the photochemical pathways involved after being excited to the S1 state, which we found to be the bright state in the FC region at higher level of theory. Ultrafast decay to the ground state has been found, facilitated by CIs involving distortions around the double bonds. We have examined whether the non-adiabatic transitions to the ground state involve twisting along one double bond or both of them. The dynamics show that the main distortions are twisting of one double bond at a time (local) followed by pyramidalization of a carbon belonging to that double bond (similar to ethylene). Among the local twisting distortions, inner carbon pyramidalization was found to be more prevalent. Only about 16% of transitions occur through a non-local or delocalized twisting of both double bonds. Several products formed on the ground state have also been observed.

Our calculations and analysis demonstrate that the internal conversion dynamics for photo-excited cc-COD are highly non-local, with trajectories exploring many different geometries at which non-adiabatic coupling drives hopping between electronic states. We looked for and were not able to find groups of trajectories that follow similar pathways. We suspect that this is related to the flexibility of the molecule and the relatively flat potential energy surface along many low frequency modes, where different initial conditions can lead to trajectories exploring quite different geometries and finding their way down to the ground state via very different points on an extended CI seam, rather than through a minimum energy pathway in which the non-adiabatic coupling and hopping is concentrated around a particular CI.

The results of this work provide insight into both the effect of conjugation on dynamics and the effect of additional low frequency degrees of freedom. The same chromophoric unit (two conjugated double bonds in this case) does not always result in the same excited state dynamics, i.e., CHD and cc-COD behave differently, while cc-COD and BD are much more similar. In addition, the low frequency degrees of freedom can also have an important effect on the dynamics by providing flexibility to the system. These observations can be useful as we try to expand our understanding of excited state dynamics and find relations between structure and dynamics.

In a separate manuscript, we will analyze time-resolved photoelectron spectroscopy measurements in light of these calculations.

See the supplementary material for additional calculations as outlined in the text.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors gratefully acknowledge the Department of Energy (DOE, Award No. DE-FG02-08ER15983 for P.C. and S.M., and DOE, Award No. DE-FG02-08ER15984 for Y.L. and T.W.) for funding. Most of the computational work was performed using the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. ACI-1548562. P.C. thanks Asst. Professor Jae Woo Park and XSEDE personnel for technical advice/help with running jobs and installation of the Bagel package.

1.
S.
Gozem
,
H. L.
Luk
,
I.
Schapiro
, and
M.
Olivucci
, “
Theory and simulation of the ultrafast double-bond isomerization of biological chromophores
,”
Chem. Rev.
117
,
13502
13565
(
2017
).
2.
C.
Dugave
and
L.
Demange
, “
Cis-trans isomerization of organic molecules and biomolecules: Implications and applications
,”
Chem. Rev.
103
,
2475
2532
(
2003
).
3.
M. A.
van der Horst
and
K. J.
Hellingwerf
, “
Photoreceptor proteins, “star actors of modern times”: A review of the functional dynamics in the structure of representative members of six different photoreceptor families
,”
Acc. Chem. Res.
37
,
13
20
(
2004
).
4.
M.
Barbatti
,
J.
Paier
, and
H.
Lischka
, “
Photochemistry of ethylene: A multireference configuration interaction investigation of the excited-state energy surfaces
,”
J. Chem. Phys.
121
,
11614
11624
(
2004
).
5.
B. R.
Brooks
and
H. F.
Schaefer
, “
Sudden polarization: Pyramidalization of twisted ethylene
,”
J. Am. Chem. Soc.
101
,
307
311
(
1979
).
6.
H.
Tao
,
B. G.
Levine
, and
T. J.
Martínez
, “
Ab initio multiple spawning dynamics using multi-state second-order perturbation theory
,”
J. Phys. Chem. A
113
,
13656
13662
(
2009
).
7.
B. G.
Levine
and
T. J.
Martínez
, “
Isomerization through conical intersections
,”
Annu. Rev. Phys. Chem.
58
,
613
634
(
2007
).
8.
T.
Mori
,
W. J.
Glover
,
M. S.
Schuurman
, and
T. J.
Martinez
, “
Role of Rydberg states in the photochemical dynamics of ethylene
,”
J. Phys. Chem. A
116
,
2808
2818
(
2012
).
9.
S. P.
Neville
,
M.
Chergui
,
A.
Stolow
, and
M. S.
Schuurman
, “
Ultrafast x-ray spectroscopy of conical intersections
,”
Phys. Rev. Lett.
120
,
243001
(
2018
).
10.
P.
Celani
,
F.
Bernardi
,
M.
Olivucci
, and
M. A.
Robb
, “
Excited-state reaction pathways for s-cis buta-1,3-diene
,”
J. Chem. Phys.
102
,
5733
5742
(
1995
).
11.
M.
Olivucci
,
I. N.
Ragazos
,
F.
Bernardi
, and
M. A.
Robb
, “
A conical intersection mechanism for the photochemistry of butadiene. A MC-SCF study
,”
J. Am. Chem. Soc.
115
,
3710
3721
(
1993
).
12.
M.
Ito
and
I.
Ohmine
, “
Nonadiabatic transition and energy relaxation dynamics in the photoisomerization of s-trans butadiene
,”
J. Chem. Phys.
106
,
3159
3173
(
1997
).
13.
C.
Nonnenberg
,
S.
Grimm
, and
I.
Frank
, “
Restricted open-shell Kohn–Sham theory for π − π* transitions. II. Simulation of photochemical reactions
,”
J. Chem. Phys.
119
,
11585
11590
(
2003
).
14.
Y.
Dou
,
B. R.
Torralva
, and
R. E.
Allen
, “
Detailed mechanism for trans-cis photoisomerization of butadiene following a femtosecond-scale laser pulse
,”
J. Phys. Chem. A
107
,
8817
8824
(
2003
).
15.
B. G.
Levine
and
T. J.
Martínez
, “
Ab initio multiple spawning dynamics of excited butadiene: Role of charge transfer
,”
J. Phys. Chem. A
113
,
12815
12824
(
2009
).
16.
W. J.
Glover
,
T.
Mori
,
M. S.
Schuurman
,
A. E.
Boguslavskiy
,
O.
Schalk
,
A.
Stolow
, and
T. J.
Martínez
, “
Excited state non-adiabatic dynamics of the smallest polyene, trans 1,3-butadiene. II. Ab initio multiple spawning simulations
,”
J. Chem. Phys.
148
,
164303
(
2018
).
17.
A. E.
Boguslavskiy
,
O.
Schalk
,
N.
Gador
,
W. J.
Glover
,
T.
Mori
,
T.
Schultz
,
M. S.
Schuurman
,
T. J.
Martínez
, and
A.
Stolow
, “
Excited state non-adiabatic dynamics of the smallest polyene, trans 1,3-butadiene. I. Time-resolved photoelectron-photoion coincidence spectroscopy
,”
J. Chem. Phys.
148
,
164302
(
2018
).
18.
P.
Celani
,
S.
Ottani
,
M.
Olivucci
,
F.
Bernardi
, and
M. A.
Robb
, “
What happens during the picosecond lifetime of 2A1 cyclohexa-1,3-diene? A CAS-SCF study of the cyclohexadiene/hexatriene photochemical interconversion
,”
J. Am. Chem. Soc.
116
,
10141
10151
(
1994
).
19.
P.
Celani
,
F.
Bernardi
,
M. A.
Robb
, and
M.
Olivucci
, “
Do photochemical ring-openings occur in the spectroscopic state? 1B2 pathways for the cyclohexadiene/hexatriene photochemical interconversion
,”
J. Phys. Chem.
100
,
19364
19366
(
1996
).
20.
M.
Garavelli
,
P.
Celani
,
M.
Fato
,
M. J.
Bearpark
,
B. R.
Smith
,
M.
Olivucci
, and
M. A.
Robb
, “
Relaxation paths from a conical intersection: The mechanism of product formation in the cyclohexadiene/hexatriene photochemical interconversion
,”
J. Phys. Chem. A
101
,
2023
2032
(
1997
).
21.
J.
Kim
,
H.
Tao
,
J. L.
White
,
V. S.
Petrović
,
T. J.
Martinez
, and
P. H.
Bucksbaum
, “
Control of 1,3-cyclohexadiene photoisomerization using light-induced conical intersections
,”
J. Phys. Chem. A
116
,
2758
2763
(
2012
).
22.
K.
Kosma
,
S. A.
Trushin
,
W.
Fuß
, and
W. E.
Schmid
, “
Cyclohexadiene ring opening observed with 13 fs resolution: Coherent oscillations confirm the reaction path
,”
Phys. Chem. Chem. Phys.
11
,
172
181
(
2009
).
23.
S.
Deb
and
P. M.
Weber
, “
The ultrafast pathway of photon-induced electrocyclic ring-opening reactions: The case of 1,3-cyclohexadiene
,”
Annu. Rev. Phys. Chem.
62
,
19
39
(
2011
).
24.
S. L.
Horton
,
Y.
Liu
,
P.
Chakraborty
,
S.
Matsika
, and
T.
Weinacht
, “
Vibrationally assisted below-threshold ionization
,”
Phys. Rev. A
95
,
063413
(
2017
).
25.
K.
Kaneshima
,
Y.
Ninota
, and
T.
Sekikawa
, “
Time-resolved high-harmonic spectroscopy of ultrafast photoisomerization dynamics
,”
Opt. Express
26
,
31039
31054
(
2018
).
26.
I.
Polyak
,
L.
Hutton
,
R.
Crespo-Otero
,
M.
Barbatti
, and
P. J.
Knowles
, “
Ultrafast photoinduced dynamics of 1,3-cyclohexadiene using XMS-CASPT2 surface hopping
,”
J. Chem. Theory Comput.
15
,
3929
3940
(
2019
).
27.
O.
Schalk
,
T.
Geng
,
T.
Thompson
,
N.
Baluyot
,
R. D.
Thomas
,
E.
Tapavicza
, and
T.
Hansson
, “
Cyclohexadiene revisited: A time-resolved photoelectron spectroscopy and ab initio study
,”
J. Phys. Chem. A
120
,
2320
2329
(
2016
).
28.
M.
Tudorovskaya
,
R. S.
Minns
, and
A.
Kirrander
, “
Effects of probe energy and competing pathways on time-resolved photoelectron spectroscopy: The ring-opening of 1,3-cyclohexadiene
,”
Phys. Chem. Chem. Phys.
20
,
17714
17726
(
2018
).
29.
O.
Schalk
,
A. E.
Boguslavskiy
,
A.
Stolow
, and
M. S.
Schuurman
, “
Through-bond interactions and the localization of excited-state dynamics
,”
J. Am. Chem. Soc.
133
,
16451
16458
(
2011
).
30.
W.
Fuß
,
S.
Panja
,
W. E.
Schmid
, and
S. A.
Trushin
, “
Competing ultrafast cis-trans isomerization and ring closure of cyclohepta-1,3-diene and cyclo-octa-1,3-diene
,”
Mol. Phys.
104
,
1133
1143
(
2006
).
31.
K.
Komori-Orisaku
,
Y.
Hirose
, and
I.
Iwakura
, “
Pulsed Nd:YAG laser-induced photoreaction of cis,cis-1,3-cyclooctadiene at 266 nm: Selective cyclization to cis-bicyclo[4.2.0]oct-7-ene
,”
Photochem. Photobiol. Sci.
16
,
146
150
(
2017
).
32.
F.
Bernardi
,
M.
Olivucci
,
I. N.
Ragazos
, and
M. A.
Robb
, “
Origin of the nonstereospecificity in the ring opening of alkyl-substituted cyclobutenes
,”
J. Am. Chem. Soc.
114
,
2752
(
1992
).
33.
Spartan’16, Wavefunction, Inc., Irvine, CA.
34.
W.
Kohn
,
A. D.
Becke
, and
R. G.
Parr
, “
Density functional theory of electronic structure
,”
J. Phys. Chem.
100
,
12974
12980
(
1996
).
35.
T.
Ziegler
, “
Approximate density functional theory as a practical tool in molecular energetics and dynamics
,”
Chem. Rev.
91
,
651
667
(
1991
).
36.
A. D.
Becke
, “
Density-functional exchange-energy approximation with correct asymptotic behavior
,”
Phys. Rev. A
38
,
3098
3100
(
1988
).
37.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the colle-salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
789
(
1988
).
38.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
, “
Accurate spin-dependent electron liquid correlation energies for local spin density calculations: A critical analysis
,”
Can. J. Phys.
58
,
1200
1211
(
1980
).
39.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
,
5648
5652
(
1993
).
40.
R.
Ditchfield
,
W. J.
Hehre
, and
J. A.
Pople
, “
Self-consistent molecular-orbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules
,”
J. Chem. Phys.
54
,
724
728
(
1971
).
41.
K.
Raghavachari
,
G. W.
Trucks
,
J. A.
Pople
, and
M.
Head-Gordon
, “
A fifth-order perturbation comparison of electron correlation theories
,”
Chem. Phys. Lett.
157
,
479
483
(
1989
).
42.
R. J.
Bartlett
,
J. D.
Watts
,
S. A.
Kucharski
, and
J.
Noga
, “
Non-iterative fifth-order triple and quadruple excitation energy corrections in correlated methods
,”
Chem. Phys. Lett.
165
,
513
522
(
1990
).
43.
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
44.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
B.
Mennucci
,
G. A.
Petersson
,
H.
Nakatsuji
,
M.
Caricato
,
X.
Li
,
H. P.
Hratchian
,
A. F.
Izmaylov
,
J.
Bloino
,
G.
Zheng
,
J. L.
Sonnenberg
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
J. A.
Montgomery
,
J. E.
Peralta
,
F.
Ogliaro
,
M.
Bearpark
,
J. J.
Heyd
,
E.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
N.
Rega
,
J. M.
Millam
,
M.
Klene
,
J. E.
Knox
,
J. B.
Cross
,
V.
Bakken
,
C.
Adamo
,
J.
Jaramillo
,
R.
Gomperts
,
R. E.
Stratmann
,
O.
Yazyev
,
A. J.
Austin
,
R.
Cammi
,
C.
Pomelli
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
V. G.
Zakrzewski
,
G. A.
Voth
,
P.
Salvador
,
J. J.
Dannenberg
,
S.
Dapprich
,
A. D.
Daniels
,
O.
Farkas
,
J. B.
Foresman
,
J. V.
Ortiz
,
J.
Cioslowski
, and
D. J.
Fox
, Gaussian 09, Revision B.01,
2009
.
45.
B. O.
Roos
,
P. R.
Taylor
, and
P. E. M.
Sigbahn
, “
A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach
,”
Chem. Phys.
48
,
157
173
(
1980
).
46.
J.
Finley
,
P.-Å.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
, “
The multi-state CASPT2 method
,”
Chem. Phys. Lett.
288
,
299
306
(
1998
).
47.
A. A.
Granovsky
, “
Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory
,”
J. Chem. Phys.
134
,
214113
(
2011
).
48.
T.
Shiozaki
,
W.
Győrffy
,
P.
Celani
, and
H.-J.
Werner
, “
Communication: Extended multi-state complete active space second-order perturbation theory: Energy and nuclear gradients
,”
J. Chem. Phys.
135
,
081106
(
2011
).
49.
B. O.
Roos
, “
The multiconfigurational (MC) self-consistent field (SCF) theory
,” in
Lecture Notes in Quantum Chemistry: European Summer School in Quantum Chemistry
, edited by
B. O.
Roos
(
Springer Berlin Heidelberg
,
Berlin, Heidelberg
,
1992
), pp.
177
254
.
50.
J. W.
Park
, “
Single-state single-reference and multistate multireference zeroth-order Hamiltonians in MS-CASPT2 and conical intersections
,”
J. Chem. Theory Comput.
15
,
3960
3973
(
2019
).
51.
H.
Lischka
,
R.
Shepard
,
I.
Shavitt
,
R. M.
Pitzer
,
M.
Dallos
,
T.
Müller
,
P. G.
Szalay
,
F. B.
Brown
,
R.
Ahlrichs
,
H. J.
Böhm
,
A.
Chang
,
D. C.
Comeau
,
R.
Gdanitz
,
H.
Dachsel
,
C.
Ehrhardt
,
M.
Ernzerhof
,
P.
Höchtl
,
S.
Irle
,
G.
Kedziora
,
T.
Kovar
,
V.
Parasuk
,
M. J. M.
Pepper
,
P.
Scharf
,
H.
Schiffer
,
M.
Schindler
,
M.
Schüler
,
M.
Seth
,
E. A.
Stahlberg
,
J.-G.
Zhao
,
S.
Yabushita
,
Z.
Zhang
,
M.
Barbatti
,
S.
Matsika
,
M.
Schuurmann
,
D. R.
Yarkony
,
S. R.
Brozell
,
E. V.
Beck
,
J.-P.
Blaudeau
,
M.
Ruckenbauer
,
B.
Sellner
,
F.
Plasser
, and
J. J.
Szymczak
, COLUMBUS, an ab initio electronic structure program, release 7.0, 2012.
52.
H.
Lischka
,
T.
Müller
,
P. G.
Szalay
,
I.
Shavitt
,
R. M.
Pitzer
, and
R.
Shepard
, “
Columbus-a program system for advanced multireference theory calculations
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
191
199
(
2011
).
53.
H.
Lischka
,
R.
Shepard
,
R. M.
Pitzer
,
I.
Shavitt
,
M.
Dallos
,
T.
Müller
,
P. G.
Szalay
,
M.
Seth
,
G. S.
Kedziora
,
S.
Yabushita
, and
Z.
Zhang
, “
High-level multireference methods in the quantum-chemistry program system COLUMBUS: Analytic MR-CISD and MR-AQCC gradients and MR-AQCC-LRT for excited states, GUGA spin-orbit CI and parallel CI density
,”
Phys. Chem. Chem. Phys.
3
,
664
673
(
2001
).
54.
BAGEL, Brilliantly Advanced General Electronic-structure Library, https://www.nubakery.org under the GNU General Public License.
55.
T.
Shiozaki
, “
BAGEL: Brilliantly advanced general electronic-structure library
,”
WIREs Comput. Mol. Sci.
8
,
e1331
(
2018
).
56.
M.
Barbatti
,
M.
Ruckenbauer
,
F.
Plasser
,
J.
Pittner
,
G.
Granucci
,
M.
Persico
, and
H.
Lischka
, “
Newton-X: A surface-hopping program for nonadiabatic molecular dynamics
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
26
33
(
2014
).
57.
M.
Barbatti
,
G.
Granucci
,
M.
Ruckenbauer
,
F.
Plasser
,
R.
Crespo-Otero
,
J.
Pittner
,
M.
Persico
, and
H.
Lischka
, NEWTON-X: A package for Newtonian dynamics close to the crossing seam, version 2, https://www.newtonx.org.
58.
R.
Shepard
, “
Geometrical energy derivative evaluation with MRCI wave functions
,”
Int. J. Quantum Chem.
31
,
33
44
(
1987
).
59.
R.
Shepard
,
H.
Lischka
,
P. G.
Szalay
,
T.
Kovar
, and
M.
Ernzerhof
, “
A general multireference configuration interaction gradient program
,”
J. Chem. Phys.
96
,
2085
2098
(
1992
).
60.
H.
Lischka
,
M.
Dallos
, and
R.
Shepard
, “
Analytic MRCI gradient for excited states: Formalism and application to the n- valence- and n-(3s,3p) Rydberg states of formaldehyde
,”
Mol. Phys.
100
,
1647
1658
(
2002
).
61.
H.
Lischka
,
M.
Dallos
,
P. G.
Szalay
,
D. R.
Yarkony
, and
R.
Shepard
, “
Analytic evaluation of nonadiabatic coupling terms at the MR-CI level. I. Formalism
,”
J. Chem. Phys.
120
,
7322
7329
(
2004
).
62.
M.
Dallos
,
H.
Lischka
,
R.
Shepard
,
D. R.
Yarkony
, and
P. G.
Szalay
, “
Analytic evaluation of nonadiabatic coupling terms at the MR-CI level. II. Minima on the crossing seam: Formaldehyde and the photodimerization of ethylene
,”
J. Chem. Phys.
120
,
7330
7339
(
2004
).
63.
J. C.
Tully
, “
Molecular dynamics with electronic transitions
,”
J. Chem. Phys.
93
,
1061
1071
(
1990
).
64.
J.
Towns
,
T.
Cockerill
,
M.
Dahan
,
I.
Foster
,
K.
Gaither
,
A.
Grimshaw
,
V.
Hazlewood
,
S.
Lathrop
,
D.
Lifka
,
G. D.
Peterson
,
R.
Roskies
,
J. R.
Scott
, and
N.
Wilkins-Diehr
, “
XSEDE: Accelerating scientific Discovery
,”
Comput. Sci. Eng.
16
,
62
74
(
2014
).
65.
C.
Zhu
,
S.
Nangia
,
A. W.
Jasper
, and
D. G.
Truhlar
, “
Coherent switching with decay of mixing: An improved treatment of electronic coherence for non-Born-Oppenheimer trajectories
,”
J. Chem. Phys.
121
,
7658
7670
(
2004
).
66.
G.
Granucci
and
M.
Persico
, “
Critical appraisal of the fewest switches algorithm for surface hopping
,”
J. Chem. Phys.
126
,
134114
(
2007
).
67.
A.
Jain
,
E.
Alguire
, and
J. E.
Subotnik
, “
An efficient, augmented surface hopping algorithm that includes decoherence for use in large-scale simulations
,”
J. Chem. Theory Comput.
12
,
5256
5268
(
2016
).
68.
J. E.
Subotnik
,
A.
Jain
,
B.
Landry
,
A.
Petit
,
W.
Ouyang
, and
N.
Bellonzi
, “
Understanding the surface hopping view of electronic transitions and decoherence
,”
Annu. Rev. Phys. Chem.
67
,
387
417
(
2016
).
69.
M. F.
Herman
, “
Nonadiabatic semiclassical scattering. I. Analysis of generalized surface hopping procedures
,”
J. Chem. Phys.
81
,
754
763
(
1984
).
70.
D. F.
Coker
and
L.
Xiao
, “
Methods for molecular dynamics with nonadiabatic transitions
,”
J. Chem. Phys.
102
,
496
510
(
1995
).
71.
G.
Käb
, “
Fewest switches adiabatic surface hopping as applied to vibrational energy relaxation
,”
J. Phys. Chem. A
110
,
3197
3215
(
2006
).
72.
A.
Carof
,
S.
Giannini
, and
J.
Blumberger
, “
Detailed balance, internal consistency, and energy conservation in fragment orbital-based surface hopping
,”
J. Chem. Phys.
147
,
214113
(
2017
).
73.
F.
Plasser
,
S.
Mai
,
M.
Fumanal
,
E.
Gindensperger
,
C.
Daniel
, and
L.
González
, “
Strong influence of decoherence corrections and momentum rescaling in surface hopping dynamics of transition metal complexes
,”
J. Chem. Theory Comput.
15
,
5031
5045
(
2019
).
74.
H.
Tao
, “
First principles molecular dynamics and control of photochemical reactions
,” Ph.D. thesis,
Stanford University
,
2011
.
75.
R. R.
Chadwick
,
D. P.
Gerrity
, and
B. S.
Hudson
, “
Resonance Raman spectroscopy of butadiene: Demonstration of a 21Ag state below the 11Bu V state
,”
Chem. Phys. Lett.
115
,
24
(
1985
).
76.
J. E.
Baldwin
and
P. A.
Leber
, “
Molecular rearrangements through thermal [1,3] carbon shifts
,”
Org. Biomol. Chem.
6
,
36
47
(
2008
).

Supplementary Material