Quantum Quality with Classical Cost: Ab Initio Nonadiabatic Dynamics Simulations Using the Mapping Approach to Surface Hopping

Nonadiabatic dynamics methods are an essential tool for investigating photochemical processes. In the context of employing first-principles electronic structure techniques, such simulations can be carried out in a practical manner using semiclassical trajectory-based methods or wave packet approaches. While all approaches applicable to first-principles simulations are necessarily approximate, it is commonly thought that wave packet approaches offer inherent advantages over their semiclassical counterparts in terms of accuracy and that this trait simply comes at a higher computational cost. Here we demonstrate that the mapping approach to surface hopping (MASH), a recently introduced trajectory-based nonadiabatic dynamics method, can be efficiently applied in tandem with ab initio electronic structure. Our results even suggest that MASH may provide more accurate results than on-the-fly wave packet techniques, all at a much lower computational cost.

Accurate computational simulations are crucial for understanding and interpreting experiments that investigate photoexcited molecular processes.Such systems are often high dimensional, involving multiple electronic potential energy surfaces and reaction channels, which makes constructing accurate parameterized models extremely challenging and time-consuming.Hence, performing ab initio simulations 'on-the-fly', such that the surfaces are computed in tandem with the propagation of the nuclear degrees of freedom, is often the only feasible option.In order to make such simulations practical, the number of electronic-structure calls per time step must be kept minimal, which requires dynamical approaches based on trajectories or localized nuclear basis functions.
Most dynamics approaches have not been specifically tested for ab initio simulations.2][3][4] However it is often unclear how far the conclusions from these simplified models can be extended to more realistic systems.Recently a test set of photoexcited molecular systems have been proposed as a benchmark for investigating the properties of different nonadiabatic dynamics algorithms using onthe-fly ab initio electronic structure. 5These three molecules -ethylene, 4-(N,N-dimethylamino)benzonitrile (DMABN), and fulvene -were initially chosen because their dynamics show similarities to the three so-called Tully models, 6 but their utility for benchmarking goes far beyond that.In particular, it is known that these systems give rise to different conical intersection topographies 5 and pathways for approaching the intersections, 7 therefore providing a rigorous test for any nonadiabatic dynamics method.This test set is currently gaining interest from the community, and it has already been the focus of a number of classical trajectory 5,8 and Gaussian wavepacket 5,7 studies.In this work, we ascertain the accuracy and utility of a novel trajectory-based dynamics technique for performing ab initio simulations, the mapping approach to surface hopping (MASH), 9 by applying it to simulate this molecular photochemistry test set and comparing our results with other more established methods.In order to determine a hierarchy in the accuracy of these approaches, we also compare with numerically exact quantum dynamics applied to linear vibronic coupling (LVC) models previously parameterized for these molecules. 7efore discussing the relaxation dynamics of the three photoexcited molecular systems, we give a brief overview of the dynamics approaches used in this work, with a particular emphasis on their similarities and differences.More information can be found in Refs.9-12.
Fewest-switches surface hopping (FSSH) 6,10 is the most popular independent-trajectory approach for simulating nonadiabatic dynamics in molecules.In FSSH, the nuclei are propagated according to (classical) Born-Oppenheimer molecular dynamics on a single surface and nonadiabatic transitions are described by stochastic changes in the 'active' surface, called 'hops'.The hopping probability is related to the rate of change of the underlying electronic wavefunction, which itself is propagated according to the associated time-dependent Schrödinger equation.One issue with altering the active surface in this way is that it can become inconsistent with the electronic wavefunction, leading to the so-called overcoherence error that is known to significantly degrade the accuracy of the obtained results.To fix this, a number of decoherence corrections have been proposed, [13][14][15][16][17][18][19][20][21] which sporadically reset the wavefunction to the current ac-tive surface.While not guaranteed, 22 it is generally accepted that decoherence corrections lead to an improvement in the accuracy of the obtained results in the majority of cases.
Despite the substantial progress that has been made in understanding many foundational aspects of FSSH, 23,24 a number of variants of the FSSH dynamics algorithm are nevertheless still widely used in the literature.The main aspect that differs between most FSSH algorithms lies in the way that the nuclear velocities are rescaled at a hop.While it is generally agreed upon that rescaling along the nonadiabatic coupling vector (NACV) is the correct thing to do, 13,[25][26][27] many other schemes are used in practical implementations of FSSH. 279][30] Additionally, an upward hop must be aborted if there is insufficient nuclear kinetic energy, often referred to as a 'frustrated hop'.The nuclear velocity along the NACV is reflected at a frustrated hop in many FSSH implementations, 13 although other suggestions have been made, [31][32][33] along with those that try to avoid frustrated hops altogether. 34,35he mapping approach to surface hopping (MASH) 9 is a recently proposed independent-trajectory approach that alleviates the problems of FSSH by utilizing the best features of mapping-based semiclassical trajectories 2,3,36,37 in a surface hopping algorithm.In many aspects the algorithm is identical to FSSH, but it contains the following key differences.In MASH, the active surface is not an additional parameter within the theory, but is uniquely determined from the electronic wavefunction.For two-state problems, this corresponds to selecting the surface for which the electronic wavefunction has the largest associated probability.In addition, the stochastic nature of hops in FSSH is replaced by a fully deterministic dynamics which guarantees that the electronic wavefunction and the active surface remain consistent.The overcoherence problem is therefore resolved in MASH without the need for ad hoc decoherence corrections.0][41][42] In fact, MASH gives a unique prescription for all aspects of the simulation algorithm, including the velocity rescaling and treatment of frustrated hops.In agreement with what many have suggested for FSSH, 13,[25][26][27] MASH determines that the velocity must be rescaled along the NACV at a hop and reflected in the case of a frustrated hop, in order for the approach to reproduce the short-time behaviour of exact quantum dynamics.In particular, the exact quantum-mechanical equation of motion for the so-called 'kinematic momentum' 43 for mode j depends on the Born-Oppenheimer and nonadiabatic contributions to the nuclear force, which are given by the two terms on the right-hand side of the following equation In addition, MASH has already been benchmarked on a range of model systems, 38,44,45 where it was regularly found to offer improvements over FSSH.
To summarize the above discussion, a hierarchy of the surface hopping approaches can be established according to their expected accuracy.Firstly, it is expected that surface hopping approaches that perform the velocity rescaling along the velocity vector (FSSH-vel) will be less accurate than those that perform the velocity rescaling along the NACV (FSSH-nacv), due to the fact that the later can describe dynamical effects arising from the nonadiabatic force.Secondly, fewest-switches surface hopping approaches that incorporate a decoherence correction (dFSSH) are expected to be more accurate than those that do not (FSSH), because it is important to impose consistency between the active surface and the electronic wavefunction.Finally, MASH is expected to be either just as accurate or more accurate than dFSSH-nacv.This hierarchy will be helpful for determining what is likely to be the correct dynamics in ab initio photochemical simulations, where numerically exact quantum dynamics is not obtainable.
A completely different approach for describing nonadiabatic transitions compared to the quantum-classical approaches described above is taken in ab initio multiple spawning (AIMS). 11,12Based on a series of approximations to full multiple spawning, 46,47 AIMS is a Gaussian wavepacket algorithm that was developed for performing onthe-fly ab initio simulations.Gaussian basis functions are propagated classically on single Born-Oppenheimer surfaces and new Gaussians are 'spawned' whenever the system enters a region of strong nonadiabatic coupling.The evolution of the Gaussian weights is then obtained by solving an approximate time-dependent Schrödinger equation spanned by the Gaussian basis.The main advantage of AIMS over the surface hopping approaches is its coupled trajectory nature, which goes beyond the independent trajectory approximation, albeit at an increased computational cost.On the other hand, given that AIMS only uses a minimal number of Gaussian basis functions, its accuracy will largely be determined by how effectively the Born-Oppenheimer forces are able to move the Gaussians into the correct regions of nuclear phase space.As a result, one way of improving upon AIMS is to instead propagate the Gaussians using equations of motion that ensure Eq. ( 1) is satisfied, as is done for example in the variational multi-configuration Gaussian (vMCG) approach. 48Given that AIMS is not a benchmark method, 49 the relative accuracy of AIMS compared to the quantumclassical approaches therefore depends on the relative severity of the independent trajectory approximation for realistic molecular simulations compared to how effectively the minimal set of Gaussian basis functions generated by AIMS spans the full support of the time-dependent wavefunction.
We now consider in more detail the relaxation dynamics of ethylene, DMABN, and fulvene, using the same electronic structure theory for each system as defined in Ref. 5. The initial conditions are taken to be of the Franck-Condon type, with the electronic system on the upper of the two considered adiabats and the nuclear system in the ground vibrational state associated with the harmonic approximation to the electronic ground state potential.MOLPRO 2012 50 and GAUSSIAN 16 51 were used for the SA-CASSCF and LR-TDDFT electronic structure calculations, while the surface hopping and AIMS dynamics were performed using the SHARC 2.0 [52][53][54] and AIMS/MOLPRO 55 codes.Files con-Figure 1: The magnitude of the effective nonadiabatic coupling, | j djvj|, and the probability distribution of the timedependent energy gap, ∆V = V+ − V−, between the two adiabatic states.These quantities are calculated for each system by averaging over the MASH trajectories.
taining the electronic structure inputs and the ground-state geometries and frequencies for each system are provided in the supplementary material.Each system displays a different type of relaxation pathway involving conical intersections.One way that this excited-state relaxation dynamics can be visualized is through the probability density of the dynamical energy gap between the two Born-Oppenheimer surfaces, as presented in the second row of Fig. 1 for the MASH trajectories.In these figures, the Franck-Condon region corresponds to a large value of ∆V , while the conical intersection seams are about ∆V ≈ 0. Also given in the first row of the same figure is the effective nonadiabatic coupling between the states (i.e., the magnitude of the dot product of the NACV with the nuclear velocity) averaged over the same trajectories.
Ethylene is known to have 'indirect' access to its conical intersections, 7 meaning that the Franck-Condon region for the S0 → S1 transition lies far away from the conical intersection seams and the initial nuclear dynamics upon photoexcitation does not directly access them.A redistribution of the vibrational energy to the appropriate modes during a few vibrational time periods is necessary before the crossing region can be accessed, giving a delayed onset of the first nonadiabatic transitions to about ≈ 25 fs.At later times, the majority of trajectories move away from the intersection region once they have reached the ground state.However the relatively slow decay in the coupling suggests that this process is relatively inefficient and that many trajectories may undergo multiple nonadiabatic transitions before remaining on the ground-state surface.
In contrast, upon a Franck-Condon type photoexcitation from the ground state to S2, DMABN has a very small energy gap between S2 and S1.As a result, the nonequilibrium dynamics in DMBAN were previously coined 'immediate' 7 because the initial nuclear wavepacket is essentially on top of the conical intersection.The small initial energy gap between S2 and S1 also means that the dynamics remain close to the conical intersection seam for relatively long times, leading to the possibility that a large number of nonadiabatic transitions take place.This is consistent with the observation that the effective nonadiabatic coupling rapidly plateaus to a non-zero value.Experimentally it is known that relaxation to the S0 state occurs at much longer times predominantly through fluorescence, and so we exclude the possibility of nonadiabatic transitions between S1 and S0 in this study.
Despite the conical intersection seams in fulvene being further away from the Frank-Condon region than in DMABN, the initial motion of the nuclear wavepacket does still allow direct access to the crossing region, 7 in contrast to ethylene.The distinct feature of fulvene is that trajectories first pass through a sloped conical intersection seam driven by a stretch in the C=CH2 moiety, 5 and then part of the wavepacket is reflected back through the crossing region at ≈ 15 fs.This system therefore provides a useful test for how well different nonadiabatic dynamics approaches can correctly describe recrossing phenomena.While a peaked conical intersection seam also exists in fulvene, driven by a twist in the C=CH2 moiety, we find that this intersection is not accessed until much later times.
We next analyze the time evolution of the electronic excited state populations, shown in the first row of Fig. 2. In order to help analyze the differences in the obtained populations from different surface hopping algorithms, the number of allowed and frustrated hops are also given in the same figure.In particular, the adiabatic populations can be exactly reproduced from the difference in the number of downward and upward hops.Despite the differing properties of the conical intersections involved in these three systems, the general trend in the results for each method is largely the same.
First, we observe that the direction in which the velocity rescaling is performed in surface hopping approaches results in a significant quantitative difference in the obtained electronic populations.For FSSH performed by rescaling along the velocity vector, dFSSH-vel deviates significantly from FSSH-vel, suggesting that the decoherence error in this case is large.In contrast, for the surface hopping results where the velocity rescaling is applied along the NACV, FSSH-nacv and dFSSH-nacv are almost identical in all cases.This can be understood from the fact that in all systems, FSSH-vel trajectories undergo a larger number of allowed hops than FSSH-nacv, making it more likely that the electronic wavefunction can become inconsistent with the active surface in the former.Decoherence corrections are therefore more necessary in FSSH-vel for reimposing consistency.
These differences between the algorithms can be further explained by considering the amount of nuclear kinetic energy available for promoting upward hops.In the case of rescaling along the velocity vector, the nuclear kinetic energy of the entire molecule is available for inducing electronic transitions, making almost all attempted hops energetically allowed.However, this behaviour is unphysical because not all of the nuclear degrees of freedom are directly coupled to the electronic transition. 57,58In contrast, the NACV rescaling direction ensures that only the kinetic energy associated with directly coupled modes is considered; this is significantly less than the nuclear kinetic energy of the entire molecule, making upward hops more likely to be energetically forbidden and therefore frustrated.The associated electronic population for the upper adiabat is therefore significantly lower for dFSSH-nacv than for dFSSH-vel.This effect has been observed in other theoretical studies, 59 including those on ethylene 60 and fulvene. 5,27n particular, this leads to qualitatively different dynamics than was previously predicted for DMABN.One of the main reasons that DMABN was previously suggested as a good candidate benchmark system was that its dynamics were expected to produce similar features to Tully's model II. 5 It was known that the adiabatic potential energy surfaces remain close in energy throughout the dynamics (as also illustrated in Fig. 1), suggesting that repeated electronic transitions between the surfaces would occur.While this is indeed observed when rescaling along the velocity vector, the observed dynamics when rescaling along the NACV instead involves a single rapid transition to the lower adiabatic state, where the system remains indefinitely.While the potential energy surfaces in DMABN do remain close together in energy relative to the kinetic energy of the entire molecule, they do not relative to the kinetic energy along the NACV.
Of the three systems, fulvene is particularly interesting because the MASH result significantly deviates from both FSSH-nacv and dFSSH-nacv.There is also a noticeable difference between these results for DMABN too, although the difference is much smaller.Figure 2 shows that this deviation between the MASH and dFSSH-nacv electronic populations is predominately due to the larger number of frustrated hops in MASH, which leads to the MASH electronic populations being slightly lower than those of dFSSH-nacv.
What is also interesting about Fig. 2 is that the electronic populations obtained by AIMS are seen to be almost identical to those obtained with dFSSH-vel.Firstly, this suggests that the independent trajectory approximation that underpins all surface hopping approaches is valid for these systems.While there are small deviations between the dFSSHvel and AIMS results around 10 fs in fulvene and towards longer times in DMABN, we note that these differences are relatively minor compared to the more significant discrepancies observed between other algorithms.More importantly, the fact that AIMS and dFSSH-vel agree so well suggests that AIMS may not be describing the effect of the nonadiabatic force, which is an effect that is also neglected in dFSSH-vel.This finding is not so surprising in the case of DMABN, where the AIMS simulation did rescale the average velocity of spawned Gaussians along the velocity vector to ensure classical energy conservation. 56However in the AIMS simulations for ethylene and fulvene, this rescaling was performed along the NACV, which at least for surface hopping algorithms is sufficient to correctly describe the nonadiabatic force.Our findings are however consistent with the observation that, unlike for surface hopping approaches, the velocity rescaling direction in AIMS is found to make almost no difference to the obtained results. 61While this point certainly requires further investigation, it does however suggest that the AIMS population dynamics may be less accurate than the best surface hopping algorithms in these cases.
In previous work, other algorithms have also been tested on these systems using the same initial conditions and electronic structure methods.In Fig. S4 in the SI, we compare AIMS and MASH to various flavours of the symmetrical quasi-classical (SQC) 62 method for the population dynamics of ethylene.The SQC results were obtained from Ref. 8. With the exception of non-gamma corrected SQC using square windows, all of the other SQC approaches give results somewhere in between AIMS and MASH.Given that the SQC results contain large statistical error, it is however hard to precisely ascertain the relative accuracies of the various approaches.While the surface hopping algorithms were performed with slightly more trajectories (1000) than the SQC approaches (500), given the observed noise in the data, we estimate that at least an order or magnitude more SQC trajectories would be needed to obtain a similar level of convergence as the surface hopping results.This highlights one of the main advantages of the surface hopping approaches (including MASH) over mean-field mapping-based approaches, in that they require significantly fewer trajectories to converge the results.
To complement the above comparison of the various dynamics approaches in the ab initio case, we also perform simulations for linear vibronic coupling (LVC) models fit to the electronic structure data for DMABN and fulvene.These LVC models have already been used to perform numerically exact quantum dynamics using the multi-configuration timedependent Hartree (MCTDH) approach, along with some vMCG calculations. 7Given that it is significantly easier to compute diabatic populations with MCTDH, and adiabatic populations with AIMS, we provide both quantities in Fig. 3.More details regarding the LVC model calculations can be found in the SI.For the diabatic populations in the DMABN model, vMCG and MCTDH are essentially indistinguishable from one another and MASH, FSSH-nacv and dFSSH-nacv produce significantly more accurate results than FSSH-vel and dFSSH-vel.While the situation is less clear cut in the fulvene model, both MASH and vMCG very accurately match the MCTDH result up to ≈ 10 fs and remain relatively close to it after that.For the adiabatic populations, the trend in the results for all of the approaches is largely the same as in the ab initio simulations.* Most importantly, the vMCG adiabatic populations agree best with those from MASH and FSSH-nacv, further suggesting that these surface hopping approaches are performing the best among all of the 'on-the-fly' approaches that we have tested.Given that the main difference between the vMCG and AIMS algorithms is how the Gaussians are propagated further suggests that this is the source of the error observed in the AIMS results.We now focus on comparing the best surface hopping approaches of FSSH-nacv, dFSSH-nacv and MASH with AIMS for nuclear observables, and we provide the FSSH-vel and dFSSH-vel results in the SI for completeness.Not all nuclear observables are particularly sensitive to the nonadiabaticity of the problem, however.For example one particularly interesting nuclear observable in the case of DMABN is the twist angle of the dimethylamino group.This is because the ground state structure is untwisted, while both of the minimum energy configurations of the S1 surface are twisted.In addition, there are two minimum energy conical intersections * Interestingly, the fulvene model unlike in the ab initio simulations seems to also provide a case for which applying decoherence corrections appears to make the FSSH results worse.
(MECIs) between S1 and S2, one of which is twisted and the other that is not. 56It is therefore interesting to ask whether the onset of twisting in the dynamics is directly connected to the nonadiabatic transition.In Fig. S3 in the SI, we give the dynamical twist angle computed using a selection of methods.The fact that the observed twist angles are within the statistical error for all methods suggests that the nonadiabatic transitions are occurring through the untwisted MECI and the observed twisting arises from the topology of the S1 surface.][65] The product yields in the photodissociation of ethylene are an example of nuclear observables that end up being more sensitive to the nonadiabaticity of the problem.In the following, we group the possible products according to the four different channels depicted in Fig. 4. 66,67 All products observed in our computational simulations match those found in the corresponding experiments, [68][69][70][71][72] with the exception of the C-C dissociation process.As commented on in previous theoretical studies of ethylene, 67 the C-C dissociation is a result of the inadequacy of the basis set used in The electronic populations of the upper adiabatic state (⟨ P+(t)⟩) and the diabatic state that coincides with this adiabat at the Franck-Condon geometry (⟨ Pdiab (t)⟩), obtained for the linear vibronic coupling (LVC) models that were constructed for DMABN and fulvene in Ref. 7. The MCTDH and vMCG results were also obtained from Ref. 7. We do not consider the analogous LVC model for ethylene, as it was found to not give rise to any electronic population transfer.
the electronic structure calculations, which gives a S0 − S1 excitation energy of 10.2 eV at the Franck-Condon geometry, which is much larger than the experimental value of 7.6 eV, 73 and most importantly significantly above the C-C bond energy of 7.7 eV.The dynamical product yields calculated for a range of dynamical methods are given in Fig. 4. Of the products considered, the yields of H elimination, C-C dissociation and ethylidene formation are largely the same for all the methods considered, at least within the statistical error.The geometries corresponding to ethylidene are known to match those of various conical intersection seams in the system. 67,74he product yield for this process therefore qualitatively describes the initial approach and subsequent exit of the conical intersection regions during the dynamics, explaining why it also qualitatively matches the average effective coupling between the Born-Oppenheimer surfaces, as given in Fig. 1.
Of particular interest is the product yield for H2 elimination, where dFSSH-nacv, MASH and AIMS all give rise to an enhanced product yield over FSSH-nacv.Such behaviour is reminiscent of the use of surface hopping to cal-culate nonadiabatic thermal rates.40][41][42] More specifically in the thermal rate problem, it was found that at short times the two-hop trajectories between the ground-state reactant and product geometries were the ones that were responsible for the incorrectly suppressed reaction rate in FSSH. 38In the H2 elimination process, the trajectories associated with an odd number of hops are important, because the reaction proceeds from the excited to the ground electronic state.As a result, we find that it is the three-and five-hop trajectories † that are the problem for FSSH.The advantage of AIMS and MASH is that the correct nonadiabatic rate is reproduced without the need for ad hoc decoherence corrections.Note that while performing AIMS simulations for long enough times to ob-Figure 4: Dynamical product yields for the major products in the photodissociation of ethylene, calculated using various dynamical methods.The width of the shading represents twice the standard error.serve the reaction is expensive, ‡ this is not the case for the independent MASH trajectories.
In conclusion, we have applied a recently proposed independent-trajectory surface hopping approach, MASH, to perform ab initio nonadiabatic dynamics simulations in molecules.We compared MASH with a set of well established methods over a series of two-state benchmark systems -ethylene, DMABN and fulvene -with the surfaces and couplings computed on-the-fly using various ab initio electronic structure methods.Both electronic and nuclear observables were considered.
Overall, MASH is likely to be the most suitable approach for performing ab initio simulations in molecules due to its accuracy and efficiency.For electronic population observables, MASH is able to correctly describe the effects arising from the nonadiabatic force, which was seen to be absent in AIMS and the most commonly used velocity rescaling scheme of FSSH.Such findings were also corroborated by the use of LVC models that were parametrized by fitting to electronic structure data for these molecules, where comparison to numerically exact quantum dynamics was possible.Photodissociation product yields can also be accurately and robustly calculated with MASH, because it solves the inconsistency problem of FSSH without the need for ad hoc decoherence corrections.
Our analysis has uncovered a potential shortcoming within AIMS and it will be interesting to ascertain whether the nonadiabatic force can be correctly incorporated into the motion of the Gaussian basis functions within the AIMS algorithm.For the systems that we have tested here, the independent trajectory approximation is seen to be valid, so that even if AIMS can be fixed, the computational simplicity of the independent-trajectory nature of MASH is still likely to be superior.An interesting question is therefore whether ‡ Because of this issue, we only run the AIMS simulations up to 200 fs.Some algorithmic advancements in AIMS have been designed to alleviate this issue, 75 but we do not directly consider them here.photochemical (or other molecular) processes can be found where the independent trajectory approximation does break down and where coupled-trajectory approaches like AIMS do offer a distinct advantage.
There are some approximations that are shared by all of the 'on-the-fly' approaches tested in this work.For example, all of them sample the initial nuclear phase-space variables from a Wigner distribution and propagate them classically.Therefore it is hard to ascertain how severe the potential issues of zero-point energy leakage and the lack of nuclear tunnelling are for these systems.MASH is guaranteed to thermalize correctly in the long-time limit with a classical nuclear bath, 76 but this is not guaranteed if some of the nuclei have a large zero-point energy.A direct dynamics version of vMCG has been used to study these systems, and it would be good to estimate the robustness of the classical nuclear approximation in these systems by comparing our results to this in future work.
One problem of comparing different methods in ab initio simulations is the difficulty in making sure that the results obtained with different code bases are compatible with one another.To this end, we have provided an extensive SI which addresses in detail the potential sources of discrepancies between different codes, as well as showing how the dynamics in SHARC 52 and the AIMS/MOLPRO package 55 can be made consistent with one another.For future benchmarking exercises, it would be nevertheless useful to have more unified code packages where the majority of the most commonly used methods are implemented, as well as electronic structure codes where all possible quantities, such as NACVs, 77 can be calculated.
The MASH algorithm used in this work is currently only applicable to two-state systems and there is a need to generalize for multi-state problems.Two distinct approaches have emerged for generalizing MASH, [78][79][80] both of which were applied to study the photochemistry of cyclobutanone in the recent JCP community challenge. 81,82We wait to see how both ideas develop in the future.
Supporting Information: The supplementary material provides all of the necessary information to reproduce the surface hopping and AIMS results in the main text.This includes the MOLPRO input and SHARC template files, the ground-state geometries and frequencies and a python script for computing the ethylene product yields.
(28) Plasser, F.; Crespo-Otero, R.; Pederzoli, M.; Pittner, J.; Lischka, H.; Barbatti is then given by where e −i H(t k ) is a matrix exponential, cα are the electronic wavefunction coefficients in the crude adiabatic basis and cα (t) = c α (t).In this work, we use n elec = 25.Once the electronic wavefunction has been successfully propagated to t + ∆t, it is then transformed back into the adiabatic basis as follows Given that the couplings in the diabatized Hamiltonian are generally less localized compared to the NACVs, we use this electronic propagation scheme even when analytic NACVs are available.
In order to rescale the velocity when analytic NACVs are not available, we use the following finite difference scheme ⟨ψ + (q)|ψ − (q + ∆xj)⟩ sgn(⟨ψ − (q)|ψ − (q + ∆xj)⟩) where j is a unit vector pointing along the Cartesian coordinate associated with nuclear degree of freedom, j.Because the NACV is only used to determine the direction of the velocity rescaling, the signs of the NACVs at different time steps do not have to be consistent.The NACV is only computed at time steps where either a hop or frustrated hop have taken place and we use ∆x = 0.001 Bohr in all of our simulations.
In order to test this finite difference scheme, j d j (q(t))v j (t) can be calculated in the small nuclear time step limit as and compared with the same quantity computed using Eq.S4.Additionally, Eq.S5 was used to S3 calculate j d j (q(t))v j (t) for DMABN in Fig. 1 of the main paper.
For both the finite difference scheme for the NACVs and the local diabatic electronic propagation scheme, the required wavefunction overlaps were computed using the associated code in SHARC.S4 In addition, a Löwdin orthogonalization is performed on all overlap matrices to ensure unitarity.
In the case of ethylene, this orthogonalization procedure is performed in the three-state space.

