The application of time-dependent wavepacket methods to reactive scattering

460

Computer Physics Communications 63 (1991) 460—481 North-Holland

The application of time-dependent wavepacket methods to reactive scattering Daniel Neuhauser

1 Department of Chemistry, Princeton University, Princeton, NJ 08544, USA

Michael Baer Department of Theoretical Physics and Applied Mathematics, Soreq Nuclear Research Center, Yavne, 70600 Israel

Richard S. Judson Sandia National Laboratories, Livermore, CA 94550, USA

and Donald J. Kouri Department of Chemistry and Department of Physics, University of Houston, Houston, TX 77204-5641, USA Received 26 March 1990

In this article, we review several methods for performing numerically-exact reactive scattering calculations using time-dependent wavepackets. The basic idea we imply is to take the multi-arrangement reactive problem and reformulate it as one or more inelastic ones. In the simplest method, we extract total reaction probabilities by calculating the flux of the wavepacket as it leaves the interaction region in the direction of the reactive arrangement. To make this practical, we use complex potentials that absorb the wavepacket before it reaches the numerical grid boundary. We describe methods that generate observables ranging from total, energy-averaged reaction probabilities up to energy- and state-resolved S-matrix elements. We also review techniques for efficiently performing the necessary inelastic wavepacket propagation.

that all numerical effort is concentrated on a

1. Introduction In the last few years, time-dependent wavepacket methods [1—42]have gained considerable importance as tools for probing surface [13—16] and inelastic scattering [10,19—25]and photodissociation [26—30]processes. The time-dependent method has some powerful advantages which we will illustrate in this article. First, it is an initialvalue problem, [qi (1) = exp( i Ht) ~p (t = 0)] so —

specific initial state of interest. The calculation of the action of the propagator, [exp( i Ht) qi ( t = 0)], can be done very accurately and efficiently using the Chebyshev expansion method devised by TalEzer and Kosloff [31—33]and generalized for use in close-coupling calculations where internal degrees of freedom are present by Mowrey, Kouri and coworkers [18—22].This method only requires repeated calculation of the action of the Hamiltothan H on a vector, without having to invert H E as is done in typical time-independent methods. The effort of calculating Hqi is espe—

—

Weizmann Fellow 1989—1990 0010-4655/91/$03.50 © 1991

—

Elsevier Science Publishers B.V. (North-Holland)

D. Neuhauser et al.

/

The application of time-dependent wavepacket methods

cially small when qi is represented on a grid, where we can use fast Fourier transforms for calculating derivatives. The scaling associated with time-dependent propagation can be very low, going as N log N or N1 + 1/D for a problem with N quatum states in D degrees of freedom. Another attractive feature is that a single wavepacket contains a wide and continuous range of energies over which we can extract S-matrix and scattering probabilities. Today, state-of-the-art reactive scattering calculations can yield numerically exact cross-sections for atom—diatom systems [44—46].However, all of the reported theoretical cross-sections have used the powerful time-independent variational [47—52] methods. Our aim has been to develop time-dependent wavepacket methods that can do these same types of reactive scattering calculations, and hopefully do them for higher energies and heavier masses. Our reasoning was that given the slow scaling of time-dependent propagation methods with problem size, time-dependent methods would have to be the most efficient available for sufficiently large problems. For this prediction to come true however, we need to solve two problems. First, we have to develop accurate methods for extracting numerically exact reactive scattering information from our wavepackets (the inelastic problem is well in hand); and second, we have to implement algorithms that are efficient enough for the slow scaling behavior to be a practical advantage i.e. to manifest itself for problems that will fit on today’s computers. The importance of this second requirement is illustrated by the fact that for collinear reactive scattering calculations, time-dependent approaches are much more expensive than the corresponding time-independent methods, but for inelastic diatom—surface scattering, wavepacket methods offer the only viable approach. The formulation of both time-independent and time-dependent reactive formalisms is difficult due to the presence of the different distinct asymptotic arrangements. In most of the time-independent methods, all of the accessible arrangements are coupled together to yield a large set of coupled equations simultaneously spanning all of the accessible configuration space [47—58].It is possi—

461

ble to devise an equivalent time-dependent formulation, but we have instead taken a different approach that is much more efficient in use of time and memory on the computer. In each of our methods, we convert the reactive problem into an essentially inelastic one. At any one time, we only propagate the wavepacket in one set of inelastic Jacobi coordinates. What makes such an approach possible is the use of complex potentials that absorb any part of the wavepacket that would physically leave the region of configuration space over which our grid is defined. In the simplest method, reaction probabilities are determined by calculating the flux of the wavefunction into each product arrangement and then absorbing the reactive portions before they run into the edge of the inelastic numerical grid. Alternatively, the total reaction probability can be determined as the complement of total elastic and inelastic scattering probabilities. We have also developed a method to obtain columns of the reactive S-matrix by propagating in the reactive arrangement and absorbing the inelastic part of the wavepacket with a cornplex potential. The conversion of the reactive problem into an inelastic one enables us to use all the machinery developed for inelastic calculation: the convenient Jacobi coordinates; a close-coupling expansion in the rotational and sometimes vibrational states; and a resolution of the scattering probabilities from a single wavepacket over a continuum of energy or momenta. We have also developed a method to obtain energy-resolved reaction probabilities from the reactive fluxes, using a time energy resolution of the wavepacket (rather than space momentum). We emphasize that our approach is exact, subject to numerical considerations, as long as the complex potential is placed far enough outside the reaction zone that its action does not interfere with the wavepacket before it is analyzed. Another essential ingredient in our method is the use of a projection-operator technique to separate the wavefunction into two parts: (1) a sum of close-coupling (vibration—rotation) states that extend to large distances in the scattering coordinate; and (2) a piece concentrated in the stronginteraction region where the vibrational motion is —~

—

462

D. Neuhauser et al.

/

The application of time-dependent wavepacket methods

represented on a grid. This lets us expend the majority of the numerical effort in the interaction zone, where the dynamics are non-trivial. The projection operator method also allows us to deal efficiently with situations in which many asymptotic channels are open those of no interest can be excluded by blocking them with a complex potential. The use of absorbing potentials (which turn reactive problems into inelastic ones) and of the projection operators are not restricted to time-dependent formulations, and indeed we have implemented them in a time-independent method similar in spirit to the wavepacket methods we will describe here. For details, see refs. [58—61,63]. In this paper we summarize our approach, developing it from simple coffinear models to full three-dimensional reactions for total angular momentum J ~ 0. Much of the material here has been presented earlier [34—40];however, here we develop it in a concrete logical structure with more consideration given to the actual implementation. To avoid tedious repetition, we often omit detailed derivations, and refer the reader to the appropriate paper. Also, some methods are presented only for collinear reactions, especially when

From section 7 onwards we concentrate on aspects that are peculiar to time-dependent methods. We first briefly review in section 7 the timedependent Chebyshev propagation method of Tal-Ezer and Kosloff [31—33]and its extension for use with systems containing internal degrees of freedom developed by Mowrey, Kouri and coworkers [18—22].The sections that follow concentrate on the extraction of scattering information over a range of energies from a single wavepacket. In section 8 we review the extraction of inelastic S-matrix elements at different momenta from the inelastically scattered portion of the wavepacket. Section 9 applies this procedure to obtain reactive S-matrix elements by propagating the wavepacket in the product arrangement. In section 10 we show how total, arrangement-selective, and vibrationally-selective reactive probabilities are obtained from a wavepacket that is solely propagated in the reagent arrangement, by calculating the reactive flux of the wavefunction. In section 11 we show how energy-resolved probabilities are obtained by a procedure based on a time energy resolution of the flux of the wavepacket. Finally, a sample calculation for F + H2 (J = 0) is

the three-dimensional equivalent only involve more indices. The balance of this paper is as follows. In the first part (sections 2—6) we concentrate on the static aspects of the reactive Schrodinger equation. In section 2, the collinear problem is formulated, including the conversion of the reactive system to an inelastic one by the use of absorbing potentials. Section 3 introduces the Projection Operator Formalism, by which the wavefunction is dissected into a mixture of (1) functions where the vibrational motion is represented on a grid and (2) close-coupling functions where the vibrational motion is represented as a superposition of asymptotic diatomic eigenfunctions. The extension from collinear to three-dimensional problems, where we use a close-coupling expansion in the rotational degree of freedom, is given in section 4 and the projection operators in three dimensions are presented in section 5. The absorbing potential is analyzed in section 6, where we derive constraints on its parameter and compare the performance of potentials with different functional forms.

presented in section 12.

—

—~

2. CoUinearjeactions

—

formulation

The collinear atom—diatom reactive problem is defined by A

