Sampling reactive regions in phase space by following the minimum dynamic path.

Understanding mechanistic aspects of reactivity lies at the heart of chemistry. Once the potential energy surface (PES) for a system of interest is known, reactions can be studied by computational means. While the minimum energy path (MEP) between two minima of the PES can give some insight into the topological changes required for a reaction to occur, it lacks dynamical information and is an unrealistic depiction of the reactive process. For a more realistic view, molecular dynamics (MD) simulations are required. However, this usually involves generating thousands of trajectories in order to sample a few reactive events and is therefore much more computationally expensive than calculating the MEP. In this work, it is shown that a "minimum dynamic path" (MDP) can be constructed, which, contrary to the MEP, provides insight into the reaction dynamics. It is shown that the underlying concepts can be extended to directly sample reactive regions in phase space. The sampling method and the MDP are demonstrated on the well-known 2-dimensional Müller-Brown PES and for a realistic 12-dimensional reactive PES for sulfurochloridic acid, a proxy molecule used to study vibrationally induced photodissociation of sulfuric acid.

In order to construct a realistic dynamical path, different approaches are necessary. For example, it is possible to formulate the task of finding a path connecting a reactant and product state as a two-point boundary value problem. 12 Such an approach has been formulated successfully in terms of a minimization problem involving the Onsager-Machlup (OM) action (which requires second derivatives) 13 or a modified target function Θ, that serves as an approximation to the OM action (involving first derivatives only 14 ), for which the trajectory is expanded in a Fourier series, as is known from Fourier path integral simulations. 15 The OM action can be thought of as a measure for the violation of Newton's equations of motion: Every Newtonian trajectory has an OM action of exactly zero. However, starting from arbitrary boundary values, paths determined by minimizing the OM action were found to not necessarily conserve energy, whereas the Θ trajectories do, but are quite similar to a MEP for the Müller-Brown PES. 14 Alternatively, molecular dynamics (MD) simulations based on accurate energy functions can be used to gain insight into dynamical processes governing a reaction. [16][17][18][19][20][21] However, for a realistic description of a chemical reaction a statistically significant number (10 4 or more) of such simulations is required. This is usually not possible with the most accurate approach which would be ab initio MD simulations at a sufficiently high level of theory. Therefore, energy functions fitted to electronic structure calculations at a relatively high level of theory (multireference CI for triatomic systems 22 or Møller-Plesset perturbation theory for larger molecules 23,24 ) have been used in the past. Running such a large number of MD trajectories is, however, computationally considerably more expensive than calculating the MEP.
For this reason, computational methods were devised to improve the sampling of such rare events. For example, in transition path sampling (TPS), once the dynamical bottleneck of a reaction -its transition state surface -is identified, reactive trajectories can be generated efficiently using Monte Carlo sampling. 25,26 Another example is milestoning, which aims to compute the time scale of complex processes with predetermined "milestones" (slices of the transition path, which are sampled with short trajectories) along a reaction coordinate. 27 Other methods to sample rare events are for example the minimum action method, 28 or the string method, 29 which is based on transition path theory. 30 Reaction rates can also be estimated from methods such as transition interface sampling, 31 transition state theory and extensions thereof, [32][33][34][35] which estimate the reactive flux through a so-called dividing surface.
In the present work the concept of a "minimum dynamic path" (MDP) is considered as an alternative approach and related to the underlying structure of phase space. The MDP corresponds to the lowest energy dynamical (following Newton's equations of motion) reactive path in phase space. Contrary to the MEP, it provides insight into the reaction dynamics and has, by definition, an OM action of zero. Once the transition state of a reaction is known, the MDP can be readily constructed with a computational effort comparable to running a single trajectory (requiring only one additional evaluation of the Hessian). The construction method can easily be extended to generate reactive initial conditions for a microcanonical ensemble of trajectories with arbitrary excess energy ∆E. It is further shown that insight obtained from the MDP is also relevant for reactive trajectories at higher energy.
First, these concepts are investigated for the well-known 2-dimensional Müller-Brown PES 36 for which exhaustive sampling is possible and serves as a validation of the results from the MDP. It is found that particular initial conditions can be prepared which lead to crossing the transition state with certainty. In a next step, the reactive dynamics for a realistic, 12-dimensional reactive PES describing the dissociation dynamics of sulfurochloridic acid 23 is investigated. This molecule is a proxy to study vibrationally induced photodissociation of sulfuric acid. 37,38 2 The Minimum Dynamic Path and Sampling Reactive Initial Conditions Since the transition state of a PES is defined as the configuration x TS with the highest potential energy V (x TS ) = E TS along the minimum energy path (MEP), every reactive trajectory with constant energy must have a total energy E tot = E pot + E kin ≥ E TS . It follows that the reactive trajectory with the lowest possible total energy E tot = E TS must pass through the TS exactly and have a kinetic energy E kin = 0 at the TS. The path this special trajectory follows through configurational space will henceforth be referred to as the minimum dynamic path (MDP). Note that the MDP and the MEP differ because a dynamical system is not only guided by forces, i.e. the gradient of the PES (which solely determines the MEP), but also keeps a "memory" of past gradients in its current momenta. In the over-damped limit, this memory is completely lost and trajectories approach the MEP (see also section S1).
Transition states are mathematically defined as first-order saddle points of the potential energy surface V (x), i.e. saddle points at which the Hessian has only one negative eigenvalue λ n with corresponding eigenvector e n . The MDP can be readily approximated by starting two trajectories at the TS with initial momenta p 0 = ± e n (where is small) and following their paths through phase space until the desired reactant or product state of the reaction is reached. Since the equations of motion are symmetric under time reversal, both paths can be combined to give the MDP.
It is also possible to extend this procedure to generate reactive initial conditions (x, p) in phase space with E tot ≥ E TS . For this, the concept of a "separating hypersurface" is introduced: Consider a PES with two minima labelled I and II and a saddle point TS separating them. Starting at an arbitrary point P in configuration space and following the direction of steepest descent, every path initiated at P will reach either I, II, or rarely the TS. The sets of points that converge to I or II form the basins of attraction for the respective minima, whereas the set of points that converges to TS forms a hypersurface separating those basins of attraction (see Figures 1 and 2A). This hypersurface must be crossed by every reactive trajectory going from I to II at some point and is referred to as the separating hypersurface.
Note that points on this hypersurface do not react equally likely to either I or II, i.e. this surface differs from the isocommittor surface. 6,7,39 It is also distinct from the concept of a dividing surface used in TPS. The dividing surface is the hypersurface with the lowest number of recrossings for which trajectories reach educt and product states equally likely. 26 The topology of the separating hypersurface on the other hand does not contain any dynamic information and depends solely on the underlying PES. It is important to point out that while the sampling method presented here does generate reactive initial conditions, they do not correspond to a thermal ensemble. In the limit of infinite sampling, rates from a microcanonical and a canonical treatment are identical, though, even for a few-particle system. 40 To sample a reactive initial condition (x, p) in phase space with energy E tot ≥ E TS , these steps are followed: 1. A point x 0 with potential energy E pot ≤ E tot is generated on the separating hypersurface as follows: Starting at configuration x = x TS + e ⊥ (see Figure 1A), where x TS is the configuration of the TS, is a small number and e ⊥ is a random direction perpendicular to e n , the gradient is followed along the direction of steepest ascent until the desired energy E pot is reached at point x 0 (see Figure 1B). The small initial displacement by e ⊥ is required because at the transition state, the direction of steepest ascent is undefined. Following the gradient ensures that the point x 0 lies on the separating hypersurface, which generally is curved.

2.
A momentum vector p 0 with random direction is drawn from an unbiased distribution and its magnitude scaled such that E tot = E kin + E pot . It should be noted that drawing momenta from an unbiased distribution and scaling their magnitude generates reactive initial conditions corresponding to a microcanonical ensemble. In order to obtain thermal reactive initial conditions the momenta would have to be drawn from a fluxweighted distribution. 41 3. A trajectory is started with initial conditions (x 0 , ±p 0 ). A short MD simulation with both initial conditions ( Figure 1C) is necessary in order to verify that the two trajectories move towards the two different basins of attraction, which is not guaranteed a priori.
4. If this requirement is met, a longer trajectory is run until time τ (which is chosen arbitrarily) starting from either (x 0 , p 0 ) or (x 0 , −p 0 ) and the final positions and momenta (x τ , p τ ) are recorded ( Figure 1D). Since the equations of motion are time reversal symmetric, a trajectory starting from the initial condition (x τ , −p τ ) is reactive and will pass the separating hypersurface after the chosen time τ . All reactive trajectories generated in this fashion form the "reactive phase space", i.e. a subset of productive initial conditions which lead to reaction (cross the separating hypersurface). Note that this sampling procedure is conceptually similar to "shooting" in TPS, 42 but here, trajectories are always started on the separating hypersurface instead of a point in configurational space reached after a random time. Figure 1: Schematic representation of the sampling procedure to generate reactive initial conditions with arbitrary energy E tot > E TS . The topology of the potential energy surface is indicated by contour lines and colours, with red tints signifying high-and blue tints low energy regions. A: The position of the transition state is x TS , I and II label the basins of attraction of the respective minima and the solid black line indicates the separating hypersurface that must be crossed by every trajectory in order to react. B: Starting from x TS + e ⊥ , the direction of steepest ascent (red arrows) is followed until the desired potential energy V (x 0 ) = E pot ≤ E tot is reached at x 0 . Note that the small displacement e ⊥ from x TS is necessary, because at the TS, the direction of steepest ascent is undefined. C: A momentum vector p 0 with random direction is drawn from a uniform distribution and scaled such that E kin + E pot = E tot . Two short trajectories (indicated by dotted lines) are started from the initial conditions (x 0 , p 0 ) (red) and (x 0 , −p 0 ) (blue) in order to confirm that both trajectories evolve towards different basins of attraction. If not, a new combination (x 0 , p 0 ) is generated. D: Either of the trajectories is followed up to time τ and its final state (x τ , p τ ) is recorded. Due to time reversal symmetry, a trajectory starting from (x τ , −p τ ) will pass the separating hypersurface after time τ and is therefore reactive.
In order to verify that the procedure described in section 2 can be applied and used to extract information about the underlying dynamics, the well-known 2-dimensional Müller-  Figure 2A).
Since the phase space corresponding to the Müller-Brown system is only 4-dimensional ((x, y, p x , p y )), it is possible to determine all regions in phase space that lead to reaction by exhaustive sampling. In order to test whether the method described in section 2 samples the same regions as such an unbiased sampling, an unambiguous definition of a reactive trajectory is required. For this purpose, the trajectory of a particle with mass m = 1 and total energy E tot = E TS + ∆E is considered to be reactive if it reaches minimum 2 (product state) within time t max = 20 starting from minimum 3 (reactant state). A trajectory is terminated when it reaches the product within t ≤ t max in which case it is "reactive" or after t max in which case it counts as "unreactive", even if it could react at a later time. Reactant-and product states are defined to be the set of points (x, y) enclosed by ellipses centered around the corresponding minima given by the parametric equations with x 0 ≈ −0.56, y 0 ≈ 1.44 and φ = π/4 defining the reactant-and x 0 ≈ 0.62, y 0 ≈ 0.028 and φ = 0 defining the product state (see blue and red ellipses in Figures 2B-F). These definitions are largely arbitrary, but needed in order to define an unambiguous reactantand product state. The equivalence between exhaustive unbiased sampling and the method described in section 2 can be tested with any arbitrary definition of states and choice of t max .
The reactive part of phase space is sampled exhaustively by generating trajectories with unbiased random initial conditions (x, y, p x , p y ) such that the total energy corresponds to If a trajectory starting from (x, y, p x , p y ) reaches the product state in time t 1 and a trajectory starting from (x, y, −p x , −p y ) reaches the reactant state in time t 2 such that t 1 + t 2 ≤ t max , the initial condition (x, y, p x , p y ) belongs to the reactive part of the phase space (see Figures 2B-F). Note that with increasing excess energy ∆E, more states on the separating hypersurface become energetically accessible, which leads to a widening of the transition region, which is also known from explicit reactive MD simulations. 43 Further, while the MDP remains representative for most of the reactive part of phase space, with increasing excess energy additional "reaction channels" open up (see Figure 2B-F) for which the time until the product state is reached can differ significantly (see Figure   S2). It is interesting to note that, depending on ∆E, the tubes emanating from either A or B do not cover the entire boundary of the ellipses but rather correspond to discrete regions in phase space (see Figure S1).
In summary, the method described in section 2 is able to efficiently sample the reactive regions in phase space and converges to the same set of initial conditions as unbiased sampling (see Figure 2). This is also the reason why the results from explicit sampling are not reported separately. Comparing the reactive part of phase space (even for large excess energy ∆E) with the MDP also shows that the MDP is representative for the reactive phase space explicitly sampled by the system, see e.g Figure 2 panels A and C. . For large ∆E, alternative "reaction channels" become accessible, but pathways resembling the MDP remain dominant. Both, exhaustive unbiased sampling and the procedure described in section 2 converge to the same reactive subspace (shown as coloured regions in panels B-F). See Figure S1 for a variant of this figure with oversaturated colours that increase visibility of reaction channels.
It should be noted that a decomposition of phase space into reactive and non-reactive subspaces has been observed in earlier work and is even possible in the absence of an imposed maximum reaction time t max due to the presence of trapped orbits. 44 4 Application to Molecular Systems: Sulfurochloridic

