On the Theoretical Determination of Photolysis Properties for Atmospheric Volatile Organic Compounds.

Volatile organic compounds (VOCs) are ubiquitous atmospheric molecules that generate a complex network of chemical reactions in the troposphere, often triggered by the absorption of sunlight. Understanding the VOC composition of the atmosphere relies on our ability to characterize all of their possible reaction pathways. When considering reactions of (transient) VOCs with sunlight, the availability of photolysis rate constants, utilized in general atmospheric models, is often out of experimental reach due to the unstable nature of these molecules. Here, we show how recent advances in computational photochemistry allow us to calculate in silico the different ingredients of a photolysis rate constant, namely, the photoabsorption cross-section and wavelength-dependent quantum yields. The rich photochemistry of tert-butyl hydroperoxide, for which experimental data are available, is employed to test our protocol and highlight the strengths and weaknesses of different levels of electronic structure and nonadiabatic molecular dynamics to study the photochemistry of (transient) VOCs.

Volatile organic compounds (VOCs) are well-known atmospheric compounds that lead to the formation of secondary pollutants, among others, and are important contributors to global warming.[1] Interestingly though, their direct interaction with sunlight is yet to be fully explored. What is known is that VOCs undergo complex networks of chemical reactions, which are typically described by ground-state rate theories.[2−4] Upon chemical reactions, primary VOCs—for example, methane, isoprene, and toluene—will be converted into a large amount of secondary, often transient VOCs. Currently, the scientific models encompassing entangled networks of VOC chemical reactions, such as the master chemical mechanism (MCM), mostly rely on the ground-state reactivity of these molecules.[5−8] Models like the MCM are used to predict the accurate composition of the atmosphere and are therefore used as a guideline for new environmental policies.

However, sunlight absorption could promote a molecule in one of its electronically excited states, potentially opening new pathways for their chemical reactivity. This observation is even more critical when one considers that numerous transient VOCs bear one or more chromophoric groups. Hence, underestimating the role of a given type of chemical reactions, like excited-state processes, may significantly diminish the predictive power of chemical mechanism models. An example of such effect is the role of photolysis in the chemical reactivity of C5-HPALD (C5-hydroperoxy aldehyde), an isoprene oxidation intermediate, which has been proposed to partly explain the order of magnitude difference between the modeled and measured OH concentrations in isoprene rich atmospheric regions.[9,10] One of the reasons behind the oversight of the VOC excited-state phenomena (e.g., photolysis) is that the reliable experiments on excited VOCs, and in particular on short-lived transient VOCs, are difficult if not impossible to conduct. As a result, experiments often rely on proxy molecules; see, for example, the case of C6-HPALD as a proxy molecule for the C5-HPALD mentioned above.[11] Because there is an urgent need to provide chemical mechanism models with pathways involving photochemical reactions, computational and theoretical chemistry may be an elegant way to fill in this gap. Here, we demonstrate that state-of-the-art quantum chemical methods and excited-state dynamics techniques are mature enough to fulfill this task and propose a reliable protocol for computing the different ingredients that compose a photolysis rate constant.

The photolysis rate constant J, a first-order decay constant, for a given process is defined as an integral over wavelength λwhere σ(λ) is the photoabsorption cross-section of a molecule that undergoes photolysis, ϕ(λ) is the photolysis quantum yield for the process of interest, and F(λ) is the solar actinic flux, or in general the flux of the light source.a In words, to know how rapidly the molecule can produce a given product from photolysis (J), one needs to gain access to the probability that the molecule will absorb a photon at a given wavelength (σ), given the availability of a certain number of photons at a given wavelength (F), and the yield of a given photolysis process once the photon at a given wavelength has been absorbed (ϕ); this process has to be scanned over the entire window of relevant wavelengths (integration). Let us stress that J corresponds to a photolysis lifetime and should not be confused with the photochemical lifetime, that is, the time for a photoexcited molecule to decay toward the ground electronic state—with or without a photolysis process—and usually ranging within a femto- to picosecond time scale. If different photolysis pathways are possible for a given VOC, each dissociation channel is assigned a different ϕ(λ), and therefore a different J. ϕ(λ) and σ(λ) values for several atmospheric molecules can be found in the literature[12] and online databases.[13] Nevertheless, data are not available for many transient VOCs. Commonly, absorption cross-sections have been deduced from chemically similar compounds, using structure–activity relationships, while photolysis quantum yields are often considered to be independent of λ, i.e., ϕ(λ) ≈ ϕconst, and estimated on the basis of proxy molecules.