A + BC (1) AB + C corresponding to inelastic and reactive scattering respectively. The three atoms A, B and C are assumed to lie on a straight line. In traditional reactive scattering formulations, non-orthogonal or hyperspherical coordinates must be employed to account for all of the open arrangements. In our formulation, however, the product arrangement AB + C is artificially closed using an absorbing potential, and we can use the very convenient Jacobi coordinates designed for purely inelastic systems. For collinear problems, these coordinates are defined by +

=

BC

II rB

_

—

rc

II

(2)

D. Neuhauser et al.

/

The application of time-dependent wavepacket methods

and MBrB+MCrC M +M

—

R= ~

B

(3\ ‘

C

,.

where h mass

~

=

-~---~

1, V is the potential, ~ is the reduced

MAMBMC ~

MA + MB +

Mc

)

1/2

(5)

=

(6) ~~

aix,

1

where r1 is the position where the complex potential is first applied, and z~sris the complex potential width. The linear form of the absorbing potential is derived in section 6. The initial wavefunction is located in the reagent region, where the colliding atom A makes a negligible to the diatom BC. In the interactionperturbation region, IN, the wavepacket is strongly distorted, and part of it reacts and then is ab-

and, as usual, the coordinates are rescaled so that only one mass appears in H: R

by IN [36,37]. The product region is eliminated using an absorbing imaginary potential —iVir(r) f r~+~r—r iV r1
where rA and MA denote the position and mass of atom A, and similarly for B and C. In Jacobi coordinates the collinear Hamiltonian reads (after separating the center of mass) 2 8 ~ H= + V( R, r), ~~4) 1 / a + —

463

-

r = —r,

sorbed by the imaginary potential. We separate RE and IN by a line at R = R~,where in practical calculations R 1 3—4 A. Our method utilizes the fact that different representations of the wavefunction are appropriate in the RE and IN regions, as discussions in the following section.

where a=

(MB +

M~)~

MBMC

i/2

(8)

A schematic picture of the potential is drawn in fig. 1. We divide space into three regions: (1) reagents, denoted by RE; (2) products, denoted by PR; and (3) a strong-interaction region, denoted

/7 =

3. Projection operators in collinear reactions The introduction of a complex potential closing off the product arrangement limits the range of the r coordinate to r1 3 A and turns the reactive problem into an inelastic one. In the asymptotic reagent region the wavefunction is aptly represented by a close-coupling superposition of the N asymptotically open vibrational states: qi(R, r, t) where

r1

qi,,

~

=

t)qi~(r),

(10)

are the eigenfunctions of the diatomic

Hamiltonian, h

(IN

R

Fig. 1. A schematic view of a collinear reactive potential in Jacobi coordinates. The product region is eliminated by placing an absorbing potential along r = r1. A close-coupling CX pansion in r is employed in the asymptotic reagent region (RE), and is supplemented by a grid function in the strong interaction zone (IN).

(11) and v(r) is the diatomic potential v(r)

lim V(R, r). R -.

(12)

464

D. Neuhauser et al.

/ The application of time-dependent wavepacket methods

In the interaction region, IN, a close-coupling representation is inexpedient, because many asymptotically closed states are required to account for the reaction. The wavefunction is better represented there on a grid in r, as the cost associated with calculating the Hamiltonian on N grid points in proportional to N log N (in the fast Fourier method [31,33]) compared with N2 in close-coupling methods. We devised a method [36,37] that couples the advantages of both grid and close-coupling descriptions, using the projection operator

and similarly .

8~(R, r,

=

iQ~

t)

iQ{H~~m(r)~m(R~t) m

=

+H~(R, r,

(19)

t))..

These equations are further manipulated to bring them into a convenient form. We define a perturbation potential

N-i

~ I4~)(~~I

(13)

U(R, r)=V(R,

(20)

r)—v(r),

n’=O

which is applied to the asymptotically open states. The wavefunction is divided into two components

and divide the Hamiltonian [eq. (4)] as H=

—

r, t) =Pqi+ Qqi =

~

2~~

~4,,(r)i~,,(R,

t) + x(R, r, t),

where Q is the complement of P, Q = 1—F,

+h(r) + U(R, r).

(21)

(14)

It is easy to see that 82/0R2 and h(r) commute with P, so their contribution vanishes in the PHQ and QHP terms [eq. (18)], as PQ = 0. The PHP terms is expanded as

(is)

(4~~ I PHP I E4),,,(r)~i,,,(R,r, m

and ~,,(R, t)=(~,,IPqi), r, t) = QI~(R,r, t).

(16) (17)

~

1

/

+ ~

82

r, t)

U,,,,,(R)i 1,,(R, r, t),

Note that since X(R, r, t) does not contain any open asymptotic states, it decays exponentially for large R and can therefore be stored on a small grid covering only the strong-reaction IN region. The close coupling functions i~,,( R, t) are nonzero in both the asymptotic-reagent and the interaction regions, RE + IN. The Schrodinger equation idqi/dt = Hqi is then replaced i~,, giving by .8~,,(R,t)

a set of coupled equations for .a =1

=

0t

and

where (~r(R) = (4,, I U(R, r)

IPHIx(R, r,

t))

I4,~).

(23)

The resulting equations are then [36,37] .

0m,~(R) /

1

Bt

1

+~

82

)~~(R~ t)

Unm(R)~im(R,

t)

m

=<~,,IPIHqi)

<4~,,IPHP+PHQIiP) (4,,(r) I PHI ~qi~(r)~~(R, + (qi,,(r)

x

(22)

+

(24)

and t))

(18)

=

Q{H~m(r)iim(R~ t) +HX},

(25)

D. Neuhauser et aL

/

The application of time-dependent wavepacket methods

465

where the operation of Q on an arbitrary vector F is defined by eqs. (17) and (18)

and z~ R is the range of the transition region. We chose this piecewise polynomial form for f since it

QF( R, r)

has relatively small gradients while beingeverywhere smooth enough (continuous second deriva-

=

F( R, r) _~~m(1)J1’(R~

r’)qi~(r’)dr’. (26)

In the Hamiltonian we calculate second derivatives by the fast2flm(R)/0R2 Fourier method Noteuse thata [eq.[33,31]. (24)], we for the transform term 8 over the long R grid (containFourier ing both the the RE + IN regions) while only the small (R, r) grid (the IN region) is used for the second derivative part in H~[eqs. (5) and (19)]. In our calculations for three-dimensional H + H 2 and F + H2 most of the numerical effort in calculating eq. (25) is associated with the H~ term. However, for heavier mass systems or for systems with more degrees of freedom, there will be many open channels, and the dominant numerical effort will be associated with the PHP close-coupling term, two withmethods a cost proportional to 2. We have derived to decrease this N cost, both of which rely on the fact that the N2 close-coupling expansion is most appropriate in the asymptotic region, RE, where the perturbation U(R) is small. In the first method [36] the close-coupling expansion is eliminating from the strong interaction region, IN, by modifying the projection operator p

—~

~ a~f(R)P=f(R)~

I~

I~

(27)

where the function f( R) vanishes in the interaction region, is equal to 1 in the asymptotic region RE, and rises smoothly in a small transition region. We have used the following form for f( R) 3 15y4 + 6y5, R R
f(R)

=

(28)

~,,

function x and the close coupling functions ii,, are restricted to the small transition region, [R1, R1 + L~R]. It can be shown that the size of this region, ~R, can be made as small as a few grid-spacings. For example, in collinear H + H2 full convergence was obtained using L~sR= 0.5 A, only 6 times as large the grid in R Weas have alsospacing derived and[36]. utilized another method to limit the cost associated with q,,,. In this method only a few of the N asymptotically open channels are included in P —

~‘

~

=

n=0

I ‘kn > (qi,, I~

N N open channels, so a complex potential V 1(R) is placed at the reagent region boundary to absorb the open states in Q. Formally, this requires adding a term QVIR(R)Q to the Hamilto—

form (25). The complex in section potential 6: to the hasH~ theterm triangular man, derived i.e. adding VIR(R) in eq. 0,

~oY~ —

R1 (29)

(30)

where N is smaller than the number of open states, N. The grid function x = Qqi does not now vanish automatically for large R, since it includes

VIR(R)=I

where R

tives) at R=R1 and R=R1+z~R. From the modified projection operator, we develop equations for x and as was done in eqs. (17)—(20). The only complication is that now PQ ~ 0. However, the derivation is still straightforward; see ref. [36]The for advantage details. of the fP form is that the closecoupling functions, ~ ( R, t) are zero in the strong-interaction region, and thus the coupling term L~,m(R)needs only be calculated (and stored) for states n and m that are close in energy. In addition, the interaction terms between the grid

R~R1, R1
(31)

where y was defined in eq. (29) and here L~R denotes the range of the complex potential.

D. Neuhauser et aL / The application of time-dependent wavepacket methods

466

In principle, we can restrict P in eq. (30) to include only one channel (corresponding to the initial diatomic state) by placing R1 far enough Out in the reagent region. However, when heavy masses are used (or in three dimensions, where diatomic rotations are included) numerical simulations indicate that a one-channel P requires very large grids for convergence (with total length R1 + z~iR 6 A), and that a better choice is to include in P a few other states that are close in energy to the initial state to account for perturbations in the near-asymptotic region, R 3—5 A. An obvious drawback to the absorption of the N N open states by the complex potential is that the inelastic scattenng information cannot be re tneved for these states in the usual method of momentum-analysis of the spatial wavefunctions sj,, ( R, T), at large times T (see section 8). In many cases, however, inelastic information in not sought only reactive probabilities are required. More importantly, we note that even with a restricted projection operator, full inelastic information can be retrieved using a time energy resolution (rather than space momentum), a method developed in section 11.

j

/

—

—

—

4 ~. . Fig. 3. A schematic view of the potential surface for three-dimensional collisions in the Jacobi coordinates (R, r, y). An absorbing complex potential along r = r1 blocks the product

arrangement.

—

—*

-9

4. Three-dimensional calculations The collinear scattering formulation is straightforwardly transformed to the complete problem of atom—diatom inelastic [19—22,25] and reactive [35,38—40]scattering in three dimensions. We use here the body-fixed three-dimensional Jacobi coordinates (R, r, y), where R and r are defined in section 2 and y is the angle between R and r (see fig. 2). Following ref. [19—22]we utilize a closecoupling expansion in the rotational degree of A~

\ \

R

3). Note that for three-dimensional reactive systems there are two possible reaction products A + BC AB + C AC + B. (32) -~

The imaginary potential blocks both product arrangements. In the main part of this section we trace the well-known derivation of the body-fixed Schrödinger equation, which is then utilized in the Projection-Operator formalism of the next section. The Schrodinger equation is derived from the Hamiltothan (see e.g. refs. [19—22,52])

1 (1 022 ~2 R1 OR2 H=——i—-—r+—--—R 2~s\ r Or

(_\ B•

freedom, y. As in the collinear case, Jacobi coordinates are very convenient for purely inelastic collisions. We again make them suitable for the reactive problem by blocking the product channel through an r-dependent absorbing potential (fig.

+ V(R, r, y),

/‘~~2 ,2 r2 R2 +—t—+— (33)

.0

Fig. 2. Jacobi coordinates for an atom—diatom system in dimensions,

~

where the operators j and 1 denote the angular momentum of the diatom BC and the orbital angular momentum of A relative to BC.

D. Neuhauser et a!.

/ The application of time-dependent wavepacket methods

The full wavefunction is characterized by six coordinates (R, r, y) and the Euler angles (0, qi, 8) which describe the orientation of the ABC triangle relative to the space-fixed coordinates. The wavefunction is decomposed as ~ 2J+1 1/2 1 JM = ~ 4 ~—qi~1(R, r, t) —

(

JM ~

)

—J

>< 0)D~(O, ~ ~) (34) where the 1/Rr factor is inserted to simplify the kinetic energy terms in R and r. J and M are the total angular momentum and its projection on a space-fixed axis, D~ denote the Wigner rotation matrix elements [65] and Y~ are spherical harmonics. In the following we eliminate the (J, M) superscripts on qi~, as they are both conserved by H. We also omit the (n 0, 10’ Q0) subscripts corresponding to the initial state. The Schrodinger equation is now converted to a set of coupled equations for ~

Ot

r, t)

~))I

