A MASH simulation of the photoexcited dynamics of cyclobutanone Open Access

In recent years, there has been a significant increase in experimental capabilities, making it possible to follow ultrafast photochemical processes in real time. Nonetheless, obtaining a clear mechanistic interpretation of molecular quantum dynamics often requires theoretical calculations to be performed in tandem. Computer simulations of photochemistry are complicated by the breakdown of the Born–Oppenheimer approximation. We thus require methods that can describe nonadiabatic dynamics involving transitions between electronic states.

In any theory, it is desirable to minimize the complexity of the description as much as possible, in order to obtain a simple intuitive picture of the key processes at play. Nonadiabatic approaches based on independent semiclassical trajectories achieve just that, of which Tully’s fewest-switches surface hopping (FSSH)¹ is the most commonly used. In addition, the favorable computational scaling of independent trajectories with system size means that a high-level of electronic-structure theory can be employed, which is crucial for making quantitative predictions of experiments. Unfortunately, FSSH has a number of well-known problems due to inconsistency and overcoherence.² 

The mapping approach to surface hopping (MASH)³ is a recently proposed alternative to the FSSH algorithm. It has a rigorous derivation based on mapping approaches^4–6 but instead of using a mean-field force, it hops between adiabatic states, similarly to FSSH. The key difference between MASH and FSSH is that MASH’s dynamics are deterministic rather than stochastic. This has an important benefit, ensuring consistency at all times between the electronic variables and the active surface. In addition, the MASH derivation uniquely determines how the momentum should be rescaled at attempted hops. One should rescale along the direction of the nonadiabatic coupling vector and reflect in the case of all forbidden hops. Although this prescription is identical to Tully’s original concept,⁷ many alternative suggestions have been made,^2,8,9 and in practice, approximations are often taken.¹⁰ Previous work has shown that MASH is often more accurate than FSSH for a range of model systems,³ and there is reason to believe that it can even be more accurate than ab initio multiple spawning (AIMS)^11,12 for photochemical problems.¹³ MASH was shown to correctly recover Marcus theory rates,¹⁴ where FSSH is known to require complicated decoherence corrections.^15,16 Additionally, unlike other mapping approaches,^4,6,17–19 MASH rigorously captures the detailed balance necessary to thermalize to the correct equilibrium distribution.²⁰ 

In this work, we simulate the photochemistry of cyclobutanone, in response to a community prediction challenge initiated by the Journal of Chemical Physics (JCP). In the experiment that we seek to predict, cyclobutanone is initially photoexcited using a 200 nm pump pulse, which is assumed to excite the molecule from the S₀ to the S₂ adiabatic electronic state. The resulting nonequilibrium dynamics are then measured using time-resolved electron diffraction. Prior to the community challenge, all previous theoretical studies of cyclobutanone have, instead, considered the dynamics starting in S₁,^21–24 meaning that the present study covers new ground. Simulating the photoexcited dynamics initialized in S₂ also poses additional theoretical challenges. S₂ is known to be a Rydberg state^25,26 and, therefore, requires a set of diffuse electronic basis functions to correctly describe the dynamics.

The original MASH formulation was limited to two-state problems and is, therefore, not directly applicable to the present study. Although a multi-state version of MASH has recently been proposed by Runeson and Manolopoulos,^27,28 we note that this does not reduce to the original two-state version and thus loses some of the key benefits of the MASH approach. This is particularly significant for photochemical applications, where the dynamics is expected to largely be a succession of effective two-state nonadiabatic transitions. We, therefore, developed a new multi-state generalization of MASH, which does recover the original version in the case that two states are uncoupled from all others. This is the method employed for the present study. This new approach will be described in detail alongside application to a series of benchmarks in a forthcoming publication.²⁹ 

The present work describes the first implementation of our new multi-state MASH method using ab initio electronic-structure methods. Due to the time constraints imposed by JCP’s challenge, we were not able to implement the most powerful version of the algorithm, and thus, we limit ourselves to studying the internal conversion between singlet states only, and employing an initial distribution obtained from a vertical transition according to the Franck–Condon principle. We note, however, that the MASH formalism can be rigorously extended to treat intersystem crossing to triplet states and to describe the excitation pulse explicitly.³⁰ Future work will test the impact of this more complete description of the nonadiabatic process. This study, therefore, serves to provide a proof of principle that MASH can be used for realistic simulations of photochemistry and can compete with more established methods, such as FSSH and AIMS.

Before describing the algorithm used for generating MASH trajectories, we first turn to the question of sampling initial conditions. Most nonadiabatic trajectory simulations are initialized using a Wigner function based on a harmonic approximation around the ground-state equilibrium geometry. Ordinarily, we would have used this standard approach within a MASH simulation. However, in the S₀ state, cyclobutanone has a low-frequency puckering mode around a C_2v geometry, which is very anharmonic and depending on the electronic-structure method used may even be predicted to be a double well with a low barrier. It is, therefore, clear that the harmonic approximation is not valid for this mode. At the level of density-functional theory (B3LYP/def2-TZVP as implemented in ORCA),³¹ we located the minimum-energy pathway between the two minima using the nudged-elastic band method³² and performed a “stream-bed walk”³³ (in Cartesian coordinates) up the other side. We call this the puckering path, and from now on, we employ mass-weighted coordinates using atomic masses of the most common isotopes, unless otherwise stated. A one-dimensional discrete variable representation (DVR) calculation³⁴ was carried out to obtain the nuclear wavefunctions along the puckering path. The predicted fundamental vibrational transition is 33 cm⁻¹ in good agreement with the experimental value of 35 cm⁻¹ from far-infrared spectroscopy.³⁵ There is also reasonable agreement with the higher excited vibrational states (supplementary material).³⁶ We learned that the experiment will slightly heat the sample to avoid condensation.³⁷ Therefore, one-dimensional positions and momenta were sampled according to the thermal Wigner function at 325 K, which can itself be evaluated in terms of the DVR wavefunctions (supplementary material).³⁸ Hessians were computed at a set of points along the puckering path and interpolated (in Cartesian coordinates). Translational and rotational modes along with the vector tangential to the path were projected out. It was found that the perpendicular frequencies are much larger than the energy-level spacing of the puckering vibrations and vary relatively slowly along the puckering path (supplementary material), which justifies our approach. Finally, the modes perpendicular to the puckering path were sampled using the thermal Wigner function within the standard harmonic approximation, and angular momentum was sampled from a classical Boltzmann distribution. More elaborate schemes for sampling Wigner distributions have been proposed, but these are not yet applicable to such large systems.³⁹ 