From the eyes of computational chemistry, the equation for photolysis rate constants (eq ) appears as an exciting playground. Over the past several years, computational photochemistry has undergone significant major developments that paved the way for an in silico evaluation of photoabsorption cross-sections σ(λ) and quantum yields ϕ(λ) from first principles.

A common theoretical approach to calculating photoabsorption spectra relies on the computation of Franck–Condon (FC) factors, along with Herzberg–Teller corrections.[14−16] This method is convenient when a small number of excited states are of interest, while the molecule possesses a single dominant conformer (see ref (17) for examples of additional strategies in this direction). Alternatively, photoabsorption spectra can be evaluated directly from the nuclear wavepacket dynamics autocorrelation function, for which various quantum[18] and semiclassical[19] simulations have been devised. Still, the latter is mainly applicable to smaller molecular systems because of its computational cost for medium-sized molecules like (transient) VOCs. Finally, the nuclear ensemble approach (NEA)[20] is a conceptually simple approach in which a photoabsorption cross-section is evaluated from vertical transitions calculated on the support of an ensemble of initial geometries. Such initial geometries are sampled from a ground-state nuclear density, either by performing long classical molecular dynamics, by using path-integral molecular dynamics, or by sampling an approximate Wigner distribution that accounts for the effect of zero point energy.[21,22] In principle, the NEA method is appropriate for the treatment of molecular systems of various sizes and with different conformers. Unlike some of the methods mentioned above, the NEA method does not yield a detailed vibronic structure of the computed photoabsorption cross-section, but it typically provides correct band shapes for multiple electronically excited states, as long as appropriate electronic structure methods are employed.[20]

On the other hand, nonadiabatic molecular dynamics methods, in which nuclear and electronic motions are coupled beyond the ubiquitous Born–Oppenheimer approximation,[54] may give a direct access to photolysis quantum yields. A plethora of methods has been developed over the decades, from wavepacket-based quantum dynamics methods[18] to trajectory-based semiclassical and mixed quantum-classical approaches.[23,24] Trajectory surface hopping (TSH), using the fewest-switches algorithm devised by Tully,[25] is a particularly efficient mixed quantum/classical method in which a swarm of classical nuclear trajectories is propagated in a set of electronic states, with the possibility of hopping between different electronic states as a result of nonadiabatic effects. Results from TSH simulations, which are efficient in practice but approximate in theory, should be justified by comparison with results of experiments and/or higher-level nonadiabatic dynamics methods such as ab initio multiple spawning.[26,27]

All methods described above also strongly rely on the quality of the underlying electronic structure method employed to compute electronic energies and other quantities required for nonadiabatic dynamics like nuclear forces or nonadiabatic couplings. For the specific case of photolysis processes, the electronic structure method should not only be tested and validated in the Franck–Condon region but also be stable and reliable until the dissociation limit. Hence, the application of computational photochemistry to atmospheric molecules is expected to challenge numerous methods commonly employed to characterize molecular excited-state properties in the Franck–Condon region, for example, LR-TDDFT (linear-response time-dependent density functional theory) or ADC(2) (algebraic diagrammatic construction up to second order), and also those used for excited-state dynamics, SA-CASSCF (state-averaged complete active space self-consistent field) or MS-CASPT2 (multistate complete active space second-order perturbation theory).

The applications of some of the advanced techniques mentioned above in the atmospheric photochemistry of VOCs are still scarce[28−30] (see also ref (31) for a recent example of mercury-based compounds). More importantly, these strategies have not been combined to determine the components of the photolysis rate constant in a coherent manner. In the following, we propose a protocol to study in silico the photolysis of (transient) VOC molecules and exemplify its performance with the study of tert-butyl hydroperoxide [tBHP (Scheme )], for which we can compare our results with available experimental data.[32] We offer a thorough benchmark of the electronic structure methods available to describe the photochemical processes of tBHP, within and beyond the FC region, and we justify our choice for the excited-state dynamics strategy employed. We show that combining NEA and TSH with XMS-CASPT2 allows for a reliable determination of (i) the photoabsorption cross-section and (ii) the wavelength-dependent quantum yields for the H, OH, and O photolysis channels of tBHP, the two key components to form photolysis rate constants given a photon flux.

