Distributed approximating function approach to time-dependent wavepacket propagation in 3-dimensions: atom-surface scattering

Computer Physics Communications

Computer Physics Communications 80 (1994) 1-16 North-Holland

Distributed approximating function approach to time-dependant wavepacket propagation in 3-dimensions: atom-surface scattering Youhong Huang 1 and D o n a l d J. Kouri 2 Department of Chemistry and Department of Physics, University of Houston, Houston, TX 77204-5641, USA

Mark Arnold, Thomas L. Marchioro II and David K. Hoffmann

3

Department of Chemistry and Ames Laboratory 4, Iowa State University, Ames, IA 50011, USA

The theoretical formalism of the distributed approximating functions (DAF) is applied to solve accurately 3D-atom-surface scattering problems. Formulated in coordinate space, the DAF approach starts from an entirely new idea: providing a "uniform" approximation everywhere to a wavepacket, and results naturally in a near-local or banded free propagator. The banded Toeplitz structure of the D A F free propagator matrix on a uniform grid makes possible the application of the most efficient codes in the matrix-vector multiplication in evolving the wavefunction of a quantum system in time, and with extremely small memory requirements. The numerical study conducted in this paper demonstrates that the DAF method outperforms the most powerful available FFT method both in CPU time and storage requirements. The DAF approach gives the same accurate results as the FFT does, and, in some cases, yields more accurate results.

I. Introduction Recently, great emphasis has been placed on developing numerical methods for solving realistic problems in quantum dynamics, using both the time-independent and time-dependent Schr6dinger equations (TISE [1] and TDSE [2]) using, e.g., the discrete variable representation (DVR) [3] and the fast Fourier transform (FFT) techniques [4-6]. Various methods have been successfully applied to atom-diatom dynamics, including scattering problems. The TDSE yields dynamic information without requiring the solution of a large linear algebraic system of equations as is required in the most popular time-independent methods. This can result in both a slow scaling of the necessary computational effort with problem size, and stability of the method when applied to large scale problems. Instead, one could imagine numerically integrating the TDSE directly by, say, evaluating the Feynman path integral from which all physically observable information can be extracted [7]. Approaches to the TDSE can be formulated in either the coordinate space or momentum space. The potential energy operator is local in the coordinate space while the kinetic energy operator is local in the momentum space. The extreme non-locality of the free particle evolution in configuration space has prevented use of the TDSE from the standard Feynman path integral formalism [7]. In particular, the Correspondence to: Y. Huang, Department of Chemistry and Department of Physics, University of Houston, Houston, TX 77204-5641, USA. 1 Supported under National Science Foundation Grants CHE89-07429 and CHE92-00598. 2 Supported in part under R.A. Welch Foundation Grant E-0608. 3 Supported in part under National Science Foundation Grant CHE92-0167. 4 The Ames Laboratory is operated for the Department of Energy by Iowa State University under Contract No 2-7405-ENG82. 0010-4655//94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved

SSDI 0 0 1 0 - 4 6 5 5 ( 9 3 ) E 0 0 9 6 - 6

2

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

free propagator e -ig~/h in quantum mechanics can be realized in the coordinate representation in the form [7] m ( x l e-ig~/hlx') = ( ~ l

] 1/2 eim(x-x')2/(2hr)'

(1)

where the eigenvector I x ) denotes a the position eigenstate,/~ = t~2/2m, is the kinetic energy operator, /~ is the momentum operator, m is the mass and r is the time step. Equation (1) gives equal relative probability over the entire space. The propagator is completely non-local due to zero uncertainty in position for the states I x ' ) and ( x I, yielding infinite uncertainty in momentum space. In a physical problem, a wavefunction after a free propagation has the form 6 ( x , t + ~') =

f

ot~


le-~/c'/nl x')