Dynamics Methods
Ab Initio Multiple Spawning (AIMS) Our AIMS simulations for ethylene and fulvene were performed using the AIMS/MOLPRO code.S7 A time step of 20 a.u. was used outside the nonadiabatic coupling region.The coupling region is entered when h j d j v j > 0.005 E h , for which the time step was reduced to 5 a.u..The value of this coupling threshold parameter was determined by reducing its value until there was no noticeable change in the obtained results.
In the coupling region, Gaussians can be spawned whenever j d j v j reaches a maximum along the trajectory.For this to occur, the overlap between the parent and child Gaussians at the spawning point must also exceed a threshold value, set for our calculations to 0.6.A minimum population of 0.01 is also required for a Gaussian to spawn a child.
The threshold for energy violation over a time step is set for ethylene and fulvene as 0.03 and 0.01 E h respectively.If this is violated, the time step is halved unless the minimum value of 1 a.u.
is reached, after which energy violation is ignored.Ad hoc features of the AIMS simulations, such as setting a decoherence time, are not used in our simulations in order to obtain the most accurate and rigorous AIMS benchmark for these systems.
When performing AIMS simulations, it is important that the initial nuclear geometry in the MOLPRO input and Geometry.datfiles exactly match for each trajectory.The initial momenta in the Geometry.datfile must be provided in atomic units.This can be obtained from the initial velocities used in SHARC by multiplying them by the appropriate atomic mass in atomic units.
For the AIMS simulations, 400 and 150 trajectories were run for ethylene and fulvene respectively.Of these, 8 trajectories failed for ethylene due to a failure in converging the CASSCF cycle.
We did not perform the AIMS calculation for DMABN ourselves and the results were instead taken S4 from Ref. S8.The details for the AIMS simulation of DMABN can therefore be found there.
For the linear vibronic coupling (LVC) models, the AIMS dynamics were performed with a modified version of the FMS90 code implemented in MOLPRO.S9 2000 AIMS trajectories were run with a times step of 2.0 a.u., reduced to 0.5 a.u. in the coupling region.The coupling region was entered when h j d j v j > 0.002 E h and a minimum population of 0.001 and an overlap threshold of 0.6 were set for spawning Gaussians.All other parameters were kept at their default value.
For AIMS simulations, the system has to be initialized in an adiabatic state, whereas for all other approaches in the LVC models, we initialized the dynamics in a diabatic state.Fig. S1 gives the MASH adiabatic populations associated with starting in an adiabatic and diabatic state, which illustrates that apart from the t = 0 value in the DMABN model, the initial conditions do not make much of a difference to the obtained populations in these LVC models.As is standard practice in AIMS, the electronic populations are calculated using the overlaps between the Gaussians S8 where N runs is the number of independent AIMS runs and N Gauss,n (t) is the number of Gaussians S5 in run n at time t.Additionally, q k,n (t) and p k,n (t) are the classical nuclear phase-space variables at which the Gaussian g k,n is centered and C k,n,+ (t) is its time-dependent weight, where the + subscript indicates that only Gaussians on the upper surface are used to compute the population of the upper adiabat.
In contrast for nuclear observables like the DMABN twist angle and the ethylene product yields, only the diagonal elements of the Gaussian basis are used.For observables Ô that are functions solely of the nuclear positions, q where the lack of the + subscript on the Gaussian weights now indicates that all Gaussians contribute to the observable, irrespective of the surface they are on.

