David Bonhommeau

Donald G Truhlar

Mixed Quantum/Classical Investigation of the Photodissociation of NH 3 ( Ã) and a Practical Method for Maintaining Zero-Point Energy in Classical Trajectories

Keywords: numbers: 31.15.xg, 31.15.xr, 31.50.-x, 33.80.Gj 1

The photodissociation dynamics of ammonia upon excitation of the out-of-plane bending mode (mode ν 2 with n 2 = 0, ..., 6 quanta of vibration) in the à electronic state is investigated by means of several mixed quantum/classical methods, and the calculated final-state properties are compared to experiments. Five mixed quantum/classical methods are tested: one mean-field approach (the CSDM method), two surfacehopping methods (the FSTU and FSTU/SD methods), and two surface-hopping methods with zero-point energy maintenance (the FSTU/SD+TRAPZ and FSTU/SD+mTRAPZ methods). We found a qualitative difference between final NH 2 internal energy distributions obtained for n 2 = 0 and n 2 > 1, as observed in experiments. Distributions obtained for n 2 = 1 present an intermediate behavior between distributions obtained for smaller and larger n 2 values. The dynamics is found to be highly electronically nonadiabatic with all these methods. NH 2 internal energy distributions may have a negative energy tail when the zero-point energy (ZPE) is not maintained throughout the dynamics. The original TRAPZ (TRAjectory Projection onto Zero-point energy orbit) method was designed to maintain zero-point energy in classical trajectories, but we find that it leads to unphysically high internal vibrational energies. The mTRAPZ (minimal TRAPZ) method, which is new in this work and provides a general method for maintaining zero-point energy in either single-surface or multisurface trajectories, does not lead to unphysical results and is much less time consuming. The effect of maintaining zero-point energy in mixed quantum/classical dynamics is discussed in terms of agreement with experimental findings. The dynamics for n 2 = 0 and n 2 = 6 are also analyzed to reveal details not available from experiment, in particular, the time required for quenching of electronic excitation and the adiabatic energy gap and geometry at the time of quenching.

I. INTRODUCTION

The photodissociation dynamics of chemical reactions can be treated classically or quantummechanically, with the accuracy depending on the details of the algorithm and the nature of the process simulated. [START_REF] Schinke | Photodissociation Dynamics[END_REF] While direct photodissociation processes are reasonably well described by mixed quantum/classical dynamics, indirect photodissociation processes, like those proceeding through a potential-energy barrier or a resonance, are more properly modeled by wave-packet dynamics. The MCTDH (Multi Configuration Time-Dependent Hartree) wave-packet algorithm has been shown to have a wide range of applicability; although this method is best suited for small systems (less than 12 degrees of freedom), [START_REF] Manthe | [END_REF][3][4] it was also successfully applied to study larger systems like pyrazine [5][6][7] and benzene. 8 Recently developed mixed quantum/classical methods also show a wide range of applicability, and their accuracy has been validated for several systems. [9][10][11][12][13][14][15][START_REF] Truhlar | Quantum Dynamics of Complex Molecular Systems[END_REF] Especially promising methods include an improved mean-field method called CSDM (Coherent Switching with Decay of Mixing), 9,10 a surface-hopping method called FSTU (Fewest Switches with Time Uncertainty), [10][11][12] and an improved FSTU method, namely the FSTU/SD (FSTU with Stochastic Decay) method. 13 These mixed quantum/classical methods combine a quantum mechanical treatment of electron dynamics with a classical mechanical treatment of nuclear motion on coupled potential-energy surfaces. Whereas older mixed quantum/classical methods for this kind of problem were less accurate than trajectory methods for single, uncoupled potential-energy surfaces, the validation tests show that the newer methods have about the same accuracy on coupled potential energies as classical molecular dynamics has for motion on single potential-energy surfaces. But they also suffer from some of the same problems, chief among which is the improper treatment of zero-point energy (ZPE). For example, a long recognized problem in single-surface trajectory calculations is that the products are sometimes produced with less than ZPE, [START_REF] Gray | [END_REF] which can lead to final translational energies higher than those that can be observed experimentally. A related, but somewhat more subtle, problem is that instantaneous vibrational energies can be lower than local ZPE requirements. While the latter are not as strict or as clear as the ZPE requirements of stationary states, they are nervertheless very important in determining the energies at which processes have an appreciable probability of passing through dynamical bottlenecks. 18,19 Several groups have developed techniques to enforce ZPE maintenance in trajectory calculations. The first treatment to fully incorporate ZPE maintenance in molecular dynamics 20,21 has several unphysical shortcomings including discontinuities in linear momenta. Problems with this method and other possible ways to enforce ZPE requirements on trajectories were analyzed carefully by Guo et al. 22 who found no existing method to be satisfactory. Lim et al. [23][24][25] suggested changing the linear momenta in the normal mode representation whenever the instantaneous vibrational energy of one mode drops below its local ZPE. This new method was called the TRAjectory Projection onto Zero-point energy orbit (TRAPZ). When applying TRAPZ to the dissociation of Al 3 , they found that the molecular fragment Al 2 was highly vibrationally excited; a large peak occurred at an energy about 1 kcal/mol (i.e., 0.043 eV) above the true ZPE. 25 This overestimation of vibrational energy was interpreted as resulting from each mode having to conserve its ZPE even if all the other modes have an energy much higher than their ZPE. Alimi et al. 26 proposed a mixed semiclassical/classical algorithm in which stiff modes were treated semiclassically and soft modes classically. This assumes that modes can be sorted into stiff and soft modes, which is not always possible for reactive polyatomic systems. Xie et al. 27 also suggested a method to overcome problems in previous methods by introducing a switching function that should eliminate couplings in the Halmitonian. The generalization of such a method to a general polyatomic system has not, however, been demonstrated, and the method was only tested on the two-mode Heinon-Heiles model. More recently Stock et al. 28,29 also discussed implications of a new semiclassical method that maintains ZPE for spin-boson models and can be applied to more complicated systems.

Although we agree with the analysis of Guo et al. 22 and their conclusions that no completely satisfactory solution to this problem is possible, we continue to encounter applications where it is unacceptable to simply ignore this problem. Therefore, in the present article we describe and apply a very general practical scheme that we have found useful. It is a modification of the TRAPZ method and is called the minimal TRAPZ (or mTRAPZ) method; we found it to give physical results. We will present this new scheme in the context of our application to coupledsurface photodissociation of ammonia, which is dominated by a conical intersection, but it is equally applicable to the kind of single-surface dynamical processes for which ZPE maintenance has usually been previously [START_REF] Gray | [END_REF][18][19][20][21] discussed.