=

(D~M(O,~

=

(D~M(0, q, 8) I Y~(

~

0) H I

~)

(Hqi ) (35) ,,~

and a straightforward calculation leads to 2 ~ 8~ 82 \ j(j + 1) (Hqi)JQ= ~—(~----~ + ‘~~) + 2p~r j2 +12 2Q2 1

L

J~(R, r)

=

(Y~(y, 0)1 V(R, r,

I

~i’) )‘Q(y,

0)). (39)

The expansion of the y coordinate in terms of spherical harmonics2 ~~2(y, enables convenient and 12 0)terms, at the expense of introducing nonlocal terms into the potential. calculation of the j We could have retained a local potential by representing the wavefunction on a grid in y, by using a set of grid functions F~(R, r, y) defined through the equation FQ(R, r, y)

~ I

r, y)Y,Q(y, 0).

(40)

On a grid, however, the ~2 kinetic term has a (removable) singularity at y = 0. 1 8 8 (41) smy8-y siny—. 8y Several ways to avoid this singularity have been •2

j

— —

—

.

—~

a cosine transform is applied to calculate the derivatives. However, this scheme still increases the maximum bandwidth of the Hamiltonian (Hm~in section 7) which increases the number of gation. Therefore, it is not needed efficient in enough in our Hamiltonian evaluations the propaatom—diatom calculations, at least up to j 50— 100. However in situations where very large values of j are needed or when three-body breakup oc—

—

2~R2

+

We also introduce the close-coupling matrix elements of the potential

advanced recently [43]. In one of these methods the grid points do not include the point y = 0, and

0)181) ~,

467

curs, these methods will provide practical alterna+ +

1

(T~qi

2 1~~1 + Tqi1~1) 2tiR ~ T’~(R,r)qi 1’~,

In to order save storage expansion. we have also examined tives the to close-coupling (36)

another method, where the potential VX(R, r), defined as

where the raising and lowering terms T ~=

[(J(J

+

1)

—

±1

Q(Q ±1))

originate in the expansion of 2J2+j2.J+j_.J_j+

,2

matrix

V~(R,r)

=

f

V(R, r, y)PA(y) d(cos y),

(37)

(42)

(38)

is J’ stored rather than the full potential matrix 1~(R,r). From V~(R,r)it is possible to generate J’~1~(R,r) when it is needed, and the number of necessary X’s is only twice the the number of

2ff

,2(Jj)

D. Neuhauser et a!. / The application of time-dependent wavepacket methods

468

f’s so the storage cost of the potential grows as Imax’ rather than as J~ax~ However, in this method the CPU time increase unacceptably due to an additional summation over A at each Hamiltonian evaluation. A third possibility is to alternate between a grid representation of tJi~(R, r, y) when V~(R, r, y)ijiv(R, r, y) is computed, and a spherical harmonics representation [qi1~(R, r)] when the kinetic term is to be calculated. The back-and-forth transformation between y and j would give rise to a j~ CPU-time cost term, but again storage would only grow linearly with Jmax• This idea has been used by both Light and coworkers in their time-independent work [67].

and ~10(R, r, t)

Q1tJiJQ(R, r,

=

(Q~

t)

1

—

Pd). (48)

By repeating the derivation of the collinear projection-operator formalism [eqs. (18), (24) and (25)], it is straightforward to show that the Schrodinger equation, eq. (35) is equivalent to O~,,JQ(R,t)

=

IHI4~,,’1’)n,,~’v’(R, t)

~

~

8t

‘

n’j’S2’

I~(R, r,

+ ~(4,,~(r) I

5. Projection operators in three dimensions. r, t)

The projection-operator formalism in three-dimensional atom—diatom collisions is analogous to the formulation in the coffinear case [35,38]. In the asymptotic reagent region, the wavefunction is again most compactly described by an expansion in a small number of vibrational states [35,38] ~Pja(R, r) ~J,,JQ(R, t)qi,,1(r), (43) —~

8t

1

><

Q1~U~(R,

=

~

r)4,,1(r)

t) + Q1~10(R,r, t),

(49)

where we defined ~

r,

(so)

(H~)1~(R,r, t)

t)

n where the qi,1(r) are now defined by 2+ v(r))4n / 1 82 + 2~r 1(r) e,,1q,,1(r). (44) —

=

and from eqs. (33)—(39)

I

(4~nj/H

I 4~n’j’) i a~ =

6nn’6jj’612Q’ j2 ~2 — 2~22

2~iR2

~

We define the projection operator P=6~~~ I~)Iqi,,

8~÷1

1(r)~~,,1(r)I
=

(45)

~

T6~~~ 2 1)

+

+8,o~T~ +

2p.R

~

(51)

I.’,

where P 1=

where

~Ii~)(4~I.

(46) =

The states (n, J) extend over all asymptotically open channels. Since P is diagonal in both ~ and J, it is easy to derive the coupled equations for pqi and Pqi and Qqi (where again Q 1 P). We define =

~njv(R,

t)=(q,,1(r)Ii~10(R, r, t))

—

(47)


IU~(R,r)

I~~’~’).

(52)

Note also that in the cross terms (QHP, PHQ) we again replace H by U, as was done to derive eqs. (24) and (25), as all other terms in H commute with P. Equations (49) constitute the required set of time-dependent linearly-coupled equations to be propagated.

D. Neuhauser et aL

/ The application of time-dependent wavepacket methods

6. The absorbing potential

469