Scheme 1

Proposed tBHP Photolysis Channels Based on Reference (32)

The tBHP molecule is a member of the alkyl hydroperoxide family, a family of ubiquitous atmospheric VOCs that also play an important role in combustion chemistry.[33] Different tBHP photolysis channels were proposed (Scheme ); however, upon the pulsed laser excitation at the low-energy tail of the absorption spectrum, only OH photolysis was measured. For the purpose of the analysis presented here, peroxide systems represent a stringent test due to their challenging electronic structure (see the Supporting Information), while their photolysis may involve the interplay of multiple dissociation channels depending on the excitation wavelength.

A crucial step for obtaining reliable σ(λ) and ϕ(λ) quantities is to carefully benchmark the level of electronic structure theory, not only in the Franck–Condon region but also beyond. For this task, we consider here LR-TDDFT using the PBE0 functional, the single-reference electronic structure method SCS-ADC(2), and the multireference XMS-CASPT2 method (see Computational Details and the Supporting Information for additional information). For the latter method, we devised an active space containing two nonbonding n orbitals as well as three pairs of bonding and antibonding σ orbitals to allow for all of the possible cleavages of O–O, O–H, and C–O bonds (orbitals selected for the active space are provided in the Figure S3).

Excitation energies obtained for all three methods are compared to those from EOM-CCSD (converged with the large basis set aug-cc-pVTZ) in Table , using a common ground-state optimized geometry (MP2/aug-cc-pVDZ). Note that peroxide bonds are notoriously hard to describe from a theoretical perspective;[34] additional details about the ground-state geometry can be found in the Supporting Information. EOM-CCSD involves a better treatment of electron correlation effects as compared to LR-TDDFT and SCS-ADC(2), while it does not depend on flexible parameters such as an active space like XMS-CASPT2. For each method, we opted for the basis sets that are sufficiently large to match the EOM-CCSD reference. We note that other high-level methods such as CC3 are even more suited as a theoretical reference,[35] if computations with large basis sets are affordable for the system of interest. Additional electronic structure methods and basis sets have been tested, and the corresponding results are listed in Table S2.

Table 1

tBHP Vertical Singlet Transitions at the MP2/aug-cc-pVDZ Ground-State Optimized Geometrya

electronic structure	n′σOO*	nσOO*	n′σOH* + Ryd
LR-TDDFT/TDA/PBE0b/def2-TZVPD	5.57	6.66	6.69
SCS-ADC(2)/def2-TZVPD	5.91	6.99	6.90
XMS(4)-CASPT2(10/8)/def2-SVPD	5.76	7.19	6.97
EOM-CCSD/aug-cc-pVTZ	5.78	6.87	6.89

Energies are in electronvolts.

LR-TDDFT within the Tamm–Dancoff approximation (TDA).

The low-lying singlet excited states of tBHP all exhibit a nσ* character, as predicted by all four methods (Table ). The transition to the S1 state is of n′σOO* character, i.e., antibonding with respect to the O–O bond (see Figure , right). The S2 and S3 states have similar energies, with a dissociative nature with respect to the O–H and O–O bonds (see Figure , right). We note that the n′σOH* state has a substantial diffuse (i.e., Rydberg) character. As such, it is crucial to employ a basis set that contains diffuse orbitals to keep the target state in the low-energy region (see Table S2). While the valence/Rydberg mixing is sometimes considered an artifact of the approximate treatment of correlation effects in π-conjugated systems,[36,37] it is a common phenomenon in many small molecules and radicals.[38] A considerable degree of mixing is observed with all of the different levels of theory tested here, which supports our assignment as n′σOH* + Ryd. Interestingly, XMS-CASPT2 and SCS-ADC(2) place this n′σOH* + Ryd state at an energy that is slightly lower than that of the nσOO* state (see Figure , right), while LR-TDDFT/PBE0 and EOM-CCSD indicate that both states have nearly the same energy. It is important to point out that the ordering of these two electronic states is strongly sensitive to the reference geometry employed (see Tables S2 and S3).

Figure 1

