We have performed trajectory surface hopping dynamics for *cis*,*cis*-1,3-cyclooctadiene to investigate the photochemical pathways involved after being excited to the *S*_{1} 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.

## I. INTRODUCTION

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 *cis*–*trans* 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.^{29} *Cis*,*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.

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 *cis*–*trans* 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 (*S*_{0} 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.

## II. COMPUTATIONAL METHODS

### A. Electronic structure calculations

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 B3LYP^{36–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-pVDZ^{43} 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) theory^{46–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.0^{51–53} and Bagel^{54,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 *S*_{1}/*S*_{0} 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 *S*_{1} and *S*_{2}, as well as minimizations on *S*_{1}, 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 *S*_{1} 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.

### B. Trajectory surface hopping dynamics

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-X^{56,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 surfaces^{58–62} (PESs) calculated on-the-fly using the Columbus 7.0 package. The trajectories were propagated starting from the *S*_{1} state which is the brighter of the two states (*S*_{1} and *S*_{2}) 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 *S*_{2}, *S*_{1}, and *S*_{0} 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 Persico^{66} has been used. In this approach, a decoherence time is defined as

where *E*_{L} and *E*_{K} are the energies of the current state, L, and any other (inactive) state, K; *E*_{kin} 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 (*C*_{K}) at every time step, Δt. This leads to the new inactive state coefficient $CK\u2032$, which is then used to rescale the current state electronic coefficient (*C*_{L}), as shown in Eq. (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.

## III. RESULTS AND DISCUSSION

### A. Conformers

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 *sp*^{3} 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 *C*_{2} 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.

### B. Vertical excitation energies and oscillator strengths

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 *C*_{2} 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 [*S*_{1}(B) and *S*_{2}(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 *S*_{1} is the bright state, while *S*_{2} 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 *S*_{1} and *S*_{2} 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}

Method . | E(S_{1}) (f)
. | E(S_{2}) (f)
. | g(S_{1})
. | g(S_{2})
. | NAC(S_{1}–S_{2})
. |
---|---|---|---|---|---|

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^{−2}) | 8.39 (0.56) | 0.37 | 0.19 | 0.47 |

SA3-CASSCF(6,6) | 6.89 (3.7 × 10^{−3}) | 8.09 (0.51) | 0.35 | 0.16 | 0.50 |

SA3-CASSCF(8,8) | 7.13 (6.8 × 10^{−3}) | 8.00 (0.46) | |||

SA3-CASSCF(10,10) | 7.11 (8.4 × 10^{−3}) | 7.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^{−3}) | 0.22 | 0.30 | 0.79 |

XMS-CASPT2/CAS(6,6) | 5.91 (0.27) | 6.72 (2.8 × 10^{−3}) | 0.22 | 0.29 | 0.79 |

XMS-CASPT2/CAS(8,8) | 6.02 (0.28) | 6.66 (5.5 × 10^{−3}) | |||

Exp | 5.44 |

Method . | E(S_{1}) (f)
. | E(S_{2}) (f)
. | g(S_{1})
. | g(S_{2})
. | NAC(S_{1}–S_{2})
. |
---|---|---|---|---|---|

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^{−2}) | 8.39 (0.56) | 0.37 | 0.19 | 0.47 |

SA3-CASSCF(6,6) | 6.89 (3.7 × 10^{−3}) | 8.09 (0.51) | 0.35 | 0.16 | 0.50 |

SA3-CASSCF(8,8) | 7.13 (6.8 × 10^{−3}) | 8.00 (0.46) | |||

SA3-CASSCF(10,10) | 7.11 (8.4 × 10^{−3}) | 7.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^{−3}) | 0.22 | 0.30 | 0.79 |

XMS-CASPT2/CAS(6,6) | 5.91 (0.27) | 6.72 (2.8 × 10^{−3}) | 0.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 *S*_{1} and *S*_{2} surfaces as well as the non-adiabatic coupling vectors. The components of these vectors are given by *g*_{A}(*S*_{I}) = ∇_{A}⟨Ψ_{I}|*H*|Ψ_{I}⟩ for the gradients and *NAC*_{A}(*S*_{I} − *S*_{J}) = ⟨Ψ_{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 *S*_{2} surface compared to XMS-CASPT2 which predicts more parallel *S*_{1} and *S*_{2} surfaces. The CASSCF method also gives the highest magnitude for the derivative coupling, indicating a very fast *S*_{1} − *S*_{2} 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}

### C. Absorption spectrum

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 *S*_{1} 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.

### D. Dynamics

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 *S*_{0} 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 *S*_{0}. 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 *S*_{0} state, the potential energy surfaces *S*_{1} and *S*_{0} 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 *S*_{0} 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 *S*_{1} (and *S*_{2}) 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.

In the TSH simulations, all trajectories were initiated on the *S*_{1} 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 *S*_{1} and an increase in the population of *S*_{2} within the first ∼14 fs because of non-adiabatic transitions to *S*_{2}. The population transfers back to *S*_{1} rather quickly leading to some increase in its population. Eventually, decay to *S*_{0} dominates leading to a rough exponential decay of *S*_{1} to *S*_{0} 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 *S*_{1} 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 spectra^{75} and calculated dynamics.^{15} The figure also shows the derivative coupling vector between *S*_{1} and *S*_{2} 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 *C*_{2} 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.

### E. Conical intersections and non-adiabatic transitions

We now focus on the non-adiabatic transitions. As expected, CIs play a crucial role in the fast dynamics we observe. A CI between *S*_{2} and *S*_{1} was found which preserves *C*_{2} 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 *S*_{1} − *S*_{0} 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 *S*_{0} minimum, while the others are within 1 eV of that energy.

Since there are two double bonds, the *S*_{1}/*S*_{0} 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 *S*_{1} − *S*_{2} 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.

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.

As we have learned from ethylene, twisting alone is not enough to make the ground state degenerate with *S*_{1}. 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.

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 *S*_{1} − *S*_{0} 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 *S*_{1} − *S*_{0} energy gap shown in Fig. 8(b).

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 *S*_{1} 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 *S*_{1} 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.

### F. Product distributions

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 *cis*–*trans* 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-shift^{29} 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.

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 *S*_{0} minimum. However, after *S*_{1} − >*S*_{0} 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.

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 *cis*–*trans* 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 *S*_{1} − >*S*_{0} 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 *cis*–*trans* isomerization, while the rest constitutes the trajectories that reverted back to cc-COD.

There were also some trajectories which failed right after the *S*_{1} − >*S*_{0} hop (failed while still in the *S*_{1}/*S*_{0} 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 *S*_{1} and *S*_{2}. 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}

## IV. CONCLUSION

In this study, we have performed trajectory surface hopping dynamics for cc-COD to investigate the photochemical pathways involved after being excited to the *S*_{1} 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.

## SUPPLEMENTARY MATERIAL

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

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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.

## REFERENCES

*ab initio*electronic structure program, release 7.0, 2012.

^{1}A

_{g}state below the 1

^{1}B

_{u}V state