A MASH simulation of the photoexcited dynamics of cyclobutanone

Joseph E. Lawrence
Department of Chemistry and Applied Biosciences, ETH Zurich, 8093 Zurich, Switzerland
Simons Center for Computational Physical Chemistry, New York University, New York, NY 10003, USA
Department of Chemistry, New York University, New York, NY 10003, USA

Imaad M. Ansari
Department of Chemistry and Applied Biosciences, ETH Zurich, 8093 Zurich, Switzerland
Jonathan R. Mannouch
Hamburg Center for Ultrafast Imaging, Universität Hamburg and the Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany
Meghna A. Manae
Kasra Asnaashari
Department of Chemistry and Applied Biosciences, ETH Zurich, 8093 Zurich, Switzerland
Aaron Kelly
Hamburg Center for Ultrafast Imaging, Universität Hamburg and the Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany
Jeremy O. Richardson
Department of Chemistry and Applied Biosciences, ETH Zurich, 8093 Zurich, Switzerland

In response to a community prediction challenge, we simulate the nonadiabatic dynamics of cyclobutanone using the mapping approach to surface hopping (MASH).

We consider the first 500 fs of relaxation following photo-excitation to the 2 state and predict the corresponding time-resolved electron-diffraction signal that will be measured by the planned experiment.
397 ab-initio trajectories were obtained on the fly with state-averaged complete active space self-consistent field (SA-CASSCF) using a (12,11) active space. To obtain an estimate of the potential systematic error 198 of the trajectories were calculated using an aug-cc-pVDZ basis set and 199 with a 6-31+G* basis set.

MASH is a recently proposed independent trajectory method for simulating nonadiabatic dynamics, originally derived for two-state problems.
As there are three relevant electronic states in this system, we used a newly developed multi-state generalisation of MASH for the simulation: the uncoupled spheres multi-state MASH method (unSMASH).
This study therefore serves both as an investigation of the photo-dissociation dynamics of cyclobutanone, and also as a demonstration of the applicability of unSMASH to ab-initio simulations.

In line with previous experimental studies
, we observe that the simulated dynamics is dominated by three sets of dissociation products, C3H6 + CO, C2H4 + C2H2O and C2H4 + CH2 + CO, and we interpret our predicted electron-diffraction signal in terms of the key features of the associated dissociation pathways.

§ INTRODUCTION

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)<cit.> is the most commonly used. In addition, the favourable 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.<cit.>

The mapping approach to surface hopping (MASH)<cit.> is a recently proposed alternative to the FSSH algorithm.
It has a rigorous derivation based on mapping approaches<cit.> 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.
Also, 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,<cit.>
many alternative suggestions have been made,<cit.>
and in practice, approximations are often taken.<cit.>
Previous work has shown that MASH is often more accurate than FSSH for a range of model systems<cit.> and there is reason to believe it can even be more accurate than ab initio multiple spawning (AIMS)<cit.> for photochemical problems.<cit.>
MASH was shown to correctly recover Marcus theory rates,<cit.> where FSSH is known to require complicated decoherence corrections.<cit.>
Additionally, unlike other mapping approaches,<cit.>
MASH rigorously captures the detailed balance necessary to thermalize to the correct equilibrium distribution.<cit.>

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 0 to the 2 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 1,<cit.> meaning that the present study covers new ground. Simulating the photoexcited dynamics initialized in 2 also poses additional theoretical challenges. 2 is known to be a Rydberg state<cit.> 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,<cit.>
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.<cit.>

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.<cit.>
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.