Now that the initial conditions for the nuclei are specified, we discuss the electronic-structure method used for the dynamical simulations. In order to capture the excited-state manifold and the bond-breaking dynamics after photoexcitation, a multireference method is required; we employ the state-averaged complete active space self-consistent field (SA-CASSCF) method in order to simultaneously describe the S₀, S₁, and S₂ states. In SA-CASSCF, the ground and excited electronic states are optimized simultaneously with a common set of orbitals but different configuration-interaction (CI) coefficients, which are constrained to form an orthonormal set. This ensures that the electronic states are orthogonal, which is particularly important when using the overlaps of the electronic wavefunction in the dynamics. The common set of orbitals also allow for efficient implementations of analytic gradients and nonadiabatic coupling vectors (NACVs). All CASSCF calculations were performed using Molpro 2023⁴⁰ and used a Slater basis with projection onto the singlet space. This is recommended as the default option within Molpro for CASSCF calculations due to increased computational efficiency over using configurational state functions.

The excitation from the ground electronic state, S₀, to the first excited state, S₁, is locally characterized by a transition from a non-bonding n orbital to an antibonding π* orbital of the carbonyl,^21,22,24,41 while the excitation to S₂ is characterized by a transition to a Rydberg 3s orbital.^25,26 In addition, previous experimental^21,22 and theoretical studies²⁴ indicate that C–C bonds are cleaved during the relaxation dynamics. It is, therefore, critical to choose a basis set that includes diffuse orbitals to accurately describe the Rydberg orbital and an active space that is able to characterize the excited state manifold while simultaneously allowing for a good description of the possible C–C bond breaking processes.

Due to the time constraints imposed by the prediction challenge, we were restricted by the size of the basis set and active space that could be used. We considered two different basis sets with diffuse functions, the Dunning aug-cc-pVDZ basis and the Pople 6-31+G* basis, and three different active spaces (detailed descriptions along with illustrations of the active spaces are provided in the supplementary material). We benchmarked vertical excitation energies (Table I of the supplementary material) and found that aug-cc-pVDZ with a (12,11) active space resulted in the best balance between accurate vertical excitation energies, computational cost, and being large enough to accurately describe the bond-breaking dynamics. Over the course of the dynamics, this (12,11) active space has sufficient flexibility to be able to describe three simultaneous bond breaking events. The calculations with aug-cc-pVDZ gave a vertical excitation energy (6.231 eV) that was closer to the experimental S₂ peak maximum (6.4 eV)²⁵ than with 6-31+G* (6.846 eV)⁴² and also much closer to the pump laser frequency of the planned experiment (6.2 eV).

On the basis of these data, unless otherwise stated, we choose to present the calculations performed using the aug-cc-pVDZ basis with a (12,11) active space in the main text. Analogous calculations using the 6-31+G* basis are, however, provided in the supplementary material, and differences between the two sets of calculations are used to help assess the sensitivity of the predicted results to the details of the electronic structure.

At the C_2v geometry, the n → π* transition is electric-dipole forbidden, while the n → 3s is electric-dipole allowed. In addition, the pump-pulse energy is on resonance with the S₂ excitation. We, therefore, utilized the Franck–Condon approximation to initialize the electronic state in S₂ (according to the MASH procedure described below). In this way, we allowed vertical transitions for the entire initial distribution and do not take account of the bandwidth of the laser pulse, as it is not obvious how to do this in a rigorous way without explicitly simulating the light field.

We optimized the structures of relevant minimum-energy conical intersections (MECIs) and crossing points (MECPs) and present their relative energies in the supplementary material. Both the S₂/S₁ MECI and the S₂/T₂ MECP are below the Franck–Condon energy, implying that they are both energetically accessible. However, the spin–orbit coupling (SOC) at the S₂/T₂ MECP is only 5 cm⁻¹, which suggests that intersystem crossing may be quite unlikely. Similar conclusions are indicated for relaxation from the S₁ state, in which the SOC is zero at the S₁/T₁ MECP. Taking cues from the study by Liu and Fang,²⁴ we identified three S₁/S₀ MECI structures, which are lower or comparable in energy to the S₁/T₁ MECP. Together, these results suggest that intersystem crossing is not important for the photodynamics of cyclobutanone. This is in agreement with previous theoretical studies of photochemistry of cyclic ketones on the S₁ state, which show that intersystem crossing only plays a role in molecules with rings of five or six carbon atoms.²³ 