In these contexts, the objective of the present study is to test the capabilities of the most recent mixed quantum/classical methods-the CSDM, FSTU and FSTU/SD-to properly model the dynamics of a real polyatomic system with a conical intersection and to determine a practical way to incorporate ZPE maintenance into this kind of dynamics.

The problem under study is the photodissociation dynamics of ammonia excited in the out-of-plane bending mode ('umbrella' mode) of its à first-excited-singlet electronic state, which has a planar equilibrium geometry. The dynamics of photoexcited ammonia has been experimentally investigated by several research groups, [30][31][32][33][34][35][36][37][38][39][40][41][42] of which we single out two. [30][31][32]43,44 Biesner et al. 30 studied the effect of excitation of the umbrella mode and found that NH 2 fragments were highly rotationally excited when n 2 = 0, where n 2 is the number of umbrella mode quanta initially in vibration, and nonadiabatic dynamics, presumed to be due to passage through or near the conical intersection, was favored. For larger n 2 values the situation became harder to interpret due to the increase in electronically and vibrationally excited products. Bach et al. 31 and Hause et al. 32 have later completed this study by examining the effect of selective excitation of each NH 3 mode on the photodissociation process. They first of all confirmed Biesner et al.'s results concerning the role of the umbrella mode. They also found that the fragmentation of NH 3 molecules excited in their antisymmetric streching mode mainly produced very slow H atoms whereas NH 3 molecules excited in a symmetric streching mode nonadiabatically dissociated. They have interpreted this surprising production of slow H atoms as a signature of the adiabatic character of the dissociation process and concluded that the dissociation was preferentially adiabatic or nonadiabatic according to which stretching mode they chose to initially excite. The present article provides the initial results of our mixed quantum/classical simulations of these dynamics. In particular, we present a study that examines the effect of exciting the umbrella mode.

The different dynamical methods (CSDM, FSTU and FSTU/SD) as well as the TRAPZ and mTRAPZ methods are presented in Section II. Mixed quantum/classical results without maintenance of ZPE are then reported in Section III. The computational efficiency and reliability of each above mentioned mixed quantum/classical method is also discussed in Section III. Results obtained with the TRAPZ and mTRAPZ methods are then collected and discussed in Section IV. Section V gives concluding remarks.

II. THEORY

II.A. Coupled potential-energy surfaces

The potential energy surfaces and their coupling matrix element are taken from previously presented fits 45 to ab initio calculations in the diabatic representation. Diagonalizing the 2 × 2 diabatic potential matrix yields the adiabatic surfaces as eigenvalues, and the nonadiabatic cou-pling vector that couples the adiabatic states is calculated from the eigenvector by equations given elsewhere. 46 The ground and first excited adiabatic electronic states of ammonia will sometimes be called adiabatic surfaces 1 and 2 respectively.

II.B. Initial conditions

The planar and most stable configuration of NH 3 in its à electronic state well is the starting point of the simulation, with the center of mass at the origin and no overall rotation. Vibrational energy is added to this configuration to generate a new configuration at a fixed total energy E tot,0

given by

E tot,0 = V 2,min + 3N -6 k=1 (n k + 1/2)hω k (1)
where V 2,min is the equilibrium potential energy of the à electronic state, the second term is the vibrational energy to be added, n k is the number of quanta of vibration for mode k, and N is 4.

In the present work n k = 0 for all modes except the umbrella mode ν 2 for which seven values of n 2 will be tested: n 2 = 0, ..., 6. We have scaled our umbrella mode frequency of 977.5 cm -1 , obtained by performing a normal mode analysis in the à electronic state well, by a factor of 0.9125 to make the starting energies calculated by Eq. ( 1) agree better with the experimental value of 892 cm -1 ; 45 this allows us to better model the energies of the initial states of the photodissociation.

All other ω k in Eq. ( 1) are unscaled values from the normal mode analysis.

Given the values of E tot,0 , we generate the initial normal mode coordinates and momenta by a combination of methods. For modes with n k = 0, we use a Wigner distribution. For other modes (in the present article, this would only be the excited umbrella mode that is n 2 > 0) we use a quasiclassical normal mode sampling (random assignment of vibrational phase).

Note that, because the vibrational modes are not separable from each other or from rotation when excited by finite amounts, populating each mode separately gives a total energy slightly different from that calculated from Eq. ( 1) and gives a total angular momentum that is slightly non-zero. These discrepancies are small and were ignored.

Once the vibrationally excited configuration is generated, atoms are rotated about the principal axes of the molecule to give a random orientation. Possible overall rotation is removed, and for this reason as well as the nonseparability mentioned in the previous paragraph, the trajectory starts at a total energy (called E α tot,0 for trajectory α) that can be slightly different from E tot,0 . The only effect of this different initial total energy on dynamical results is to slightly change the theoretical excess energy but this is not critical (we will briefly discuss it in Sec. III). The fixed-energy dynamics is then performed in the adiabatic representation by means of several mixed quantum/classical methods that are described in Sec. II D.

II.C. Excess energy

From the knowledge of E tot,0 and the dissociation energy of the ground electronic state we can deduce the average energy available for fragmentation; this is called the excess energy, denoted E exc . This information is crucial to compare to experiments. Knowing the relative kinetic energy E rel (H) of departing H atoms (ionized for detection purposes)

E rel (H) = µ red 2 V 2 R ( 2
)
where

µ red = m(NH 2 )m(H)/m(NH 3 ) is the reduced mass and V R = v(H) -v com (NH 2 )
is the velocity of the departing H atom relative to the center of mass of NH 2 , and the excess energy E exc as calculated from the initial state of the experiment, experimentalists can directly deduce the NH 2

internal energy E int (NH 2 ) according to

E int (NH 2 ) = E exc -E rel (H) (3)
and consequently shed some light on the photodissociation mechanism. Table I collects theoretical excess energies for n 2 = 0, ..., 6 as calculated from the initial conditions of Eq. (1) with scaled ω 2 , the experimental excess energies of Biesner et al. 30 for n 2 = 0, ..., 6, the experimental excess energies of Bach et al. 31 for n 2 = 0, ..., 4, and the experimental excess energy of Hause et al. 32 for n 2 = 0. Theoretical excess energies are slightly higher than experimental ones, due in part to the fact that the dissociation energy used by experimentalists is not exactly the same as the dissociation energy of our potential-energy surface, but this should not strongly affect the results and final conclusions. Note that the classical excess energy would be 0.519 eV higher since

E ZP E [NH 2 ( X)] ≈ 0.
519 eV in the ground electronic state at dissociation, but we define E exc with the quantum zero of energy for comparison to experiment (see Fig. 1 for a schematic energy diagram of adiabatic surfaces).

II.D. Dynamical methods

