State-to-state reaction probabilities from the time-independent wavepacket reactant-product decoupling equations: application to the three-dimensional H + H2 reaction (for J = 0)

State-to-state reaction probabilities from the time-independent wavepacket reactant-product decoupling equations: application to the three-dimensional H + H2 reaction (for J = 0)

29 August 1997 CHEMICAL PHYSICS LETTERS ELSEVIER Chemical Physics Letters 275 (1997) 173-180 State-to-state reaction probabilities from the time-in...

494KB Sizes 0 Downloads 22 Views

29 August 1997

CHEMICAL PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 275 (1997) 173-180

State-to-state reaction probabilities from the time-independent wavepacket reactant-product decoupling equations: application to the three-dimensional H + H 2 reaction (for J = O) Stuart C. Althorpe

a,l,

Donald J. Kouri

a,2, David

K. H o f f m a n b

a Department of Chemistry andDepartment of Physics, University of Houston, Houston, TX 77204-5641, USA b Department of Chemtstry and Ames Laborato~ 3, Iowa State University, Ames, IA 50011, USA Received 16 April 1997; in final form 27 May 1997

Abstract The reactant-product decoupling (RPD) equations are a rigorous formulation of quantum reactive scattering recently introduced by Peng and Zhang. Solving the RPD equations yields state-to-state reaction probabilities from a set of decoupled wavepacket propagations, each of which is confined to an arrangement of the reaction. Previously, we derived the time-independent wavepacket version of the RPD equations, and developed an efficient Chebyshev-based propagator for solving them. In this Letter we apply our method to 3D reactive scattering, presenting state-to-state reaction probabilities for the H + H 2 reaction (for J = 0). © 1997 Elsevier Science B.V.

I. Introduction

wavefunction for reactive scattering, X(t), is partitioned according to

Peng and Zhang [ 1-4] 4 have recently introduced a new set of equations for calculating state-to-state reaction probabilities, whereby the time-dependent

X ( t ) = X r ( t ) + ~-~Xp(t). p

(1)

The arrangement functions, Xr(t) and Xp(t), satisfy the reactant-product decoupling (RPD) equations,

~Xr(t) i h - - - - - ~ = H x r ( t) - i ~,, Vpr Xr( t) ,

(2)

P

i Supported under National Science Foundation Grant CHE9403416. 2 Supported in part under R.A. Welch FoundationGrant E-0608. Partial support from the Petroleum Research Fund, administered by the AmericanChemical Society, is also acknowledged. 3 The Ames Laboratory is operated for the Department of Energy by Iowa State University under Contract No. 2-7405ENG82. 4 Several generalizations of the basic approach are given in Ref. [2].

ih aXp(t) Ot =Hxt,(t)+iVprXr(t),foreach

p ,

(3)