Our approximation of neglecting the triplet states in the dynamics of cyclobutanone is further tested by selecting ten trajectories from the final aug-cc-pVDZ set, along which we simulated the electronic dynamics including three triplet states (i.e., a six-state SA-CASSCF) with a (12,11) active space, starting on the S₂ state. The resulting electronic dynamics showed less than 0.6% population transfer to the triplet manifold. Although this simple test is not completely reliable, especially considering the small number of trajectories considered, it, nonetheless, lends further weight to justify our neglect of the triplet states. Further details of these calculations and plots of the triplet populations over time for a couple of representative trajectories are given in the supplementary material.

We now turn to our choice of dynamics method, MASH. For clarity, here, we give a brief discussion of the important features of the method, highlighting key differences in FSSH. For a more detailed discussion, see the supplementary material as well as Refs. 3, 14, and 29. Before discussing the multi-state generalization, it is instructive to introduce the key ideas behind the original two-state implementation of MASH.³ In this case, there are two key differences to FSSH: first, how the active surface is determined, and second, how the initial electronic variables (wavefunction coefficients) are chosen. Unlike FSSH, in MASH, the active state is obtained deterministically. In the two-state case, the active state is determined by the sign of the z component of the Bloch sphere, $Sz(2,1)=|c2|2−|c1|2$⁠, such that when $Sz(2,1)>0$⁠, the active state is n = 2, and when $Sz(2,1)<0$⁠, the active state is n = 1. That the dynamics can be fully deterministic might at first seem surprising, particularly given that it is the stochastic nature of the hopping in FSSH that allows it to capture wavepacket bifurcation. However, the stochastic nature of the FSSH hops is effectively replaced in MASH by an initial sampling of the wavefunction coefficients. For example, to initialize a system in a pure state on adiabat 2, in FSSH one chooses $Sz(2,1)=1$ and hence $Sx(2,1)=Sy(2,1)=0$⁠. However, in MASH $Sz(2,1)$ is instead sampled from the probability density $ρ2(Sz(2,1))=2h(Sz(2,1))|Sz(2,1)|$ (where h(x) is the Heaviside step function), with $Sx(2,1)$ and $Sy(2,1)$ chosen uniformly from the corresponding circle on the Bloch sphere. It is due to this ensemble that MASH is able to describe wavepacket bifurcation.

MASH has been shown to offer a number of formal improvements over FSSH. First, unlike FSSH, it can be rigorously derived as a short-time approximation to the quantum–classical Liouville equation (QCLE), meaning that it can, in principle, be systematically improved toward this limit. Second, and perhaps most importantly, MASH does not suffer from the inconsistency error of FSSH. This is because the deterministic nature of the MASH algorithm means that the electronic variables are always directly related to the current active state. In contrast, in FSSH, the wavefunction coefficients can become completely inconsistent with the active state. These methodological improvements are expected to lead to more reliable predictions for no extra computational cost. In fact, for a series of different model systems, MASH has been shown to be as accurate or more accurate than FSSH at reproducing quantum-mechanical benchmark results.^3,14 Decoherence corrections can be derived rigorously for MASH,³ although in the vast majority of cases, the dynamics are already accurate enough without them.¹³ 

While the original MASH method was derived for two-state systems only, we have recently proposed an N-state generalization of MASH ideally suited for simulating photochemical systems, which we call the uncoupled-spheres multi-state MASH method (unSMASH).²⁹ The unSMASH method generalizes the original two-state MASH by treating possible transitions between pairs of adiabatic states independently. This is done by introducing N − 1 independent effective two-state Bloch spheres between the current active state, n, and each of the other states: $S(n,j)$ for j ≠ n, j = 1, …, N. Each sphere then evolves as it would in the original two-state MASH for the truncated electronic space consisting of the two corresponding adiabatic states. Attempted hops occur when one of the $Sz(n,j)$ changes sign. As in the two-state theory, the hops are accepted or rejected according to whether there is enough kinetic energy in the direction of the NACV between the active state and the possible new state. The component of the momentum along the NACV is then either rescaled to conserve energy in the case of allowed hops or reflected in the case of rejected (frustrated) hops. The unSMASH method is a rigorous short-time approximation to the QCLE when there is only coherence between one pair of adiabatic states at a given time. It is, therefore, well suited to photochemical problems involving a series of successive separate transitions between adiabatic surfaces, as one would expect in the photochemical relaxation of a typical organic molecule, such as cyclobutanone.

The integrator used to evolve the unSMASH equations of motion is closely related to those suggested previously for FSSH.^43–46 A full mathematical description of the integrator is given in the supplementary material; here, we simply give some of its key features. The basic structure of the integrator involves first propagating the nuclear positions and momenta from t to t + δt using velocity Verlet, before the electronic variables are evolved from t to t + δt using a unitary operator based on information calculated at geometries q(t) and q(t + δt). Finally, any attempted hops are treated, along with their associated momentum rescalings. As with FSSH propagation schemes, one must contend with the fact that close to conical intersections, the NACV is very sharply peaked. This means that algorithms that rely solely on the NACV can require arbitrarily small time steps in order to capture the corresponding electronic transition. For this reason, it is common to use an effective time-averaged nonadiabatic coupling that can be calculated from the overlap between the adiabatic wavefunctions at successive time steps rather than to compute the NACVs explicitly.^43–46 This can be implemented for MASH in the same way as it is for FSSH. However, in the present work, we found that the calculation of NACVs was significantly less computationally expensive than computing the overlaps. For this reason, the integrator we employ uses a mixed scheme, only calculating overlaps when the adiabatic surfaces come close together (|ΔV| < 2000 cm⁻¹), and otherwise using the NACVs. We note that using the NACVs rather than the overlaps has an additional advantage when combined with CASSCF, as it means that discontinuities in the active space do not lead to spurious nonadiabatic transitions. This is similar to the advantage of calculating the forces analytically rather than by finite difference, which can lead to unphysical sudden large changes in the momenta when encountering a discontinuity in the energy.