§ METHODS

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 0 state, cyclobutanone has a low-frequency puckering mode around a 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),<cit.> we located the minimum-energy pathway between the two minima using the nudged-elastic band method <cit.> and performed a “stream-bed walk”<cit.> (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 <cit.> was carried out to obtain the nuclear wavefunctions along the puckering path.
The predicted fundamental vibrational transition is 33 cm$^{-1}$ in good agreement with the experimental value of 35 cm$^{-1}$ from far-infrared spectroscopy.<cit.>
There is also reasonable agreement with the higher excited vibrational states (supplementary material).
[Such good agreement was not achieved when using a simple scan of the potential along the normal mode corresponding to the puckering motion.]
We learned that the experiment will slightly heat the sample to avoid condensation.<cit.>
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).
[We note that a classical Boltzmann sampling along the puckering path would probably have been sufficient due to the fact that the thermal energy is sufficiently larger than the energy-level spacing of the puckering vibrations.  However, our approach is more generally applicable at any temperature and is consistent with the quantum-mechanical treatment of the perpendicular modes.]
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.<cit.>

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 photo-excitation, 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 $\S0$, $\S1$ and $\S2$ states.

In SA-CASSCF, the ground and electronic excited 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 (NACV).
All CASSCF calculations were performed using Molpro 2023,<cit.> 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, 0, to the first excited state, 1, is locally characterized by a transition from a non-bonding $n$ orbital to an antibonding $\pi^{*}$ orbital of the carbonyl<cit.> while the excitation to $\S2$ is characterized by a transition to a Rydberg $3s$ orbital.<cit.>
Additionally, previous experimental<cit.> and theoretical studies<cit.> 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 information).
We benchmarked vertical excitation energies (Table I in the supplementary information) and found that the 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 aug-cc-pVDZ calculations gave a vertical excitation energy (6.231) that was closer to the experimental 2 peak maximum  (6.4)<cit.> than with 6-31+G* (6.846)[This is true even after accounting for the width of the rather wide absorption peak observed in <cit.>, with a full width half maximum of around 0.3.] and also much closer to the pump laser frequency of the planned experiment (6.2).

On the basis of this 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 details of the electronic structure.

At the 2v geometry, the $n \to \pi^{*}$ transition is electric-dipole forbidden while the $n \to 3s$ is electric-dipole allowed.
In addition, the pump-pulse energy is on resonance with the 2 excitation.