where H is the Hamiltonian for the reactive system, and where, for each value of p, - i V p r is an artificial partitioning potential, which is placed at the start of the pth product channel, just after the 'point-ofno-return' (i.e. after the point beyond which the system, to an acceptable level of accuracy, cannot

0009-2614/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PH S0009-2614(97)00744-6

174

S.C. Althorpe et al. / Chemical Physics Letters 275 (1997) 173-180

retum to the reactant channel). Eqs. (1)-(3) are exact, since the sum of Eqs. (2) and (3) is the time-dependent Schr~Sdinger equation for X(t). The first RPD equation (Eq. (2)) is the equation which is now applied routinely in time-dependent reactive scattering whenever one wishes to calculate either (definite energy) total reaction probabilities, or vibrationally resolved reaction probabilities [5-7]. The 'product partitioning potentials', - iVrr, absorb Xr(t) where it is travelling away from the interaction region into the product channels, with the result that Xr(t) is exactly equal to X(t) in the 'reactant arrangement' of the reaction (which we here define to include both the reactant (entrance) channel and the strong interaction region). The original authors pointed out that Eq. (2) is applicable to all chemical reactions for which it is possible to distinguish experimentally between reactants and products. The new feature of the RPD equations is the second set of equations (i.e. Eq. (3) for each value of p), in which the absorbed part of Xr(t) appears in the t-dependent source term, iVpr Xr(t). The latter is located after the point-of-no-retum at the start of the pth product channel, and is travelling away from the reactant arrangement, down the product channel. Each solution of Eq. (3), Xp(t), is thus confined to the pth product arrangement of the reaction, where it is exactly equal to X(t). Wavepacket methods for solving the RPD equations are likely to prove the most efficient methods which have yet been devised for calculating (exact) state-to-state reaction probabilities. Such methods will yield the exact wavefunction by means of a set of decoupled propagations, each of which is confined to one arrangement of the reaction, and which therefore scales better than n 2 (where n is the number of basis functions used in the one arrangement). Moreover, solving the RPD equations ameliorates the notorious 'coordinate problem' of reactive scattering, since one can solve each RPD equation in its own set of arrangement coordinates (provided one fits the source terms, iVprxr(t), of Eq. (3) from r-arrangement to p-arrangement coordinates). To our knowledge there are currently two efficient wavepacket methods for solving the RPD equations. The first of these was introduced by Peng and Zhang [1-3], and is based on the modified Cayley propagator [8]. The second method, which was de-

veloped by the authors of the present Letter, uses an extended form [9] of the Chebyshev propagator [10] to solve the time-independent wavepacket (TIW) version of the RPD equations [9,11] (or TIW-RPD equations). The latter equations are derived from the TIW equation for reactive scattering in a manner analogous to Eqs. (1)-(3). Information and references on the TIW equation - which is applicable to all branches of scattering - may be found in Refs. [12-15]. The object of this Letter is to report the first application of our Chebyshev-based method to full 3D reactive scattering. We summarise the TIW-RPD equations and the extended Chebyshev propagator in Section 2, and discuss their application to the 3D A + BC reaction (for J = 0) in Section 3. In Section 4, we test the method by applying it to the 3D H + H 2 reaction (for J = 0). Section 5 concludes the Letter.

2. The time-independent wavepacket reactantproduct decoupling (TIW-RPD) equations We shall assume that the reader is familiar with the time-dependent wavepacket (TDW) approach to scattering, in which an initial wavepacket located at a finite distance from the scattering potential is propagated by the time-evolution operator until it has exited from the scattering region. The time-independent wavepacket (TIW) equation [12-15] is derived by applying a partial Fourier transform to the TDW equation. We shall here be considering the integral form of the TIW equation [13],

=

i

G+-(E)x(O),

(4)

in which G+-(E) is the causal/anti-causal Green's function. It can be shown [13] that, in the region between the initial wavepacket and the scattering potential, the solutions, X-+(E), are proportional to the causal/anti-causal Lippmann-Schwinger solutions (of the time-independent SchriSdinger equation). In Refs. [9,11] we showed that, when the TIW equation, Eq. (4), is applied to reactive scattering, the solution, X +(E), can be partitioned according to, +-(E) = ~r+-(E) q- E ~ p + - ( E ) , p

(5)

S.C. Althorpe et a l . / Chemical Physics Letters 275 (1997) 1 7 3 - 1 8 0

where the functions, ~ r ± ( E ) and +p-+(E), satisfy the time-independent wavepacket reactant product decoupling (TIW-RPD) equations,

and T~p0 = ~Tpr 0 , r/pl = 2 H n o r m r / p 0

i ~r -+( E ) = - ~ a ;

(E)/1((0),

(6)

175

r/pn = 2 nnorm r/pn

+ ~prl

,

l -- r / p n - 2 ~- O'prn -- O'prn- 2 .

(12) +r,+-( E ) = - G + - ( E ) F e ) ( E ) + f ( E ) , foreach p .