In total, we sampled 200 sets of initial positions and momenta from the Wigner function at 325 K. In each case, the initial active state was S₂, and the initial spheres were randomly sampled as described above. These initial samples were used to launch 200 separate unSMASH trajectories (using the aug-cc-pVDZ basis) with a time step of 0.5 fs. Of these initial 200, a total of 198 trajectories were run for the desired total time of 500 fs. For 123 of these trajectories, spin contamination in the excited-state manifold meant that the SCF step in the SA-CASSCF cycle did not converge within the maximum number of allowed iterations, which was chosen as 160. This typically only occurred at later times after the molecule had already reached S₀ and had already undergone the primary dissociation step, with the majority (77) having already reached at least 200 fs before the SA-CASSCF failed to converge. Rather than discarding these trajectories (which would bias the results), we elected to finish those that had already reached S₀ by running them for the remaining time on the ground state. This was done using state-specific CASSCF (SS-CASSCF) with the same basis and active space size as for the three-state SA-CASSCF calculations. Note that the initial momentum and position for the SS-CASSCF part of the trajectory were simply taken from the end of the previous convergent SA-CASSCF step, and the trajectory continued for the remaining time using velocity Verlet integration. The two trajectories that could not be completed were those that crashed due to spin contamination while on an excited electronic state; these trajectories are excluded from all analysis that follows.

In Fig. 2, we present the electron-diffraction signal predicted by our MASH simulation, both before and after convolution. From these results, one immediately obtains a qualitative picture of the dynamics after photoexcitation. In the unconvolved signal, after only 50 fs, the predicted electron-diffraction signal shows a significant positive peak at around $3.25Å$⁠, along with a corresponding negative peak in the region 1.25–2.5 $Å$⁠. At 100 fs, the positive peak has broadened and shifted toward $4Å$⁠, while the negative peak has deepened. This trend carries on until around 350 fs, after which point the negative peak approaches a steady state and the positive peak has broadened and shifted to such large distances as to become almost invisible. Although in the convolved signal the locations and heights of the peaks are somewhat modified, the same qualitative behavior can be observed. This behavior is clearly indicative of a rapid dissociation following the photoexcitation of cyclobutanone. The dissociation leads to a depletion of short bond distances due to bond breaking, with a corresponding increase at continually larger and larger distances as the resulting fragments move apart. In contrast, a simple ring-opening reaction, without dissociation, would result in a persistent positive signal between $4Å<r<6Å$⁠, as seen, for example, in the photoinduced ring-opening of cyclohexadiene.⁵⁰ From these results, we can ascertain that the majority of the dissociation occurs within the first 250–300 fs, with the onset of dissociation occurring very rapidly at around 50 fs.

Although the electron-diffraction signal contains a large amount of information and gives an immediate qualitative picture of the nuclear dynamics, it is, nevertheless, hard to immediately extract detailed mechanistic information from the signal alone. This is why, despite the ever-increasing resolution of experiments, molecular simulations are an important tool in understanding complex photochemical processes. In the following, we consider what the additional information available from our molecular simulations tells us about the dissociation process before returning to discuss how signatures of these features could be observed in the experimental electron-diffraction signal.

We first consider what our simulation predicts about the electronic dynamics after excitation. Figure 3 shows the average population on each adiabatic state as a function of time. From this we can clearly see that the system undergoes rapid electronic relaxation, with a half-life on S₂ of about 50 fs. It appears that the system primarily undergoes a sequential transition, first from S₂ to S₁ and then from S₁ to S₀. The resulting half-life for the combined excited-state manifold (S₂ + S₁) is predicted to be about 100 fs, with about 90% of the molecules having relaxed to S₀ by 250 fs. This timescale matches closely the dynamics seen in the electron-diffraction signal, indicating that the energy released into the nuclear degrees of freedom by the electronic relaxation leads rapidly to dissociation.

To obtain a clearer quantitative picture of the photodissociation process, it is helpful to analyze the trajectories according to the fragments formed in the dissociation process. We define a molecular fragment as a series of atoms that form a connected graph, where the nodes of the graph correspond to atoms and the edges of the graph indicate that the distance between the corresponding pair of atoms is within a cutoff distance of $r<2.0Å$⁠.⁵¹  Table I shows the total yield of all observed molecular fragments at 500 fs. We note that, due to the presence of secondary dissociation processes that occur after 500 fs, this is not likely to be equivalent to the final product distribution. We can see that the most abundant product is carbon monoxide (CO), closely followed by cyclopropane/propene (C₃H₆). There are also significant numbers of ethene (C₂H₄) as well as ketene (C₂H₂O) and the highly reactive methene (CH₂), along with a handful of other fragments that are only observed in a small number of trajectories. Most notable of these are the seven remaining C₄H₆O molecules that have not yet dissociated, of which three have already undergone ring-opening, three remain as cyclobutanone, and one has undergone a rearrangement to form cyclopropanal. We note that, based on the analysis that follows, we expect the majority of them to eventually dissociate to give C₃H₆ + CO.