Fewest-Switches Surface Hopping (FSSH)
In FSSH, S10 the nuclei are propagated according to Newton's equations of motion on the active Born-Oppenheimer surface and the electronic wavefunction is propagated according to the timedependent Schrödinger equation [Eq.S2b].The velocity-verlet scheme is used for the nuclear dynamics with ∆t = 0.25 fs.In order to initialize the electronic subsystem in the excited state, the initial wavefunction coefficients are given by c − (0) = 0 and c + (0) = 1.
In order to describe nonadiabatic transitions, the active surface is changed stochastically with the following probability where Re[•] returns the real part of the quantity and any negative probabilities are set to zero.
Electronic population observables are constructed in terms of the active state, rather than the electronic wavefunction coefficients.
At an attempted hop, the kinetic energy along the NACV must first be compared with the electronic transition energy.A hop from state α to β is energetically allowed if

S6
in which case the velocity along the NACV is rescaled according to In the event that Eq.S9 is not satisfied, the hop is frustrated and the nuclear velocity along the NACV is reflected as follows This treatment of frustrated hops in FSSH required a modification to the SHARC code in order to ensure that the nuclear velocity is always reflected at a frustrated hop.
To add decoherence into FSSH simulations, we use the energy-based decoherence scheme de- where E kin is the total nuclear kinetic energy and C = 0.1 E h is the decoherence parameter.
All FSSH simulations were performed using SHARC 2.0.S6 For ethylene, DMABN and fulvene, 1000, 600 and 600 trajectories were run respectively to ensure well converged results.However far fewer trajectories would be needed for a good qualitative description of the observables.For ethylene only 5 trajectories or fewer failed as a result of being unable to converge the CASSCF cycle, whereas for fulvene it was 3 trajectories or fewer.
For the LVC models, the surface hopping approaches were initialized in a diabatic state, to allow direct comparison with the MCTDH and vMCG results.Ref. S13 explains how to compute diabatic populations with FSSH.We choose to use a very large number of trajectories (100000) to essentially remove any statistical error within the results, which is easy to do in model calculations.
However, far fewer trajectories could have been used in practice to qualitatively reproduce the ) where V λ is the Born-Oppenheimer surface for state |ψ λ ⟩, Pλ = |ψ λ ⟩ ⟨ψ λ | is the associated electronic population operator, d (µ,λ) j is the NACV between states |ψµ⟩ and |ψ λ ⟩ and σ(µ,λ) x = |ψµ⟩ ⟨ψ λ | + |ψ λ ⟩ ⟨ψµ| is the associated electronic coherence operator.More details associated with this formula are given in the supporting information (SI).MASH is constructed to describe the Born-Oppenheimer (adiabatic) force through the active surface and the nonadiabatic force through the velocity rescaling performed along the NACV.