(7) The additional Green's function, G ~ ( E ) , is associated with the reactive system plus the set of energydependent partitioning potentials, F p ) ( E ) , each of which is located at the start of the pth product arrangement, just after the point-of-no-return. In Ref. [9], we showed that the TIW-RPD equations may conveniently be solved by using an extended version of the familiar Chebyshev propagator. To do this, Fp)(E), must be assigned the form, Fpr+ ( E ) = T-iV,~AHsin ¢5,

(8)

where ~b is determined by cos qb = ( E - H ) / A H , and the scaling parameters, A H and H, are chosen so that the spectrum of the scaled Hamiltonian, H .... = ( H - H ) / A H , is confined to the interval [ - 1 , 1 ] . The solutions, ~r+-(E) and ~p+-(E) can then be expanded as, !

U

E e - i"'~r/r. 2 v A H s i n ~b n=0

~f(E)

1

(9)

N

~ . ± ( E ) = 2 ~ A H s i n ~ b E e-i"'C'r/p., n=O

(10)

where the energy-independent basis functions, r/~n and pp., are generated from the recursion relations, 1 r/r0

- - x ( O )

1 + l/tot

1 + Vtot

2//norm r/r 0

M ~r+-(E) 2 r r k H s i n & Ej e-i"J6r/r,j,

(13)

where the subscript, n j, signifies that the sum is over every Mth value of n only, and the parameter, M, is decreased until suitably converged results are obtained. In this approximation, the r/rn functions must still be generated by Eq. (11); the r/pn functions, however, are now generated by the modified recursion relation, rlpn

for t ' / : nj = 2 HnormVlpn- l - Dpn- 2 + O'prn forn - 2 = nj = 2Hnormrlpn- I - rlpn-2 - O'ern--2 otherwise. =2Hnormrlpn_l-rlpn- 2

(14) This recursion contains only 1 / M t h the number of source terms, O'prn, a s does Eq. (12), and thereby reduces the labour required to fit the ~rpr. terms from r-arrangement to p-arrangement coordinates (by a factor of M).

'

1 r/r I =

In Eq. (11) the potential, Vtot, is given by, Vtot = ~.Vpr, and serves to damp the r/r. basis functions to zero in the region of the partitioning potentials, Fp+(E) In Eq. (12), the source terms, O'pr., are given by, o'pr. = Vprr/r., and are localised in the region of Vpr (i.e. just after the point-of-no-return at the start of the pth product channel). The basis functions, r/rn and r/p., are therefore confined to the r-arrangement and to the p-arrangement respectively. Instead of computing the above recursions (which are exact in the limit that N ~ oo), we exploit the near linear dependency of the r/r . functions, and make the approximation [9,11 ],

'

1

r/r2

- -

r/r,,

1 - - [ 2 H 1 + Vlot

1 + V~o,

[2Hnormr/rl

-- 2r/r0],

3. A p p l i c a t i o n o f the T I W - R P D 3 D A + B C r e a c t i o n (for J = 0) ....

r/rn--'--( 1 --

e q u a t i o n s to the

gtot) r/rn- 2] ,

(11)

For the 3D reaction, A + BC, there are a total of 3 TIW-RPD equations, each corresponding to either

S.C. Althorpe et al. / Chemical Physics Letters 275 (1997) 173-180

176

the r-arrangement (A + BC), or to one of the parrangements (AB + C or AC + B). To solve the r-arrangement equation, Eq. (6), we use the rarrangement Jacobi coordinates, ( R r, r r, Or), of which R r is the BC bond length, r~ is the distance between A and the BC centre of mass, and Or is the angle between R r and r r. In these coordinates, the J = 0 Hamiltonian has the familiar form [16], 1

h 2

Hr

~2

h2

a 2R~ 2tz~ R r R, h2

2 m r r~ ~r 2

rr

h2

+ ~ +

1 a2

2mrr-~

j2 + V( R r, r~, Or) ,