Acid
In order to study the MDP for a concrete molecular system, the dissociation of sulfurochlo- In order to verify that the umbrella motion -found to be an important degree of freedom in the MDP -is also relevant in reactive trajectories with excess energy, 100 initial conditions that lead to dissociation within 175 fs were run with an excess energy of ∆E = 5 kcal/mol, see Section 2. MD simulations using a custom code were run with the velocity Verlet in-tegrator 47 and a time step of 0.1 fs for a total of 2000 time steps. All these trajectories follow a similar motion compared to the MDP prior to the elimination reaction (see Figure   3): Averaging all reactive trajectories exhibits the same oscillatory behaviour as is observed for the MDP. For short time-scales (∼ 0.2 ps) prior to the reaction, it appears that reactive trajectories all exhibit a similar "concerted" motion, at least for this particular case, and that this coupling can be visualized and analyzed. Of course, the details of this observation depend on the underlying PES.
For quantifying key aspects of the underlying motion, the total energy of the MDP trajectory and the 100 reactive trajectories is decomposed into normal mode 48 contributions. Table 1 reports the harmonic frequencies ω of all 12 normal modes corresponding to internal degrees of freedom together with their associated motion.
which ensures that the true potential energy E pot (x) of configuration x is divided exactly among the normal modes, such that i E i pot = E pot (x).
For the kinetic energy, the velocities v derived from the momenta p are transformed to a vector w in normal mode space in the same way x is transformed to q. Note that due to the way the transformation into normal modes is defined, the entries w i of w correspond to the "momentum" of each normal mode i divided by the square root of the associated effective mass. The kinetic energy E i kin of normal mode i is defined accordingly as where E kin (p) is the true kinetic energy according to the momenta p, which ensures that it is divided exactly among the normal modes.
The total energy of a normal mode is then simply the sum of its kinetic and potential energy.
It should be noted that the method described above does not guarantee E tot = i E i kin +E i pot because of rovibrational coupling, which leads to some vibrational energy being unaccounted for in the transformation to the Eckart frame. Nonetheless, the total energy is conserved approximately and will fluctuate around a constant mean with a typical amplitude of < 0.1 kcal/mol, which is sufficient for the present purpose. Also, it should be noted that the normal mode decomposition assumes a harmonic PES and is thus only strictly valid close to the equilibrium geometry. When the normal mode decomposition is performed in highly anharmonic regions, results can get distorted and should therefore always be considered with care. Still, the decomposition provides a quantitative comparison of the relative importance of normal modes between different trajectories. Note that other methods to decompose the energy into normal mode contributions are possible. 50,51 Some of this energy is transferred to mode 16 during the next 100 fs. Between -75 fs and -50 fs before the reaction, modes 14 and 16 start to lose energy, whereas mode 18 becomes excited. Around t = −30 fs modes 13 and 18 contain by far the largest fraction of the total energy (33% and 30% for the MDP, 34% and 20% for the averaged results). Note that while it is difficult to directly correlate timepoints between MDP and the averaged results due to a difference in kinetic energy of up to ∆E = 5 kcal/mol, both results show similar trends and the dynamics of the energy flow between modes is comparable. Figure 4 shows the normal mode energy decomposition for different times prior to the reaction for the MDP and an average for an ensemble of 100 reactive trajectories. While the MDP and the ensemble statistics differ slightly, they both follow qualitatively similar trends.
Further, the analysis confirms that mode 13, which corresponds to an "umbrella motion" of the SO 3 moiety, and mode 18, which corresponds to the OH stretch vibration, are excited prior to the reaction, a trend which could already be deducted from a visual inspection of the MDP. Thus, the molecular picture that arises from this analysis suggests that energy flows from other degrees of freedom into these modes, which promote dissociation. This is consistent with a previous study which demonstrated that vibrationally induced photodissociation can be promoted via vibrational energy redistribution by exciting the -OH stretching motion. 23

