12016
J. Phys. Chem. A 2010, 114, 12016–12020
Controlling Product Selection in the Photodissociation of Formaldehyde: Direct Quantum Dynamics from the S1 Barrier Marta Arau´jo,† Benjamin Lasorne,*,‡ Alexandre L. Magalha˜es,† Michael J. Bearpark,§ and Michael A. Robb§ REQUIMTE, Departamento Quimica e Bioquimica, Faculdade de Ciencias, UniVersidade do Porto, 4169-007 Porto, Portugal, Department of Chemistry, Imperial College London, London SW7 2AZ, United Kingdom, and Institut Charles Gerhardt Montpellier, UMR 5253, CNRS-UM2-ENSCM-UM1, CTMM, UniVersite´ Montpellier 2, CC 1501, Place Euge`ne Bataillon, 34095 Montpellier, France ReceiVed: October 5, 2010
Controlling the selectivity between H2+CO and H+HCO in the S1/S0 nonadiabatic photodissociation of formaldehyde has been investigated using direct quantum dynamics. Simulations started from the S1 transition state have suggested that a key feature for controlling the branching ratio of ground-state products is the size of the momentum given to the wavepacket along the transition vector. Our results show that letting the wavepacket fall down from the barrier to the conical intersection with no initial momentum leads to H2+CO, while extra momentum toward products favors the formation of H+HCO through the same conical intersection. Quantum dynamics results are interpreted in semiclassical terms with the aid of a Mulliken-like analysis of the final population distribution among both products and the reactant on each electronic state. Introduction 1
Despite having been thoroughly studied for decades, recent literature shows that there is a renewed interest in formaldehyde photochemistry.2-7 Scheme 1 intends to be a summary of all mechanisms that have been proposed so far for the dissociation of H2CO to H+HCO or H2+CO in the ground singlet state. They are presented in ascending order of photoexcitation energy, and involve either intersystem crossing (triplet-singlet crossing, T1/S0 X),internalconversionthrougharealcrossing(singlet-singlet, S1/S0 conical intersection (CoIn)), or an S1/S0 Fermi golden rule (FGR) type of radiationless decay to S0 at finite energy gap. This work focuses on the higher energy end of the scheme (highlighted in blue), where dynamics are determined by the bottleneck in the excited state, S1 TS.4,7 S1 TS is relatively high in energy compared to the FC point, which implies that FCf S1 TS f C1 CoIn is not a significant pathway when low excitation energies are used experimentally.6 However, the probability of reaching the region of the S1 TS from FC is not zero and should increase with added excess energy. Our aim in this work was to verify that product selection could be achieved once the barrier had been reached. We performed direct quantum dynamics computations, which confirm that giving extra momentum to S1 TS with various magnitudes and signs along the transition vector (TV) direction controls the ground-state product branching ratio H2+CO: H+HCO. Computational Details All excited-state stationary points (S0 min, S1 min, and S1 TS), S1/S0 conical intersections (minimum-energy C1 CoIn and Cs CoIn, a saddle point in the intersection space) and the S1/S0 * To whom correspondence should be addressed. E-mail: benjamin.
[email protected]. † Universidade do Porto. ‡ Imperial College London. § Universite´ Montpellier 2.
SCHEME 1: Mechanisms Proposed in the Literature for the Dissociation of Formaldehyde, Following Photoexcitation to the S1(n-π*) Statea
a
The mechanism studied in this paper is shown in blue.
seam MEP (minimum-energy path connecting Cs CoIn to C1 CoIn within the intersection space8) represented in Figure 1 have been computed and characterized as described in ref 4. Energetics is recalled in Table 1. Quantum trajectories have been computed using the direct dynamics variational multiconfiguration Gaussian (DD-vMCG) wavepacket method9 implemented in a development version of the Heidelberg MCTDH package.10 Potential energy surfaces have been calculated on-the-fly along those coupled quantum trajectories via an interface with a development version of the GAUSSIAN program,11 using the complete active space selfconsistent field (CASSCF) method with a reduced SA-2CAS(8,7) active space7 (SA-2 standing for S1/S0 state-averaged orbitals) and a 6-31G* basis set to compute analytical gradients and Hessians. The S1 and S0 coupled electronic states have been represented in a diabatic picture, as explained in ref 7, using FC and C1 CoIn as references for the diabatization. In all
10.1021/jp109549r 2010 American Chemical Society Published on Web 10/27/2010
Controlling Product Selection in H2CO Photodissociation
J. Phys. Chem. A, Vol. 114, No. 45, 2010 12017
Figure 1. 4-GBF case: quantum trajectories plotted in the frame of the HCH angle and the difference between the two CH bond lengths (gray dotted line: from FC with zero momentum; red line: from S1 TS with zero momentum; blue line: from S1 TS with extra momentum corresponding to an extra kinetic energy of 0.1 eV along the TV toward product; green line: from S1 TS with extra momentum corresponding to an extra kinetic energy of 0.1 eV along the TV toward reactant). Top right insert: potential energy profiles of the same trajectories (an identical color scheme is used), as a function of the HCH angle. Black solid line: extrapolated direction of the TV. Black dashed line: S1/S0 seam minimum-energy path connecting Cs CoIn to C1 CoIn.4 H11 and H22 label the two diabatic potential energy profiles (diagonal elements of the diabatic potential energy matrix), whereas H12 labels the potential-like coupling profile between both diabatic states (off-diagonal element of the diabatic potential energy matrix).
TABLE 1: Energies of the Most Significant Critical Points Given in kcal mol-1 above S0 mina SA-2-CASSCF(8,7)/6-31G* CASSCF(10,9)/6-31G* S0 S0 min S1 min S1 TS S0/S1 C1 CoIn S0/S1 Cs CoIn
(-113.963265) 23.6 82.5 106.2 141.1
S1
S0
S1
104.1 91.0 113.6 106.2 141.1
(-113.996199) 35.6 81.1b 107.4b 144.8b
103.4 88.5 116.0b 107.4b 144.8b
a Absolute energies in hartree in brackets. b Optimized with S1/S0 state-averaged orbitals and no CP-MCSCF (coupled perturbed multiconfiguration self-consistent field) contribution to the gradient (zero correction at CoIn but not at TS).
simulations, the global wavepacket has been expanded in a timedependent basis set of up to eight Gaussian basis functions (GBF), making it possible for the wavepacket to split into multiple components after crossing the S1/S0 seam of intersection. Some simulations were started either from the FC point or from S1 TS with no excess energy. In those cases, momentum developed from an initial zero value as the wavepacket evolved. Note that a wavepacket on a constant potential energy surface will spread symmetrically. In classical terms, the corresponding Wigner distribution contains initial conditions with nonzero velocities equally distributed with respect to both signs. A wavepacket centered on a barrier will thus spread on both sides even faster than on a flat surface, with a preference for the side with the most negative curvature if the barrier is not symmetrical. Other simulations were started from S1 TS where a nonzero mean momentum was given to the initial wavepacket. In all cases, the Mulliken-type population analysis proposed in ref 12 has been used to calculate the gross Gaussian population (GGP). This quantity directly gives a measure of the contribution of a given GBF to the global density of probability in each electronic state and can be seen as the weight of the quantum trajectory followed by the center of each GBF. Depending on the region where each GBF terminates its course,
the final value of the GGP can be used to provide the probability to either regenerate the reactant or form H2+CO or H+HCO in each electronic state. Results and Discussion In this section we show with direct quantum dynamics simulations that giving extra momentum to S1 TS with various magnitudes and signs along the TV direction can alter the ground-state product branching ratio H2+CO:H+HCO. The most representative results obtained using a basis set of 4 GBF and of 8 GBF are shown in Figures 1 and 3, respectively. The smaller basis set was used mainly for testing mechanistic ideas, whereas 8 GBF were used to make it possible for the wavepacket to split into multiple components after crossing the S1/S0 seam of the intersection. In those figures, quantum trajectories (i.e., trajectories followed by the centers of the GBF used to expand the wavepacket) are plotted in the frame of two geometrical parameters that, when combined, allow for a clear distinction between reactants and each type of products: an internal coordinate that directly conditions the proximity of the two hydrogens, the HCH angle, and the difference between the two CH bond lengths, which reflects the possible symmetry breaking of the system. This distinction is further illustrated by the potential energy profiles (plotted against the HCH angle) of all trajectories included in Figure 1 (see insert), where three wells (reactants/H2+CO/ H+HCO) can be noted. Progress from “reactant region” into “product region” is characterized by angle tightening. At the same time, the difference between bond lengths goes to zero if molecular products are formed (since the two CH bonds stretch to infinity at a similar rate), whereas it keeps increasing when dissociation to H+HCO takes place. Before discussing results originating from S1 TS, let us verify whether a wavepacket originating from the FC point with no excess energy can reach the S1 transition barrier. The gray dotted line in Figure 1 represents the trajectory of the basis function
12018
J. Phys. Chem. A, Vol. 114, No. 45, 2010
Figure 2. (a) The S1 TS TV (pointing toward products). (b) The derivative coupling vector for C1 CoIn.
that passes the closest to the S1 TS. The corresponding GGP is 13%. Although S1 TS is relatively high in energy compared to the FC point, after several oscillations for 297 fs around S1 min, the trajectory reaches a geometry that differs by an average 4.5% relative to the S1 TS internal coordinates and by 0.4% relative to its excited-state potential energy. The calculation level used in this work underestimates the energy of the barrier (10 kcal mol-1) compared to more accurate descriptions (13 kcal mol-1 with a CAS(10,9)/6-31G* calculation;4 19 kcal mol-1 with a high-level MRCISD(Q)/AV5Z calculation6). This makes passing the S1 transition barrier happen faster in our simulation for artificial reasons. It is thus to be understood that the time scale shown here is not to be compared to experiments when low excitation energies are used.6 However, the probability of reaching the region of the S1 TS from FC is not zero and should increase with added excess energy. Now, let us focus on cases originating from S1 TS. Each trajectory in Figure 1 was chosen as the most representative out of four for each case. From S1 TS with extra momentum along the TV toward reactant, all basis functions regenerate the reactant (total S1:S0 ) 94%:6%). The trajectory represented (green line) corresponds to the largest final GGP equal to 39% of the global wavepacket (almost exclusively on S1). In the case of the wavepacket originating from S1 TS with zero momentum, the trajectory shown (red line) has the largest weight. The final GGP (summed over the two states) is equal to 44% of the global wavepacket (S1:S0 ) 11%:33%). This is the only trajectory that does not describe regeneration of the reactant. From S1 TS with extra momentum along the TV toward product, three basis functions regenerate the reactant. The fourth trajectory (blue line) corresponds to 23% of the global wavepacket (S1:S0 ) 1%:22%). Movies for each of the trajectories represented in Figure 1 are available in the Supporting Information in QuickTime format. As expected (by analogy with MEPs starting from S1 TS, which have the same driving momentum: the TV), when starting from S1 TS with no extra momentum, part of the wavepacket dissociates to molecular products in the ground state;following decay at the S1/S0 seam close to the minimum-energy point, C1 CoIn;and part goes back to the S1 basin region. Starting from S1 TS, but now adding extra momentum along the TV in either direction, showed that giving a large momentum in the direction represented in Figure 2 (panel a) leads to H+HCO formation (whereas adding a large momentum in the opposite direction leads back to the minimum region in S1, as expected). Both results differ from the no-momentum case described above, in which, although the wavepacket naturally acquired momentum along the TV, that momentum was redistributed to the left-hand side of Figure 1 when the wavepacket crosses the seam (i.e., to the molecular products region, where the difference between the two CH bond lengths does not exceed 1 Å).
Arau´jo et al. For comparison, we also plotted the direction of the TV (Figure 2, panel a) in Figure 1 (solid black line). This line corresponds to a rectilinear displacement of the nuclei (a normal mode), although it appears curved in Figure 1 because it is represented in the space formed by two internal coordinates that are curvilinear. This pathway would be followed by the system if it started from the S1 barrier with a very large momentum parallel to TV (such that the potential energy could be neglected in front of the kinetic energy). Competition between two different momenta thus seems to control the H2+CO:H+HCO branching ratio: the momentum given initially and the momentum acquired when the wavepacket is deflected when crossing the seam. This can be confirmed by analyzing the evolution of the nonadiabatic coupling terms after a large adiabatic population is transferred at the S1/S0 seam of intersection. For the Gaussian function that represents the formation of H2+CO, the absolute value of H12 (potential-like coupling between both diabatic states) increases significantly (see Figure 1). This means that the corresponding quantum trajectory has been deflected along the derivative coupling vector (direction of the gradient of H12), represented in Figure 2 (panel b), when crossing from S1 to S0. In contrast, when dissociation to H+HCO takes place, H12 is roughly preserved after the crossing (see Figure 1). This effect is well-known in nonadiabatic quantum dynamics.13 In fact, it is used in the semiclassical trajectory surface hopping method,14 for which the energy loss is compensated by extra momentum given along the direction of the derivative coupling when the trajectory hops down at finite energy gap. In quantum terms, a wavepacket starting on S1 and arriving at a conical intersection will split into two components, one on each electronic state. The ground-state component can thus spread along the sides of the lower cone due to the shape of the potential energy surface. When the wavepacket is fast, the original momentum dominates, and it does not spread as much as a slow wavepacket. Therefore, for this system, giving momentum along the TV in the products’ direction seems to act against the nonadiabatic deflection that needs to occur at the conical intersection to produce H2+CO. This change of direction should be prevented if the wavepacket still has a large momentum along the TV (leading to H+HCO) when arriving in the seam region. If so, then a possible intermediate regime in favor of H2+CO might be obtained if a small (negative or positive) momentum were given along TV. We confirmed this prediction with a more systematic study that covered a wider range of values for the extra momentum. Note that the case with 0.1 eV of excess energy (blue line) could be surprising in that directed and large momentum leads to minor product formation in a range of 23% while the major product corresponds to the regeneration of the reactant, which requires an inversion of the momentum. This can be due to two factors: first, the excess energy is relatively small (about 2 kcal mol-1), which means that a significant part of the wavepacket can still spread back toward the well around S1 min, while the rest will move forward, following the extra momentum; second, part of the S1-forward-moving wavepacket is deflected back toward the reactant when crossing the seam, because the wavepacket will spread on the S0 surface all over the lower sheet of the double cone (aforementioned momentum redistribution). Using 4 GBF resulted in the reactant being regenerated plus either one or the other product being formed in each simulation (more in agreement with what could be expected from semiclas-
Controlling Product Selection in H2CO Photodissociation
J. Phys. Chem. A, Vol. 114, No. 45, 2010 12019
Figure 3. 8-GBF-case: quantum trajectories plotted in the frame of the HCH angle and the difference between the two CH bond lengths (red lines: dissociation to CO+H2; blue lines: dissociation to HCO+H; green lines: regeneration of reactant). Black crosses: diabatic S1/S0 intersections. Each panel shows, for a given initial geometry/initial momentum pair, all 8 trajectories followed by individual components of the global wavepacket. The magnitude of the extra momentum is given in terms of the corresponding extra kinetic energy.
Figure 4. Final distribution of population among geometries and electronic states for each global wavepacket (basis set: 8 GBF). Green arrows: initial momentum given along the TV toward reactant; blue arrows: initial momentum (of increasing magnitude from left to right) given along the TV toward product. The magnitude of the extra momentum is given in terms of the corresponding extra kinetic energy.
sical dynamics), whereas using 8 GBF allowed for the wavepacket to split and give a finite branching ratio between different products and reactant, as can be seen in Figures 3 and 4. Furthermore, by comparing the final populations of similar cases obtained using 4 GBF and those obtained using 8 GBF, we did not find a quantitative agreement; but results do agree qualitatively (in terms of which product dominates, as a function of initial conditions). Therefore, simulations of the first type can be regarded as computationally cheaper explorations aimed at easing the interpretation of the effect of giving an initial momentum to the wavepacket; using 8 GBF gives betterconverged results (by construction, due to the variational character of the method), but confirms the same trend. In Figure 3, the case with no extra momentum shows three trajectories describing the formation of H2+CO (in red), while two such trajectories are observed in the case with an extra momentum corresponding to 0.1 eV. This information is not sufficient for a comparison of the two cases because there is no
indication of the relative weights of the trajectories in the global wavepacket or of how the population is distributed between S1 and S0. To this end, Figure 4 shows the final distribution of population among geometries and electronic states for selected trajectories starting from S1 TS and giving extra momentum, as shown in the bottom of the figure. For each trajectory represented there, all eight components of the wavepacket have been considered in determining the total population for each final electronic state/final geometry pair. The central case is the S1 TS-starting, no-momentum case. To the right, there is a growing intensity of the momentum being added along the TV (as shown in Figure 2, panel a), making dissociation to H+HCO ever more likely and ultimately decreasing H2+CO formation. To the left, there is an addition of momentum along -TV; as a result, more reactant is regenerated. In all four cases, part of the global wavepacket follows a trajectory that traces a path toward products. As expected, when the momentum given initially is directed toward reactant, there is a delay in population transfer at the S1/S0 seam; the opposite occurs if momentum is directed toward products. From S1 TS to reactant, with 0.1 eV of excess energy, the only crossing occurs at 12.8 fs. From S1 TS with no initial kick, the first crossing occurs at 4.9 fs and with 0.1 and 0.2 eV of excess energy at 3.7 and 3.0 fs, respectively. As a final remark, some trajectories represented in Figure 3 are diverted back to the reactant region after hitting the S1/S0 seam, a typically quantum behavior that corresponds to the wavepacket splitting into different components when the nonadiabatic transition occurs. With negative momentum (directed from S1 TS to reactant), the difference between the two CH bond lengths stays around zero (see Figure 3), and regeneration of S0 H2CO can be understood in terms of FGR decay (with a low efficiency when looking at the final S1:S0 population ratio in Figure 4, as expected). With positive momentum (from S1 TS to products), a fraction of S0 H2CO is formed after the system has passed in the region of the seam. This means a nonadiabatic, direct S1/S0 decay around the S1/S0 conical intersection, with formation of predissociated H · · · CHO that recombines either into H2CO (see green trajectories in Figure 3, where the difference between the two CH bond lengths can reach about 1
12020
J. Phys. Chem. A, Vol. 114, No. 45, 2010
Å) or, with more efficiency, dissociates into products (see red and blue trajectories in Figure 3). Conclusions We showed with direct quantum dynamics simulations that the H2+CO:H+HCO branching ratio in the S1/S0 nonadiabatic photodissociation of formaldehyde is controlled by the direction and size of the mean momentum of the wavepacket when it crosses the seam of conical intersection. The dominant product results from a competition between the momentum developed on the S1 potential energy surface and the momentum acquired when the wavepacket is deflected along the derivative coupling vector when crossing from S1 to S0. To prove this, we started the simulations from the S1 transition state and varied the size of the initial mean momentum given along the direction of the TV toward products. Our results show that letting the wavepacket fall down from the barrier to the conical intersection with no initial momentum leads to H2+CO, while extra momentum toward products favors the formation of H+HCO through the same seam of conical intersection. This study is not intended to simulate any practical experiment. Rather, we have suggested a “dynamics mechanism” for controlling motion on a potential energy surface at a quantum level. Acknowledgment. Support from the European Social Fund and FCT-Portugal (Grant No. SFRH/BD/36197/2007) and the EPSRC-U.K. (Grant No. EP/F028296/1) is acknowledged. All calculations were carried out using the Imperial College High Performance Computing Service: http://www3.imperial.ac.uk/ ict/services/teachingandresearchservices/highperformancecomputing. Supporting Information Available: Movies for each of the trajectories represented in Figure 1 in QuickTime format. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Yeung, E. S.; Moore, C. B. Photochemistry of single vibronic levels of formaldehyde. J. Chem. Phys. 1973, 58, 3988–3998. (2) Townsend, D.; Lahankar, S. A.; Lee, S. K.; Chambreau, S. D.; Suits, A. G.; Zhang, X.; Rheinecker, J.; Harding, L. B.; Bowman, J. M. The roaming atom: Straying from the reaction path in formaldehyde decomposition. Science 2004, 306, 1158–1161.
Arau´jo et al. (3) Yin, H. M.; Kable, S. H.; Zhang, X.; Bowman, J. M. Signatures of H2CO photodissociation from two electronic states. Science 2006, 311, 1443–1446. (4) Araujo, M.; Lasorne, B.; Bearpark, M. J.; Robb, M. A. The photochemistry of formaldehyde: Internal conversion and molecular dissociation in a single step. J. Phys. Chem. A 2008, 112, 7489–7491. (5) Maeda, S.; Ohno, K.; Morokuma, K. Automated global mapping of minimal energy points on seams of crossing by the anharmonic downward distortion following method: A case study of H2CO. J. Phys. Chem. A 2009, 113, 1704–1710. (6) Zhang, P.; Maeda, S.; Morokuma, K.; Braams, B. J. Photochemical reactions of the low-lying excited states of formaldehyde: T1/S0 intersystem crossings, characteristics of the S1 and T1 potential energy surfaces, and a global T1 potential energy surface. J. Chem. Phys. 2009, 130, 114304–10. (7) Araujo, M.; Lasorne, B.; Magalhaes, A. L.; Worth, G. A.; Bearpark, M. J.; Robb, M. A. The molecular dissociation of formaldehyde at medium photoexcitation energies: A quantum chemistry and direct quantum dynamics study. J. Chem. Phys. 2009, 131, 144301–8. (8) Sicilia, F.; Blancafort, L.; Bearpark, M. J.; Robb, M. A. New algorithms for optimizing and linking conical intersection points. J. Chem. Theory Comput. 2008, 4, 257–266. (9) Worth, G. A.; Robb, M. A.; Lasorne, B. Solving the time-dependent Schro¨dinger equation for nuclear motion in one step: direct dynamics of non-adiabatic systems. Mol. Phys. 2008, 106, 2077–2091. (10) Worth, G. A.; Beck, M. H.; Ja¨ckle, A.; Burghardt, I.; Lasorne, B.; Meyer, H. D. The MCTDH Package, Development Version 9.0; University of Heidelberg, Heidelberg, Germany, 2008 (see http://www.pci.uniheidelberg.de/tc/usr/mctdh/). (11) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Scalmani, G.; Mennucci, B.; Barone, V.; Petersson, G. A.; Caricato, M.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Li, X.; Hratchian, H. P.; Peralta, J. E.; Izmaylov, A. F.; Kudin, K. N.; Heyd, J. J.; Brothers, E.; Staroverov, V. N.; Zheng, G.; Kobayashi, R.; Normand, J.; Sonnenberg, J. L.; Ogliaro, F.; Bearpark, M.; Parandekar, P. V.; Ferguson, G. A.; Mayhall, N. J.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Burant, J. C.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Chen, W.; Wong, M. W.; Pople, J. A. Gaussian DeVelopment Version, Revision G.01; Gaussian, Inc., Wallingford CT, 2007. (12) Allan, C. S. M.; Lasorne, B.; Worth, G. A.; Robb, M. A. A straightforward method of analysis for direct quantum dynamics: Application to the photochemistry of a model cyanine. J. Phys. Chem. A 2010, 114, 8713–8729. (13) Desouter-Lecomte, M.; Dehareng, D.; Leyhnihant, B.; Praet, M. T.; Lorquet, A. J.; Lorquet, J. C. Nonadiabatic unimolecular reactions of polyatomic molecules. J. Phys. Chem. 1985, 89, 214–222. (14) Tully, J. C.; Preston, R. K. Trajectory surface hopping approach to nonadiabatic molecular collisions: The reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562–572.
JP109549R