Electronic structure of tBHP at and beyond the Franck–Condon region. Rigid scan along the O–O bond of tBHP starting from a ground-state optimized geometry at the MP2/aug-cc-pVDZ level of theory (left panel). Comparison of the electronic energies obtained with XMS(4)-CASPT2(10/8)/def2-SVPD (black), LR-TDDFT/TDA/PBE0/def2-TZVPD (palatinate), and SCS-ADC(2)/def2-TZVPD (gray) for the four lowest electronic states, S0 (squares), S1 (triangles), S2 (inverted triangles), and S3 (diamonds). Main orbital contributions characterizing the lowest excited singlet states based on XMS(4)-CASPT2(10/8)/def2-SVPD at the Franck–Condon point (right panel).

We also assessed the behavior of the different electronic structure methods along the dissociation profile of the O–O bond (left panel of Figure ). While XMS-CASPT2 exhibits a proper behavior with a clear dissociation limit, LR-TDDFT and SCS-ADC(2) fail to reproduce the proper topography of the resulting potential energy curves as soon as the molecule is distorted away from the Franck–Condon region. We further note that the S3 state exhibits a strong double-excitation character after 1.7 Å at the XMS-CASPT2 level of theory, a type of electronic state that both LR-TDDFT and SCS-ADC(2) are unlikely to reproduce accurately.[39,40] The inadequacy of the latter two methods, both considered single-reference methods, is also reflected by the increase in the extent of multireference character along the O–O bond cleavage. The D1 diagnostic, a measure of the ground-state multireference character, quickly soars beyond the recommended limit of 0.04[41] for an O–O bond length after 1.8 Å and reaches values of >0.3 (Figure S2). A scan along the O–H bond is proposed in the Supporting Information (Figure S1) and shows similar trends between the different methods as observed here for the O–O bond and agrees with earlier work on MeOOH.[34] On the basis of all of the results presented in this section, it is clear that an adequate description of the photolysis pathways for tBHP requires the use of XMS-CASPT2 (combined with a def2-SVPD basis set).

Having determined an adequate level of theory for the electronic structure of tBHP, we can now proceed with the calculation of the first quantity of interest, the photoabsorption cross-section σ(λ). Figure compares the σ(λ) computed with the NEA (see Computational Details and the Supporting Information for more information) to the ultraviolet cross-section values measured at several discrete wavelengths[32] relative to the infrared cross-sections. The dashed line through the experimental points serves as a guide for the eye. Note, however, that the authors of ref (32) measured σ(λ) above 210 nm via diode array spectroscopy and scaled their spectrum to optimize the agreement with the discrete σ values obtained at 228 and 254 nm.

Figure 2

Calculated and experimental photoabsorption cross-sections for tBHP. The inset shows the 500 molecular geometries sampled from a harmonic Wigner distribution (MP2/aug-cc-pVDZ) that were used to calculate the photoabsorption cross-section, while numerated energy windows are indicated in different colors. Experimental data were obtained from ref (32). This photoabsorption cross-section incorporates only transitions toward the first three (adiabatic) electronic excited states.

The theoretical photoabsorption cross-section obtained with XMS-CASPT2 is in excellent agreement with the available experimental data, correctly predicting the offset and the increase in the absorption intensity. It is important to note that the photoabsorption spectra presented here incorporate only transitions toward the three lowest excited electronic states; as such, progressing toward shorter wavelengths would require the inclusion of more excited electronic states. An overlap of all of the 500 nuclear configurations sampled within the NEA is provided in the inset of Figure and shows the diversity of geometries sampled with this strategy. The NEA spectra computed with LR-TDDFT/TDA/PBE0 and SCS-ADC(2), on the same set of 500 nuclear configurations, agree less well with experiment (see Figure S4); the cross-sections are overestimated as a result of the shift in excitation energies and larger oscillator strengths. Overall, the NEA combined with XMS-CASPT2 allows for an accurate simulation of the tBHP photoabsorption cross-section.

On the basis of the calculated σ(λ), we can now progress toward the photodissociation of tBHP. In the following, we show how we can use the nonadiabatic molecular dynamics method TSH to unravel the wavelength-dependent quantum yields associated with the photochemistry of tBHP. “Wavelength-dependent” here means that the molecule shall be excited at different wavelengths, and following its nonradiative decay toward the ground electronic state, one can determine the quantum yield for the different photolysis products emerging from each excitation wavelength. It is crucial to note at this stage that the wavelength dependence of photolysis quantum yields is often neglected for transient VOCs, due to the challenges associated with determining such quantities at different excitation energies. The idea of the computational protocol presented here is to offer a strategy for recovering this critical information.