Reactive Trajectories versus Vibrational Energy Relaxation
In order to test whether the insights about the reaction dynamics gathered from analysis of the MDP are applicable in a more general context, the differences between non-reactive and dissociative trajectories of HSO 3 Cl after OH-stretch overtone excitation were studied using the normal mode decomposition scheme described earlier. The reactive MD simulations (with the exception of using the velocity Verlet integrator 47 ), were carried out along the same lines as in the previous study. 23 The change of the integrator was necessary to allow a meaningful normal mode energy decomposition analysis.
Individual trajectories were started from a geometry optimized structure of HSO 3 Cl. The system was heated to 300 K. The equations of motion were propagated using the leapfrog  Brown PES the initial conditions do not decompose phase space into two types of trajectories that could be distinguished after initial preparation (here vibrational excitation), at least when energy content in the participating modes is used to differentiate between them. This is consistent with the intuitive notion that the fate of a trajectory -whether it leads to reaction or not -is decided in the phase immediately before bond breaking occurs, see also

Conclusion
The concept of an MDP was introduced, which is related to the MEP but includes dynamical effects due to inertia. It was shown that the MDP resembles the dynamics of an ensemble of reactive trajectories, thus providing valuable insight into the reaction dynamics of a system of interest with a single trajectory. Further, a method was described which allows direct sampling of reactive phase space, making the generation of reactive initial conditions for MD simulations more efficient. The techniques were demonstrated for the 2-dimensional Müller- Brown model system and a more realistic 12-dimensional reactive PES of sulfurochloridic acid. Conclusions drawn from the MDP about which modes promote dissociation in HSO 3 Cl also hold for OH-stretch overtone induced photodissociation, which was confirmed by energy decomposition analysis for reactive and non-reactive trajectories of HSO 3 Cl.
The techniques discussed here can also be applied to ab initio molecular dynamics simulations in the gas phase. Because for an initial assessment of a particular reactive trajectory between a given reactant and product state only one MDP simulation is required, it will be possible to use high-level electronic structure theory methods (MP2 or even CCSD(T) depending on system size) to obtain information about the participating degrees of freedom.
Similarly, applications to condensed phase systems, e.g. to investigate competitive ligand rebinding 53 or bimolecular reactions in solution 54 should be possible to determine relevant degrees of freedom accompanying the reaction. Once such active degrees of freedom are known, they could be used to drive a chemical system by controlled excitation of the corresponding motions.