We therefore utilized the Franck–Condon approximation to initialize the electronic state in 2 (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 (MECI) and crossing points (MECP) and present their relative energies in the supplementary material.

Both the 2/1 MECI and the 2/2 MECP are below the Franck–Condon energy, implying that they are both energetically accessible.  However, the spin-orbit coupling (SOC) at the S$_2$/T$_2$ MECP is only 5 cm$^{-1}$ which suggests that intersystem crossing may be quite unlikely. Similar conclusions are indicated for relaxation from the S$_1$ state, in which the SOC is zero at the 1/1 MECP.
Taking cues from the study by Liu and Fang,<cit.> we identified three S$_1$/S$_0$ MECI structures, which are lower or comparable in energy to the S$_1$/T$_1$ 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 1 state, which show that intersystem crossing only plays a role in molecules with rings of 5 or 6 carbon atoms.

Our approximation of neglecting the triplet states in the dynamics of cyclobutanone is further tested by selecting 10 trajectories from the final
set, along which we simulated the electronic dynamics including 3 triplet states (i.e., a 6-state SA-CASSCF) with a (12,11) active space, starting on the 2 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, we give here a brief discussion of the important features of the method, highlighting key differences to FSSH.
For a more detailed discussion see the supplementary material as well as  Refs. ,  and . Before discussing the multi-state generalisation it is instructive to introduce the key ideas behind the original two-state implementation of MASH.<cit.> 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, $S^{(2,1)}_z=|c_2|^2-|c_1|^2$, such that when $S^{(2,1)}_z>0$ the active state is $\activestate=2$, and when $S^{(2,1)}_z<0$ the active state is $\activestate=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 $S^{(2,1)}_z=1$ and hence $S^{(2,1)}_x=S^{(2,1)}_y=0$.  However, in MASH $S^{(2,1)}_z$ is instead sampled from the probability density $\rho_2(S^{(2,1)}_z)=2 h(S^{(2,1)}_z) |S^{(2,1)}_z|$ (where $h(x)$ is the Heaviside step function), with $S^{(2,1)}_x$ and $S^{(2,1)}_y$ 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. Firstly, 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 towards this limit. Secondly, 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.<cit.>
Decoherence corrections can be derived rigorously for MASH,<cit.> although in the vast majority of cases, the dynamics are already accurate enough without them.<cit.>

While the original MASH method was derived for two-state systems only, we have recently proposed an $N$-state generalisation of MASH ideally suited for simulating photochemical systems, which we call the uncoupled-spheres multi-state MASH method (unSMASH).<cit.>

The unSMASH method generalises 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 and each of the other states: $\bm{S}^{(n,j)}$ for $j\neq n$, $j=1,\dots,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 $S_z^{(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.<cit.>  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+\delta t$ using velocity Verlet, before the electronic variables are evolved from $t$ to $t+\delta t$ using a unitary operator based on information calculated at geometries $\bm{q}(t)$ and $\bm{q}(t+\delta 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.<cit.> 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 ($|\Delta V|<2000\,{\rm cm}^{-1}$), 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.

Calculated steady-state PDF (in $\AA^{-2}$) for the initial distribution.  Also shown in orange is a histogram of the atom pair distribution (carbons and oxygens only).
The inset shows the atom pair distances (in $\AA$) for the $\C{2v}$ geometry.

§ RESULTS

Simulated ultrafast electron-diffraction results.
The panels on the left show the change in the probability density function relative to the initial configuration.
The panels on the right show the same data convolved with a 160 fs (FWHM) Gaussian to simulate the instrument response function.
Blue is loss, red is gain, with equally spaced contour levels showing the height of the $\Delta \mathrm{PDF}(r)$ signal relative to the maximum peak height in the steady-state PDF.

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 2, 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 0, 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 0 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, 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 2 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.

To compare our simulated results with the results of the proposed experiment, the final set of 198 trajectories (of length 500 fs) were used to generate electron-diffraction signals. This was done based on elastic-scattering
calculations within the
independent-atom model <cit.> using the ELSEPA program<cit.> as described in the supplementary material.
Note that it would in principle be possible to go beyond this approximation, in the spirit of Ref. , using the ab initio two-electron density of the MASH active state, provided by the CASSCF calculations.
However, in this work, we assume the independent-atom model to be sufficient.
The electron-diffraction signal was transformed (including a Gaussian damping function)<cit.> to obtain the atomic pair distribution functions (PDF).
Then $\Delta$PDF was defined as the difference between the PDF at time $t$ and the steady-state PDF as shown in Fig. <ref>.
As experimental electron-diffraction results are typically presented with arbitrary units,<cit.>
all results are given relative to the maximum peak height in the steady-state PDF, i.e.

\begin{equation} \label{PDF}
\Delta \mathrm{PDF}(r,t) =\frac{ \mathrm{PDF}(r,t)-\mathrm{PDF}_{\rm SS}(r)}{\max(\mathrm{PDF}_{\rm SS}(r))}
\end{equation}

Finally, the results were convoluted with a Gaussian (160 fs FWHM) to simulate the instrument response function (accounting for both the width of the pump pulse as well as the detector).

In Fig. <ref>, we present the 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 photo-excitation.

In the unconvolved signal, after only 50 fs the predicted electron-diffraction signal shows a significant positive peak at around $3.25\,\AA$, along with a corresponding negative peak in the region 1.25–2.5 $\AA$.
At 100 fs the positive peak has broadened and shifted towards $4\,\AA$ 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 behaviour can be observed. This behaviour 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\,\AA<r<6\,\AA$, as seen for example in the photo-induced ring-opening of cyclohexadiene.<cit.> 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.

§.§ Electronic dynamics

We first consider what our simulation predicts about the electronic dynamics after excitation. Figure <ref> 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 2 of about 50 fs. It appears that the system primarily undergoes a sequential transition, first from 2 to 1 and then from 1 to 0. The resulting half life for the combined excited-state manifold (2 + 1) is predicted to be about 100 fs, with about 90% of the molecules having relaxed to 0 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.

§.§ Reaction products

To obtain a clearer quantitative picture of the photodissociation process it is helpful to analyse 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 nodes of the graph correspond to atoms and edges of the graph
indicate that the
distance between the corresponding pair of atoms
is within a cutoff distance of $r<2.0\,\AA$.[Note the same cutoff was used for all pairs of atoms]  Table <ref> 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 (C3H6). There are also significant numbers of ethene (C2H4) as well as ketene (C2H2O) and the highly reactive methene (CH2), along with a handful of other fragments that are only observed in a small number of trajectories. Most notable of these are the 7 remaining  C4H6O molecules that have not yet dissociated, of which 3 have already undergone ring-opening, 3 remain as cyclobutanone and 1 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 C3H6 + CO.

Average (unconvolved) adiabatic populations as a function of time for the 198 trajectories. Shaded region shows an approximate 95% confidence interval (the Wilson score interval).<cit.>

Further insight can be gained by grouping the trajectories according to the molecular fragments they produce. Table <ref> shows the three most common sets of products present at 500 fs, along with the frequency at which they are observed. That these three reaction products should dominate is not a surprise based on previous theoretical and experimental work<cit.> such as that of Trentelman et al.<cit.> who rationalised their

measurements on the CO produced
after photo-excitation of cyclobutanone with 193 nm light in terms of these 3 possible reaction products. Following this earlier work, it is helpful to distinguish reactions according to whether the  primary dissociation event produces ethene (C2H4), labelled the C2 pathway, or cyclopropane/propene (C3H6), labelled the C3 pathway.
The C3 pathway is the simplest,

C4H6O →C3H6 + CO,

involving the cleavage of two carbon–carbon bonds to form carbon monoxide and a C3H6 diradical, which typically rapidly forms highly vibrationally excited cyclopropane.
In a small number of trajectories the C3H6 radical was observed to undergo a rearrangement to form the more stable propene, and we note that on longer timescales one would expect the excited cyclopropane to also undergo this rearrangement.

The C2 pathway is more complicated. In their work Trentelman et al. considered this to consist of a primary dissociation step,

C4H6O →C2H4 + C2H2O,

forming ethene and ketene (ethenone), with a possible secondary dissociation step,

C2H2O →CH2 +CO,

in which ketene dissociates to form carbon monoxide and methene. Here, however, we also consider the possibility of a third process in which both (R<ref>) and (R<ref>) occur in a single primary dissociation step,

C4H6O →C2H4 + CH2 +CO.

Total product yields  500 fs after initial excitation. Note the total number of initial C4H6O molecules is 198. Fragments are identified by using a cutoff radius of $2\AA$.

Fragment     Count     Yield

CO     158     80%

C3H6     141     71%

C2H4     45     23%

C2H2O     30     15%

CH2     15     7.6%

C4H6O     7     3.5%

H     5     2.5%

C4H5O     3     1.5%

C3H5     2     1.0%

Main reaction products at 500 fs. Reaction products are identified by using a cutoff radius of $2\AA$.

I     II     III

Products     C3H6 + CO     C2H4 + C2H2O      C2H4 + CH2 + CO

Count     141 (71%)      30 (15%)     15 (7.6%)

§.§ Time-dependent fragment formation

To obtain a more detailed understanding of these dissociation pathways,
in Fig. <ref> we plot the time-dependent yield of each of the 5 major reaction products. It is instructive to first consider the yields of C2H4 and C3H6. 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 C2H4 and C3H6 are greater than 90% of their final value. It is notable that although the onset of formation of C3H6 begins slightly earlier than C2H4, the timescale associated with the formation of C2H4 is significantly shorter. In fact, while the yield of C2H4 is essentially constant between 200 and 500 fs there is a notable $50\%$ increase in the yield of C3H6 over this timescale.
This can be understood by noting that (as shown in the supplementary material Fig. S24), trajectories which form C2H4 stay on 2 for slightly longer but reach 0 earlier than those which form C3H6.

The rapid formation of C2H4 is thus associated with a sudden successive relaxation from 2 to 1 and then 1 to 0, releasing a large amount of energy and resulting in a rapid (almost concerted) bond breaking. While C3H6 can also form in this way, we see that there is another slower formation mechanism. This slower mechanism involves trajectories becoming temporarily trapped on 1 without dissociating; when they eventually relax to 0 they then predominantly form C3H6 rather than C2H4. It seems likely that this is because the formation of C2H4 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 2 to 1.

Average (unconvolved) fragment yields, for the 5 most common fragments, as a function of time for the 198 trajectories. Shaded region shows an approximate 95% confidence interval (the Wilson score interval).<cit.>

Returning to Fig. <ref> we note that, in contrast to C2H4, the yield of CH2 continues to increase beyond 200 fs.

Specifically, the yield of CH2 increases from around 4.5%
at $t=200$ fs to around 7.6% at $t=500$ fs corresponding to an approximately $70$% increase in population. This can be understood as arising from a secondary dissociation of ketene, as defined by (R<ref>).

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 CH2 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 (R<ref>).

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

Minimum value of the simulated electron-diffraction signal given as a change in probability density function relative to the initial configuration.
Dashed lines show the unconvolved results from the 198 trajectories, and solid lines show the convolved results using a 160 fs (FWHM) Gaussian to simulate the instrument response function.
Note the height of the $\Delta \mathrm{PDF}(r)$ signal in all cases is given relative to the maximum peak height in the steady-state PDF.

Weighted average over the three dominant reaction products of the simulated ultrafast electron-diffraction results at 500 fs. $f_1$ corresponds to the fraction of products I and $f_2$ to the fraction of products II used in the weighted average, with the remaining fraction corresponding to products III.
Results are shown here  for the signal convoluted with with a 160 fs (FWHM) Gaussian to simulate the instrument response function. Dashed line shows the prediction from the aug-cc-pVDZ trajectories.

§.§ C3 vs C2 ratio

These observations, along with directly visualising the trajectories, confirm
that the dominant reaction pathways are those given in (R<ref>) to (R<ref>), 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. [We note that in 2 cases we observe that C2H4 is formed via dissociation of a very short lived C3H6 radical. However, given that in both cases CO is still in close proximity, and given that it would be difficult to distinguish these trajectories experimentally we feel it is reasonable to include these trajectories with the C2 pathway, effectively considering them as being described by (R<ref>).]
The trajectories not in groups I, II or III can be categorised 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 7 remaining undissociated C4H6O,  the fact that the yield of C3H6 in Fig. <ref> still shows a significant positive slope at 500 fs makes it seem likely that the majority will eventually dissociate according to (R<ref>).
For this reason we additionally categorise these trajectories as following the C3 pathway.

On this basis, of the trajectories that follow either the C3 or C2 pathways, we observe that $77\%$ dissociate via the C3 pathway, and can estimate the influence of statistical error using a 95% Wilson score interval<cit.> 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. <ref>.

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.<cit.> who obtained a value of 1.2. This experiment used a full arc zinc lamp rather than laser excitation, observing 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<cit.> 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 nm and 248 nm, which correspond to excitation to 1 (centered at 280 nm) rather than 2. We additionally note that, if one takes the C3/C2 ratios for the lowest two wavelengths reported by Denschlag and Lee (253.7 nm 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 1 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.

§.§ How will the different reaction pathways influence the experimental signal?

Analysing 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 analyse how sensitive the planned experiment is to the relative fractions of reactions that follow each pathway.

In Fig. <ref> we show the minimum value of the $\Delta \mathrm{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\,\AA$. This corresponds to the distance between the oxygen and the two -carbons in cyclobutanone (Fig. <ref>), 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 also by the fact that the distance between the carbons and the four hydrogens on the neighbouring carbons in cyclopropane is about 2.5 $\AA$ (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 C3H6 and C2H4 curves in Fig. <ref>. 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 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 [PDF].[To obtain this quantity from the experimental measurements one would of course need to take into account the fraction of molecules excited by the pump pulse.]

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 $\Delta\mathrm{PDF}(r)$ signal at 500 fs one can additionally make use of the lineshape to help distinguish the reaction products present. Figure <ref> 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_1$ and $f_2$ of products I and products II respectively, with the remaining contribution from products III ($f_1+f_2+f_3=1$). [Note the individual signals that are being averaged can be seen in Figs. S25-S27 in 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\lesssim f_1 \lesssim 0.85$. Note 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\lesssim f_1 \lesssim 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 lineshape at 500 fs alone could be used to distinguish different reaction products. Figure S28 shows the possible signals at 500 fs normalised 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_1$, of reaction products I, is the depth of the shoulder at about 1.8 $\AA$. Although the changes are relatively subtle compared to the absolute differences of Fig. <ref>, we see that, with a small enough relative error, it should nevertheless be possible to distinguish between different fractions, $f_1$, 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.

§.§ Potential systematic errors

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:

* electronic-structure theory

* initial conditions/treatment of the pulse

* dynamics (including zero-point energy leakage)

* 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 therefore 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 dynamical correlation. Along the path that leads from the Franck-Condon geometry to the 2/1 MECI and the 1 minimum, CASPT2 energies are very similar to CASSCF. This is also observed along the path connecting the 1 minimum to the first 1/0 MECI. However, the barriers to access the second and third 1/0 MECIs increase by 0.7 eV and 0.8 eV, respectively, at the CASPT2 level of theory. Since, the third 1/0 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 1 window) implies a strong wavelength dependence when exciting to 2, then 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 1 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 to 4.6 for aug-cc-pVDZ vs. 0.75 to 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 0 to 2 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 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,<cit.> 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 supplementary material) is difficult to distinguish by eye from the signal in Fig. <ref> or the 6-31+G* result (also in 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. <ref> and additional figures in 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.<cit.> These measurements were used to estimate the timescale of 2 decay, and found a timescale of about 740 fs. If correct, this would indicate that the timescale for 2 decay obtained in the current study is too fast. We note that our scans between the Franck–Condon point and the 2/1 MECI (shown in supporting information Figs. S12 and S13) 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 2 decay.

§ CONCLUSION

In conclusion, we have used molecular simulations to investigate the dynamics of cyclobutanone after photo-excitation to the 2 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. Firstly, 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.<cit.> 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.<cit.> Secondly, 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.<cit.> 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.<cit.>

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. Firstly, 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 an forthcoming publication.<cit.>  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 there would be time for sufficient population transfer to the triplet manifold to significantly influence the dynamics.
Secondly, 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 photo-excitation.<cit.> The Franck–Condon approach we have used, taking the initial distribution as the ground state nuclear Wigner distribution placed on the 2 state and mimicking the finite width of the exciting laser pulse by convolution with a Gaussian, is a standard way of initialising a simulation in photo-excited systems,<cit.> 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.<cit.>

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.<cit.> 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).<cit.> There has therefore been significant interest in the development of a nonadiabatic version of RPMD.<cit.> 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,<cit.> 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.<cit.> 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 initialising 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 thermalise to the correct distribution.<cit.>

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 minimises 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.<cit.>
Although this would significantly increase the cost of an on-the-fly simulation, by  exploiting modern machine-learning techniques<cit.> 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.

§ SUPPLEMENTARY MATERIAL

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

§ ACKNOWLEDGEMENTS

JEL was supported by an Independent Postdoctoral Fellowship at the Simons Center for Computational Physical Chemistry, under a grant from the Simons Foundation (839534, MT).
JRM and AK are supported by the Cluster of Excellence “CUI: Advanced Imaging of Matter” of the Deutsche Forschungsgemeinschaft (DFG) – EXC 2056 – project ID 390715994.

§ AUTHOR DECLARATIONS

§.§ Conflict of interest

The authors have no conflicts to disclose.

§.§ Author Contributions

Please add/modify the description of your contribution as appropriate

I would say let's include a detailed description in addition to using the somewhat limited Credit taxonomy:
Data curation
Formal analysis
Funding acquisition
Project administration
Writing - original draft
Writing - review & editing
(lead, equal or supporting)

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), conceptualisation (lead), project management (lead).

§ DATA AVAILABILITY

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