(15) where /-£r, mr are the appropriate reduced masses, and j2 is the usual angular momentum operator in Or. We represent H r on a (3D) grid of discrete points, (R~r, r i, 07), such that each rl~,(R r, r r, 0n) function is expanded as, 1

Tlrn( g r , rr, Or) = ~ r r r

F,f~.( R~, rr)lOr~},

(16)

k

where the kets, 107), denote a set of Gauss-Legendre DVR functions [17] in 0r. The radial functions, f~,(R,, r,), are further expanded in terms of the (Ri,, r~) grid representation (which will be specified in the next section). The functions, "r/,,(R,, r,, Or), are of course confined to the r-arrangement by the partitioning potentials, Fp~(05), one of which is placed at the start of either p-arrangement. To solve the p-arrangement TIW-RPD equations (namely Eq. (7) for p = 1, 2) we use the p-arrangement Jacobi coordinates, (Rp, re, 0p), in terms of which the Hamiltonian, Hp, resembles Eq. (15) (with a p in place of each r). We represent H r on a 3D grid of discrete points such that e a c h Ylpn(R p, rp, Op) function is expanded analogously to ~),n(Rr, r r, Or) in Eq. (16). In general, a given grid point in the p-arrangement coordinates is not coincident with a grid point in the r-arrangement coordinates, and hence the source terms, O-p~,, of Eq. (14) must be fitted from the r-arrangement coordinates (in terms of which

they are initially obtained) to the p-arrangement coordinates. We use the simple fitting scheme,

O'prn(Rip, r j,

Ok ) = g A( Ripr, Rlr) A( r Jr, rrm ) lmn

Xa(O2r,O~")~,..(R'r, rrm,Or") (17) in which the notation, (Ripr, r pJr , ~,kr), signifies that . the p-arrangement grid point, (R~, rpj, 0p~), is expressed in terms of the r-arrangement coordinates, ( R r, r r, Or). The symbol, A, denotes an approximation to the (1D) Dirac delta function in each rarrangement coordinate; we will give the expressions used for each A in the next section. Note that the general fitting scheme of Eq. (17) is applicable whenever both H r and lip are represented by 3D grids of discrete points; (in non-grid representations, the (1 D) A matrices will be coupled together into a 3D fitting matrix). The advantage of Eq. (17) over similar fitting schemes (that have been used in calculations which did not solve the RPD equations) is that one need sum over only those grid points for which o-p~n is non-zero, and which are therefore confined to the region of the partitioning potential, Fpr(05). Thus the application of Eq. (17) will scale as only N [ N F, with Ni and N1e denoting the number of grid points that overlap Fpr(05) in the r and p grids. As we mentioned above, further efficiency may be obtained by employing the approximation of Eq. (13), which requires that one fit a total of only N / M of the O'prn terms. Having calculated the arrangement functions, ~:r+(E) and ~p(E), we extract the state-to-state reaction probabilities (and inelastic probabilities) by means of the TIW S-matrix Kohn variational formula [131,

ih2 k~okv S;[~°r ( E ) =

2

+

/z(2"rr) Btor(k.o)BtA( kv

)*

X(X,,za(O)]G+(E)lx~,otor(O)}.

(18)

The action of G + ( E ) on I X~010r(0)} is here a symbolic representation of the calculation of ~(E), with A denoting either r (for the inelastic S-matrix) or p

177

S.C. Althorpe et al. / Chemical Physics Letters 275 (1997) 173-180 ^+

which g,ja(R]0) is non-zero, and where hF(k ,,R) is the Ricatti-Hankel function which satisfies,