To obtain a more detailed understanding of these dissociation pathways, in Fig. 4, we plot the time-dependent yield of each of the five major reaction products. It is instructive to first consider the yields of C₂H₄ and C₃H₆. As was observed in the electron-diffraction signal, the onset of dissociation occurs at around 50 fs, where the number of observed fragments begins to increase sharply, and the vast majority of the primary dissociation is over by around 300 fs, where the yields of both C₂H₄ and C₃H₆ are greater than 90% of their final value. It is notable that although the onset of formation of C₃H₆ begins slightly earlier than C₂H₄, the timescale associated with the formation of C₂H₄ is significantly shorter. In fact, while the yield of C₂H₄ is essentially constant between 200 and 500 fs, there is a notable 50% increase in the yield of C₃H₆ over this timescale. This can be understood by noting that (as shown in Fig. S24 of the supplementary material), trajectories that form C₂H₄ stay on S₂ for slightly longer but reach S₀ earlier than those that form C₃H₆. The rapid formation of C₂H₄ is thus associated with a sudden successive relaxation from S₂ to S₁ and then S₁ to S₀, releasing a large amount of energy and resulting in a rapid (almost concerted) bond breaking. While C₃H₆ can also form in this way, we see that there is another slower formation mechanism. This slower mechanism involves trajectories becoming temporarily trapped on S₁ without dissociating; when they eventually relax to S₀, they predominantly form C₃H₆ rather than C₂H₄. It seems likely that this is because the formation of C₂H₄ requires a greater amount of energy in the corresponding carbon–carbon bond stretch and that this slower pathway results in a more even distribution of the energy released from the relaxation from S₂ to S₁.

Returning to Fig. 4 we note that, in contrast to C₂H₄, the yield of CH₂ continues to increase beyond 200 fs. Specifically, the yield of CH₂ increases from around 4.5% at t = 200 fs to around 7.6% at t = 500 fs corresponding to an ∼70% increase in population. This can be understood as arising from a secondary dissociation of ketene, as defined by (R3a). This is confirmed by the concomitant decrease in the yield of ketene over the same period. However, we note that the appearance of methene before 100 fs indicates that secondary dissociation is not the only pathway to its formation. If it were, one would expect the rate of formation of CH₂ to be proportional to the population of ketene. We, therefore, conclude that a significant fraction of methene is formed via the primary dissociation reaction shown in (R3b).

Although there are differences in the relative fractions of each product formed, the 6-31+G* calculations show a similar qualitative behavior to these aug-cc-pVDZ calculations. In Subsection III D, we will discuss these quantitative differences further.

These observations, along with directly visualizing the trajectories, confirm that the dominant reaction pathways are those given in (R1)–(R3b) and allow us to make a prediction of the relative yields of the C3 and C2 pathways (the C3/C2 ratio). To do so, we choose in all cases to identify reaction products I as belonging to the C3 pathway and reaction products II and III to the C2 pathway.⁵⁵ The trajectories not in groups I, II, or III can be categorized as those that have undergone hydrogen dissociation and those that have not yet dissociated. We choose to exclude those reactions that have undergone hydrogen dissociation from the discussion of the C3/C2 ratio. For the seven remaining undissociated C₄H₆O, the fact that the yield of C₃H₆ in Fig. 4 still shows a significant positive slope at 500 fs makes it seem likely that the majority will eventually dissociate according to (R1). For this reason, we additionally categorize these trajectories as following the C3 pathway. On this basis, of the trajectories that follow either the C3 or C2 pathway, we observe that 77% dissociate via the C3 pathway and can estimate the influence of the statistical error using a 95% Wilson score interval⁵² to give an lower and upper bound to our prediction of (70%, 82%). The corresponding C3/C2 ratio is found to be 3.3, and propagating the 95% Wilson score interval gives upper and lower bounds on the statistical error as (2.35,4.6).

Following the same analysis of the 6-31+G* product yields (given in the supplementary material), one arrives at a somewhat different C3/C2 ratio. There, we find the fraction of all C3 and C2 trajectories that dissociate via C3 to be 0.5 with a 95% Wilson score confidence interval of (43%, 57%) corresponding to a C3/C2 ratio of 1 with the statistical confidence interval at 95% of (0.75,1.35). The two sets of simulations, therefore, show a statistically significant difference in this quantity, and we will return to discuss how this influences our estimation of the systematic error in our predicted electron-diffraction signal in Sec. III F.

As far as we are aware, there does not exist a definitive experimental measurement of the C3/C2 ratio at 200 nm. The only attempt at a direct measurement that we could find in the literature was that of Shortridge et al.⁵⁶ who obtained a value of 1.2. This experiment used a full arc zinc lamp rather than laser excitation, observing that only the 202.6 nm line showed appreciable absorption, and had a buffer gas of cyclopropane that somewhat complicated the interpretation of their results. In later work, Trentelman et al. chose to ignore this result and, instead, extrapolated the measurements of Denschlag and Lee⁵³ to obtain a ratio of 1.3 at 193 nm. It is, however, questionable whether this extrapolation is valid, given that it was based on excitation at wavelengths between 318 and 248 nm, which correspond to excitation to S₁ (centered at 280 nm) rather than S₂. We additionally note that, if one takes the C3/C2 ratios for the lowest two wavelengths reported by Denschlag and Lee (253.7 and 248 nm), then a linear extrapolation in energy results in a significantly higher C3/C2 ratio of 2.5, in much closer agreement with the aug-cc-pVDZ ratio prediction. In any case what is clear from the existing experimental results is that within the S₁ band, the C3/C2 ratio has a strong energetic dependence. Given this, systematic errors in the electronic structure and approximations used in the initial conditions might be expected to have a profound effect on the ratio seen in the simulations, which would certainly be consistent with the observed difference between the 6-31+G* and aug-cc-pVDZ results.