Let us start by defining different energy windows on the calculated photoabsorption cross-section of tBHP (see color-coded regions in Figure ), centered around 5.0 eV (window 1), 5.5 eV (2), 6.0 eV (3), 6.5 eV (4), 7.0 eV (5), and 7.5 eV (6).b It is clear that the number of energy windows chosen here and their width (0.3 eV) are arbitrary; the choice is essentially guided by the balance between the accuracy desired in the calculation of the wavelength dependence of ϕ(λ) and the overall computational efficiency of the protocol. The sampling can be systematically improved by considering more initial conditions from the Wigner distribution, which allows setting narrower windows. However, the results appear to be rather insensitive to the precise parameters of sampling windows when the excited electronic states considered are mostly dissociative in nature and, as such, do not possess a clear vibronic structure (as is the case here with tBHP).c

Upon photoexcitation, tBHP exhibits a surprisingly rich photochemistry. The swarm of trajectories initiated from the energy windows described above and propagated using the TSH method shows that the system undergoes ultrafast photolysis, typically in <30 fs. Interestingly, the photolysis of tBHP follows very different routes depending on the excitation window, and Figure depicts this diversity by showing how specific bonds of tBHP evolve over time along >500 TSH trajectories. All trajectories (lines) start from the localized distribution around the equilibrium geometry (black dots at the bottom left corner of Figure ) and undergo bond cleavage. However, the fate of the system during the excited-state dynamics strongly depends on the initial excitation. For the low-energy windows (gray, blue, and green lines in Figure ), OH dissociation is a dominant photolysis channel; the trajectories all show a ballistic stretch of the O–OH bond (horizontal axis in Figure ). This observation is fully consistent with the n′σOO* character observed earlier for the lowest excited state. Going toward higher energies (i.e., shorter wavelengths), the H abstraction channel starts to gain importance, as shown in Figure with more trajectories starting to evolve along the O–H coordinate (vertical axis). This change in the photolysis mechanism is caused by the higher excitation energy, reaching the state with a strong n′σOH* character. When trying to probe even higher energy windows, we find less common mechanisms (see the diagonal of Figure ). Among these mechanisms, one observes the concerted dissociation of both H and O radicals, as well as oxygen dissociation leading to the formation of tert-butyl hydroxide as a photoproduct. Interestingly, this O dissociation occurs along with the roaming motion of the H atom of the peroxo moiety (see the Movie in the Supporting Information). Roaming processes are well-known in ground-state chemistry as they allow one to bypass transition-state theory[45] and have occasionally been discussed in the context of excited-state dynamics.[46] Importantly, O dissociation has been observed for the parent molecule CH3OOH[47] (see below). The Supporting Information provides the interested reader with a representative trajectory for all of the processes described in this paragraph, highlighting the interplay between electronic states as well as the strict conservation of total energy in the TSH dynamics.

Figure 3

Analysis of the photochemical processes observed during the excited-state dynamics of tBHP. Projection of the 541 (Table S4) trajectories on two critical bond lengths, O–O and O–H (see the structure in the inset for an identification of these coordinates). Each TSH trajectory is color-coded according to its initial excitation window (windows are highlighted in Figure ). Inserted molecular geometries represent snapshots of the TSH dynamics highlighting four possible mechanisms H, O + H, O, and OH photodissociation (left to right, respectively).