(for the reactive S-matrix); I X~0t0,(0)) is the initial wavepacket (X(0) of Eq. (6)), and (X, la(0)l is a 'test function' onto which is projected ~:~-(E). We take the test functions and the initial wavepacket to have the form,

lira h~(k~R) = e -+i(k"R-,=/2) ,

with k~ being the appropriate momentum.

X~oz,,a(Ra, re, 0el0) = g~0z0a(R,10)4~o,o(re)/3,0(0a), (19)

4. Numerical application of the TIW-RPD equations to the H + H 2 reaction (for J = 0)

where ~b.do(%) is the vibrational wavefunction for the diatom, Pt(0a) is a (normalised) Legendre polynomial, and g.otoa(R, lO) is the Gaussian wavepacket,

We have tested the theory described above by calculating the state-to-state reaction probabilities for the 3D H + H 2 reaction (for J = 0) on the LSTH potential surface [18]. The results are illustrated in Figs. 1 and 2, where they are compared with the results of a previous calculation by Zhang et al. [ 16]. In calculating these results, we took full advantage of the symmetry of the reaction, solving only one of the (two) p-arrangement TIW-RPD equations, and including only half of the Gauss-Legendre DVR points in 0r We represented the R r, r r and R e coordinates in the discrete distributed approximating

g.,,,,,a(R.dO) 1 ((Ra-R°)2) - vRl~/~7~- exp -2R~

" e'k"'R'"

(20)

in Eq. (18) are given by The coefficients, Bt~(k,.), + [13], 1 oR2^+

B,~(k~) = -~---£wjR ' hF(k,,R)g,4.~( RlO) dR,

(21)

where R~ ~
i

i

t

i

i

r

/'/"

0.09 0.08

/

~

,

0.07

,

(22)

g ~

i

"'",

/

',

]

" ~

"

-

0.06

=~ ~

. ............... 7~'.:.< ............ [] .........

/

0.05

13..

0.04

,,' /

. -..

,; / 0.03

~'"

/ /

0.02

," /

0.01

,'..' ,I.

."

'~ r~. . . . . . . . . . . m ....

.,

"-"

....

',x

"

'... '

,

"• 'x . . . . . . . . .

'

" ..........

--.X--

0 0.6

-"

0.7

0.8

0.9

1

1.1

i 1.2

1.3

1,4

Total Energy (eV)

Fig. l. Calculated state-to-state reaction probabilities, P ( v 0 = 0, Jo = 0 -..* v = O, j ) for the 3D H + H 2 reaction (for J = 0). The curves are the TIW-RPD results of this work for j = 0 (solid line), j = 1 (dashed line), j = 2 (dotted line), and j = 3 (chained line). The points are the results of Zhang et al. taken from Ref. [16].

178

S.C. Althorpe et al, / Chemical Physics Letters 275 (1997) 1 7 3 - 1 8 0 0.025

. . . .

/'] ..... .... ,,,

0.02

,S' \

0.015

',

,

n

ir I

0.01

/i

', o

///

',,

J

/ //

+ ,.

/,

0.005

0.8

0.9

1

1.1

1.2

1.3

1.4

Total Energy (eV) Fig. 2. Same as Fig. I for P ( v o = O, Jo = 0 --~ ~, = 1, j ) .

functional (DAF) representation [19], in which the grid points are distributed at equally spaced intervals, and the action of the kinetic energy operator is computed very efficiently by means of a fast convolution algorithm [20]. In the rp coordinate, however, we used the DVR [17] generated from the diatom vibrational wavefunctions, c~oo(re), as the latter are a very good basis set in this coordinate. The dimensions of all the grids are given in Table 1 together with other numerical parameters used in the calculation. An advantage of using the discrete DAF representation in (R r, rr), is that it enabled us to represent each of the fitting matrices, A ( R i p r , Rtr) and A(r~,, rr") , of Eq. (17) by a continuous DAF [19]. The latter acts as an accurate fitting matrix because it provides a "well-tempered approximation" [19] to the Dirac delta function, which is as accurate off the grid as on the grid. We represented the A(0~,, Or") matrix by expanding the Dirac delta function in terms of Legendre polynomials. In addition to the two partitioning potentials (one of which was placed at the start of each product

channel), we also placed absorbing potentials along the R r a n d R p coordinates at the grid boundaries (in order to prevent non-physical reflections). All the

Table 1 Values (in atomic units) taken by parameters in the reactant and product calculations Parameter

Reactant

Product

grid dimensions

Rmin Rmax rmi. rmax

1.0 11.0 0.2 6.5

3.5 13.5 0.2 3.5

grid points

NR N,

72 48 12

72 6 11

No

initial wavepacket

test functions

R0

6.5

-

Rt kny

0.5 - 7.0

-

Ro Rt kay AH = / ~

6.5 0.35 7.0 0.6

8.0 0.35 7.0 0.6

S.C. Althorpe et al. / Chemical Physics Letters 2 75 (1997) 173-180

partitioning/absorbing potentials took the simple form,

V~b~(R A, r~, Oh)

( R _RI )3 = 0.4 R2 = 0

RI

for Rx > R 1 otherwise, (23)

where A indicates the appropriate substitution of p or r. The interval, [R 1, R2], was set to [3.5,6.5] au in the partitioning potentials, [8.0, 11.0] au in the r-arrangement absorbing potential, and [9.5, 13.5] au in the p-arrangement absorbing potential. Finally, we note that adequate convergence (to better than 1%) was obtained with N (the highest value of n in Eqs. (9) and (10)) set to 2000, and M (see Eq. (13)) set to 16. The (small) differences between our results and the results of Zhang et al. [16] in Figs. 1 and 2 can be attributed mostly to the absorbing/partitioning potentials, and (probably in the case of the discrepancy at 1.1 eV) to small convergence errors in the results of Zhang. Solution of the TIW-RPD equations has proved to be an accurate and efficient method for calculating the state-to-state reaction probabilities of the 3D H + H 2 reaction.

5. Conclusions In this Letter we have shown how to solve the TIW-RPD equations for the 3D A + BC reaction (for J = 0) using our extended Chebyshev propagator, and have successfully tested the method by calculating state-to-state reaction probabilities for the (3D) H + H 2 reaction. We are now in a position to apply the method to more challenging reactions, where it is likely to prove much more efficient than most of the other wavepacket methods currently in use (which do not make use of the RPD equations). For example, we intend to calculate state-to-state reaction probabilities for the (fully 6D) H + H 2 0 and H 2 + OH reactions, systems for which state-to-state calculations have only very recently become possible [21-23] because of the large basis sets required. Clearly much smaller basis sets will be required in the equivalent TIW-RPD calculations, since each

179

RPD equation can be solved in its own optimised basis set, A further advantage of the TIW-RPD method is that, although it requires a coordinate transformation from reactant to (each set of) product Jacobi coordinates, the transformation need be applied only in the region of the partitioning potential. (We note that if hyperspherical coordinates had been used, the latter transformation would still be necessary, together with an additional transformation for each final state analysis. This is because the hyperangles appropriate to each arrangement are different, and because the physical hyperplanes on which the final states are defined are not orthogonal to the hyper-radius.) In future work, we shall be testing further developments of the RPD equations. For example, in a recent paper [2], it was shown that the RPD equations can be further partitioned into sets of vibrationally decoupled equations, and that the RPD equations also provide a very elegant formulation of molecular photodissociation. We also intend to compare the present TIW-RPD method with Peng and Zhang's time-dependent method [1,3] of solving the RPD equations. We note here that the TIW-RPD method exactly parallels Peng and Zhang's method and takes full advantage of the form of the RPD equations. For this reason, we expect the TIW-RPD method to be more efficient, since recent computational tests by Kroes and Neuhauser [24] have shown that TIW-based approaches are more efficient than the corresponding time-dependent wavepacket approaches. It is likely, however, that both methods will find wide application, given the importance of the RPD equations to the future of reactive scattering.

Acknowledgements We thank Dr. Wei Zhu and Professor J.H.Z. Zhang for useful discussions.

References [1] T. Peng, J.Z.H. Zhang, J. Chem. Phys. 105 (1996) 6072. [2] D.J. Kouri, D.K. Hoffman, T. Peng, J.H.Z. Zhang, Chem. Phys. Lett. 262 (1996) 519.

180

S.C. Althorpe et al. // Chemical Physics Letters 275 (1997) 173-180

[3] W. Zhu, T. Peng, J.Z.H. Zhang, L Chem. Phys. 106 (1997) 1742. [4] J. Dai, J.Z.H. Zhang, J. Chem. Soc. Faraday Trans. 93 (1997) 699. [5] D. Neuhauser, M. Baer, J. Chem. Phys. 90 (1989) 4351. [6] D. Neuhauser, M. Baer, R.S. Judson, D.J. Kouri, J. Chem. Phys. 90 (1989) 5882. [7] D. Neuhauser, M. Baer, R.S. Judson, D.J. Kouri, Comput. Phys. Commun. 63 (1990) 460. [8] R.S. Judson, D.B. McGarrah, O.A. Sharafeddin, D.J. Kouri, D.K. Hoffman, J. Chem. Phys. 94 (1991) 3577. [9] S.C. Althorpe, D.J. Kouri, D.K. Hoffman, J. Chem. Phys. 106 (1997) 000. [10] H. Tal-Ezer, R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [11] S.C. Althorpe, D.J. Kouri, D.K. Hoffman, J.Z.H. Zhang, J. Chem. Soc. Faraday Trans. 93 (1997) 703. [12] W. Zhu, Y. Huang, G.A. Parker, D.J. Kouri, D.K. Hoffman, J. Phys. Chem. 98 (1994) 12516. [13] D.J. Kouri, D.K. Hoffman, Few-Body Syst. 18 (1995) 203; D.J. Kouri, Y. Huang, W. Zhu, D.K. Hoffman, J. Chem. Phys. 100 (1994) 3662. [14] D.J. Kouri, Y. Huang, D.K. Hoffman, in: Dynamics of

[15] [16] [17] [18]

[19]

[20] [21] [22] [23] [24]

Molecules and Chemical Reactions, R.E. Wyatt, J.Z.H. Zhang (Eds.), Marcel Dekker, New York, 1996, pp. 307-321. S.C. Althorpe, D.J. Kouri, D.K. Hoffman, N. Moiseyev, Chem. Phys. 217 (1997) 289. J.Z.H. Zhang, D.J. Kouri, K. Haug, D.W. Schwenke, Y. Shima, D.G. Truhlar, J. Chem. Phys. 88 (1988) 2492. Z. Bacic, J.C. Light, Annu. Rev. Phys. Chem. 40 (1989) 469. B. Liu, J. Chem. Phys. 58 (1973) 1925; P. Siegbahn, B. Liu, J. Chem. Phys. 68 (1978) 2457; D.G. Trnhlar, C.J. Horowitz, J. Chem. Phys. 68 (1978) 2466; J. Chem. Phys. 71 (1979) 1514(E). D.K. Hoffman, T.L. Marchioro II, M. Arnold, Y. Huang, W. Zhu, D.J. Kouri, J. Math. Chem. 20 (1996) 117, and references therein. Y. Huang, D.J. Kouri, M. Arnold, T.L. Marchioro II, D.K. Hoffman, Comput. Phys. Commun. 80 (1994) 1. D.H. Zhang, J.C. Light, J. Chem. Phys. 105 (1996) 1291. W. Zhu, J. Dai, J.Z.H. Zhang, D.H. Zhang, J. Chem. Phys. 105 (1996) 4881. R.T Pack, G.A. Parker, J. Chem. Phys. 87 (1987) 3888. G.J. Kroes, D. Neuhauser, J. Chem. Phys. 105 (1996) 8690.