Three dynamical methods were applied to the photodissociation dynamics of ammonia. These methods have been detailed elsewhere 9,[11][12][13] and we will briefly recall differences between these methods and older ones. Our discussion here assumes that the dynamics are calculated in the adiabatic representation, which is the case in the present article, although the methods are more general and may be applied in any representation.

The CSDM (Coherent Switching with Decay of Mixing) method is a mean-field approach in which the average surface governing the nuclear dynamics decays toward an adiabatic surface, with the decay proceeding faster as the system gets farther from a region of strong coupling. The first-order decay time used to control the demixing of the average surface to a quantized state is defined as

τ = h ∆V c + E 0 T s ( 4
)
where ∆V is the difference between the adiabatic electronic energies, c and E 0 are two parameters, and T s is the kinetic energy associated with the component of the momentum where energy is being added or removed as the trajectory demixes. This demixing algorithm eliminates the main problem encountered in earlier mean-field methods, namely the subjective interpretation of the number of fragments on each adiabatic surface. In the rest of this work we will present results obtained for c = 1, E 0 = 0.1 (the original default values of CSDM) and also with c = 10, E 0 = 0.1 and with c = 1, E 0 = 0.001. (Note that c is unitless, and E 0 is always given in hartrees.) For brevity, CSDM simulations that use these three sets of parameters will be respectively called CSDM(1), CSDM(2), and CSDM(3) for the sake of simplicity.