Figure 2 :
Figure2: The electronic population of the upper adiabatic state (⟨ P+(t)⟩) for ethylene, DMABN and fulvene, computed for a wide variety of different methods.The width of the shading represents twice the standard error, which is less than the line width if not visible.The AIMS result for DMABN was obtained from Ref.56.For the surface hopping approaches, the average number of upward hops (n hops ↑ ) downward hops (n hops ↓ ) and frustrated hops (n frust ) are also given.

Figure 3 :
Figure3: The electronic populations of the upper adiabatic state (⟨ P+(t)⟩) and the diabatic state that coincides with this adiabat at the Franck-Condon geometry (⟨ Pdiab (t)⟩), obtained for the linear vibronic coupling (LVC) models that were constructed for DMABN and fulvene in Ref.7.The MCTDH and vMCG results were also obtained from Ref. 7. We do not consider the analogous LVC model for ethylene, as it was found to not give rise to any electronic population transfer.

Figure S1 :
Figure S1: Figure showing the MASH dynamical populations of the upper adiabatic state in the LVC models, when starting in either an adiabatic or diabatic state.The MASH results when starting in an initial diabatic state are the same as given in Fig. 3 of the main paper.
scribed in Ref. S11,S12.For a trajectory with current active surface β, the electronic coefficients are updated at each time step as follows

Figure S2 :
Figure S2: A flowchart defining the various products formed after the photoexcitation of ethylene.The first row of the diagram determines whether any C-H bonds have been cleaved, the second whether the C-C bond has been broken, the third whether any hydrogen atoms have migrated to another carbon and the fourth looks for any hydrogen bridging atoms.

Figure. S2 Figure S3 :
Figure.S2 contains a flowchart that shows how we determined the products formed at a given time step of a trajectory.This was based on the scheme used in Ref. S16 and was tested by analyzing a sufficiently large random selection of trajectories by hand.The bond-length criteria used in our product specifications are very different from the equilibrium bond lengths in ethylene, because ethylene is highly vibrationally excited once it relaxes back to the electronic ground state after the initial electronic photoexcitation.
Is max(r C-H i )> 2.4 Bohr?