qJ(x', t) d x ' ,

(2)

The matrix elements ( x l e x p ( - iI£r/h) I x'), with increasingly large I x - x ' I, due to increasingly large momentum, add up destructively and their effects cancel out during the evolution in time. This numerical cancellation, however, is tremendously costly and makes such calculations impossible in practice [7]. To solve the problem, the F F T approach involves discretizing the phase space associated with the position and momentum such that the uncertainties in both representations are finite. The irrelevant high momenta are, thereby, eliminated [7g]. The realization of the F F T method makes use of the orthogonal basis functions, {ei2=ik/N;k=--(ln--1),..

.,In ~ },

(3)

for j = 0, 1. . . . . N - 1. The correspondences between the indices j, k and the physical variables, position xj and momentum Pk, are given by L

x~ =j-~,

2"rr

Pk = hk--£-.

(4)

This is equivalent to defining a square window [-Pmax, Pm~x] in the discretized momentum space, where

Pma~= h2avN/L and L is the length of the 1D physical domain. The F F T algorithm enables the calculation of the action of kinetic energy operator on the wavefunction with the computational effort scaling as @(Nlog2N). However, the global character of the FFT, by which the calculation of the value of the wavefunction at a grid point involves all the values on the grid, limits the method for the calculations of high dimensional time-dependent problems. The global nature of the F F T can be viewed as a fundamental consequence of the use of exp(ipxx/h) or exp(i2~jk/N) as the basis, plus the fact that an interpolation condition holds for the discrete representation of the wavepacket. This condition requires that the discretized Fourier representation of the wavepacket be exact on the coordinate grid points. This leads to the Fourier representation of e x p ( - i I ~ / h ) which is a full matrix, hence the "global property" of the FFT. Furthermore, the global character also limits the F F T with respect to performing efficient calculations on distributed memory massively parallel computers, due to the communication time required. Communication time and memory requirements are two of the most important criteria for a competitive numerical method in the coming era of massively parallel distributed memory machines [8-13]. In a recent series of papers, we have developed an entirely new theoretical formalism, that of distributed approximating functions (DAF), which leads to a novel, powerful approach for accurately evolving the wavefunction of a quantum system in time [14-23]. To avoid the global property attendant to eq. (1) and FFTs, and make the calculation as efficient as possible, the D A F has been constructed to

Y. Huang et aL / DAF approach to time-dependent wavepacketpropagation in 3D

3

provide a "uniform" approximation everywhere to the wavepacket, as opposed to giving it exactly on the grid points. This relaxation of the interpolation property leads to the introduction of an alternative "basis" comprised of " H e r m i t e functions" (the products of Hermite polynomials times their generating function). This leads to a much more localized (in coordinate space) approximate discretized DAF propagator. This banded operator is obtained at the cost of no longer having a unitary propagator. As a consequence, one may propagate safely only for a finite number of time steps. Said differently, the DAF-propagator is constructed to be used only for a specific class of wavepackets (the "DAF-class") comprised of those which can, to a predetermined accuracy, be approximated as an ( M + 1)th degree polynomial under the D A F envelope [14-17]. Thus, in place of a unitary DAF, one has a DAF which is the identity for ( M + 1)th degree polynomials. Since no L2-wavepacket is ever exactly a polynomial, the DAF propagation has an approximate nature. One may propagate only for a finite number of time steps; the parameters of the DAF can be chosen such that propagation can be as long as is required for a given system. It is also interesting to note that in Cartesian coordinates the DAF-free propagator is a convolution, and is diagonalized by a discrete Fourier transform. Thus, in " m o m e n t u m space", the DAF propagator is diagonal, and due to the occurrence of the Gaussian generator of the Hermite polynomials, it is concentrated in a limited region of momentum space. The Hermite polynomials tend to broaden the DAF in momentum space, but the Gaussian ultimately forces it to zero (actually quite rapidly, but smoothly, after a certain point) as the momentum is made larger. Of course, the choice of the Hermite functions to represent the DAFs was originally guided by the desire that the action of the free evolution operator on the D A F representation of the wavepacket should be known analytically. As a consequence, since D A F is diagonal in momentum space, one may view the DAF also as "prefiltering" the propagator with a "window" function in the momentum representation, giving rise to a highly banded Toeplitz matrix structure for the free propagator both in continuous and discretized coordinate spaces. This special structure of the D A F matrix makes possible the application of the most efficient codes available for the matrix-vector multiplication, with extremely small amount of memory requirements. The bandedness or near-locality of the DAF makes it highly suitable for massively parallel, distributed memory machines because of the markedly lower need for interprocessor communication. Clearly, the successes enjoyed in recent years by the F F T and DVR methods provide a practical measure against which any new approach must be judged. Our objective here is to compare the DAF and F F T approaches. The comparisons we present were carried out on a vector processing supercomputer. As noted above, due to the lower interprocessor communication requirements, we expect the DAF to be even more efficient relative to the F F T on massively parallel supercomputers. This paper is organized as follow: The formalism of the DAF and of the atom-surface scattering problem will be presented in section 2. In section 3, results of the numerical studies are given, which include a fast algorithm for the banded DAF matrix-vector direct multiplication, and the numerical results of two examples of atom-surface scattering by both the F F T and the DAF approaches are shown. In section 4 we present our conclusions.

2. Formalism

2.1. Distributed approximation function (DAF) We have mentioned in the introduction that the DAF introduces a "window" function, ~,/2 (k~2(0) p2)n g(M)( p ) = E e-~2(°)P2/z, n=O n!

(5)

4

Y. Huang et a L / DAF approach to time-dependent wavepacketpropagation in 3D

in its "eigenspace" to prefilter out the uninteresting high momentum components. The smooth, Gaussian decay reflects the non-unitary approximate nature of the DAF as the identity for M t h degree polynomials. The parameters, M, o-(0) and coordinate spacing Ax (which for discretized DAFs are not independent) must be chosen to ensure that the DAF can be applied as many times as needed for the particular application. The filtered wavefunction is given by

o(M)( p ) = g(M)( p ) tO(p)

(6)

in the momentum space. The width of the filter is controlled entirely through the parameters 0.(0) and M. In the limits of 0.(0) tending to zero, or M tending to infinity, the width is infinite, g(m)(p) = 1, and all momenta pass through the filter. The width of the filter is determined by the DAF-class of the wavepacket, tO(p, 0) and the potential energy distribution of the system. After a single free propagation step, the resultant wavefunction and the D A F function are

to(M)(p, t+'c) = p ( U ) ( p , r ) to(p, t)

(7)

and M/2

(10.2(0) p2) n

p(M)(p, ¢) = E

n!

n=O

e-~2(~)P2/2,

(8)

respectively, where

o'2('c) = 0.:(0) + izh/m

(9)

is the time dependent "width" of the Gaussian in coordinate space. That the DAF propagator is non-unitary is reflected in the fact that the modulus of its eigenvalues, I ff(M)(p, Z) I, ranges smoothly from zero to one. The eigenvalue is one only at p = 0. This is in stark contrast to the Fourier-based approach where the modulus of the eigenvalues is 1 for Pmax --< P -- Pmax" To obtain the corresponding coordinate representation form of eq. (7), given by

o(n)(x, t + z ) = l

F(n)(x-x'lz)

(lO)

to(x', t) dx',

-

we note that

F(M)(x -- X' [ 'r) =

e-ip(~-~')h d p

1 - 0.(0) exp

( (x-x') 2 20.:(r)

1

n=0

(2

(X--X')

n!:)-1/2H2n v 0.(z)

'

(11) which is the time-dependent D A F function in the coordinate space. H2n is a Hermite polynomial [14-17]. Equation (10) is a convolution fitting. In the DAF approach, the uncontrolled oscillations are circumvented by the choice of the finite value of M and non-zero value of a(0) according to the physical energy distribution in the system and the initial wavefunction. In the limit as 0.(0) goes to zero or M to infinity, the time evolved D A F function eq. (11) becomes the free propagator in eq. (1) and the fitting becomes exact. The free propagator is short-ranged in I x - x ' l due to the Gaussian in eq. (11), which reflects the smooth filtering out of excessively high momenta. The calculation of the quantum dynamics in the coordinate space is facilitated by the near-local property of the DAF function and the local

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

5

property of the potential energy operator. The DAF approach in continuous space can be performed by Monte Carlo integration of the real time Feynman path integral [16]. In this paper we focus on a uniform grid quadrature representation of the convolution in eq. (10). A Gauss-class quadrature DAF representation will be presented in a subsequent paper. In the uniform grid or trapezoidal quadrature representation, the integral (10) is approximated by a sum on a uniform grid, namely ~b(M)(j, t + ~') =

~F~M)(j -- kit) ~O(k,t),

(12)

k

where the discretized DAF propagator,

F(M)(j -

k I t ) = ----(-~-exp

M/2

2o.2(, )

o.(O) 1

X n=0E( 4./ to.(T) ]

1(2.n!2)_1/2 ( g2n

) "~c~'--('~

,

(13)

is a banded symmetric Toeplitz matrix. Here, A is the x-grid spacing. In calculations, we define a matrix Fj(fft) by

Fj~M)={F(M'(j-k[r)

f°rlj-k[
(14)

otherwise, The value of w is based on the error which can be tolerated. The bandwidth of the DAF matrix is thus defined by 2w + 1. The width of the banded matrix is controlled by the parameters M, o./A and r. For a DAF-class wavefunction, the approximating scheme is asymptotic in the sense that an increase of the value of M beyond an optimal value leads to a deterioration of the degree of the accuracy [18]. The larger the value of o.(0)/A is, the larger is the optimal value M and the higher the degree of the accuracy. The asymptotic nature of the approximation results from the uniform grid representation. The bandedness of the DAF matrix leads to extremely small memory storage requirements, and it localizes boundary effects. In comparison, the FFT approach accurately represents a wavefunction whose momentum is bounded by [ Pmax [ in terms of the discrete Fourier series, where Pmax = ,rr/A. However, the non-locality of the Fourier series (arising from the interpolation constraint and the structure of the Fourier basis) results in the delocalization of the boundary effect, and in large storage requirements. For a uniform quadrature grid, the DAF parameters, o.(0) and M, are related to each other in an approximately linear fashion [24], and are thus to be viewed as a parameter-pair, constrained by the requirement that all the eigenvalues of the DAF matrix be less than or equal to one. The Toeplitz structure of the DAF matrix enables one to obtain an analytic solution of the eigenvalues for the infinite matrix [25]. Another important feature of the application of the discretized DAF approach to the TDSE is that the overall error in wavepacket propagation is linearly dependent on the number of the free propagation steps [25]. We are focusing on the method of the symmetric split operator [26] in this paper, so that g,(t + r) = e-i/4~/h~b(t) = e -ip~/zh e -ig~/h e-iP~/2h~b(t)

(15)

in general and =

k

-)

t),

(16)

6

Y. Huang et a L / DAF approach to time-dependent wavepacket propagation in 3D

in the DAF representation. The generalization to higher dimensions in the Cartesian coordinate is straightforward.

2.2. Atom-surface scattering The T D S E has been solved to obtain scattering information in atom-surface and molecule-surface scattering [2]. One of the merits of the method is that its computational effort scales slowly with the number of coupled diffraction channels, which can easily be many hundreds of states, and standard methods for time-independent coupled channel (CC) [1,2] calculations become formidable for such cases. Also, results for a range of energies are easily obtained with little additional effort in T D S E approaches. For a perfect surface, the structure of the surface is periodic. The system is described by a 3D Cartesian coordinate system. The periodic boundary conditions in x and y are highly suitable for the F F T method. In the D A F approach, however, the periodic condition can also be easily introduced by the following treatment.

Periodic DAF As it is shown in eq. (14), the D A F matrix can be viewed as a vector, namely,

Fj(kM)={Fo(M)(j--k["r) forlj-kl<_w,

(17)

otherwise, where the vector, F(M)(jlT), is given in eq. (13) To impose the periodic boundary condition on the D A F matrix, the so called "wrap-around" [27] technique is employed to construct a periodic DAF vector, PF ~M), namely p F ( M ) ( j I ~ ) =F(M)(jI~),

for j = - w , . . . , w,

pF(M)(N+jI~) =F(M)(jI~),

for j = - 1 , . . . , - w ,

p F ( M ) ( - N + j l z ) =F(M)(jlr),

for j = 1 . . . . . w,

(18)

pF(M)(j[z) = O,

otherwise,

(19)

where N is the size of the grid and 2w + 1 is the bandwidth. The resultant periodic DAF matrix has the structure given by

PF~(~M)=pF(M)(j--kIr),

for j = 0 . . . . . N - l ,

k=0 ..... N-l,

(20)

and is a circulant matrix. The circulant matrix-vector multiplication can be performed with computational effort scaling as ~e(N log N ) [25]. It is noted that, for an imperfect surface, F F T methods require consideration of the boundary effects whereas the original D A F method doesn't require any additional modification. The Cartesian coordinates are chosen such that the surface is parallel to the x-y plane. An incident atom approaches normal to the surface and parallel to the z axis. The wavepacket describing the atom, after a propagation in time, is given by the form

o(M)( I, j, k) = E e-iV"'J'~)¢/zhPF~tYx)( ~) PFj~¥')(z) l~i'k'

XFk(~,D('r) e -iV(l',j',k')¢ /2h ~b(/',

j', k').

(21)

The forms of F~Mz)('r), PF ~Mx) and PF fM~) are given in eqs. (14) and (20), respectively. The scattering information is extracted out from the time-dependent wavepacket by the following final state analysis [22].

Y. Huanget al. / DAF approachto time-dependentwavepacketpropagationin 3D

7

2.3. Analysis of scattering in the atom-surface scattering problems Atomic units are used in the analysis, i.e., h = 1, e = 1, m e = 1. The analysis is carried out by two transforms, i.e. coordinate to momentum and time to energy [22]. The Hamiltonian of a typical atom-surface scattering has the form 1 IQ = - - - z a + V ( x , y, z). 2/z

(22)

For a perfect surface, the potential V(x, y, z) is periodic in x and y. The time-dependent, fixed energy Schr/Sdinger equation is H 6 = i-~p.

(23)

The incident atom approaches normal to the surface and the initial state is 1 exp( - ikav z ) exp ( O(00kavf xyz 10) = al ~----~-

(z )- 2zav)2 ~2 ,

(24)

where a is the lattice constant. Here, k~v is the average momentum along z, Zav is the average z-position of the projectile, and ,~ is the Gaussian width-parameter. The initial wavepacket is constructed to have negligible overlap with the interaction potential V(x, y, z). A coordinate-to-momentum transform then gives not only the range of energies of the initial system but also the amplitude of the corresponding fixed energy states. For the initial wavepacket eq. (24), the transform can be carried out analytically, namely,

O(OOkavlXyzlO)=f~_odklexp(-ikz ) A ( k ) ,

(25)

where the square lattice with the lattice constant a is assumed and the amplitude of the momentum component, A(k), is obtained by a Fourier transform and has the form,

A ( k ) = lx/~-(rr)-3/4 e x p ( - ½~2(k - k a v ) 2) exp(i(k - kay) Zav),

(26)

From the basic theory of wavepackets and time dependent scattering [28], we also know that after propagating for time t, the wavepacket is given by ~(00kavl xyz It) = L ~ d k A( k ) ~ + (00k I xyz ) e x p ( - i(k2/2tz)t ),

(27)

where O+(00k I xyz) is the solution of Lippmann-Schwinger equation. In the free propagation region, the solution of Lippmann-Schwinger equation has the form 1

O+(00klxyz) = - e x p ( - i k z z ) + a

E

Se(OOImn)Ymn(X, Y) eik'nz,

(28)

,

(29)

m,n=O

where the function Ymn(X, y) is

gmn(X, y )

1 =

exp i

a

x+

a

y

8

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

and the total energy of the system is 2jz

2U + 2/~

(m2 +

nz)"

(30)

The system is highly degenerate and yields hundreds of open channels easily. Care must be taken in using eqs. (25) and (27). The same amplitude A(k) occurs in both equations provided that the initial wavepacket has no overlap with the interacting potential V(x, y, z) in t < 0 and A(- I kl) is essentially zero. The final state analysis is performed by taking the time to energy transform of eq. (27), namely, 1

2~rfodt exp(iE/) 0(0, 0, ka~lx, Y, zlt) -2~fo

=

dtfo

dk~ A(k~) O+(O'O'kzlx' Y'z)exp -i - ~ - E

A(kz) ,+(0, 0, kzlx, y,

t

-e)

(31)

In going from the integral over dk z to dr/, we have essentially changed from momentum to kinetic energy as the integration variable. The wavepacket in the energy domain is then given by

q,e(O,O, kavlX, y,z)= 1-~f ~dt 2"rr 0

exp(iE/)q,(0,0,

kavlX, y, z[t)

= -~A(kz) q~+(O,O, k z Ix, y, z),

(32)

where the change of the variable k z to ~ = k2/2~ is made in the integration. Projecting the state ~0E(0, 0, kay I x, y, z) onto the orthonormal function Ym,,(x, Y)' we are able to obtain a simple formula for the state-to-state transition amplitude information, i.e.

fodX fodY qJe(O,O, k,v Ix, y, zs)Ymn(X, y) =

~zA(kz)(e-ikzZ~mo~no + SE(O, 0 l m , n)eik""zs),

(33)

where

fofoY, n,,(x, y) Ym,,,,(x, y) dx dy = (~mm,~nn,

(34)

is used in obtaining eq. (33). Finally, the transition probabilities (the absolute square of the S-matrix elements) are given by

kmn ISe(0, Poo;mn(E) = ---~-0

0 l m , n)

12.

In the following section, we use the methods given in this section to obtain the numerical results.

(35)

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

9

3. Numerical study 3.1. Algorithm In the coordinate space, the potential energy operator is local and the D A F free propagation operator is near-local as shown in eqs. (11) and (14). The most costly step involved in calculating the time propagation of the wavefunction is the matrix-vector multiplication. To explore the banded structure of the D A F matrix, eq (14), and to take advantage of the fact that the DAF matrix-vector multiplication is a discrete convolution eq. (12), we have adopted two algorithms designed for such special structures.

Banded matrix-vector multiplication Since the D A F matrix in eq. (14) is a banded symmetric Toeplitz matrix in the uniform grid representation, it can be represented by a vector of length equal to the one half of the bandwidth of the D A F matrix, namely b ( j ) = 6(M)(j),

for j = 0 , . . . , w .

(36)

The algorithm for i = 1 . . . . . N y(i) = b(O) x x(i) end for j = 1. . . . . w for i = j , . . . , N y(i) = b ( j ) x x ( i - j + 1) + y ( i ) end fori=l,...,N-j y(i) = b(j) X x(i + j) + y ( i ) end end maximizes the length of the inner loop and keeps the vector pipeline as full as possible. The algorithm is easily generalized to cases with dimensions greater than one.

DAF / Fast convolution algorithm The D A F matrix-vector multiplication involved in the calculation of propagated wave-function in eq. (12) in the uniform grid representation is equivalent to a discrete convolution calculation. In this respect, the D A F function corresponds to a response function while the wavefunction to be propagated can be viewed as a signal. The well-known overlap-save algorithm [5] is used to partition the wavefunction into segments. The minimum length of each segment is the bandwidth, 2w + 1, of the DAF matrix. Each segment is padded at the ends by the values of the neighboring segments. The length of each pad is one half of the bandwidth. The convolution on each segment is computed by an F F T algorithm. The computational effort thus scales as @(Nlog2ms), where m s >_ 4w + 1. The output in the unpadded region of each segment is the desired result. The procedure may be detailed according to the following: 1. Perform an F F T for the D A F matrix (14) in a "wrap-around" format once (an "overhead step"). 2. Partition the wavefunction into segments with length _> 2nw + 1. 3. Pad each segment on both ends with values from neighboring segments with lengths > nw (and with zeros at the ends of the grid). 4. Perform an F F T of each padded segments ( > 4nw + 1).

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

10

5. Multiply the values of each padded segment with the convolved values of the "wrap-around" DAF vector. 6. Perform another F F T on each padded segment. 7. Extract out the output as the unpadded region of each transformed, padded segment. The result is the propagated wavefunction. This procedure not only takes advantage of the bandedness of the DAF matrix but also uses the very powerful fast convolution algorithm to do the matrix multiplication involved in the time-dependent wavefunction propagation. It is noted that the D A F / Convolution method is different from the standard F F T treatment of applying the free propagator both with regard to the floating point operations and storage requirements. For the D A F / C o n v o l u t i o n method, the memory required is that of one segment. As we demonstrate below, the D A F / C o n v o l u t i o n outperforms the widely used F F T methods both with respect to CPU time (which scales as •(Nlog2m s) vs. O(Nlog2N)) and memory. We discuss the memory issue further in the subsection reporting numerical results.

3.2. Numerical results The examples considered are He and Ne scattered off LiF(001) and W(ll0). Atomic units are employed in the calculation, i.e., h = 1, e = 1 and m e = 1. The results of the CC calculation reported by A.T. Yinnon et al. [29] are used as the accepted values of various transitions. The potential of the system has the form given by

V(x, y, z)

e x p [ a ( z o -z)]{exp[a(z o - z ) ] - 2.0} - 2/30 e e x p [ 2 a ( z 0 - z ) ] [cos(2vrx/a) + cos(2~ry/a)],

=D e

(37)

which somewhat oversimplifies the problem (in that, e.g. no lattice vibration is included) but provides an ideal test for T D S E methods. In particular, it is ideally suited to the F F T method due to the periodicity of the potential in x and y. Thus, if the D A F / c o n v o l u t i o n method is shown to be more efficient than standard FFTs for this problem, it is essentially certain that this will be true also for any other scattering problem. To reduce the grids, an absorbing potential is introduced in the z coordinate in a form of a linear ramp given by Z --Z a

Vi(x, y, z ) = v0i

•

(38)

Zmax I Za

It is required that there be negligible overlap between the physical potential in eq. (37) and the absorbing potential eq. (38), namely,

V(x, y, z)Vi(x, y, z) = O.

(39)

Table 1 Parameters of the incident atoms, the potential energy and the absorbing potential (abs. pot.). Parameters

He/LiF(001)

Ne/W(ll0)

Eto t

2.2842 × 10-3 5.792 5.37 0.582 0.1 2.8 × 10-4 5.0 5.0 x 10-3

2.3245 x 10 -3 13.066 5.18 0.582 0.01 6.3 X 10-4 4.0 5.0 × 10-3

kz a fl De z0 abs. pot. Voi

11

Y. H u a n g et al. / D A F approach to time-dependent w a v e p a c k e t propagation in 3 D

Table 2 P a r a m e t e r s of the grid, t h e w a v e p a c k e t , the D A F , the n u m b e r of t i m e steps n t for N e / W ( 1 1 0 ) with t i m e step ~" = 800. Wavepacket I

W a v e p a c k e t II

= 0.4 Za~ = 13.0 ka~ = 13.066

~ = 2.0 Zav = 28.0 k a v = 13.066

Grid

DAF1

DAF2

nt

Grid

DAF3

nt

Xmin = 0 . 0 Xmax = a N x = 32 Ymin = 0 . 0 Ymax = a ivy = 32 Zmin = 0.0 Zmax = 32.0 N z = 256 z a = 17.0 ~) Z s = 15.0 b)

O-x(0) = 1.65Ax M x / 2 = 12 w x = 15 try(0) = 1.65Ay M y / 2 = 12 nx r = 15 trz(0) = 2.75za z M z / 2 = 36 w z = 20

o'x(0) = 1.65za x M x / 2 = 12 w x = 15 % ( 0 ) = 1.65zay M y / 2 = 12 w r = 15 ~rz(0) = 2.6A z M z / 2 = 32 w z = 20

500

Xmin = 0.0

cry(0) = 1.65A~ M x / 2 = 12 w x = 15 try(0) = 1.65A r M y / 2 = 12 wy = 15 o'z(0) = 2.75A z M z / 2 = 37 w z = 20

1200

Xmax = a N x = 32

Ymin = 0.0 Ym~x = a Ny = 32 Zmin = -- 5.0 Zmax = 55.0 N~ = 512 z~ = 40.0 a) Z s = 39.0 b)

a) Za is the p o s i t i o n w h e r e the a b s o r b i n g p o t e n t i a l starts. It e n d s at Zmaxb) Zs is the p o s i t i o n w h e r e the analysis of s c a t t e r i n g is done.

The

parameters

absorbing

corresponding For

the

CRAY-YMP two the

z 0,/3

parameters FFT

and

are provided

given

herein

limit determined

are Some

by the

1. F o r

grid and

an

the

optimized,

A single

for all the cases.

a together

in table

for the

calculation,

is employed.

methods

segments

D e, a ,

potential

CPU

used.

the energies

case, two bandwidths

DAF

are listed

assembler

is used

For

with each

the

throughout

2 and

algorithm we

given

atoms,

and

the

3SCI

For the DAF

W a v e p a c k e t II

,~ = 0.8 Z~v = 16.0 kay = 5.792

~ = 2.0 z~v = 28.0 k a y = 5.792

Grid

the

grid

the number in section

library

into

two

only

of partitions 3.1, we pad

nt

Grid

DAF2

nt

Xmin = 0.0

O'x(0) = 1 . 6 5 A x

600

Xmin = 0.0

g x/2

trx(0) = 1.65A x M s / 2 = 12 wx = 15 try(0) = 1.65Ay M r / 2 = 12 wy = 15 o-z(0) = 2.5Az M ~ / 2 = 20 wz = 25

1200

Xma x

a

N x = 32 Ymin = 0 . 0

Ymax = a Ny = 32 Zmin = 0.0

Zma~ = 37.0 N z = 256 za = 22.0 z s = 20.0

= 12

wx = 15 O'y(0) = 1.65Ay M y / 2 = 12 wy = 15 O'z(0) = 2.75A z M z / 2 = 36 w z = 20

X m ax = a

N x = 32 Ymin = 0 . 0 Ymax = a

Ny=32 Zmin = --5.0 Zmax = 55.0

N z = 512

z a = 40.0 z s = 39.0

on the

DAF1

=

The

approach,

Table 3 P a r a m e t e r s of the grid, the w a v e p a c k e t , the D A F , the n u m b e r of t i m e steps n t for H e / L i F ( 0 0 1 ) w i t h t i m e step ~" = 375. Wavepacket I

the

for calculations.

3. from

partition

by increasing

the procedure

incident

are used

the calculations.

is possible

Following

in tables

FFT

DAF/Convolution,

slight improvement

DAF-bandwidth.

3D

of the normal

to

each

12

Y. Huang et al. / D A F approach to time-dependent wavepacket propagation in 3D

Table 4 Transition probabilities for N e / W (110) at Eto t = 2.3245 X 10 -3, ]~ = 0.4. Final state mn

DAF1

DAF1/Convolution

DAF2/Convolution

FFT

CC

00 10 20 30 11 21 31 22 32

0.4046 (0.7%) 0.1106 (0.6%) 0.0061 (0.0%) 0.00014 (0.0%) 0.0298 (0.3%) 0.0016 (5.8%) 3.5 X 10-5 (2.8%) 0.0001 (0.0%) 1.8 )< 10- 6

0.4008 (0.2%) 0.1098 (0.1%) 0.0061 (0.0%) 0.00014 (0.0%) 0.0296 (1.0%) 0.0016 (5.8%) 3.5 X 10-5 (2.8%) 0.0001 (0.0%) 1.8 x 10 - 6

0.4098 (2.0%) 0.1107 (0.7%) 0.0061 (0.0%) 0.00013 (7.0%) 0.0297 (0.7%) 0.0016 (5.8%) 3.3 X 10-5 (8.3%) 0.0001 (0.0%) 1.7 x 10- 6

0.4046 (0.7%) 0.1106 (0.6%) 0.0061 (0.0%) 0.00014 (0.0%) 0.0297 (0.7%) 0.0016 (5.8%) 3.5)<10-5(2.8%) 0.0001 (0.0%) 1.8 )< 10 - 6

0.4016 0.1099 0.0061 0.00014 0.0299 0.0017 3.6X 10 -5 0.0001 _

APE a):

1.3%

1.2%

3.0%

1.3%

CPU (s):

2213.5

1184.5

1184.5

1216.0

a) APE is the average percentage error

segment by half of the bandwidth of the banded DAF matrix. The 3D FFT is applied to each segment. For the direct-DAF, the algorithm in section 3.1 is used for the calculation of the z dimension of the wavefunction. For x and y dimensions, it is still possible to keep the innermost loop as the "z-loop" such that vector pipelining is used most efficiently. The results for the two initial wackpackets for N e / W ( l l 0 ) are shown in tables 4 and 5, obtained by both the FFT and DAF methods. The CC results are listed as the accepted values. The percentage errors are provided in parentheses. The same thing was done for He/LiF(001) and the results are listed in tables 6 and 7. As mentioned previously, to have a stable propagation, there is a relation between o-(0) and M which must be satisfied for a given grid. Such a pair of parameters are referred to as a stable parameter-pair. However, the relation between the choice of the stable parameter-pair and the degree of the accuracy of the fitting equation (12) for a given grid and time step r is empirical at this stage. Empirically, the degree of the accuracy in a practical calculation is determined by propagating the initial wavepacket forward one time step r first, and then propagating backward another time step. The average error, 2e, of the two step propagation is determined by comparing the propagated wavepacket to the initial wavepacket. The final error of the total propagation scales as n e [25], where n is the number propagation steps. Since a typical scattering problem involves a number of propagation steps of order of 103, the

Table 5 Transition probabilities for N e / W (110) at Eto t = 2.3245 x 10 -3, ~ = 2.0. Final state mn

DAF3/Convolution

FFT

CC

00 10 20 11 21 22

0.3954 (1.5%) 0.1090 (0.8%) 0.0061 (0.0%) 0.0297 (0.7%) 0.0016 (5.8%) 0.0001 (0.0%)

0.4080 (1.6%) 0.1104 (0.5%) 0.0061 (0.0%) 0.0298 (0.3%) 0.0016 (5.8%) 0.0001 (0.0%)

0.4016 O.1099 0.0061 0.0299 0.0017 0.0001

APE:

1.5%

1.4%

CPU (s):

4599.9

5929.9

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

13

Table 6 Transition probabilities for He/LiF(001) at Eto t = 2.2842 × 10-3, ~ = 0.8. Final state

DAF1

DAF1/Convolution

FFT

CC

00 10 20 30 11 21 31 22 32 33

0.02086 (0.3%) 0.00992 (0.8%) 0.003581 (3.1%) 0.01802 (3.0%) 0.01935 (1.0%) 0.00916 (0.5%) 0.01495 (0.7%) 0.06810 (3.4%) 0.01928 (4.8%) 0.00264

0.02093 (0.7%) 0.00993 (0.9%) 0.003588 (2.9%) 0.01806 (2.7%) 0.01935 (1.0%) 0.00918 (0.8%) 0.01495 (0.7%) 0.06828 (3.2%) 0.01920 (4.3%) 0.00264

0.02083 (0.2%) 0.00992 (0.8%) 0.003590 (2.9%) 0.01799 (3.1%) 0.01938 (0.9%) 0.00918 (0.8%) 0.01492 (0.8%) 0.06810 (3.4%) 0.01923 (4.5%) 0.00264

0.02079 0.00984 0.03696 0.01857 0.01955 0.00911 0.001505 0.07053 0.01840 0.00138 a)

APE:

2.0%

1.9%

1.9%

CPU (s):

2597.5

1425.5

1456.1

mn

a) We have done calculations for 00 to 33 transition in various parameters. The results of the DAF and the FFF methods are in agreement.

one-time-step-error, E, is required to be on the order of 10 -6 or less. The results of two choices of the stable parameter-pair, o-(0) and M, listed in table 4 illustrate the performance of the DAF's. Regarding the storage requirements for best performance in the calculations, 4 N x N y N z real words of storage must be provided as the workspace for the 3D FFT method. For the partitioned D A F / Convolution, the storage for the workspace is 4nln2n 3 real words, ni, (i = 1, 2, 3), is the size of the segment in the x , y , z dimensions, respectively. For the simple vector pipelining-DAF approach, 2nln2n 3 real words of storage are required for the work space. Here, nl, n 2 and n 3 are the bandwidths of the DAF matrices in x , y , z , respectively. Updating the value on the grid is done almost locally in the DAF method while it involves all the values on the grid in the FFT method. The purpose of using wider initial wavepackets in tables 5 and 7 is to show the advantage of the overlap-save technique in terms of CPU time for relatively large grids.

Table 7 Transition probabilities for He/LiF(001) at Eto t = 2.2842 × 10 -3, ~ = 2.0. Final state

DAF2

DAF2/Convolution

FFT

CC

00 10 20 30 11 21 31 22 32 33

0.2046 (1.6%) 0.01009 (2.5%) 0.03717 (0.6%) 0.01811 (2.5%) 0.01907 (2.5%) 0.00919 (0.9%) 0.01458 (3.1%) 0.06869 (2.6%) 0.01811 (1.6%) 0.00236

0.02089 (0.5%) 0.01005 (2.1%) 0.003716 (0.5%) 0.01799 (3.1%) 0.01912 (2.3%) 0.00919 (0.9%) 0.01457 (3.2%) 0.06865 (2.7%) 0.01891 (2.8%) 0.00234

0.02086 (0.3%) 0.01005 (2.1%) 0.003717 (0.6%) 0.01803 (2.9%) 0.01910 (2.3%) 0.00919 (0.9%) 0.01460 (3.0%) 0.06868 (2.6%) 0.01896 (3.0%) 0.00234

0.02079 0.00984 0.03696 0.01857 0.01955 0.00911 0.001505 0.07053 0.01840 0.00138

APE:

2.0%

2.0%

2.0%

CPU (s):

11708.8

4683.6

6165.5

mn

14

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

3.3. Discussion of the results The numerical results of the two examples of atom-surface scattering problems illustrate the properties of the DAF method in Cartesian coordinates and the power of the method in all the respects: accuracy, CPU time and storage requirements. The results are all compared to those obtained by the straight forward FFT method on the same grid, for a problem which is especially suited to the FFT method, namely the problem of an atom scattering off a perfectly periodic, corrugated surface. In table 4, we see that the accuracy of the DAF approach is increased for Ne/W(001) atom-surface scattering as the parameter-pair {o,(0), M} is tuned more finely from the stable value {2.6, 32} to {2.75, 36}. The average percentage error is smaller than for the FFT. The CPU time in the DAF/Convolution calculation is only a little bit faster than the FFT calculation since the ratio of the overlap to each segment (40/125) is not small enough to lead to a large advantage. Although the straight DAF calculation is the slowest of all, it is the most general algorithm which can be applied to any general banded matrix-vector multiplication. The non-uniform grid DAF approach [25] and the non-Cartesian DAF approach [23] both result in general banded matrices. In these cases, the FFT method and DAF/Convolution method are not applicable. Nevertheless, it is noted that the CPU time taken by the direct-DAF approach is only about twice as much as required by the most powerful 3D FFT method. Furthermore, the direct-DAF approach requires the smallest storage space of all the methods, while the FFT needs the most. Table 5 shows the calculations for wider initial wavepackets. Consequently, the grid is larger than the previous one. The DAF/Convolution with 2 segments outperforms the FFT method noticeably by the ratio 1.0/1.35, for the CPU time, which is better than the theoretical ratio logz28/log229= 1.0/1.15. The reason may be related to the machine hardware in processing a large memory calculation. All of the calculations are done by allocating the storage required for the best performance. Again, the direct-DAF calculation only doubles the CPU time compared to that taken by the FFT calculation, but it involves much lower storage requirements. The accuracy obtained by the DAF method with the chosen parameters is about the same as that of the FFT method. In tables 6 and 7 for H e / L i F , we again found that the DAF method yields accurate results with faster speed than the FFT method. The same parameter-pair {tr(0), M} is used in the example since this pair gave accurate results for k = 13.066 in Ne/W, which is higher than k = 5.79 in He/LiF. 4. Conclusion

The numerical study conducted in this paper demonstrates the power of the DAF method in solving quantum dynamical problems. In the DAF method, the "wrap-around" technique is applied to the DAF functions on the surface to treat the periodicity. The study shows that the DAF method outperforms the most powerful available FFT implementation both in CPU time and storage requirements, even though the examples chosen are especially suited to the FFT algorithm. The DAF approach, in some cases, also yields more accurate results than the FFT does. The degree of accuracy in the DAF is determined by the choice of the stable DAF parameter-pair, o-(0) and M. The basic, banded-DAF also performs well, requiring only about twice the CPU time of the DAF/Convolution and 3D FFT methods. However, unlike them, the simple banded DAF continues to be applicable for non-uniform grids and non-Cartesian coordinates. Of course, it also requires the lowest storage of the three methods. The next test for the DAF will be reactive scattering. Based on the successful application of the DAF method in solving 1D, 2D and 3D inelastic scattering problems with relatively low CPU requirements, low memory requirements, and the near-locality of the DAF, we expect that the DAF method will prove extremely powerful in solving general quantum dynamical problems, not only on vector supercomputers, but also on distributed memory massively parallel machines.

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

15

References [1] For the hyperspherical approach, see: (a) A. Kuppermann and P.B. Hipes, J. Chem. Phys. 84 (1986) 5962. (b) J. Linderberg, Int. J. Quantum Chem., Quantum Chem. Symp., 19 (1986) 467. (c) R.T. Pack and G.A. Parker, J. Chem. Phys. 85 (1987) 5252. (d) J.M. Launay and B. Lepetit Chem. Phys. Lett. 144 (1988) 346. For the non-variational amplitude density approach, see: (e) J.Z.H. Zhang, D.J. Kouri, K. Haug, D.W. Schwenke, Y. Shima and D.G. Truhlar, J. Chem. Phys. 88 (1988) 2492. (f) M. Baer and Y. Shima, Phys. Rev. A 35 (1987) 5252. For the variational amplitude density approach, see: (g) D.W. Schwenke, K. Haug, M. Zhao, D.G. Truhlar, Y. Sun, J.Z.H. Zhang and D.J. Kouri, J. Phys. Chem. 92 (1988) 3202. For the S-matrix Kohn variational method, see: (h) J.Z.H. Zhang, and W.H. Miller, J. Chem. Phys. 88 (1988) 449. For the arrangement decoupling absorbing potential method, see: (i) D. Neuhauser and M.J. Baer, J. Chem. Phys. 88 (1988) 2858. For the log-derivative Kohn variational method, see: (j) D.E. Manolopoulos and R.E. Wyatt, J. Chem. Phys. 91 (1989) 6096. For the "scattered-wave" variational approach, see: (k) Y. Sun, D.J. Kouri, D.G. Truhlar, and D.W. Schwenke, Phys. Rev. A 41 (1990) 4857. [2] See, e.g., I.C. Kulander, ed., Time Dependent Methods for Quantum Dynamics, Thematic Iss., Comput. Phys. Comm. 63 (1991), and references therein. [3] J.V. Lill, G.A. Parker and J.C. Light, Chem. Phys. Lett. 89 (1982) 483. J.C. Light, J.P. Hamilton, and J.V. Lill, Chem. Phys. Lett. 82 (1985) 1400. [4] J.W. Cooley and J.W. Tukey, Math. Comput. 19 (1965) 297. [5] H.J. Nussbaumer, Fast Fourier Transform and Convolution Algorithms (Springer, New York, 1982). [6] M.D. Feit and J.A. Fleck, J. Chem. Phys. 78 (1982) 301. D. Kosloff and R. Kosloff, J. Chem. 79 (1983) 1823. [7] See, e.g.: R.P. Feynman, Rev. Mod. Phys, 20 (1948) 367. For attempts to avoid the problems of these cancellations, see, e.g.: (a) J.D. Doll and D.L. Freeman, Science 234 (1986) 1356. (b) D. Thirumalai, E.J. Bruskin and B.J. Berne, J. Chem. Phys. 79 (1983) 5029, 5063; 81 (1984) 2512. (c) W.H. Miller, S.D. Schwartz and J.W. Tromp, J. Chem. Phys. 79 (1983) 4889. (d) K. Yamashita and W.H. Miller, J. Chem. Phys. 82 (1985) 5474. (e) R.D. Coalson, J. Chem. Phys. 83 (1985) 688. (f) N.E. Hendrikson and E.J. Heller, Chem. Phys. Lett. 148 (1988) 277. (g) N. Makri, Chem. Phys. Lett. 159 (1989) 489. (h) N. Makri and W.H. Miller, J. Chem. Phys. 90 (1989) 904. (i) A. Sethia, S. Sanyal and Y. Singh, J. Chem. Phys. 93 (1990) 7268. (j) T.L. Marchioro, J. Math. Phys. 30 (1990) 2935. (k) T.L. Marchioro and T.L. Beck, J. Chem. Phys. 96 (1992) 2966. [8] See, e.g.: (a) C.D. Polychronopoulos, Parallel Programming and Compilers (Kluwer, Dordrecht, 1988). (b) J. Dongarra, I. Duff, P. Gaffney and S. McKee, Vector and Parallel Computing (Wiley, New York, 1989). (c) G.F. Carey, Parallel Supercomputing: Methods, Algorithms and Applications (Wiley, New York, 1989). (d) R. Suaya and G. Birtwistle, VLSI and Parallel Computation (Morgan Kaufmann, San Mateo, CA, 1990). [9] D. Chasman, R.J. Silbey and M. Eisenberg, Theor. Chim. Acta 79 (1991) 175. [10] D.K. Hoffman, O.A. Sharafeddin, D.J. Kouri, M. Carter, N. Nayar and J. Gustafson, Theor. Chim. Acta 79 (1991) 297. [11] C.C. Martson, G.G. Balint-Kurti and R.N. Dixon, Theor. Chim. Acta 79 (1991) 313. [12] N. Nayar, D.K. Hoffman, D.J. Kouri, J.L. Gustafson, and G.M. Prabhu, Proc. Parallel Comput. Process. '91 (1991). [13] A. Lagana, Comput. Phys. Comm. 70 (1992) 223. [14] D.K. Hoffman, N. Nayar, O.A. Sharafeddin and D.J. Kouri, J. Phys. Chem. 95 (1991) 299. [/5] D.K. Hoffman and D.J. Kouri J. Phys. Chem. 96 (1992) 1179. [16] D.J. Kouri, X. Ma, W: Zhu, B.M. Pettitt and D.K. Hoffman, J. Phys. Chem. 96 (1992) 9622. [17] D.J. Kouri and D.K. Hoffman, J. Phys. Chem. 96 (1992) 9631. [18] D.K. Hoffman, M. Arnold and D.J. Kouri, J. Phys. Chem. 96 (1992) 9637. [19] D.K. Hoffman, M. Arnold and D.J. Kouri, J. Phys. Chem. 97 (1993) 1110.

16 [20] [21] [22] [23] [24] [25] [26]

Y. Huang et al. / DAF approach to time-dependent wavepacket propagation in 3D

N. Nayar, D.K. Hoffman, X. Ma. and D.J. Kouri, J. Phys. Chem., in press. D.K. Hoffman, M. Arnold, W. Zhu and D.J. Kouri, J. Chem. Phys., in press. D.J. Kouri, M. Arnold and D.K. Hoffman, Chem. Phys. Lett. 203 (1993) 166. D.K. Hoffman and D.J. Kouri, J. Phys. Chem., in press. M. Arnold, D.K. Hoffman and D.J. Kouri, to be published. M. Arnold, T.L. Marchioro, D.K. Hoffman, Y. Huang, and D.J. Kouri, to be published. M.D. Feit and J.A. Fleck, J. Chem. Phys. 80 (1983) 2578. P. De Vries, in: NATO ASI series Vol. B171 (1988) p. 113. A.D. Bandrauk, H. Shen, Chem. Phys. Lett. 176 (1991) 428. [27] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes (Cambridge Univ. Press, Cambridge 1986). [28] M.L. Goldberger, and K.M. Watson, Collision Theory (Krieger, Huntington, NY, 1975). [29] A.T. Yinnon, S. Bosanac, R.B. Gerber and J.N. Murrell, Chem. Phys. Lett. 58 (1978) 364.