For the wavepacket not to be reflected at any point before the end of the grid at x 0, the semiclassical momentum p(x) [2f1(E + iV1(x))]~’2,should vary as little as possible over one cycle of the wavefunction =

In our formulation an essential ingredient is the use of an absorbing imaginary potential in the product channels, which reduces the reactive problem to an inelastic one. The imaginary potential is not chosen arbitrarily or by parameter-fitting; rather, we use two quantitative criteria that together ensure that the reacting portion of the wavepacket will be fully damped by the potential; any imaginary potential satisfying these criteria is appropriate. The first requirement is that the potential should vary smoothly (in a sense that we will discuss below); otherwise, a wavepacket impinging on the optical potential will reflect from it without absorption. In addition, it must have a sufficiently large magnitude, so that a wavepacket transversing it will be fully damped. In this section we briefly review the quantitative implications of the no-reflection and full-absorption constraints, and discuss how to choose a potential in actual calculations. The derivation is based on a one-dimensional model, as drawn in fig. 4 (see also ref. [34]). A negative imaginary potential, iV1(x) with a typical magnitude V10 ~ located in the region x E [0, Z~x].A reflecting wall is present at the point x end 0 which in actual application to the of thecorresponds grid. A wavepacket with a typical energy E is impinging from the right, x> ~x. —

=

=

~p(~v) A ~ 1

(53)

where 1 A

(54)

~.

—

It then follows (ref. [34]) that reflection at a point x is negligible if

I 8V

1/Ox I

~