One can finally deduce the in silico wavelength-dependent quantum yield for each of the different photolysis channels observed for tBHP by compiling the results of the TSH dynamics for the selected energy windows (Figure ). For the lowest-energy window, centered at around 5.0 eV (248 nm), the computed OH quantum yield amounts to (0.95 ± 0.03), whereas the H quantum yield equals (0.05 ± 0.03). Baasandorj et al.[32] measured the tBHP ϕOH at 248 nm with pulsed laser photolysis and obtained values of (1.04 ± 0.07) relative to the ϕOH of HNO3 and (1.11 ± 0.10) relative to the ϕOH of H2O2. As for a molecule like tBHP, ϕOH can hardly exceed 1, and the authors conclude that the actual value is close to unity; the quantum yields for the other photolysis channels were not investigated but were expected to be minor. However, in an earlier study[47] of the photolysis of a parent molecule, CH3OOH, the OH, H, and O quantum yields at 248 nm were found to be (1.00 ± 0.18), (0.038 ± 0.007), and <0.007, respectively. A small amount of H fragmentation at 248 nm indicates the proximity of the O–H antibonding excited state, which has modest but non-negligible intensity in the low-energy spectrum, and is in agreement with our simulations. As a further matter of clarification, we note that the energy bands of dissociative states are intrinsically spread due to the nuclear motion present even at the ground-state zero point vibrational energy. Bond stretching can substantially alter the vertical transition energies of an antibonding state as compared to the transitions at the ground-state minimum (Table ). It is therefore key when simulating the photodynamics of such molecules to account for such non-Condon effects, as they are likely to alter the ordering of the electronic-state characters (or their degree of mixing). Overall, the simulated quantum yields at 248 nm are in fairly good agreement with experimental evidence. More importantly, the simulation protocol described here allows us to straightforwardly compute and predict the wavelength dependence of ϕ(λ) for tBHP. Despite the fact that the low-wavelength region may be out of the atmospherically relevant “actinic range”, ϕ(λ) was calculated here as a proof of principle. In the higher-energy part of the spectrum, H dissociation starts to gain importance, with ϕH ≈ ϕOH at 165 nm. This can be attributed to the higher-lying nσOO* and n′σOH* antibonding excited states, both being close in energy. Below 180 nm, O atoms start to appear, due to O + H dissociation but also due to O dissociation (with formation of tert-butyl hydroxide). These mechanisms occur as a result of the interplay between the nσOO* and nσOH* states: the system undergoes transitions through the manifold of excited states, which breaks both O–H and O–O bonds. Occasionally, the leaving H atom may be captured in the orbit of the remaining O atom on the parent molecule, leading to the formation of tert-butyl hydroxide. Finally, the release of OOH fragments (Scheme ) was not observed in the investigated wavelength range.

Figure 4

Calculated wavelength-dependent photolysis quantum yields for tBHP. Experimental data from refs (32) and (47) are shown in the inset.

To summarize, we propose here a general protocol for evaluating the key constituents of the photolysis rate constant J: the photoabsorption cross-section and the wavelength-dependent photolysis quantum yields. We show how one can evaluate the sought-for ingredients—and, consequently, the corresponding photolysis rate constants for any given functional form of the photon flux—by combining a state-of-the-art electronic structure method (XMS-CASPT2) and nonadiabatic molecular dynamics (trajectory surface hopping). We benchmark our protocol with the photolysis of tBHP, a compound for which reliable experimental data are available. We carefully compare the different strategies available for the level of electronic structure theory, as well as the excited-state dynamics, highlighting the shortcomings of commonly employed (single-reference) methods in offering a proper description of the potential energy surfaces beyond the Franck–Condon region for the different photolysis pathways suffered by tBHP. The results obtained for tBHP are in solid agreement with available experimental data. We note that the strategy described in this work can straightforwardly be extended to account for additional chemical reactivity coming from hot ground-state dynamics following the nonradiative decay, or reactions involving intersystem crossing processes. Hence, the protocol proposed here paves the way for the prediction of photolysis rate constants for more complex, multichromophoric transient VOCs, for which experimental data are currently unavailable.

Computational Details

A detailed Computational Details section is available in the Supporting Information. We summarize here only the most important information.

All XMS-CASPT2[48] (extended multistate complete active space second-order perturbation theory) computations were conducted with the BAGEL 1.2.0[49] program package. A single-state single-reference (SS-SR) contraction scheme[50] was used, while the vertical shift was set to 0.5 au to avoid problems with intruder states.

Photoabsorption cross-sections were calculated with the NEA[20] based on a harmonic Wigner distribution; we assume that the molecule is in the ground vibrational state for all vibrational modes. Five hundred initial conditions were sampled, and the electronic transitions (three lowest excited singlet states) were computed with different electronic structure methods (Figure S4) and broadened with Lorentzians using a phenomenological broadening of 0.05 eV.

TSH[25] nonadiabatic dynamics simulations were performed with the SHARC 2.1 code[51,52] coupled to BAGEL. Four electronic states were considered. Initial conditions were randomly sampled from the Wigner distribution described above, based on narrow energy windows in the photoabsorption cross-section. The standard deviations of photolysis quantum yields were estimated following the method of Persico and Granucci.[53] TSH dynamics was also validated by comparisons with ab initio multiple spawning (Supporting Information).