The FSTU (Fewest Switches with Time Uncertainty) method is a trajectory surface-hopping method that improves the way that classical forbidden hops are taken into account as compared to what is done in the standard TFS (Tully's Fewest Switches) algorithm. 47 In trajectory surface hopping calculations carried out in the adiabatic representation the classical dynamics of the nuclei is governed by one adiabatic surface at a time. To take into account the multisurface character of the problem, a hopping probability is estimated that allows or prevents electronic transitions between surfaces. An electronic transition, also called a hop, can be forbidden at time t 0 if the potential energy of the target surface (surface after hopping) is greater than the conserved total energy. These hops are called classically forbidden hops or frustrated hops. In this case, standard TFS methods omit the hop even if the electronic equations of motion call for it, but the FSTU method looks backward (t h < t 0 ) and forward (t h > t 0 ) in time for a possible hopping time t h that is different from t 0 . In particular, FSTU searches for a phase-space region where a hop would be energetically allowed. The electronic transition is allowed at t h = t 0 if a hopping point is reachable within the Heisenberg interval of time uncertainty, that is when

|t 0 -t h | ≤ h 2∆E (5)
where ∆E is the energy that would need to be borrowed at t 0 to allow a hop.

FSTU/SD (FSTU with Stochastic Decay) is an improved FSTU algorithm in which some decoherence is included. 13 The decoherence time in the FSTU/SD method is modeled as the time for phenomenological decay of the off-diagonal elements of the electronic density matrix due to imagined nuclear wave packets propagating in the different electronic states at different velocities, which finally leads to their dephasing and the separation of their centers. 48 All FSTU and FSTU/SD calculations employ the gradV criterion 49 for adjusting momentum after a frustated hop.

II.E. The TRAPZ and mTRAPZ methods

II.E.1. General description

Many methods have been designed to incorporate quantum effects and especially ZPE corrections in classical trajectory simulations. We focussed our attention on one of these methods, namely the TRAPZ (TRAjectory Projection onto Zero-point energy orbit) method. This method was detailed by McCormack et al. 25 for use in the Cartesian frame, and we only recall here its essential features before presenting differences with the mTRAPZ method. We adopt the following notation: vectors with 3N components, where N is the number of atoms, will be bolded, and 3N × 3N matrices will be put between brackets.

The key element of the TRAPZ method is to perform an instantaneous normal mode analysis at each time step of the dynamics. At the current time t 0 , positions and linear momenta are first mass-scaled according to

x3(i-1)+j = m i µ x 3(i-1)+j (6a) p3(i-1)+j = µ m i p 3(i-1)+j (6b)
where j ∈ {x, y, z} and i ∈ {1, ..., N}, x and p are the Cartesian atomic positions and linear momenta, x and p are the mass-scaled positions and linear momenta, m i is the mass of atom i, and µ is a constant mass taken as 1 amu.

The projected Hessian [K P ] at time t 0 is then determined in mass-scaled Cartesian coordinates from the knowledge of the Hessian [K] and a projector 50,[START_REF] Truhlar | Theory of Chemical Reaction Dynamics[END_REF] [P] by

[K P ] = ([I] -[P])[K]([I] -[P]). (7)
For a nonlinear molecule with i ∈ {1, ..., N}, α, β, γ ∈ {x, y, z}, the projector is given by

P iγ,i ′ γ ′ (t 0 ) = m -1 N H 3 (m i m i ′ ) 1/2 δ γγ ′ (8) + α,α ′ ,β,β ′ ǫ αβγ x iβ [I x 0 (t 0 ) -1 ] αα ′ ǫ α ′ β ′ γ ′ x i ′ β ′
where [I x 0 (t 0 )] is the 3 × 3 inertia matrix estimated at time t 0 and ǫ αβγ is the totally antisymmetric unit pseudotensor of rank 3

ǫ αβγ =          0 if αβγ are not all distinct 1 if αβγ is an even permutation of xyz -1 if αβγ is an odd permutation of xyz (9)
The first term of the projector projects onto the three directions locally corresponding to overall translation, and the second projects onto the three directions corresponding to overall rotation.

The diagonalization of the projected Hessian deduced from Eq. ( 7) then supplies 3N -6 vibrational modes, represented by the eigenvectors L k (k ∈ {1, ..., 3N -6}), and eigenvalues κ k of Eq. ( 7); the instantaneous vibrational frequencies Ω k at time t 0 are given by

Ω k (t 0 ) = κ k (t 0 ) µ (10)
Mass-scaled linear momenta are then transformed into instantaneous normal-mode momenta

P k (k ∈ {1, ..., 3N -6}) by P k = 3N i=1 L i k (t 0 ) pi (11)
where L i k (t 0 ) are the 3N components of the eigenvector L k at time t 0 . The instantaneous vibrational energy of mode k at time t 0 is calculated according to

E k (t 0 ) = 1 2µ P 2 k + D k (t 0 ) Ω k (t 0 ) 2 ( 12
)
where D is the mass-scaled vector of the first derivatives of the potential-energy in instantaneous normal mode coordinates. Note that an instantaneous normal mode analysis provides at most 3N -6 vibrational modes with real non-zero frequencies since some frequencies can be imaginary at some geometries. The number of real non-zero frequencies is therefore 3N -q(t 0 ) with q(t 0 ) ≥

6.

The TRAPZ method consists in comparing E k (t 0 ) for these 3N -q(t 0 ) modes to the local ZPE of mode k, given by

E ZPE,k (t 0 ) = 1 2 hΩ k (t 0 ) (13)
whenever the instantaneous vibrational energy of mode k drops below the local ZPE, the linear momentum P k is adjusted to compel ZPE maintenance

P ′ k = sign(P k ) µ hΩ k (t 0 ) - D k (t 0 ) Ω k (t 0 ) 2 (14)
Assuming that n(t 0 ) modes violate the ZPE conservation, the total energy conservation is ensured by modifying the 3N -q(t 0 ) -n(t 0 ) other modes according to

P ′ k = β P k , with β = 3N -q(t 0 ) k=1 P 2 k - n(t 0 ) k=1 P ′2 k 3N -q(t 0 ) k=n(t 0 )+1 P 2 k (15)
For the q(t 0 ) -6 modes with imaginary frequencies Ω k (t 0 ), we have P ′ k = P k . The momenta are finally backtransformed into Cartesian coordinates according to

p ′ 3(i-1)+j = m i µ p′ 3(i-1)+j (16)
where the 3N components of p′ are given by

p′ l = 3N -q(t 0 ) k=1 L l k (t 0 ) P ′ k + 3N k=3N -q(t 0 )+1 3N j=1 L l k (t 0 ) L j k (t 0 ) pj . ( 17
)
For the present application the second term of Eq. ( 17) is zero because the center of mass is at rest and there is no overall rotation.

This procedure is applied at each time step in the original TRAPZ method. McCormack et al. 25 noticed that the internal energy of Al 2 , the product of the dissociation of Al 3 , that is obtained with this algorithm was peaked at an energy significantly higher than the true ZPE. This raises a problem, namely that constraining all instantaneous mode energies to equal or exceed their local ZPE adds too much vibrational energy at the expense of translational energy. Consequently, the fragmentation may slow down and lead to unphysical results. The mTRAPZ method partly alleviates this shortcoming by only applying the linear momenta transformation when the total instantaneous vibrational energy is below the total local ZPE:

TRAPZ criterion: E k (t 0 ) ≤ E ZPE,k (t 0 ) (18a) mTRAPZ criterion: 3N -q(t 0 ) k=1 E k (t 0 ) ≤ 3N -q(t 0 ) k=1 E ZPE,k (t 0 ) (18b)
This criterion ensures that the system's vibrational energy never drops below the local ZPE while limiting the impact of the original TRAPZ method that tends to pump out away too much vibrational energy.

II.E.2. Technical shortcomings of TRAPZ and mTRAPZ methods

The TRAPZ and mTRAPZ methods have three minor shortcomings. First, the limited accuracy of eigenvectors deduced from the diagonalization of the projected Hessian necessarily leads to a limited accuracy of the transformed momenta. An accumulation of small errors over the course of a trajectory, naturally worse with TRAPZ than mTRAPZ, may lead to a steady loss of the centerof-mass location and linear momentum conservation. Knowing that the projector in Eq. ( 8) is only defined in the center-of-mass frame and that the results are very sensitive to small deviations of the center of mass from the origin, we decided to eliminate the problem by re-initializing the center-of-mass location to the origin and re-zeroing the center-of-mass linear momenta just before applying Eqs. ( 14)- (17) in either TRAPZ or mTRAPZ algorithms.

The second issue arises because the projected Hessian is calculated in Cartesian Coordinates which are not globally suitable for describing rotational eigenvectors. [START_REF] Murry | [END_REF] Curvilinear coordinates would be more appropriate, but are impractically complicated to implement for real polyatomic systems. In the light of this problem, TRAPZ and mTRAPZ were not applied at some identified geometries where the projected Hessian was not properly estimated. This is most commonly encountered when the system is close to the conical intersection or in a region corresponding to the repulsive part of the potential.

Finally, a few trajectories ended with a total angular momentum greater than 0.1h due to accumulated inaccuracies in various numerical steps when the TRAPZ or mTRAPZ methods are applied; these are discarded with no appreciable effect on the results.

II.F. Observables

Several observables were analyzed with special attention to the final NH 2 internal energy which is the key result directly compared to experiments. The final NH 2 internal energies are interpolated by using a Legendre moment method 53 which allows one to generate smooth NH 2 internal energy distributions that can be directly compared to experiments. In the present work, the order of normalized Legendre polynomials only exceeds 18 when using the TRAPZ method due to the very atypical shape of the obtained distribution in Sec. IV A. For a given simulation, distributions are however not sensitive to the maximum order of these polynomials provided this order is not too small (data would be badly interpolated) nor too large (noisy oscillations would appear). This method will also be applied for plotting N-H distance distributions but will not be applied for plotting adiabatic energy gaps in Sec. IV B 2.

The nonplanarity of ammonia has been estimated by computing a nonplanarity angle denoted α np . Since planar NH 3 molecules have three θ i angles (i ∈ {1, 2, 3}) whose sum equals 360 deg, α np has been defined in degrees as follows

α np = 360 - 3 i=1 θ i . ( 19
)
This angle is zero when the molecule is planar and positive otherwise.

II.G. Computational details

The coupled-state mixed quantum/classical dynamics calculations were carried out with the ANT08 code, [START_REF] Ant | Adiabatic and Nonadiabatic Trajectories[END_REF] which was updated to include the TRAPZ and mTRAPZ algorithms for the present study. The equations of motion were integrated using a spacialized algorithm 46 presented previously, with a step size of 1 fs. Except where specified otherwise, each reported simulation corresponds to running five batches of 1000 trajectories. These trajectories stop either after an integration time of 100 ps, which is the simulation time limit, or when a hydrogen atom is detected at 10 Å (at least) from the nitrogen atom of the remaining NH 2 fragment during 50 successive time steps; the final time of the trajectories will be denoted t f .

The main experimental result obtained from the study of the photodissociation of ammonia [30][31][32] is the NH 2 internal energy. We computed the NH 2 internal energy by Eq. ( 3). Results obtained by performing mixed quantum/classical simulations with the CSDM(1), CSDM(2), CSDM(3), FSTU, and FSTU/SD methods for n 2 = 0 are reported in Fig. 2 where they are compared to available experimental distributions. 30,32 The first striking point, as shown in Fig. 2(a), is that distributions are quite similar for whatever mixed quantum/classical method is used; in particular, the energy distribution is bimodal with a main peak at 0 eV and another one at about 0.8 eV.

These distributions are in qualitative agreement with the experimental distributions of Biesner This important problem will be discussed in Sec. IV.

We also note, as discussed in Sec. II C, that the excess energies of the theoretical distributions can be slightly higher than theoretical excess energies specified in Table I. Furthermore, as discussed in Sec. II B, a single trajectory may have an energy slightly different than the nominal one. We have checked that using E α tot,0 (different for each trajectory α) instead of E tot,0 has no significant impact on NH 2 internal energy distributions.

Note that we have plotted FSTU/SD distributions but not CSDM and FSTU ones in the comparison to experiments. This choice was made for three reasons: (i) FSTU/SD will be shown in Sec. III C to lead to less frustrated hops and to be less time consuming than FSTU, and we will consequently carry out calculations with ZPE maintenance only for this method in Sec. IV, (ii) we wanted to avoid overloading the figures, and (iii) Fig. 2(a) shows that CSDM, FSTU, and FSTU/SD calculations yield similar dynamics.

When the umbrella mode is excited (n 2 > 0) the general shape of the distribution becomes unimodal with one peak that steadily moves toward higher energies (as n 2 is increased) as shown in Fig. 3 for FSTU/SD simulations. (A comparison of results obtained for CSDM(1), CSDM(2), CSDM(3), FSTU, and FSTU/SD simulations is presented in Fig. S1 of supporting information.)

The same behavior is observed experimentally in Fig. 3 but the experimental peak is more pronounced close to the maximum excess energy which means that the hydrogen atoms are on average colder in the experiment than in our mixed quantum/classical simulations. This is a signature of the classical overestimation of the excess energy. Also note that the experimental distribution obtained for n 2 = 2 has a peak at low energies. This peak was assigned 30 to very low rotational states of ammonia. We do not observe such a peak in our mixed quantum/classical simulations.

III.B. Dissociation mechanism

Whatever the mixed quantum/classical method and the number of quanta of umbrella-mode vibration, the dissociation was found to be mainly electronically nonadiabatic as found by experimentalists. [30][31][32] The proportion of trajectories that dissociate on adiabatic surfaces 1 and 2 are reported in Table II. In all cases, at least 92.6% of the trajectories end in the ground electronic state (adiabatic state 1), with the average over all calculations for various initial states and methods being 97.0%.

Among the trajectories that dissociate on surface 1 or 2, some of them can be considered as directly dissociating. Within the framework of the two surface-hopping methods tested here (FSTU and FSTU/SD), we consider that a trajectory starting in the first excited state (surface 2) directly dissociates into the ground electronic state (surface 1) when there is only one nonadiabatic transition or hop before dissociation, and a trajectory starting in the first excited state directly dissociates into the excited electronic state when there is no hop before dissociation. In this last case, the dissociation can be considered purely adiabatic. We found that at least 65% of trajectories that dissociate in the ground electronic state are direct but at most 34% of trajectories that dissociate in the excited state are (see Table S1 of supporting information). Most of trajectories that are found to dissociate in the excited state thus undergo two or more hops before dissociating. The FSTU and FSTU/SD methods are in qualitative agreement with each other, and we will now focus our attention on the efficiency and reliability of all the mixed quantum/classical methods we have tested on ammonia.

A fundamental source of error in TFS-like methods is the occurrence of classically forbidden hops. A great number of classically forbidden hops can lead to discrepancies when comparing to quantum or experimental results. 10,13,46,[START_REF] Hack | [END_REF][56][57] Table III compares the fraction of accepted hops for the FSTU and FSTU/SD methods. More than 50% of attempted hops are accepted within the FSTU/SD calculations, but this number drops to less than 10% when the FSTU method is used.

Note that the percentage of accepted hops is also about 50% when ZPE maintenance is included in FSTU/SD. The FSTU/SD method should be more accurate than the FSTU method in general 13 even if there is no noticeable difference in final results for ammonia. The FSTU/SD method is also about 1.6 times less time consuming than the FSTU method which provides another reason to prefer this recently developed method over FSTU.

III.C.2. The CSDM method

The three versions of the CSDM method that are tested here lead to similar results. However, the choice of parameters can influence the computational cost of the method. The average computational cost of 5000 trajectories run on the à electronic state of ammonia with the CSDM (1) method is about 36 CPU hours on a single 1.3 GHz IBM Power 4 microprocessor, this number has to be multiplicated by about 1.2 with the CSDM(2) method and by 2 with the CSDM(3) method.

Although the CSDM(1) method is about 1.4 times less time consuming than the FSTU/SD method, this method will not be applied in Sec. IV because the TRAPZ and mTRAPZ methods are not designed to ensure the ZPE maintenance on average surfaces but rather on uncoupled surfaces, and trajectory surface hopping methods such as FSTU and FSTU/SD provide more natural algorithms for combining with TRAPZ and mTRAPZ.

The goal 23 of the TRAPZ method is to maintain the ZPE, and it was designed to remove the instabilities introduced by previous methods for doing this. We have consequently applied the TRAPZ method during a FSTU/SD simulation for n 2 = 0, 1, and 6. Distributions obtained in Fig. 4 reveal that the TRAPZ calculations fail to reproduce even the qualitative character of the experimental distributions. The distributions produced when using the TRAPZ algorithm are peaked at energies close to the maximum excess energy, and, have exactly the same shape for n 2 = 0, 1, and 6. The large amount of final vibrational energy is not physical and is due to the fact that the TRAPZ method tends to pump translational energy into NH 2 vibrational energy to an excessive degree. The consequence of this unphysical behavior is that fragmenting H atoms are very cold. Note that we have only run 2000 trajectories for each simulation in Fig. 4 since simulations including TRAPZ are four to thirty times more time consuming than the original FSTU/SD simulations, with the precise timing ratio depending on the umbrella mode excitation.

IV.B. Simulations with the mTRAPZ method

IV.B.1. NH 2 internal energies

As discussed in Sec. II E, the mTRAPZ method is based on the same methodology as the TRAPZ method but is less constraining due to a more tolerant criterion for modifying momenta (see Eq. (18b)). The distribution obtained after performing a FSTU/SD+mTRAPZ simulation for n 2 = 0 is reported in Fig. 5. The main peak is at about 0.4 eV from the minimum NH 2 internal energy and the qualitative agreement with experiment is better than that obtained for FSTU/SD without ZPE maintenance since the ZPE maintenance prevents negative NH 2 internal energies.

We compare the partitioning of internal and translational energy in Table IV. For n 2 = 0, 49% of the excess energy flows into translation in our simulation, and the experiment finds a value of 53 ± 2%. The partitioning of the energy is very satisfactory with an average deviation from experiment of less than 3% (The results labeled mTRAPZ* in this table will be explained below).

Figure 6 presents mTRAPZ energy distributions obtained for n 2 > 0 and compares them to distributions obtained with neither TRAPZ nor mTRAPZ. The result obtained with the mTRAPZ algorithm for n 2 = 6 in Fig. 6 is completely different from the TRAPZ result in Fig. 4. The TRAPZ algorithm changes the momenta to increase some mode energies after most time steps even when the initial n 2 is 6; this shows that a high initial vibrational energy does not prevent the energies of individual modes from dropping below their local ZPE. However, we found that in the same kind of situation, the mTRAPZ algorithm will rarely change the momenta since it is less likely that the total instantaneous vibrational energy drops below the total local ZPE. FSTU/SD+mTRAPZ distributions consequently tend to resemble FSTU/SD simulations for high initial vibrational energies.

Although mTRAPZ distributions obtained for n 2 ≥ 2 both maintain the ZPE and are qualitatively close to mixed quantum/classical distributions reported in Fig. 3, the mTRAPZ distribution obtained for n 2 = 1 maintains the ZPE but disagrees with experiment. We also tested another TRAPZ version that is less constraining than mTRAPZ and that improves the results obtained for n 2 = 1. This method consists in applying the TRAPZ method when and only when the instantaneous vibrational energy is below the product ZPE, using E ZPE [NH 2 ( X)] = 0.519 eV as the zero-point requirement at all times (see Fig. 1). We call this method mTRAPZ* to keep in mind that it is a simplification of the mTRAPZ method. Results obtained with the mTRAPZ* method are reported in Fig. 7 for 0 ≤ n 2 ≤ 3 (which are the quantum numbers for which the effect of ZPE maintenance is most important as shown in Fig. 5 and Fig. 6). The ZPE is maintained as in mTRAPZ but distributions are smoother and closer to experimental results. However, mTRAPZ* is only suitable here since the ZPEs of NH 2 ( Ã) and NH 2 ( X) are quite close to each other with that for NH 2 ( X) lower and with both product values considerably lower than the reactant value.

It would be difficult to apply if several types of molecular fragments could be produced by dissociation on different adiabatic surfaces and would seem not to enforce the intended physics if the local ZPE throughout the dynamics was often lower than the product ZPE. In the present cases, though, the partitioning of energy in Table IV has the same level of accuracy for mTRAPZ and mTRAPZ*. In particular, neither result is fully satisfactory for n 2 = 1, but the agreement with experiment for n 2 > 1 is better with mTRAPZ*. The mTRAPZ and mTRAPZ* methods are also comparable in terms of computational cost; they are 2 to 3 times more time consuming than FSTU/SD simulations without ZPE maintenance.

IV.B.2. Dissociation dynamics

The dynamics is still found to be highly electronically nonadiabatic with mTRAPZ. In fact, it is even more electronically nonadiabatic than without ZPE maintenance; in particular, all the trajectories dissociate in the ground electronic state for n 2 = 0, 1, and 2. This is an improvement compared to the mixed quantum/classical results reported in Table II since the adiabatic dissociation is energetically forbidden for n 2 = 0, 1, and 2 in the experiment. (For n 2 = 2 the adiabatic dissociation is not energetically forbidden in our simulation, but the energy difference between the average initial total energy E tot,0 and the quantum energy of the à electronic state at dissociation is only about 0.08 eV.) The mTRAPZ method therefore allows one to overcome one of the most troublesome limitations of classical mechanics by preventing dissociations that are energetically forbidden in the quantum world.

For 3 ≤ n 2 ≤ 6, the amount of adiabatic dissociation in the FSTU/SD+mTRAPZ simulation never exceeds 0.1% whereas the experimentalists have found that 80 ± 2% of the excess energy flows into internal energy for n 2 = 6 and, of that 80%, 54 ± 4% corresponds to internal energy of rovibrationally excited NH 2 ( X) and 26 ±4% to internal energy of electronically excited NH 2 ( Ã). 30 Putting 26% of the excess energy into electronic excitation energy, coupled with the nature of the vibrational energy distributions of product molecules formed in state 2 (see Fig. 2b in Ref. 30) implies that about 28 to 30% of the amino radicals are produced in the à state when n 2 = 6. Thus, our simulations underestimate the amount of electronic adiabaticity for systems with high n 2 . This discrepancy with experiment is not apparently due to a weakness in the mTRAPZ algorithm since the amount of adiabatic dissociation does not exceed 7% even in various simulations without ZPE maintenance (see Table II).

The percentage of trajectories with N hop allowed hops (N hop ≥ 0) is reported in Table V for the FSTU, FSTU/SD, FSTU/SD+TRAPZ, FSTU/SD+mTRAPZ, and FSTU/SD+mTRAPZ* calculations when n 2 = 0 and n 2 = 6. Knowing that the dynamics starts on the à electronic state, an even number of hops leads to adiabatic dissociation, and an odd number of hops leads to nonadiabatic dissociation. All the simulations reveal a preference for nonadiabatic dissociation and the adiabatic dissociation almost disappears when ZPE maintenance is ensured. When n 2 increases, we expect a decrease in the number of trajectories that dissociate directly (i.e., with N hop = 1) into the ground electronic state (see Sec. III B for the definition of direct trajectories) due to the larger probability for upward hops. The number of such direct trajectories indeed decreases from 91.7% to 69.3% for the FSTU/SD+mTRAPZ calculations, but this decrease does not increase the number of trajectories with an even number of hops although the adiabatic dissociation channel is energetically allowed for n 2 ≥ 2 in our simulations. We can also notice that the number of trajectories with an odd number of hops is much closer to FSTU and FSTU/SD results for FSTU/SD+mTRAPZ and FSTU/SD+mTRAPZ* simulations than it is for FSTU/SD+TRAPZ simulations, which is another indication of the preferability of mTRAPZ as compared to TRAPZ.

In order to better understand the physical origin of the results in Table V, we performed further analysis of the FSTU/SD+mTRAPZ trajectories for the subset of trajectories with N hop = 3 for n 2 = 6. We found that the average time to the first hop is 39 fs; the average time between the first and the second hop is 151 fs; and the average time between the second and the third hop is 19 fs. Thus, the first downward hop occurs rather quickly and, when the system hops up to the excited state, it tends to return to the ground state still more quickly than for the first downward hop. (The first downward hop is delayed by the time required to reach the region where the surfaces are strongly coupled.) This subset of trajectories with N hop = 3 has, on average, 2.1 frustrated hops, which are (as always) frustrated upward hops since downward hops are never frustrated. These characteristics tend to produce products on the lower surface.

The final time of the dynamics, denoted t f , as explained in Sec. II G, is larger when the umbrella mode is more excited; average values of t f are t f = 686 ± 1171 fs for n 2 = 0 and t f = 1008 ± 4430 fs for n 2 = 6. (The fact that the standard deviations exceed the mean is an indication that the distributions are very skewed.) Note that two trajectories with ZPE maintenance (mTRAPZ method) were found to reach the maximum time limit for n 2 = 6; they were removed and, therefore, the averages of t f are not affected by this limit. Unlike t f , the time at the last downward hop, denoted t hop , is smaller when the umbrella mode is more excited; average values and standard deviations of t hop are t hop = 135±339 fs for n 2 = 0 and t hop = 110 ± 205 fs for n 2 = 6. Comparing the average times of the last downward hop shows that a higher initial excitation of the umbrella mode allows the system to more rapidly deexcite. (The t hop distributions obtained for n 2 = 0 and n 2 = 6 are plotted in Fig. S2 of supporting information.)

The typical adiabatic energy gap at the last downward hop is ∼ 0.3 eV for n 2 = 0 and ∼ 0.4 eV for n 2 = 6, which indicates that most of the electronic transitions occur in a region close to the conical intersection (see Table S2 of supporting information). Figure 8 shows the adiabatic energy gap distribution at the last downward hop. Distributions obtained for n 2 = 0 and n 2 = 6 are peaked at 0 eV although the latter distribution spreads to larger energies which demonstrates that more electronic transitions occur in regions far from the conical intersection when n 2 increases.

To shed some light on the important role of the conical intersection as a function of n 2 values, we have studied ammonia geometries at the last downward hop. First of all, the minimum and maximum N-H distance distributions are plotted in Fig. 9 for n 2 = 0 and n 2 = 6. The minimum N-H distance provides information about the typical N-H bond length in the incipient NH 2 fragment, and the maximum N-H distance is the dissociation coordinate. The maximum N-H distance is peaked at about 2 Å when n 2 = 0 which corresponds to the conical intersection on the present set of coupled potential-energy surfaces, 45 and the average maximum N-H distance is 2.05 ± 0.21 Å.

The average number of allowed hops by trajectory is 1.2 ± 0.8 which confirms that most of the trajectories are direct (see Table V). When n 2 increases a second peak appears close to 3 Å, and the average maximum N-H distance increases to 2.46 ± 0.56 Å, which shows that some transitions occur farther from the conical intersection. This bimodal distribution is also accompanied by a drop in the number of direct trajectories to 69.3% (see Table V). For n 2 = 6, the average number of allowed hops is about 1.6±1.4 when the maximum N-H distance is below 2.5 Å (first peak) and increases to 2.5 ± 2.5 (second peak) when the maximum N-H distance exceeds 2.5 Å. Therefore, the larger is the number of hops, the larger is the maximum N-H distance, on average. N-H distances are howewer not sufficient to determine the proximity to a conical intersection geometry; one must also consider H-N-H angles. We have therefore studied the nonplanarity angle α np defined in Sec. II F. According to Ref. 45, the lowest-energy conical intersection is planar for ammonia; the coefficient of planarity at this critical geometry is therefore zero. In our simulations, the average α np at the start of the trajectory equals 2.1 ± 1.5 deg and 17.6 ± 17.1 deg for n 2 = 0 and n 2 = 6 respectively. At the last downward hop, these values increase to 5.0 ± 4.3 deg and 36.3 ± 35.8 deg, respectively. Thus we see that the last downward hop occurs for near-planar geometries when the initial n 2 = 0. When n 2 = 6, we find that hops occur, on average, farther from the conical intersection, which may be expected since exciting the out-of-plane bending mode favors nonplanar geometries. We may compare the value of 36 deg to value of α np = 32 deg for a regular tetrahedron; this shows that, on average, the geometry of a downward hop when the initial n 2 = 6 is closer to tetrahedral than to planar. We note, however, that many ammonia molecules are close to planar at the last downward hop for n 2 = 6; the distribution of α np values at the last downward hop is peaked at α np ≈ 0 deg and is very broad for n 2 = 6.

V. CONCLUSION

We have studied the photodissociation of ammonia in the ground vibrational state of its à excited electronic state and also after excitation of 1 to 6 quanta of the umbrella mode of the à state. We report several interesting findings:

1. The coupled-state dynamics are very similar with three different mixed quantum/classical methods: coherent switching with decay of mixing (CSDM), fewest switches with time uncertainty (FSTU), and FSTU with stochastic decay (FSTU/SD); and for CSDM the results do not depend strongly on the parameters of the decoherence time. This, plus the successful validation of these mixed quantum/classical trajectory methods in previous studies, [9][10][11][12][13][14][15]58 gives us confidence in the essential correctness of all three methods. Production runs were then made with FSTU/SD.

2. We developed a satisfactory method, called minimal trajectory projection onto zero-point energy orbit (mTRAPZ), for maintaining zero-point vibrational energy in the mixed quantum/classical trajectories, whereas uncorrected trajectories give unphysical final states with less than zero-point energy, and the original TRAPZ algorithm gives highly vibrationally excited products in severe disagreement with experiment.

3. The FSTU/SD+mTRAPZ mixed quantum/classical calculations primarily produce dissociation products in the ground electronic state, in good agreement with experiment, although the amount of electronically adiabatic dissociation seems underestimated for high n 2 . The reason for this underestimation is not understood yet.

4. The FSTU/SD+mTRAPZ calculations place 50% (n 2 = 0) to 71% (n 2 = 6) of the excess energy of ground-electronic-state product in internal energy (with the rest in translational energy), where n 2 is the initial vibrational quantum number of the umbrella mode. This is in agreement with experiment, where this percentage rises from 47% at n 2 = 0 to 80% at n 2 = 6.

5. The qualitative success of the FSTU/SD+mTRAPZ calculations in accounting for the experimental observables of the ground electronic-state products gives us some confidence in the predictions for unobserved quantities. These include the average time of quenching, 135 fs (with a large positive tail on the distribution) for n 2 = 0 and 110 fs (again with a long positive tail) for n 2 = 6, and the average adiabatic gap at quenching, 0.3 eV for n 2 = 0 and 0.4 eV for n 2 = 6. We find that the out-of-plane bend angle at quenching is correlated with the initial n 2 . Even though increasing the initial vibrational excitation opens the adiabatic dissociation channel, the system is found to reach strong coupling regions more rapidly and to dissociate electronically nonadiabatically despite possible upward hops. The last downward hop is finally found to occur on average for ammonia geometries close to the conical intersection, at least for n 2 = 0. In particular, the average bond length of the breaking N-H bond is 2.05 Šfor n 2 = 0 and 2.46 Šfor n 2 = 6, as compared to 1.99 Šfor the lowestenergy conical intersection (versus 1.04 Šfor the equilibrium geometry in the à electronic state of ammonia).

The present study has shed some light on the photodissociation dynamics of ammonia upon electronic and vibrational excitation and its possible treatment by means of mixed quantum/classical methods. The present study only concerned the out-of-plane bending mode, and the effect of exciting symmetric and antisymmetric stretches is under study in our group. The state selectivity of the adiabatic and nonadiabatic dissociation pathways will be investigated to compare with most recent experiments of Hause et al. 32 in order to more completely understand the dissociation mechanisms of photoexcited ammonia, a prototype system for studying the state-selected dynamics of a system with a conical intersection. obtained by FSTU/SD simulations (red curve) compared to the experimental distribution of Hause et al. 32 (black curve). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

FIG. 3: (Color online). NH 2 internal energy distributions at the end of FSTU/SD simulations without ZPE maintenance (red curves) compared to experimental results of Biesner et al. 30 (black curves) for 1 ≤ n 2 ≤ 6. The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison. 32 (black curve) and FSTU/SD simulations (red curve). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

FIG. 6: (Color online). NH 2 internal energy distributions at the end of FSTU/SD+mTRAPZ simulations for n 2 > 0 (green curves) compared to experiments (black curves) 30 and FSTU/SD simulations (red curves). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

FIG. 7: (Color online). NH 2 internal energy distributions at the end of FSTU/SD+mTRAPZ* simulations for 0 ≤ n 2 ≤ 3 (green curves) compared to experiments (black curves) 30 and FSTU/SD simulations (red curves). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison. Each curve is normalized in such a way that the maximum of the distribution is unity. compared to the experimental distribution of Hause et al. 32 (black curve). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison. The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison. 32 (black curve) and FSTU/SD simulations (red curve). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

37

FIG. 6: (Color online). NH 2 internal energy distributions at the end of FSTU/SD+mTRAPZ simulations for n 2 > 0 (green curves) compared to experiments (black curves) 30 and FSTU/SD simulations (red curves).

The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

  et al. and Hause et al. (see Fig. 2(b) and Fig.2(c)). However they have a tail in the negative energy range, and this behavior is unphysical. We can understand the origin of this tail in terms of classical dynamics: there is no constraint in these mixed quantum/classical calculations to ensure ZPE maintenance, and the departing H atom may take away more than the quantum excess energy.

FIG. 1 :

 1 FIG.1:(Color online). Energy as a function of the largest N-H bond distances for the two lowest adiabatic surfaces based on fitted adiabatic surfaces from Ref.45. The three N-H distances are denoted r i (i ∈ {1, 2, 3}). The minimum of the à electronic state is planar withr 1 = r 2 = r 3 = 1.039 Å.The saddle point is planar with r 1 = 1.289 Å and r 2 = r 3 = 1.040 Å, and the lowest-energy conical intersection is also planar. The minimum of the X state is tetrahedral with r 1 = r 2 = r 3 = 1.016 Å. The zero of energy is the minimum of the X ground electronic state. Reactant (i.e., NH 3 ) zero-point and vibrational excitation energies are indicated in red and product (i.e., NH 2 ) zero-point energies in green. Note that the diagram is only schematic (not to scale) between 1 eV and 5 eV and beyond 2 Å to make relevant quantities more discernable.

FIG. 4 :

 4 FIG.4:(Color online). NH 2 internal energy distributions at the end of FSTU/SD+TRAPZ simulations for n 2 = 0 (black curve), n 2 = 1 (red curve) and n 2 = 6 (green curve). Each curve is normalized in such a way that the maximum of the distribution is unity.

FIG. 8 :

 8 FIG. 8: (Color online). Adiabatic energy gap distributions at the last downward hop for n 2 = 0 (black bars) and n 2 = 6 (red bars) from FSTU/SD+mTRAPZ simulations. Each curve is normalized in such a way that the maximum of the distribution is unity.

FIG. 1 : 33 FIG. 2 :

 1332 FIG. 1: (Color online). Energy as a function of the largest N-H bond distances for the two lowest adiabatic surfaces based on fitted adiabatic surfaces from Ref. 45. The three N-H distances are denoted r i (i ∈ {1, 2, 3}). The minimum of the à electronic state is planar with r 1 = r 2 = r 3 = 1.039 Å. The saddle point is planar with r 1 = 1.289 Å and r 2 = r 3 = 1.040 Å, and the lowest-energy conical intersection is also planar. The minimum of the X state is tetrahedral with r 1 = r 2 = r 3 = 1.016 Å. The zero of energyis the minimum of the X ground electronic state. Reactant (i.e., NH 3 ) zero-point and vibrational excitation energies are indicated in red and product (i.e., NH 2 ) zero-point energies in green. Note that the diagram is only schematic (not to scale) between 1 eV and 5 eV and beyond 2 Å to make relevant quantities more discernable.

34

 34

FIG. 3 :

 3 FIG.3:(Color online). NH 2 internal energy distributions at the end of FSTU/SD simulations without ZPE maintenance (red curves) compared to experimental results of Biesner et al.30 (black curves) for 1 ≤ n 2 ≤ 6.

6 FIG. 4 :

 64 FIG.4:(Color online). NH 2 internal energy distributions at the end of FSTU/SD+TRAPZ simulations for n 2 = 0 (black curve), n 2 = 1 (red curve) and n 2 = 6 (green curve). Each curve is normalized in such a way that the maximum of the distribution is unity.

38

 38

FIG. 7 :FIG. 8 :

 78 FIG.7:(Color online). NH 2 internal energy distributions at the end of FSTU/SD+mTRAPZ* simulations for 0 ≤ n 2 ≤ 3 (green curves) compared to experiments (black curves)30 and FSTU/SD simulations (red curves). The maxima of the experimental curves is normalized to 1, and the maxima of the theoretical curves to 0.5 for ease of comparison.

TABLE V :

 V Percentage (%) of trajectories with N hop allowed hops for the FSTU, FSTU/SD, FSTU/SD+TRAPZ, FSTU/SD+mTRAPZ, and FSTU/SD+mTRAPZ* methods for n 2 = 0 and n 2 = 6. n 2 N hop FSTU a FSTU/SD a FSTU/SD+TRAPZ FSTU/SD+mTRAPZ FSTU/SD+mTRAPZ*

	0	0.	< 0.1	0.	0.	0.
	1	91.4	92.5	89.4	91.7	91.7
	2 < 0.1	0.1	0.	0.	0.
	0 3	6.4	5.5	9.4	6.5	6.4
	4 < 0.1	< 0.1	0.	0.	0.
	5	1.6	1.4	1.0	1.5	1.4
	≥ 6 0.6	0.5	0.2	0.3	0.4
	0	1.4	1.7	0.	0.	0.
	1	61.6	64.2	81.9	69.3	69.4
	2	2.0	1.7	0.	0.	0.
	6 3	20.4	20.2	13.1	19.8	20.4
	4	0.9	0.8	0.	0.	0.
	5	7.8	6.3	3.5	7.1	6.7
	≥ 6 5.9	5.0	1.6	3.7	3.6

a without ZPE maintenance.