(E

+ V1 (x

)~

3/2

%/~i’.

(55)

We turn this equation to a constraint on the typical magnitude of V10 by approximating ~ ~ > (56) Ox t~sx and by imposing eq. (55) at x ~x, where the imaginary potential vanishes. Then it follows that 3~2. (57) V~0~ ~ ~xE The second requirement on the potential is that a wavepacket transversmg the whole region [0, ~x] —~

=

~

Fig. 4. Schematic view of the one-dimensional absorbing complex potential

—

iV

1(X) with a typical magnitude V10 and placed in the 2 —A~y3 dependence = (x region [0, ax]. The potential is constructed to damp wavepackets impinging from x > ~x. reflecting wall where (or theyend of—the grid) is placed at x = 0. Solid, and dashed lines present respectively potentials with a y and 2y

D. Neuhauser et a!. / The application of time-dependent wavepacket methods

470

I

(x )2 I

the asymptotic IN channel, where the complex potential is placed on the long one-dimensional grids for q~(R,1) (see section 3). However, along

exp(—2 z~tV10). (58)

(59)

the vibrational coordinate r, the potential should be made as narrow as possible; for example, in total reactivity calculations for H + H2 and F + H2, we find that the start of the complex potential can be placed as close as r1 2 A, so even a complex potential with a range of 1 A will increase the length of the r grid by 50%.and Weintherefore employ a short-range potential, that

where v is the average velocity of the wavepacket. It follows from eqs. (58) and (59) that for significant damping, the potential must satisfy

case neither constraint (no-reflection and full-absorption) is absolutely fulfilled. We then must optimize the shape of the potential by requiring that 8V1/Ox is minimal while maximizing f V, ( x)

(60)

dx [see eqs. (55) and (58)]. This leads to a triangular-shaped potential, as was used in eq. (31)

will be strongly damped. The damping of is approximately equal to

dtl

exp{_2fvi(i(t))

—

~‘

Here, ~(t) denotes the position of the wavepacket center at a time t, and z~tis the time it takes the wavepacket to traverse [0, ~x] back and forth 2 ~x v

v

=

(~) 2E

1/2 ‘

~ Vio.

Equations (57) and (60) constitute the two sought constraints on the height and width of the imaginary potential. Note that the value of V10, determined by the geometrical average of eqs. (57) and (60), is equal to the typical energy of the wavepacket V10 Etypicat. (61) The width of the imaginary potential, ~x, is determined by the range of energies in the wavepacket. Specially, a wavepacket containing energies in the range [E~, Ema.j [and Etypjcai (E~nEmax)1~2] will be absorbed if E1~4 max
—

V1(x)

=

{

~

0
0,

L~X~ X,

z~x x —

=

(63) We have also examined potentials that rise less steeply than linear at the point where the complex potential starts. For instance, the linear dependence on y in eq. (63) can functional be replacedform by a of the or 2, the potential is 2y2 which has the same area higher term. With 2/3y3 y —

0.0 ________________________________

—

—20

For example, in simulations performed on reactive F+ H 2, typical minimum and maximum kinetic energy of the reacting fragments were E,~ 0.02 eV, E~ax 2 eV, so an optical potential width of r ~ 1 A is required. [40] We now come to the question of the shape of the potential. If there is no significant numerical cost in taking a large ~x, both constraints in eqs. (57) and (60) can be easily satisfied. In thatat case the potential should rise very smoothly the point x ~ix where the potential is first applied,

.0

/ _______________________________

—

—

=

so that very little reflection occurs. This is the case, for example, for absorption in the far ends of

~ —

0

—0.5

0.0

0.5

1.0

LOG ( E) Fig. 5. A log—log display of the absorption a function deof 2 — ~~y3as dependence, energy for potentials with y and 2y noted respectively by a solid and dashed lines for 1~sx= 0.5 A and dotted and dot-dashed lines for ~x = 1 A. parameters are V 10 = 0.3 eV; and ~ts= mH/sfi. Note that for a small value of t~xthe linear form is slightly less effective at energies near V10, but is better at very low energies.

D. Neuhauser et a!. / The application of time-dependent wavepacket methods

471

as the linear potential, a maximal derivative at x 0 but a zero derivative at x Lxx. For each of these potentials, we determined the absorption probability as a a function of energy by integrating the time-independent Schrodinger equation (d2i~i/dx2 2p~(V1(x) E)i~i)subject to the con-

The initial wavepacket should be positioned sufficiently far out into the asymptotic region that there is negligible overlap between the wavepacket and the interaction potential. We typically choose

dition qi(x 0) 0. Figure 5 illustrates how the absorption probabilities rise at very low and very high energies; the advantage of the the linear potential for a small t~xis evident. Our analysis is not detailed enough to produce a potential that is optimally absorbing in a global sense. However in practice, the guidelines given in this section usually provide a potential that adequately meets the needs of our calculations.

We now briefly describe the Chebyshev propa-

=

=

=

—

=

=

R0

3A + 8a.

=

(65)

gation scheme [31—33,18—22].In this method, the Schrodinger equation d ~p j~

=

—

i H4’

(66)

for a hermitian Hamiltonian H is formally integrated to give (67)

7. Initial conditions and propagation

The propagator is expanded as

In the previous sections we described how to convert the reactive Hamiltonian into an complex inelastic one, and showed how to efficiently represent the resulting inelastic wavefunctions. The forms of the Hamiltonian and wavefunctions that we arrived at are useful for both time-dependent and time-independent treatments [58—61,63]. In this section, we begin to concentrate on the timedependent aspects of the calculation. Here we comment on the choice of initial wavefunction and then briefly review the Chebyshev expansion of the propagator utilized for calculating the time-dependent wavefunction, eqs. (24), (25) and (49). These formulas will arise again in section 11 when we describe the time energy resolution of the wavepackets. In collinear calculations, the initial wavefunction is simply a Gaussian wavepacket in R times a vibrational eigenfunction

ehhltqi(T)

—~

r, t=O) =‘si~0(R,t=0)4~0(r) 2

=

exp(

—

4a2 (R—R0)

+

=

e~”~ ~ (2

—

8~~)

n=0

x (— i) ~ ~

)T (H_k)

~ (T), (68)

where J~and 7, are the Bessel and Chebyshev functions [64], and (H

+~

)/2,

z~iH

=

(Hmax

—

H~0)/2, (69)

where Hmax and H~0 are upper and lower boundaries on the spectrum of H which can be easily estimated [20—22,32,33].The argument of the Chebyshev polynomials is a rescaled Hamiltoman Xm (H H)/z~H,constructed to have eigenvalues within the region (—1, 1), assuring that the Chebyshev recursion formula —

7~(X)~p=21~1(x)5p— i~2(x)~p

(70)

is stable. In eq. (68), the infinite sum in n is truncated at a value that can be estimated from

ikR)~~ 0(r)~

(64)

where R0, a and k are respectively the center, width and central momentum of the translational wavepacket, and n0 is the initial vibrational level,

the values of ~ and Hmax before the propagation has begun, thereby freeing us from the need to perform repeated convergence checks. In a typical calculation, on the order of 10000 terms are needed for the complete propagation.

472

D. Neuhauser et a!.

/

The application of time-dependent wavepacket methods

The Chebyshev propagation is also applicable to the projection-operator formalism; for example, in collinear calculations, the wavefunction ( R, r, I) is replaced by {~(R, r, t), s~ (R, t)} and H is replaced by the projection-operator Hamiltonian, eqs. (24) and (25) (or eq. (49) in three dimensions). The complex absorbing potential could in principle present a difficulty, since for a complex argument the Chebyshev polynomials are not bounded in the region [—1, 1]. Therefore the expansion of the decaying propagator exp[ i( H

the corresponding eigenstate of the asymptotic reagent Hamiltonian, P2/2.t + h ( r):

iji

—

~ S~(E)=

(‘I’flE(R, r)

I 4’(R,

r, T)),

(73)

where jk R

‘I’~(R, r)

=

e

‘

k~as i/211 (E

/~(r),

—

(74)

—

iV 1) t] may involve delicate cancellation between large terms. However, we find that the expansion is still stable and accurate if ~ H is increased slightly (by about 5%) from its original value, eq. (69), and the time-steps used are not too large, I

~ 10/V~,

(71)

where V10 is a typical value of the complex potential. The main advantage of the Chebyshev expansion is that for sufficiently large time steps t, the number of required Hamiltoman operations in the expansion [eq. (68)] is ‘~max— t ~H.

(72)

Since in practical applications we chose V1 0.5_i eV while L~His at least 10 eV, it follows that t~H~ 100, which is large enough to make the Chebyshev expansion efficient, —

8. Energy resolution — inelastic probabilities Time-dependent wavepackets contain a range of energies. This fact allows simultaneous determination of scattering information at many energies from a single wavepacket. For inelastic scattering, Mowrey and Kouri developed a simple method to obtain energy resolved S-matrix elements [18—22],and we outline their formulation for collinear calculations. The wavepacket i~(R,r, I) is propagated until a large time T such that the scattering process is complete and the inelasticaily scattered part of the wavefunction is receeding away into the asymptotic reagent zone. For a total energy E, the amplitude for transition to a final state n is given by the overlap of ( R, r, I) with iji

and n labels the vibrational state, and aE gives the weight that an energy E has in the initial wavefunction aE

=

f e~0’~,j~0(R, t 0) dR,

(75)

=

where =

~I21i(E

—

,,

).

(76)

°

Note that in three dimensions the expression for aE is complicated by the fact that the eI~~R should be replaced by the eigenfunction of the asymptotic Hamiltoman, in this case a spherical Bessel or Hankel function. A stable method for performing the necessary integrals is outlined in ref. [20]. Direct application of this method requires large grids along the R coordinate, since by the time that the slower components of the wavepacket emerge from the interaction region, the faster components have already advanced far into the RE asymptotic region. This difficulty is remedied by repeatedly extracting that part of the wavefunction that emerges from the interaction region (ref. [19,28,29,36,39]). Specifically, at a series of times t~,the wavefunction is dissected as ~Ji(R, r,

t~)—sf(R)4(R, +

(1

—

r,

t~)

f( R)) iJ~(R,

r,

t,),

(77)

where as in section 5, f rises from f 0 in the interaction region to f I asymptotically. The fiji part is resolved into individual energies and =

=

D. Neuhauser et aL

/

The app!icalion of time-dependent wavepacket methods

473

specific vibrational states, and its contribution to the S-matrix is accumulated

from reagent Jacobi coordinates to product Jacobi coordinates; for collinear collisions

S~0...~1(E) 0 ~ e~1E(~I~flOE(R, r) I fiji(R, r,

P= II rA and

t,)).

—

rB

I

(80)

(78) After each dissection, the propagation is continued with (1 f ) ip replacing i,b. In initial applications, f was taken as a step function [19]; Jackson, Metiu and Heather [28,29] realized that a smooth f is needed to avoid spurious oscillations at low energies and Neuhauser and Baer introduced the polynomial form [eq. (28)], which is found to be most economical, leading to small grid sizes [36]. It is important to realize the different requirement on the f employed here and in section 5. For direct energy-resolution as employed here, f must be different from 0 only in the fully-asymptotic region, and the range, ~ R, of the transition region must be wide enough to ensure that the application of fiji does not contain any incoming momenta. Specifically,

R

=

II ~

MBrB+MArAII MB + MA

—

~,

(81)

—

(79)

~ ~

where ~ is the minimum momenta in the wavepacket. In section 8, where f was used as a part of the projection operator formalism, the requirement spacings) on f wassince muchthere less fiji stringent few grid and (1(~R~ —f)iji awere dynamically coupled.

9. Reactive analysis arrangement

—

where a “hat” denotes the product arrangement, C + AB. The wavefunction is interpolated onto the product arrangement coordinates grid, ~p(R, r, ~) 1p(.k, P, ~) (82) and an absorbing optical potential is set on the reagent side —~

V1(r) V1(P). (83) The transformed wavefunction is then propagated forward in time in the product arrangement. The inelastically scattered portion is absorbed by V1 (P), while the reactive part moves out into the product —~

asymptotic region, and an inelastic-type analysis (as described in the previous section) yields the reactive S-matrix elements. The interpolation of the wavefunction to product coordinates must be done very accurately. Currently we are using Fourier-type interpolation, [38,39] which for collinear reactions reads: t’I’(P, p, i), (84) ~ P, ~) P,p ~ e~”~~” =

where ‘I’(P, p, ~) is the momentum-space representation of the wavepacket in the reagent coordinates, evaluated by a fast Fourier transform

propagation in the product ‘I’(P,

~‘

i)

=

NN R

In our “inelastic” approach the product arrangements are blocked by a complex potential, and it is therefore impossible to apply directly the asymptotic S-matrix analysis for the reactive channels. We have therefore developed two different approaches to obtain reactive information. In one approach [39], the wavepacket is initially propagated as described, in the reagent coordinates. At a time ~ when the wavefunction is fully localized inside the interaction region, we switch

~ T

e

J~r)~p(Rr,

Rr

(85) Here, NR and N denote the number of grid points in R and r. In three dimensions the interpolation is more complicated. This is because the bodyframe wavefunction is written in terms of the projection f~on the body-frame z-axis, which changes when going from one arrangement to another. States of projection Q will transform into a linear combination of states ~. Also, the high

D. Neuhauser et aL / The app!ication of time-dependent wavepacket methods

474

dimensionality of the transformation (four old degrees of freedom transformed to four new ones) makes a one-step transformation prohibitively expensive. We had to factor the transformation into a series of steps to yield an affordable algorithm; see ref. [38,39] for details. In our studies of threedimensional H + H2 and F + H2, the total interpolation time takes between 5% and 30% of the total propagation time, depending on the details of the grid parameters. An essential requirement of this basic scheme is that the wavefunction must be completely contamed in the region of overlap between the inelastic and reactive grids. In practice, this forces the grids to be relatively large, which drives up the cost of the propagation. This is an especially difficult requirement for very exoergic systems (e.g. F + H2), where, unless large grids are used (RIN 5 A), we find that the “front” part of the wavefunction reacts and moves quickly away into the product asymptotic region before the tail of wavepacket has even reached the reaction zone. One solution to this problem is similar in spirit to the dissection of the asymptotic wavefunction for final-state analysis described in section 8, eq. (77). The wavefunction should be dissected a few times into an interpolated product-part iji~(t~)and and a remaining reagent part, iji2(t1) iji1(11)

=giji(t1)

iji2(t1)

=

(1 —g)iji(11),

pletely description of this method will be published separately.

10. Reactive probabilities from fluxes The methods we have been describing for propagation of the wavefunction in the reactive arrangement are appropriate when full S-matrix information is sought. However, often we will be satisfied with more limited information, such as rotationally summed or total reaction probabilities. In this section we describe a complimentary reactive analysis method in which the time-dependent wavefunction is propagated solely in the reagent Jacobi coordinates, and reaction probabilities are calculated by analyzing the flux of the wavefunction in the P direction. Note that besides elimination of the cost of interpolation and the large grids needed, propagation in the reagent channel takes only half as many rotational channels for symmetric systems A + B2 in three dimensions due to the j 0, ±2selection rule. Total reactivity ~R is obtained by integrating the flux along a line r r1. For collinear reactions ~

=

=

—

g(ti)(e_~~huijii(ti)+ e_hhhhiji2(tj))

r=r~ dR

dt,

88

~R

=

dR dt.

(89)

~ Imf ~ ji7t~( R, r, i

Oiji

Q(R

r

Or

(87)

and similarly for ~P2. Another variation on this scheme is to use a reactive grid that is large enough to contain the initial, 0 inelastic wavepacket. In this case, the interpolation can be performed without having to do any inelastic propagation. Our initial trial runs indicate that this scheme will not require much larger grids than the others we have mentioned and will be simpler to use numerically. A com=

t)

and in three dimensions

X =

ImJiji*(R, r, Oiji(R, r, t) Br

(86)

where the function g(R, r) smoothly decreases from 1 in RE to 0 in PR. The two parts of the wavefunction can then be propagated separately on the corresponding grids, until the next dissection time, when we mix them and dissect again:

=

In the projection-operator formalism, iji(R, r, t) can be replaced by X(R, r, t), as the asymptotic open channels are exponentially small at the fluxline r r1. For nonsymmetric systems in three dimensions we can also calculate the partial reactive cross-seetion into each of the two product arrangements (C + AB, B + AC). If the wavefunction emerges with an angle ‘1’ > ~ at r r1 we determine that it reacted into the C + AB arrangement, otherwise it =

=

D. Neuhauser et al.

/

The app!ication of time-dependent wavepacket methods

reacted into B + AC. More formally, an order function D~(y)is defined

where

D~(y)

ljiflfJfOf(R)

{

=

~

y

0,

y>IT,

<~T

(90)

=

BC

—~

Imf ~

~.

(96) —

B + AC)

~ I~7~ ij.’7~( R,

r,

so direct evaluation of state-resolved fluxes would require very large IN grid. However, transitions between different vibrational states cease at much shorter distances (R 2.5 A) so we can use IN

t)

17 ii’

—

Oij,~( R, r, Or

x

t)4~111(P) dP.

Individual state-to-state transitions converge in three dimensions only at large distances, .k 4 A,

and it then follows that ‘~R(A +

f iji1~(R,r,

=

475

t) I I

I

dR dt,

(91)

Ir=ri

grids which are not much larger than the grids required to obtain total reaction cross-sections. Further details are given in ref. [36].

where 2~Tf~

=

0)Y~,’17(y,0)1~,(y)d(cos Y)

11. Reactive fluxes (92)

with similar equations holding for the B + CA arrangement. Individual state-to-state reaction probabilities are also available in a flux form. In collinear reactions

(.k)

=

~-

ImJ~p~’ (R)

o~p,,f( ~) d t,

BR

(93)

where

~

t)

f 4’(R,

r, t)4~~(P) dP,

(94)

where in eq. (94) R and r, the inelastic coordinates, are now taken to be functions of I and P, the reactive coordinates. Note that the vibrationally resolved probabilities require interpolation of onto the product grid (k, ~), as in section 9. However, it is only required at a few values of I at which the flux is to be calculated, as opposed to the product coordinate propagation method of section 9 where the entire wavefunction had to be interpolated. The three-dimensional analogs of eqs. (93) and (94) are iji

(I) ~ Imf’P~J.~Q1 (

flfJ~f

=

~

(I) ‘PflfJf ~1

d t, (95)

—

resolution to individual en-

ergies In the previous section we described how to calculate transition probabilities by integrating over time the outgoing flux of a wavepacket. This method yields an average of the energy-dependent probabilities over the initial energy bandwidth of the wavepacket; thus, these fluxes do not yield direct energy-resolved reaction probabilities, unless we use wavepackets that are narrow in momentum space, which are correspondingly broad in configuration space. Even then, we only get information at one energy per propagation. In order to calculate fluxes corresponding to many energies from a single, narrow (in R) wavepacket, we have developed another procedure. As in ref. [15] we calculate the energy-dependent scattering wavefunctions by a Fourier analysis (from time to energy) of the propagating wave packet

I ‘P (E))

1 =

~fI~(t)>

&Etdt,

(97)

where the proportionality constant aE [see eq. (75)] gives the weight of a component with an energy E inthe the energy-resolved initial translationalwavefunctions wavepacket. From I i/i(E)), the fluxes are obtained using the formulas in section 10, except that now the integral over time in the fluxes is dropped, and instead they are divided by the translational velocity; for

D. Neuhauser et a!. / The app!ication of time-dependent wavepacket methods

476

example, the expression for the collinear total-reaction probability [eq. (88)] is here: Pr(E)

=

~

energy in each Chebyshev step, where N denotes the total number of grid points and basis functions in ‘P. This is to be compared with the 2N log

ImJip*(R, r1, E)

~1~

B ‘P (R, r1, E) X

—

eq. (99) only requires the extra work of accumulating iji(X, E). This takes order N operations per

dR,

(98)

rr~

N cost for each Hamiltonian operation. Thus, for a moderate number of energies to be analyzed (NE 20), the energy resolution does not signifi—

with similar equations for the state-to-state and three-dimensional probabilities. The resolution of the wavepacket to different energies must be done very accurately. In ref. [15] a split-operator formulation was adopted, but we are now using a more accurate procedure, which is based on the Chebyshev propagation [62].First we note that for the fluxes the wavepacket is sought solely within the IN region, where the wavepacket vanishes at time t <0 (t 0 denotes the starting propagation time, in which incoming wavefunction is placed in the asymptotic reagent part). Thus, in eq. (97) we replace the t oc lower limit by t 0. The integral is then divided into sections of length ~T, where L~Tis the length of each Chebyshev propagation period and the expansion of eq. (68) is inserted

cantly increase the propagation time. The storage poses a more serious difficulty here, since for each energy an entire copy of wavefunction must be stored. The required analysis time and storage space are made negligibly small when only total (or arrangement-specific) reaction probabilities are sought, since then only the single slice of the wavefunction at r r1 is calculated. We also need to calculate the spatial derivative of I ‘P(E)) at r

aEI’P(E))=~JdtexP[iE(m 0 XI ‘P(t + m L~T)> ~T+t)J LiT dt exp[iE(m i~T+t)

The calculation B r and of ~Hthe derivative I’P(R,,i,m~T) is performed (101) here efficiently accurately by the Fourier. transform; for a general vector F(R, r), BF(R, r) ~G(r1 r)F(R, r1), (102)

=

=

=

0

~ exp[iEm L~T]~C(E,n)7,

(2

~

eIE~~~T~C(E, n)

‘H

—

—

H

—

r=r~

where (99)

G(r1—r)=~(—ip)exp[ip(r1—r)].

(103)

Note that since only values at one point r r1 are sought, the computation of the gradient in eq. =

where the coefficients (— i)

=

r

H—H ~H )‘Pm~fl~

=

‘P (R, r, E) r

=

—iHtJ I ‘P(m ST))

C ( E, n)

=

—

=

=

=

~,,

)j LiT exp

Xj (L~Ht)dt

. [i(

E

—

—

H) t 1 (100)

are prepared by numerical integration. Since the action of the Chebyshev polynomials is already performed to carry out time-dependent propagation, [see section 7, eq. (68)], the evaluation of

(102) takes N operations, vs. 2N log N in the fast Fourier transform required for the Hamiltoman evaluation. . As another application of flux energy-resolution, we note that it allows a direct evaluation of inelastic probabilities and even S-matrix elements by calculating ‘P (R, r, E) along the line R R1 in the asymptotic reagent region without the need =

D. Neuhauser et aL

to use the section 8.

fiji

/ The app!ication of time-dependent wavepacket methods

accumulation method described in

12. Sample results

—

F + H2, (I = 0)

Our time-dependent reactive formalism was applied to obtain reaction probabilities and rate constants in collinear H + H2, [36,37] reaction probabilities [35,38], and S-matrix elements [39] for the three-dimensional H + H2, (J 0) reactions. Currently we are applying this method to F + H2, for both the J 0 case and for higher Js [40,41]. In this section we review sample results for total reactivity in F + H2 (n0 0, j0 0, J 0), obtained with the S5A potential of refs. [67] The simulations reported here concentrated on total reactivities, for which very small grids are appropriate. Vibrationally-resolved reaction probabilities for the F + H2 reaction were also obtained in this approach, and show good agreement with time-independent results; see ref. [40] for details. For the reactive probability, the xJ17(R, r, I) Q‘P grid extends over the range

I

1.0

0.8

-

0.6

-

I’

477 I

I

I’’

2.20

2.30

-

02

-

=

=

=

=

R E [0.0, R1]

j

=

0, 2, 4,

.

. .

=

[0.0, 3.0],

,32,

-

1.60

1.70

1.80

1.90

2.00

2.10

2.40

Energy (eV)

=

=

r E [0.2, 3.0],

00

(104)

where all distances are in A and we utilize the symmetry of the diatomic H2 to restrict i to even

Fig. 6. Total reaction probabilities for F+ H2 (J = 0) on an SSA potential. Time-dependent (indicated by dots) are superimposed on the results ofresults the time-independent hyperspherical-propagation of ref. [56,571.

[0,15]. The parameters of the initial Gaussian wavefunction [eq. (64)] were a 0.25, R0 6. We ran two wavepackets centered respectively at translational energies of 0.2 eV and 0.4 eV, and for each of these wavepackets we analyzed the reaction probabilities at 20 energies. (Without ad=

=

ditional effort we could have obtained reaction fluxes over a continuum of energies see section ii). One check the that accuracy of the calculations is based on theonfact probabilities at several —

values. The number of grid points in r and R are respectively, 32 and 45. Note again that propagation is performed only in the F + H2 arrangement; total reaction probabilities are obtained from the flux of the wavepacket in the r coordinate, and no transformation to the H + FH arrangement is required. The parameters of the complex potential [see eq. (9)] are

energies between 0.2 and 0.4 eV were calculated in both runs and gave essentially identical answers. The propagation time steps for the Chebyshev

V10=0.4eV,

calculations of Kress, Bacic, Parker and Pack [56,57]. The results are in satisfactory agreement less than 2% difference. Our two runs took 10 minutes each on a Cray-2. Due to the slow scaling of the time-dependent runs as a function of the total angular momentum, J, we are confident of being able to calculate total integral reactive cross-sections with this method. Runs with up to

[r1, r1

+ ~r]

=

[2.0, 3.0].

(los)

The close coupling functions ~ ( R) include all open vibration—rotation states with energies that are at most 0.8 eV above the H2 ground state, i.e. (n, j) (0, 0), (0, 2), . . . , (0, 10), and (1, 0), . . . ,(1, 6). The grid for ~nj( R) covers the region R C =

expansion [eq. (68)] were z~T=5 x i0’~ s, and the wavepacket was run until a time Ttinai 3 X 10-13 ~. Figure 6 shows the reaction probabilities obtained from our method compared against results of time-independent hyperspherical propagation =

—

478

D. Neuhauser et aL / The application of time-dependent wavepacket methods

20 have already been performed and will be reported separately.

J

=

13. Conclusions In this last section, we show how our methods emerge into an overall scheme for improving the efficiency of reactive scattering calculations. We also mention some other promising ideas that we have not discussed. Our goal is to perform reactive scattering calculations in three dimensions and ultimately to get cross-sections for systems of practical interest. This is such a big problem that the only way to make progress is to break it up into smaller, more manageable pieces. An obvious example is the decomposition of the total wavefunction into states of definite total angular momentum J which are totally decoupled. This turns the problem of calculating cross-sections into a set of smaller but independent problems. Another decoupling that we take advantage of in the time-dependent framework is between different initial states. Only those of interest need to be propagated, which offers an advantage over noniterative time-independent methods were all final and all initial states are rigorously coupled together. Our basic idea has been to look for other, approximate, but physically reasonable decoupling schemes that can be implemented in such a way that the numerically exact nature of the calculations is retained. The principle one we have described here is the decomposition of the multiarrangement reactive scattering problem into a set of approximately decoupled inelastic problems. This decoupling arises because the wavepackets we use are localized in space at any given time the wavepacket only lives in a few localized regions of configuration space. Once the pieces of the total wavepacket reach different asymptotic arrangements, they are rigorously decoupled. While they are together in the strong interaction region, they are coupled but can be conveniently and efficiently represented in terms of any one of the asymptotic arrangement Jacobi coordinate sets. To turn this idea into a practical computational method requires several things. First and fore—

most, we must be able to treat parts of the wavepacket that leave our localized computational grid in such a way that the portion remaining retains its correct behavior. Our main tool here is the complex potential that absorbs any part of the wavepacket just before it leaves the computational grid. An alternative method is to use a cut-off function, the f of section 8. Periodically, the wavepacket is divided into two pieces, one close to the edge of the grid and the other in the interior. The edge piece is then thrown away and the interior piece is propagated for another time step. In certam applications, this method works as well as the complex potential. A second requirements is that we have an efficient propagation scheme within the inelastic region. Our basic tool is the Chebyshev propagation method [31—33,18—22]that has the advantage of being stable, easy to program and efficient. Additionally, we choose the representation of the wavefunction that is most appropriate in each particular region and use the projection operator method to couple the different representations and regions together. There are of course other representations and other propagation schemes available that may have advantages over the ones we discuss. In particular, the split operator method of Feit and Fleck [7] and the frozen Gaussian basis set method of Heller and co-workers [4] are both elegant and potentially useful in our overall computational scheme. Both of these can be used with time-dependent potentials that can arise either from a physical cause (e.g. when light-molecule interactions are allowed) or a computational cause (as occurs in certain semiclassical coordinate decoupling methods.) The Chebyshev method may lose its great advantage in efficiency in these cases. The hyperspherical coordinate representation [53— 57] could also be used. A related technique that take advantage of the localized nature of the wavepacket to reduce computational expense is the floating grid method [12,26], where the grid moves along with the wavepacket. In the principle, one could use several grids simultaneously, spawning off a new one each time the wavepacket bifurcates. This would involve some complicated book-keeping but the reduction in the overall size of the grids might make it worthwhile.

D. Neuhauser et aL

/ The app!ication of time-dependent wavepacket methods

The third requirement for our scheme to work is an accurate method for extracting final-state information from the wavepackets. We do this either by manipulating the wavefunction itself to get S-matrix elements or by first calculating the flux of the wavepacket across some boundary and then extracting probabilities but no phase information. The advantage of the S-matrix method is that we can ultimately produce state, energy and angularly resolved cross-sections. The flux methods have the advantage of being somewhat more robust numerically because small errors that creep into the phase of the wavefunction cancel in the final analysis step. In the S-matrix method, these phase errors translate directly into errors in the S-matrix. Finally, we will mention one other class of approximate decoupling schemes that the time-dependent method allows us to take advantage of. First, under typical chemical conditions, sudden changes in rotational level with large L~fare unlikely during a short interval of time only transitions of a few j levels are seen. This arises because the potential matrix is strongly banded in j, i.e. the magnitude of the j—j’ coupling elements drops off rapidly away from the diagonal. We can exploit this fact by neglecting terms in U~when j and j’ differ by more than say 8 or 10 in certain regions of (R, r) space. To maintain accuracy, we must also use small time steps, but the effective time step in the Chebyshev method is already sufficiently small. For a system with 100 open .1 states, this approximation can reduce the computational expense by a factor of 10. A second, similar idea exploits the fact that the size of the basis set needed while the wavepacket is (bnefly) in the strong interaction region is larger than is needed either before or after. Lemoine and Corey [17,33]use this fact by having the size of the basis set be a dynamical parameter. When more basis functions are needed, they are added, and later removed when they are no longer needed. In this method, the savings are greatest during the first portion of the propagation when the wavepacket has little probability outside of the initial state. In our approach, such a separation can be effected, e.g. by initially including only the close-coupling functions, and only incorporating the Q’P part —

479

once the wavefunction arrives in the strong-interaction region. What we have presented here is by no means the last word in exact time-dependent reactive scattering. The first three-dimensional results for J 0 were published barely two years ago, which should be contrasted with the time-independent methods that first produced cross-sections more than 15 years ago. We anticipate that in a few years, our methods will be but one of a family of practical time-dependent schemes. =

Acknowledgements D.J.K. and R.S.J. acknowledge partial support by NASA—Ames Research Center under grant NAG 2-503. D.J.K. also acknowledges support from the R.A. Welch Foundation under grant E-608. D.N. acknowledges partial support by the Office of Naval Research. R.S.J. acknowledges partial support by the Departement of Energy under contract DE-ACO4-76DP00789.

References [1] J. Mazur and R.J. Rubin, J. Chem. Phys. 31(1959)1395. [21E.A. McCullough and R.E. Wyatt, J. Chem. Phys. 54 (1971) 3578, 3592. [3] Ch. Zubert, T. Kamal and L. Zulike Chem. Phys. Lett. 36 (1975) 396. [4] E.J. Heller, J. Chem. Phys. 62 (1975) 1544; 65 (1976) 979; 68 (1978) 3891; 75 (1981) 2923. [51E. Kellerhals, N. Sathymurthy and L.M. Raff, J. Chem. Phys. 64 (1976) 818. P.M. Agrawal and L.M. Raff, J. Chem. Phys. 74 (1981) 5076. [6] K.C. Kulander, J. Chem. Phys. 69 (1978) 5064. J.C. Gray, GA. Fraser, D.G. Truhlar and K.C. Kulander, J. Chem. Phys. 73 (1980) 5726. A.E. Orel and K.C. Kulander, Chem. Phys. Lett. 146 [7] M.D. Feit and J.A. Fleck, J. Chem. Phys. 79 (1983) 301. [8] Z.H. Zhang and D.J. Kouri, Phys. Rev. A 34 (1986) 2687. [91T. Joseph and J. Mans, Mol. Phys. 58 (1986) 1149. (10] S. Thareja and N. Sathymurthy, J. Phys. Chem. 91(1987) 1970. [11] D. Huber and EJ. Heller, J. Chem. Phys. 87 (1987) 5302 89 (1988) 4752 D. Huber, E.J. Heller and R.G. Littlejohn, J. Chem. Phys. 89 (1988) 2003.

480

D. Neuhauser et aL

/ The app!ication of time-dependent wavepacket methods

D. Huber, S. Ling, D.G. Imre and E.J. Heller, J. Chem. Phys. 90 (1989) 7317. [12] V. Mohan and N. Sathymurthy, Comput. Phys. Rep. 7 (1988) 213. [13] G. Drolshagen and E.J. Heller, J. Chem. Phys. 79 (1983) 2072. [14] R.C. Mowrey, F. Bowen and Di. Kouri, J. Chem. Phys. 86 (1987) 2441. [151 R. Viswanathan, S. Shi, E. Villalonga and H. Rabitz, J. Chem. Phys. 91(1989) 2333. [16] F. Bowen, R.C. Mowrey and Di. Kouri, J. Chem. Phys. to be submitted. [17] D. Lemoine and G.C. Corey, J. Chem. Phys., in press. [18] R.C. Mowrey and D.J. Kouri, Chem. Phys. Lett. 119 (1985) 285. [19] R.C. Mowrey and D.J. Kouri, J. Chem. Phys. 84 (1986) 6466. [20] Y. Sun, R.C. Mowrey and D.J. Koun, J. Chem. Phys. 87 (1987) 339. [21] Y. Sun and D.J. Kouri, J. Chem. Phys. 89 (1988) 2958. [22] Y. Sun, R.S. Judson and Di. Kouri, J. Chem. Phys. 90 (1989) 241. [23] D. Hoffmann, 0. Sharafedin, R.S. Judson and D.J. Kouri, J. Chem. Phys. 92 (1990) 4167. [24] 0. Sharafedin, R.S. Judson, D.J. Kouri, and D. Hoffmann, J. Chem. Phys. 93 (1990) 5580. [25] R.S. Judson, C-K. Zhang, D.J. Kouri and R.L. Jaffe, in preparation. [261 C. Leforestier, Chem. Phys. 87 (1984) 241; and in: Theory of Chemical Reaction Dynamics, ed. D.C. Clary (Reidel, Dordrecht, 1986) pp. 235—246. [27] 5. Sawada, R. Heather, B. Jackson and H. Metiu, J. Chem. Phys. 83 (1985) 3009; B. Jackson and H. Metiu, J. Chem. Phys. 85 (1986) 4192. [281 B. Jackson and H. Metiu, J. Chem. Phys. 86 (1987) 1026. [29] R. Heather and H. Metiu, J. Chem. Phys. 86 (1987) 5009. [30] F. le Quere and C. Leforestier, J. Chem. Phys. 92 (1990) 247. [31] R. Kosloff and D. Kosloff, J. Chem. Phys. 79 (1983) 1823; J. Comput. Phys. 52 (1983) 35. [32] H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [33] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [34] D. Neuhauser and M. Baer, J. Chem. Phys. 90 (1989) 4351. [35] D. Neuhauser, M. Baer, R. Judson and D.J. Kouri, J. Chem. Phys. 90 (1989) 5882. [36] D. Neuhauser and M. Baer, J. Chem. Phys. 91 (1989) 4651. [37] D. Neuhauser and M. Baer, J. Phys. Chem. 93 (1989) 2872. [38] D. Neuhauser, M. Baer, R. Judson and Di. Kouri, J. Chem. Phys. 93(1990) 312. [39] R. Judson, Di. Kouri, D. Neuhauser and M. Baer, Phys. Rev. A 42 (1990) 351. [40] D. Neuhauser, M. Baer, R. Judson and D.J. Kouri, Chem. Phys. Lett. 169 (1990) 372.

[41] D. Neuhauser, M. Baer, R. Judson and Di. Kouri, J. Chem. Phys. to be submitted. [42] D. Neuhauser, R.S. Judson, R.L. Jaffe, M. Baer and Di. Kouri, Chem. Phys. Lett., in press. [43] M.D. Feit and J.A. Fleck, Opt. Lett. 14 (1989) 662. R. Ratowsky, J. Comput. Phys., in press. B.R. Johnson, J. Chem. Phys. 92 (1990) 574. CE. Dateo, V. Engel, R. Almeida and H. Metiu, Comput. Phys. Commun. 63 (1991) 435 (this volume). [44] M. Mladenovic, M. Zhao, D.G. Truhiar, D.W. Schwenke, Y. Sun and Di. Koun, J. Phys. Chem. 92 (1988) 7035. N. Blais, M. Zhao, D.G. Truhlar, D.W. Schwenke and D.J. Kouri, Chem. Phys. Lett., in press. M. Zhao, D.G. Truhiar, D.W. Schwenke and D.J. Kouri, J. Phys. Chem., in press. [45] J.Z.H. Zhang and W.H. Miller, J. Chem. Phys. 88 (1989) 4549; J. Chem. Phys. 89 (1989) 4454; Chem. Phys. Lett. 159 (1989) 124; J. Chem. Phys. 91 (1989) 1528. [46] R. Wyatt and D. Manolopolous, Chem. Phys. Lett. 59 (1989) 129. [47] A. Kuppermann and PG. Hipes, J. Chem. Phys. 84 (1986) 5962. P.G. Hipes and A. Kuppermann, Chem. Phys. Lett. 133 (1987) 1. [48] (a) K. Haug, D.W. Schwenke, Y. Shima, D.G. Truhiar, J.Z.H. Zhang and D.J. Kouri, J. Phys. Chem. 90 (1986) 6757. (b) J.Z.H. Zhang, Di. Kouri, K. Haug, D.W. Schwenke, Y. Shima and D.G. Truhlar, J. Chem. Phys. 88 (1988) 2492. [49] F. Webster and iC. Light, J. Chem. Phys. 85 (1986) 4744; 90 (1989) 265 and 300. [50] D.W. Schwenke, K. Haug, D.G. Truhlar, Y. Sun, J.Z.H. Zhang and D.J. Kouri, J. Phys. Chem. 91 (1987) 6080. D.W. Schwenke, K. Haug, M. Zhao, D.G. Truhiar, Y. Sun, J.Z.H. Zhang and D.J. Kouri, J. Chem. Phys. 92 (1988) 3202. [51] J.Z.H. Zhang and W.H. Miller, Chem. Phys. Lett. 140 (1987) 329; J. Chem. Phys. 88 (1988) 449. [52] M. Baer, Theory of Chemical Reaction Dynamics, Vol. I (CRC, Boca-Raton, Fl, 1985). M. Baer and Y. Shima, Phys. Rev. A 35 (1987) 5252. M. Baer, J. Phys. Chem. 91 (1987) 5846. D. Neuhauser and M. Baer, J. Chem. Phys. 88 (1988) 2858. M. Baer, J. Chem. Phys. 90 (1989) 3043; Phys. Rep. 178 (1989) 99. [53] J. Linderberg and B. Vessel, Int. J. Quant. Chem. 31 (1987) 65. J. Linderberg, Padkjaer, Y. Ohm and B. Vessel, J. Chem. Phys. 90 (1989) 6254. [54] G.A. Parker, R.T. Pack, B.J. Archer and RB. Walker, Chem. Phys. Lett. 137 (1987) 564. R.T. Pack and G.A. Parker, J. Chem. Phys. 87 (1987) 3888. [55] J.M. Launay and B. Lepetit, Chem. Phys. Lett. 144 (1988) 346.

D. Neuhauser et aL

/ The application of time-dependent wavepacket methods

[56] J.D. Kress, Z. Bacic, GA. Parker and R.T Pack, Chem. Phys. Lett. 157 1989) 484. [57] Z. Bacic, J.D. Kress, GA. Parker and R.T Pack, J. Chem. Phys. 92 (1990) 2344. [58] D. Neuhauser and M. Baer, J. Chem. Phys. 92 (1990) 3419. [59] D. Neuhauser and M. Baer, J. Phys. Chem. 94 (1990) 186. [60] D. Neuhauser, M. Baer, and D.J. Kouri, J. Chem. Phys. 92 (1990) 2499. [61] D. Neuhauser, C.Y. Ng and M. Baer, Chem. Phys. Lett., in press. [62] D. Neuhauser, J. Chem. Phys. 92 (1990) 264. [63] M. Baer, D. Neuhauser and Y. Oreg, Trans. Faraday Soc., in press.

481

[64] M. Abramowitz and IA. Stegun, eds., Handbook of Mathematical Functions (National Bureau of Standards DC, Washington, 1964). [65] D.M. Brink and G.R. Stachier, Angular Momentum (Oxford Univ. Press, Oxford, 1971). [66] iC. Light, I.P. Hamilton and J.V. Lull, J. Chem. Phys. 82 (1985) 1400. [67] F.B. Brown, R. Steckler, D.W. Schwenke, D.G. Truhiar and B.C. Garrett, J. Chem. Phys. 82 (1985) 188. R. Steckler, D.G. Truhlar and B.C. Garrett, J. Chem. Phys. 82 (1985) 5499.