Analyzing our calculated trajectories has allowed us to give a detailed prediction of what happens during the dissociation. However, in order for our predictions to be testable, we need to connect them to observable features of the planned experimental signal. Hence, in the following, we return to consider how key features of the different reaction pathways could be observed in the experimental electron-diffraction signal. To do this, we begin by calculating three hypothetical electron-diffraction signals corresponding to the trajectories that produce reaction products I, II, and III. The full signals are given in the supplementary material. In the following, we connect key features of these signals to the total signal and analyze how sensitive the planned experiment is to the relative fractions of reactions that follow each pathway.

In Fig. 5, we show the minimum value of the ΔPDF curve as a function of time for each of the reactions alongside the total signal. The minimum value occurs in all reactions and at all times close to $2.5Å$⁠. This corresponds to the distance between the oxygen and the two β-carbons in cyclobutanone (Fig. 1) and also to the distance between the oxygen and the β-carbon in ketene. Hence, this minimum is closely associated with the formation of CO. We note, however, that drawing a direct connection between the depth and the amount of CO is complicated somewhat, both by the presence of two rather than just one β-carbons in cyclobutanone and by the fact that the distance between the carbons and the four hydrogens on the neighboring carbons in cyclopropane is about 2.5 $Å$ (although the lighter H atoms have a weaker signal). Nevertheless, comparing the curves for products I and products II, we can see a clear parallel with the C₃H₆ and C₂H₄ curves in Fig. 4. The curve for products II has clearly plateaued by around 350 fs, while products I continues to decrease gradually. Using this connection, we can attribute the notably shorter timescale of the curve for products II vs products I to the “concerted” vs “sequential” nature of the dissociation reactions. In a similar manner, the slope at long time for products III can be understood as arising from the secondary dissociation of ketene to form methene and carbon monoxide.

Although these qualitative differences in the slopes are interesting, the most notable difference between the curves is their magnitude at 500 fs. This would, therefore, appear a good way of using the experimental result to judge the accuracy of our predicted product yields (and the C3/C2 ratio). We remind the reader that in order to compare the results shown here to the future experimental results, one should bear in mind that the results here are given relative to the maximum peak height of the steady-state PDF [Eq. (1)].⁵⁷ If there were only two dominant sets of reaction products, the depth of the minimum alone would allow one to get a good estimate of the ratio of products, which could be done by fitting a weighted average of the corresponding peak depths. Since there are three dominant reaction products, determining the relative ratios of all three from the electron-diffraction signal is not so straightforward. However, by considering the full ΔPDF(r) signal at 500 fs, one can additionally make use of the line shape to help distinguish the reaction products present. Figure 6 does exactly that. Alongside the signal obtained from our MASH simulation (dashed line), it shows the hypothetical signal that would be obtained for different fractions f₁ and f₂ of products I and products II, respectively, with the remaining contribution from products III (f₁ + f₂ + f₃ = 1). (Note that the individual signals that are being averaged can be seen in Figs. S25–S27 of the supplementary material.) This serves as a relatively simple testable quantitative prediction that one can use to assess the relative fractions of these key reaction products. While it is not highly sensitive, there is only a relatively small part of the parameter space, which is consistent with the aug-cc-pVDZ results, i.e., 0.65 ≲ f₁ ≲ 0.85. Note that the equivalent plot for the 6-31+G* calculations is given in the supplementary material and shows that in this case, only the region 0.35 ≲ f₁ ≲ 0.65 would be consistent with the predicted result. This demonstrates that, although the predicted electron-diffraction signal is not radically different in the two cases, with a low enough experimental error, one would be able to distinguish between them.

In case it is difficult to accurately determine the depth from the experimental signal, in the supplementary material, we also consider how the line shape at 500 fs alone could be used to distinguish different reaction products. Figure S28 shows the possible signals at 500 fs normalized so that the depth at their minimum is −1. From this, one can see that the most notable change to the line shape with varying fractions, f₁, of reaction products I, is the depth of the shoulder at about 1.8 $Å$⁠. Although the changes are relatively subtle compared to the absolute differences of Fig. 6, we see that, with a small enough relative error, it should, nevertheless, be possible to distinguish between different fractions, f₁, of reaction products I and, hence, different C3/C2 ratios. This could, therefore, also be used as a direct test of our prediction of the product ratio.

Finally, having focused on how sensitive the experiment would need to be to distinguish between different possible reaction pathways, it is natural to ask how confident we are in our predicted signal and how large our systematic errors are likely to be. There are four potential sources of error as follows:

1. electronic-structure theory

2. initial conditions/treatment of the pulse

3. dynamics (including zero-point energy leakage)

4. calculation of the electron-diffraction signal.

The established sensitivity of the C3/C2 ratio on the excitation wavelength means that we expect that the results will be particularly sensitive to errors in the electronic-structure theory and initial conditions. Given the accuracy of MASH in previous benchmark tests and the established nature of the independent atom approximation, we focus here on assessing the errors due to the electronic-structure theory and initial conditions. We have, therefore, performed a number of additional calculations for which additional figures are given in the supplementary material.

