401
Chemical Physics 111 (1987) 401-408 North-Holland, Amsterdam
PROTON-HYDROGEN COLLISIONS: TRAJECTORY SURFACE HOPPING CALCULATIONS Ch. SCHLIER,
U. NOWOTNY
Fakultiir ftir Physik,
Universitiit Freihurg, D 7800 Freiburg, FRG
Received
25 September
AND NEW MEASUREMENTS
and E. TELOY
1986
The trajectory surface hopping (TSH) method is used to compute product cross sections for the reactions H + + D,, D+ + D, and D+ + H,. They are compared with (mostly) integral cross sections obtained with a guided beam apparatus. The agreement is excellent in view of the absolute character of this comparison. Some systematic deviations of the computation can be traced to deficiencies of the potential (DIM), and of the TSH method.
1.Introduction Fifteen years ago Tully and Preston [l] * invented the trajectory surface hopping (TSH) method for the simulation of collisions (reactive or non-reactive) involving two or more potential energy surfaces. This method extends the quasiclassical trajectory (QCT) method to cases, where the “hop” from one potential energy surface to the other can be localized at “seams”, where a pair of surfaces is strongly coupled. In many cases these seams are avoided crossings, which exist abundantly in reactive molecular open shell systems. It is important to note that in molecular (i.e. at least u-i-atomic) systems the surface hopping picture will be a valid one, whenever energetic and other constraints allow the trajectories to reach the seam with more than marginal velocity. Since the coupling (which in many cases can be expressed by the gap width between the pair of adiabatic surfaces) will in general vary strongly in configuration space, the bundle of incoming trajectories will at a fairly well defined place in configuration space be divided between mostly adiabatic (P,, = 0.0) and mostly diabatic ( Phop =
* Actually the TSI method Dfiren [2].
was indepently
also invented
by
0301-0104/87/!§03.50 0 Elsevier Science Publishers (North-Holland
Physics
Publishing
Division)
1.0) behaviour. In many cases (including the present one), a change in the coupling causes only a small displacement of this location, and this has no large consequence for the trajectory bundle. It follows that (in contrast to common belief) it is not important to know the surface coupling very accurately in order to produce good TSH calculations. This insensitivity has already been noted in ref. [l]. In spite of its potential benefits, not too many TSH calculations have been performed until now. We mention the systems H+ + H, and its isotopes [3-51, At-+ + H, [6], Ne + Hz [7], alkali-halogen systems [S], rare gas-I; [9], CH + 0 [lo], C + 0, [ll], Na* + H, [12], K + 0, [13], and Cl- + H, [14]. One reason for this small number is, of course, the lack of good potentials. Without reasonable potential energy surfaces the results of TSH calculations are sometimes of as little value as those of some QCT calculations. Another reason may be the (often unnecessary, see above) fear, that in addition to the surfaces, the couplings might be needed from the quantum chemist. Finally, there is up to now little knowledge about the errors involved in the TSH method. A valid assessment of the TSH method must use a computation founded on accurate surfaces, and must compare it to precise experiments, since a comparison to a 3D quantum mechanical calculation B.V.
402
Ch. Schlier et ul. / Trcjectoty
sur~uce hoppwg m proton - hwlrogen collisions
is still unfeasible. The system, for which this assessment can be most easily done is H: and its isotopes. The present paper is a first step to such an assessment. It is still hampered by the use of the DIM potential (31 of Hc. Better potentials for the whole configuration space will be available soon [15], and the computation can then be repeated. Nevertheless, it was thought worthwhile to publish our results now, since they include experimental data, which hitherto exist only in unpublished form [16,17]. Another piece of comparison concerned with vibrationally resolved differential cross sections at 20 eV (c.m.) will be published separately [18].
2. Experimental method The experiments were performed with an improved version of the guided beam apparatus described in refs. [19-211. In short, the guided beam method uses inhomogeneous rf fields (here in the shape of an octopole) to guide the primary ion beam through the scattering chamber, and to collect all secondary products and guide them to the detector. It is especially suited for the determination of integral cross sections, since a true 4~ collection of products can be obtained. (The situation is thus much better than in neutral beam scattering.) Intensities are high enough to make the statistical errors of integral cross sections negligible. The largest single systematic error which remains is that of the pressure in the scattering chamber. Having compared three baratron manometers, we believe that this error is no more than 5%. The next important systematic error in our energy range, the loss of primary beam by large angle scattering, resulting in transverse energies which cannot be contained by the guiding rf field, was carefully monitored. It should be negligible below = 5 eV c.m., and at most a few percent at the largest energies shown. So we believe that the absolute scale of our experiments can be marked with a _+lO% systematic error, while the relative errors of the curves are much smaller. Measurements for different isotopes were referenced
against each other, so again the relative error of the curves in relation to others should be only a few percent. The only differential cross section shown in this paper (fig. 5 below) was measured in another apparatus, which has also been described [20,21]. Here, no absolute scale for da/dS2 is available. The main error is statistical since not only is the data very detailed but also most of the H+ + D, + HD+ intensity goes to much larger angles than the 5O shown.
3. Computational
method
Our TSH program, which has been documented in a thesis [22], has been written in close similarity to the prescriptions of ref. [l]. Like its authors we use the DIM potential for Hz with Morse functions as input [3], the data of which are presented in table 1. The lowest surface has a deep equilateral well (Ri = 0.92 A, D = 4.92 eV below H+ + Hz). This surface is coupled to the next higher one only at a seam in the asymptotic channels. This seam was fixed at r = 1.32 A, where r is the smallest of the three internuclear distances. The gap width between the two lowest adiabatic surfaces decreases exponentially with R, the atom-diatom distance. A section through the seam region is shown in fig. 1. In practice, all seam crossings for R < 3 A were completely adiabatic (and not counted as such in the program). This means that the processes involving “strong interaction” (rearrangement and vibrational energy transfer) are almost completely decoupled from surface hopping. Surface hopping occurs either in the entrance channel (if the vibrational excitation of the target diatom is large), or in the exit chan-
Table 1 Morse parameters for the three potentials which serve as an input to the DIM procedure for H:. Note that H: (*X,) is an “anti-Morse” potential Potential
0, (W
R, (A)
B, (A-‘)
H,(‘Z,) H; (*Es) H; (‘2,)
4.7462 2.7942 17.102
0.7413 1.0600 0.1821
1.9736 1.3594 1.6456
Ch. Schlier et al. / Trajecton, surfuce hopping in prototl
1
0
I
2 I
R [Al
3
v kV -1
-2
-3
-4
Fig. 1 ciection of the two lowest adiabatic potentials showing the coupling “seam” near r( = R,,) =1.32 collinear arrangement.
of H: i\ in a
nel, if vibrational excitation has been high enough (u’ > 4 for H,, u’ > 6 for D2). The hopping probability, P,,, was computed by the same quasi-Landau-Zener formula as in ref. [l]. Hops which were neither completely adiabatic nor completely diabatic, and thus contributed to the branching between both surfaces, were largely limited to R = 4.5 * 1.0 A. This explains the enormous robustness of the cross sections against changes in the coupling parameter. Unpublished computations by Schupp [23] showed that a change of the Landau-Zener coupling by a factor 2.5 in both directions hardly affects the integral cross sections reported below. (The changes were less than 25% for H+ + D, at 2.5 and 4.5 eV, and less than 10% at 7 eV.) This can be understood by computing the corresponding shift in the main hopping region, which turns out to be from 4.5 to 3.8 (5.2) A when the coupling is
-
hydrogen collisions
403
increased (decreased). Since the density of the trajectories at the seam is quite uniform, such a shift means only that other trajectories hop than before, not another number of them. For each energy and each mass combination 2000 trajectories were computed. Vibrational and rotational initial energies were fixed by setting ui = 0 and ji = 1 (the most probable j which is populated at room temperature). All other initial conditions were set by the usual Monte Carlo selection including the impact parameter with weight P(b) = b db. The maximum impact parameter, sufficient to catch all trajectories which lead to reaction and/or charge transfer, was between 1.7 and 0.9 A depending on the energy, and somewhat on the masses. Trajectory branches were followed up to a branching depth of 4, and for 0.05 < Phop < 0.95. This is called the “ants” procedure in ref. [l]. Otherwise only one branch, randomly selected according to P,, was continued (the “ant-eater” method). Above 4.7 eV dissociation of the product molecule is also possible, and plays an increasing role if the energy becomes larger. Dissociated molecules were assigned equally to the possible ionic and atomic product channels. The dissociation energy (4.74 eV) of the molecular products seems also to be the absolute upper energy limit for the formation of long-lived complexes. (This is an exact result for collinear complexes.) We have discussed complex formation in the system HD: in a former paper [24]. In general, the energy range where complex formation plays the predominant role, and that one where surface hopping is prominent are well separated (below and above E = 2.5 ev). Nevertheless there is a non-negligible number of long-lived complexes well in the surface hopping range, among which many hop more than 16 times. As before, we define here a long-lived complex in a mass-independent way by counting the number of “minimum exchanges”, i.e. changes of the identity of R,, = min(R,,, R,,, R,,), e.g. from R,, to R AB’ A pair of minimum exchanges can be thought of as “one vibration” of the complex, and we have shown elsewhere [25] that demanding N > 8 minimum exchanges is equivalent to other definitions of long-lived and chaotic trajectories. Com-
Ch. Schlier et 01. / Trujecto~
404
surface hopptng itt proton - h.vdrogen collisions
Table 2 Complex formation cross section (defined with more than eight minimum exchanges, see text) for different mass combinations. In addition to the translational energy, the internal energy corresponding to u = 0, j =l (0.27 eV for H,, 0.19 for D2) was present Translational
Complex formation cross section (10m3 A*)
energy (eV)
H++D,
D++D,
D+ +H,
2.5 3.0 3.5 4.0 4.5
705 f 329 + 162 * 31+ 5rt
9091~64 520 + 43 217 + 27 58&13 5& 4
1254*73 627 + 49 237+28 90*17 2+ 2
53 30 21 9 3
Table 3 shows that our results are virtually identical to those of Krenos et al. [4] at the three energies where they are reported.
4. Results Integral
cross sections
for the reactions
04
H++D*+D++HD, +D+HD+,
(lb)
--j H + D; ;
(lc)
(2)
D++D,+D+D;; plex formation cross sections for E = 2.5-4.5 eV are shown in table 2. Some of the long-lived trajectories have not been finished in the calculation even allowing thirty minimum exchanges and 400 fs of flight time. Though these were never more than 3% of all trajectories, they must not altogether be neglected, but distributed into the available channels. We assumed that complexes which are stopped without any surface hop would remain on the lower surface, while those with one (and that meant in each case: several) surface hops have a 50 : 50 percent chance to reach each surface. Within one surface each channel was considered to be equally probable. This procedure was checked by looking at the branching of trajectories with fifteen to thirty minimum exchanges into the different channels, and was found to be valid. Since some of the channels to be corrected are low-probability ones, the correction may amount to much more than 3% for those channels. The largest correction was 24% (HD+ from H+ + D, at 2.5 eV).
(34
D++H,+H++HD,
w
--+H+HD+;
are presented in figs. 2-4. We discuss first fig. 2, showing the three experimentally observable channels from H + + D,. The measured data have already been published [26]. We note the following: (1) The overall agreement (on an absolute scale!) is excellent even if some systematic differences are visible. (2) The agreement is almost perfect for the rearrangement channel (D+) below 4.5 eV. Above this energy, more and more D+ comes from dissociation, mostly from highly excited DC. This contribution is 20% at 5 eV, 74% at 7 eV, 91% at 9 eV, and 96% at 12 eV. Even then, the agreement between experiment and calculation is very good in view of the errors in the D: channel. (3) The computed threshold behaviour of the charge exchange channels (HD+ and D: with and without rearrangement) is displaced by 0.4 eV compared to the experimental values (AH’ = 1.83
Table 3 Comparison of TSH Results of Krenos et al. [4] and this work. Units are 10-r .&’ Energy
Ref.
D+
HD+
D:
[41 this work
31 *4 31.85 1.8
20 +3 21.9 + 1.5
40 *4 36.15 1.9
_
]41 this work ]41 this work
9 +1 11.5+1.3 5 +1 6.5 +0.9
20 &2 18.3 * 1.7 14 *2 9.9+ 1.1
73 *3 61.9* 3.1 72 +4 63.1& 2.8
7 8.3 16 18.3
DCL
HL
(eV) 4.0 5.5 7.0
_ _ +1 5 1.1 ?2 + 1.5
2 *1 2.4kO.6 3 *1 3.6? 0.7
Ch. Schlier
et al. /
Trujectoty
surfucr
hopping
in proton-hydrogen
405
collisions
Fig. 2. Integral cross sections for the reaction Hf + D, + products. Experimental points as indicated in the figure, their error is discussed in the text. Computed points (0, a, 0) have been partially connected by a thin line to guide the eye. Three typical statistical errors of the combination are shown, also (arrows) the thresholds for surface hop and dissociation. Note that the comparison between experiment and theory is on an absolute scale
eV including zero-point energies), and by 0.3 eV compared to the potentials used (cf. table 1). This is a definite shortcoming of the TSH procedure as it is implemented. The situation can be understood from fig. 1. It shows a section through the potential energy surfaces at the seam. The program demands that a surface hopping is always “vertical”. It is therefore allowed in this conformation only if the transverse total energy (i.e. the total energy minus the kinetic energy of the motion parallel to the seam) exceeds ‘2.36 eV. This is
much more than the energy difference between both valleys (1.96 eV), or the energy needed to reach the seam (2.06 eV) on the lower adiabatic surface. It is easily conceivable, that a quantum calculation would allow some “leakage” sideways across the seam, which would lower the energy requirements. Actually, assuming a linear threshold law, the deconvolution of the experimental data of fig. 1 for the charge exchange process H+ + D, -+ H +
2
Fig. 3. Same as fig. 2 for the reaction D+ + D,. charge exchange channel is observable.
Only
the
3
4
5
6
7
6
E,,[eVI
Fig. 4. Same as fig. 2 for the reaction Dt + H,. The channel leading to H, + is not observable. Note the different ordinate scales of figs. 2 to 4.
406
Ch. Schlier et al. / Trajectory surf&
D: yields a threshold of 1.80 + 0.05 eV, while the accepted value is 1.868 eV [27]. This shows that the real charge transfer reaction has no additional activation barrier within our error limits. (4) Between 2.5 and 6 eV the HDf channel shows a doubly peaked structure with an intermediate minimum near 4 eV. As Ochs and Teloy [26] remarked, the D+ channel shows a small hump near that energy, and it turns out that the sum of D+ and HD+, i.e. the rearrangement channel irrespective of the charge distribution, is very smooth and structureless. A similar structure can be seen as an S-shaped shoulder near 3 eV in the DC channel, and similarly in all other isotopic variants (figs. 3 and 4). The TSH calculation shows all these structures, albeit displaced on the energy scale (no wonder in view of what has been said above), and with some distortion. So it is clear that it is not an interference phenomenon, which leads to this “wavy” behaviour. A closer look on the trajectories as they cross the seam shows that the intermediate minimum in the D’ channel at 4 eV is correlated with a reorientation of the vibrational phase of the product molecule with respect to the effective beginning of the seam at R = 3.5 A. (5) The only severe systematic deviation of the computation from the experiments for reaction (1) is the falloff of the “charge transfer” channel (lc) at energies above 6 eV. The experimental curve descends much faster than the computed. one. Since the amount of intensity in this channel is determined by the amount of inelasticity transferred to the D, diatom in the collision, we tend to believe, that the responsibility for the wrong behaviour of the calculation lies in the faulty repulsive core of the DIM potential. This is a point which must be further investigated. (6) For the isotopic variants D+ + D, (fig. 3) and D’ + H, (fig. 4) the number of measurable channels is limited by the identity of particles, and the mass equality of Hz and D+. Almost everything which has been said before can be repeated here. An additional systematic deviation exists in reaction (3); however, the only case where the overall disagreement is definitely outside the error limits is channel (3b). Experiment and computation show large kinetic
hopping in proton - h.vdrogen collisions
isotope effects. Three pairs of reactions can be compared: Ht + D, + D+ (la) with D+ + H, -+ H+ (3a), H+ + D, -+ HD+ (lb) with Dt + H, -+ HD+ (3b), and the sum of Ht + D, -+ HD+ and + D: ((lb) + (1~)) with D+ + D, + DC (2). The general tendency is that the cross sections increase from H’ + D, to D+ + D,, and from there to Di + H,. This can be understood as a pure mass effect. Mahan [28] has introduced the model of a collinear sudden collision of A + BC. In this case the fraction of translational energy transferred to BC vibration can be computed to be
A&i, f=.=
trans
4m,m,mcM
= sin2(2P),
(mA + mBj2(mg + mc >’ (4)
where M=m,+m,+mc, and /3 is the usual skewing angle. For our three cases the f values are 5/9, 3/4, S/9. So, even if we have neither collinear collisions nor hard sphere potentials, a tendency to transfer more of the initial energy to the products going from reaction (1) to (2) and (3) can be expected. Assuming that higher vibrational transfer on an incorrect repulsive core leads also to larger errors in this transfer can perhaps explain, that the largest errors of our computations occur for D+ + H,. (7) More specific cross sections are, of course, included in the TSH calculation. Differential cross sections in specific vibrational states for H+ + H, at 20 eV have been measured at Gottingen, and compared with TSH calculations in a separate paper [18]. Many doubly differential cross sections d20/d52 dE’ of the system Hf + D, and some of H+ + H, have been measured by Gerlich [29] (see also ref. [21]), but most of these are in the energy range below the surface hopping threshold. For the others the comparison of experiment and TSH calculation is possible. It is, however, more difficult and more ambiguous than for integral cross sections, since (a) experimental resolution must be considered, (b) the experimental cross sections are not determined on an absolute scale, and (c) the greater detail increases the statistical error of the computation (or the computational work to reduce it) enormously. So we show here only one comparison (fig. 5),
Ch. Schlrer et ul. / Trajectoq sur$we hopping in proton - h_vdrogencoi~isions I
I
1
0
I
I
I
t
t
2
3
4
_:*_. I
E,, CaVl
Fig. 5. One example of a differential cross section: da/dSE at 5” (tab) versus the center of mass energy. The weak forward component of H+ + Dz + HD’+ + D exists only in the first of the two peaks (at = 3.5 eV not at = 5 eV), which the integral cross section shows in fig. 2.
the differential cross section do/dO for HC + D, -+ HD+ at 8,, = 5* as a function of collision energy. The experiment shows that only in the first of the two peaks showing up in the integral cross section (cf. fig. 2) there is fo~ard-scattered HD+. The calculation agrees with this, albeit with large statistical errors. It gives, moreover, an absolute scale (note the smallness), and tells us that only a small fraction (10%) of HI)’ is scattered forward, while most of it goes to large angles. A further assessment of the TSH method by monitoring differential cross sections should involve new calculations, and experiments specifically aimed at this comparison.
5. Conclusions The foregoing discussion shows that TSH calculations are well apt to simulate chemical reactions involving two or more potential energy surfaces, If the surfaces are correct, and one avoids effects dominated by quantum mechanics (e.g., rainbows) an overall accuracy of l&20% (or perhaps better) can be expected. In the implementations used so far, surface hopping at threshold must be considered a quantum effect, and will not be reproduced correctly by the TSH method. Further research should be done to improve this situation. In addition, it will be worthwhile to repeat our calculations as soon as complete and
407
accurate potential energy surfaces are available, in order to find out, how much of the small discrepancies remaining in this paper are due to the potentials.
We thank our collaborators G. Ckhs, G. Stasek, and W. Beyer, who measured the integral cross sections of H+ + D,, D” + H, and Dt + D,, respectively. D, Gerlich contributed much to the success of the guided beam method, and we thank him for many discussions. Support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References [l] J.C. TulIy and R.K. Preston, J. Chem. Phys. 55 (1971) 562. [Z] R. DSren, J. Phys. B 6 (1973) 1801. [3] R.K. Preston and J.C. Tully, J. Chem. Phys. 54 (1971) 4297. (41 J.R. Krenos, R.K. Preston and R. Wolfgang, J. Chem. Phys. 60 (1974) 1634. [5] R.K. Preston and R.J. Cross, J. Chem. Phys. 59 (1973) 3616. [6J S. Chapman and R.K. Preston, J. Chem. Phys. 60 (1974) 650; S. Chapman, J. Chem. Phys. 82 (1985) 4033. [7] P.J. Ktmtz, J. Kendrick and W.N. Whitton, Chem. Phys. 38 (1979) 147. [8] C.W.A. Evers and A.E. de Vries, Chem. Phys. 15 (1976) 201; C. Evers, Chem. Phys. 21 (1977) 355; J.A. Aten, C.W.A. Evers, A.E. de Vties and J. Los, Chem. Phys. 23 (1977) 125; C.W.4. Evers, A.E. de Vries and J. Los, Chem. Phys. 29 (1978) 399; C.W.A. Evers, Chem. Phys. 30 f1978) 27. PI B. Garetz, M. Rubinson- and J.I. Steinfeld, Chem. Phys. Letters 28 (1974) 120; B.A. Garetz, L.L. Poufson and J.I. Steinfeid, Chem. Phys. 9 (1975) 385. II01 M. McGregor and R.S. Rerry, J. Phys. B 6 (1973) 181. 1111S.R. Kinnersly and J.N. Murretl, Mol. Phys. 33 (1977) 1479. I121 N.C. Blais and D.G. Truhlar, J. Chem. Phys. 79 (1983) 1334. [I31 G. Parlant, M. Schriider and S. Goursaud, Chem. Phys. 75 (1983) 175. U4] M. Sizun, E.A. Gislason and 6. Parland, Chem. Phys. 107 (1986) 311.
408
Ch. Schlier et al. / Trujector3, surface hopping in proton - h_vdrogen collisions
[15] W. Meyer, P. Botschwina and P. Burton, J. Chem. Phys. 84 (1986) 891; private communication. [16] G. Stasek, Diploma Thesis, Freiburg, FRG (1975). [17] W. Beyer, Diploma Thesis, Freiburg, FRG (1976). [18] G. Niedner, M. NOB. J.P. Toennies and Ch.G. Schlier, to be published. [19] E. Teloy and D. Gerlich, Chem. Phys. 4 (1974) 417. [20] E. Teloy, in: Invited Papers, 10th ICPEAC, ed. G. Watel (North-Holland, Amsterdam, 1978) p. 591; [21] D. Gerlich, in: Electronic and atomic collisions, eds. D.C. Lorents, W.E. Meyerhof and J.R. Peterson (Elsevier. Amsterdam, 1986).
1221 U. Nowotny, Thesis, Freiburg, FRG (1985). [23] W. Schupp, private communication. [24] D. Gerhch, U. Nowotny, Ch. Schlier and E. Teloy, Chem. Phys. 47 (1980) 245. [25] Ch. Schlier, in: Energy storage and distribution in molecules, ed. J. Hinze (Plenum Press, New York, 1983). [26] G. Ochs and E. Teloy, J. Chem. Phys. 61 (1974) 4930 L. [27] K.P. Huber and G. Herrberg. Molecular spectra and molecular structure, Vol. 4, Constants of diatomic molecules (Van Nostrand-Reinhold, New York, 1979). [28] B.H. Mahan, J. Chem. Educ. 41 (1974) 308. [29] D. Gerlich, Dissertation, Freiburg, FRG (1977).