Ideally to confirm the accuracy of our prediction, we would perform the same simulations at a higher level of electronic-structure theory, e.g., CASPT2 with the same (12,11) active space. Given the time constraints of the challenge, this is not possible. However, we have calculated CASPT2 energies with the aug-cc-pVDZ basis along linear interpolated (in internal coordinates) paths (Fig. S12), to examine the effect of dynamic correlation. Along the path that leads from the Franck–Condon geometry to the S₂/S₁ MECI and the S₁ minimum, CASPT2 energies are very similar to CASSCF. This is also observed along the path connecting the S₁ minimum to the first S₁/S₀ MECI. However, the barriers to access the second and third S₁/S₀ MECIs increase by 0.7 and 0.8 eV, respectively, at the CASPT2 level of theory. Since the third S₁/S₀ MECI seems to lead to the C2 products (Fig. S20), it is likely that repeating our aug-cc-pVDZ simulations at the CASPT2 level would further increase the C3/C2 ratio.

Another thing that might influence the C3/C2 ratio is the initial conditions. Hence, to assess the extent to which the initial conditions may influence the predicted electron-diffraction signal, we considered how the vertical excitation energy of each initially sampled geometry correlated with the final reaction products. Under the assumption that the strong wavelength dependence observed in the product ratio at longer wavelengths (within the S₁ window) implies a strong wavelength dependence when exciting to S₂, one would expect the initial excitation energy in our simulation to correlate with the final product. To test this hypothesis, we performed a two-sample z-test on the mean of the vertical excitation energy for products I and products II + III. This resulted in a difference in the means that was not statistically significant, with a p-value of 0.55. This goes against the hypothesis stated above. This implies that the C3/C2 ratio is not so sensitive at these shorter wavelengths, probably because the large excess of energy makes all conical intersections accessible, in contrast to excitation to S₁ where shorter wavelengths may open new channels, significantly altering the product distribution. Hence, it appears from these tests at least that the initial conditions are perhaps not such a large source of systematic error. Of course, there are other features of the initial conditions, such as the dependence of transition-dipole moments, which we have not included but could affect the product ratios.

Our most direct information about possible systematic errors comes from the difference between the aug-cc-pVDZ calculations and the 6-31+G* calculations. As discussed above, when comparing the dynamics within a particular channel, both sets of calculations give very similar results. However, where they differ is in their predicted C3/C2 ratio (statistical intervals giving 2.35–4.6 for aug-cc-pVDZ vs 0.75–1.35 for 6-31+G*). We have shown that this difference could, in principle, be observed in the electron-diffraction signal, but the overall difference is small. If one were forced to choose between the two without considering prior photochemical experiments, one would say that aug-cc-pVDZ should be preferred since it is a slightly larger basis than 6-31+G*. It was for this reason, in addition to the more accurate S₀ to S₂ vertical excitation energy that we chose to focus on the aug-cc-pVDZ results in the main paper. We note, however, that the aug-cc-pVDZ calculations are not expected to be so much more accurate than the 6-31+G* calculations that we should discount them entirely. Given that prior experimental results indicate that the C3/C2 ratio should be around 1.2–1.3,^53,54,56 one might well suggest that the most accurate prediction could be obtained by taking an average of the two sets of calculations. We note, however, that the electron-diffraction signal that results from averaging (shown in the supplementary material) is difficult to distinguish by eye from the signal in Fig. 2 or the 6-31+G* result (also in the supplementary material). Hence we can conclude that, while there is some uncertainty in our predicted C3/C2 ratio, it seems reasonable to assume that the fraction of trajectories following the C3 pathway is between about 0.5 and 0.85. Overall, given the sensitivity of the electron-diffraction signal observed in Fig. 6 and additional figures of the supplementary material, we can, therefore, be confident that our results are likely to accurately reproduce the timescales and other key features of the planned experimental signal.

The only caveat to this is that, since writing our original manuscript, we have become aware that there do exist previous measurements of time-resolved photoelectron and mass spectra of cyclobutanone after excitation with a 200 nm pulse.⁵⁸ These measurements were used to estimate the timescale of S₂ decay and found a timescale of about 740 fs. If correct, this would indicate that the S₂ decay predicted by the current study is too fast. We note that our scans between the Franck–Condon point and the S₂/S₁ MECI (shown in Figs. S12 and S13 of the supplementary material) do show a significant barrier at the CASPT2 level that is not present in CASSCF, which could contribute to such an error in the timescale of S₂ decay.

In conclusion, we have used molecular simulations to investigate the dynamics of cyclobutanone after photoexcitation to the S₂ state. We have focused here on details of the dynamics most relevant to JCP’s community challenge. In particular, we have made quantitative predictions of the electron-diffraction signal that will be observed in the planned experiment and explained this signal in terms of the main fragmentation reactions observed in our simulation.

In addition to understanding the dissociation of cyclobutanone, this study also serves as the first application of the newly proposed multi-state MASH method, called unSMASH, to an ab initio simulation. By choosing to use MASH for our simulations instead of the more commonly used FSSH, we expect to have gained the following advantages. First, our simulations require no ad hoc decoherence corrections. This is because the deterministic dynamics of MASH means that the electronic variables always remain consistent with the active surface, therefore fixing the inconsistency problem of FSSH.^3,14 This is likely to be important for accurately describing the photodissociation process, where it is known that the inconsistency error can lead to suppressed product yields.¹³ Second, MASH is able to correctly describe the effects of the nonadiabatic force through its uniquely determined momentum rescaling algorithm. This is crucial for accurately describing the electronic population dynamics, which is known to be particularly sensitive to how the momentum rescaling is performed.^13,59 While FSSH can, in principle, describe this effect by rescaling the momenta along the NACV at a hop and reflecting in the case of a frustrated hop, in practice, it is common to use an isotropic momentum rescaling in ab initio FSSH simulations.¹⁰ 

While we expect the present study to have captured the most important details of the system, there are a number of areas in which the simulation could be improved. First, we could include intersystem crossing between the singlet manifold and the triplet manifold. This is something that can be rigorously achieved within the MASH framework, as will be described in a forthcoming publication.³⁰ We note, however, that given the rapid nature of the dissociation observed in the present study and the weak spin–orbit coupling, it is unlikely that there would be time for sufficient population transfer to the triplet manifold to significantly influence the dynamics. Second, we may seek to improve the accuracy of our initial conditions. The branching ratio between C3 and C2 channels may depend sensitively on the amount of energy imparted to the system by the photoexcitation.^53,54 The Franck–Condon approach we have used, taking the initial distribution as the ground-state nuclear Wigner distribution placed on the S₂ state and mimicking the finite width of the exciting laser pulse by convolution with a Gaussian, is a standard way of initializing a simulation in photoexcited systems;⁵⁰ however, it is not without approximation. It will, therefore, be interesting in future studies to explore the sensitivity of the results to the initial conditions used, or even explicitly simulate the pulse within the MASH dynamics.³⁰ 

We note, however, that despite the possible improvements that could be made to the initial conditions, at present, they would all rely on the use of a Wigner-transformed initial density. This has the advantage that it introduces zero-point energy into the initial distribution. However, this is not without its limitations. Specifically, since the underlying nuclear dynamics is classical, the zero-point energy will eventually (due to anharmonicity) become evenly distributed among the molecular degrees of freedom. This so-called zero-point energy leakage is a well-known problem.⁶⁰ In condensed-phase problems in thermal equilibrium, this issue has been largely solved by imaginary-time path-integral methods, such as ring-polymer molecular dynamics (RPMD).^61–63 There has, therefore, been significant interest in the development of a nonadiabatic version of RPMD.^64–70 However, for photochemical problems that are inherently far from the linear-response regime, where imaginary-time path-integral methods have been shown to give very accurate results,^71–75 it is unclear whether a nonadiabatic RPMD would be the final solution to this problem. This is because RPMD effectively assumes a rapid decoherence of vibrational modes that may introduce additional errors in such low-pressure gas-phase systems.^76,77 The search for an optimal method for such problems, therefore, continues. In future studies, it would be interesting to assess the importance of zero-point energy leakage in this system by initializing all modes from a classical Boltzmann distribution and comparing the resulting simulations. If the classical nuclear limit is valid, then one has the additional advantage that MASH has the correct detailed balance to thermalize to the correct distribution.³ 

The aspect that will probably have the most significant impact on our dynamics is the accuracy of the electronic-structure theory. In the present study, we have chosen a level of theory that minimizes the cost while still being sufficient to allow the simulation to describe the key qualitative features of the electronic subspace. We used two different basis sets and found that although the qualitative description of the electron-diffraction pattern is similar, the C3/C2 product ratio is significantly different. In future work, it will, therefore, be interesting to investigate the system using even more accurate electronic-structure theory methods, such as using larger basis sets and including dynamic correlation via CASPT2, MRCI, or coupled cluster theory.⁷⁸ Although this would significantly increase the cost of an on-the-fly simulation, by exploiting modern machine-learning techniques,^79,80 one may hope to make such calculations tractable.

There are more advantages of developing a machine-learned model than just making the use of high accuracy electronic-structure theory achievable. In particular, having trained a model, it would then be comparatively inexpensive to perform a systematic analysis of the sensitivity to different initial conditions. Furthermore, it would give the opportunity to use more expensive dynamics methods that are impractical for on-the-fly calculations. Such studies would complement the aims of the JCP challenge, by providing an objective comparison of the accuracy of different dynamics methods and the approximations they make and helping to push the boundaries of accuracy in excited-state simulations. Given the right potential-energy surfaces and couplings, we believe that MASH can be a very powerful simulation tool for obtaining reliable predictions of photochemistry.

See the supplementary material for additional numerical results referred to in the main text.

J.E.L. was supported by an Independent Postdoctoral Fellowship at the Simons Center for Computational Physical Chemistry, under a grant from the Simons Foundation (839534, MT). J.R.M. and A.K. are supported by the Cluster of Excellence “CUI: Advanced Imaging of Matter” of the Deutsche Forschungsgemeinschaft (DFG)—EXC 2056—Project ID 390715994.

The authors have no conflicts to disclose.

J.E.L. MASH dynamics code development (lead), data analysis (lead), writing manuscript first draft (lead), running trajectories (lead), code testing (equal), development of interface between MASH code and Molpro (supporting). I.M.A. selection of electronic-structure method (lead), development of interface between MASH code and Molpro (equal), initial conditions (supporting), triplet simulations (supporting), writing manuscript first draft (supporting), code testing (equal), analysis of accuracy of electronic-structure theory (supporting). J.R.M. code testing (equal), writing manuscript first draft (supporting). M.A.M. analysis of accuracy of electronic-structure theory (lead), development of interface between MASH code and Molpro (equal), writing manuscript first draft (supporting). K.A. triplet simulations (lead), writing manuscript first draft (supporting). A.K. writing manuscript first draft (supporting). J.O.R. initial conditions (lead), calculation of electron-diffraction signal (lead), writing manuscript first draft (supporting), conceptualization (lead), project management (lead).

Joseph E. Lawrence: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Imaad M. Ansari: Investigation (supporting); Software (equal); Validation (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Jonathan R. Mannouch: Conceptualization (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Meghna A. Manae: Investigation (supporting); Software (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting). Kasra Asnaashari: Validation (equal); Writing – original draft (supporting). Aaron Kelly: Writing – original draft (supporting). Jeremy O. Richardson: Conceptualization (lead); Formal analysis (equal); Investigation (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.