J. Phys. Chem. 1992, 96,9622-9630
9622
ARTICLES Monte Carlo Evaluation of ReaCTime Feynman Path Integrals for Quantal Many-Body Dynamics: Distrtbuted Approximating Functions and Gaussian Sampling Donald J. Kouri,* Wei Zhu, Xin Ma, Department of Chemistry and Department of Physics, University of Houston, Houston, Texas 77204-5641
B. Montgomery Pettitt, Department of Chemistry, University of Houston, Houston, Texas 77204-5641
and David K.Hoffman Department of Chemistry and Ames (Received: September 4, 1991)
Iowa State University, Ames, Iowa 5001 1
In this paper we report the initial steps in the development of a Monte Carlo method for evaluation of real-time Feynman path integrals for many-particle dynamics. The approach leads to Gaussian importance sampling for real-time dynamics in which the sampling function is short ranged due to the Occurrence of Gaussian factors. These Gaussian factors result from the use of a generalization of our new discrete distributed approximating functions (DDAFs) to continuous distributed approximating functions (CDAFs) so as to replace the exact coordinate representationfreeparticle propagator by a “CDAF-class, free-particle propagator” which is highly banded. The envelope of the CDAF-class free propagator is the product of a “bare Gaussian”, exp[-(x’- ~)~$(0)/(2u~(O) + h2r2/m2)], with a “shape polynomial” in (x’- x ) ~where , a(0) is a width parameter at zero time (associated with the description of the wavepacket in terms of Hermite functions), T is the time step (7 = r / N , where r is the total propagation time), and x and x’are any two configurations of the system. The bare Gaussians are used for Monte Carlo integration of a path integral for the survival probability of a Gaussian wavepacket in a Morse potential. The approach appears promising for real-time quantum Monte Carlo studies based on the timedependent SchrMingerequation, the time-dependent von Neumann equation, and related equations.
I. Introduction In recent years there has been great interest in the use of time-dependent quantum mechanics to describe a variety of dynamical systems. A good indication of the range of methods and problems of interest is provided by the many articles in a recent volume of Computer Physics Communications devoted solely to this subject.’ Of particular interest have been methods which involve the evaluation of the action of the evolution operator F ( t ) = exp(-iHt/h) (1) where H is the usual (time independent) Hamiltonian and t is the time of evolution. The time evolution of both the wave function (in the Schrbdinger picture) and the density operator (for an isolated or closed system) can be obtained in principle if one has an algorithm for evaluating the effect of F(t).I Furthermore, by means of integral equation methods, more general systems can be treated rigorously,2-’ including timedependent interactions and open system^.^ The integral equation approach enables one to separate the actions of various noncommuting operators, say a reference H, and a disturbance Hd,comprising the Hamiltonian, H (e.g., the kinetic energy K and potential energy V, operators). Various possible choices for H, and Hdcan be made,3 leading to corresponding choices of F,(t) = exp(-iH,t/h) (2) Several alternatives for evaluating eq 1 have received attention. The Chebychev expansion method of Tal-Ezer and Koslofa has proved very powerful for a variety of problemsOgTwo popular alternative methods are the so-called “split operator” approach ‘Am- Laboratory is ope-rated for the U S . Department of Energy by Iowa State University under Contract No. 2-7405-ENG82.
and the more accurate “symmetric split operator” The Chebychev method is distinctive because it yields accurate results for an arbitrarily long time t if H is time independent over the entire interval, provided only that sufficient terms are included in the Chebychev expansion.* The split operator approach relies on the neglect of higher order commutators in expressing exp(-iHt/h) aslo exp(-iHt/h) = exp(-iKt/h) exp(-iVt/h) (3) where the error is of order t2. The symmetric split operator approach is given exp(-iHt/h) = exp(-iVt/2h) exp(-iKt/ h ) exp(-iVt/2h) (4)
or alternatively exp(-iHt / h ) exp(-iKt / 2 h) exp(-i Vr / h ) exp(-iKt / 2 h ) (5) where two different forms are possible because of the fact that one may interpret either K as the reference Hamiltonian (eq 4) or V as the reference Hamiltonian (eq 5).3 Clearly, the split operator expressions require t to be sufficiently small in order for them to be accurate. In most applications of the above methods, the action of the kinetic energy operator has been carried out in the momentum representation and of the potential energy has been carried out in the coordinate representation. Fast Fourier transforms (FFTs)have been used to transform between these two representations.I0J2 This essentially limits such approaches to systems in which no more than three degrees of freedom are treated by wave packet propagation. The integral equation method^^-^ make use of the fact.that the exact solution of the Schrbdinger equation, for any time interval
0022-3654/92/2096-9622%03.00/0 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 24, 1992 9623
Real Time Feynman Path Integrals and choice of H,(time independent) and Hd (which cun be time dependent), can be expressed as I+(t)) = (-i/h)jldt'exp[-iHr(t 1-7
- t?/h]HdI+(t?)
+
exp(-iHrT/h)l+(t - 7 ) ) (6) Then various approximations to the integral over t' yields different algorithms; e.g., the use of the trapezoidal rule and H,= K,Hd = V yields the kinetic referenced modified Cayley method4s5 (KRMC):
I+W)
= (1 + (i7/2h)V)-' exp(-iKT/h)(l - (i7/2h)V)I+(t - 7 ) ) (7)
It is obvious that eq 7 is identical to the kinetic referenced symmetric split operator method (KRSSO) up to T ~ However, . it is not necessarily true that the KRSSO is more accurate than the KRMC, even through the former contains terms to infinite order in i7V/2h and the evolution operator of eq 4, in contrast to eq 7,is unitary. This is because the coefficients of the higher order terms in the KRSSO are not correct. An obvious example is afforded by the case when the potential time dependent, V(t). Then the operator exp(-iHt/h) is no longer the correct evolution operator, while eq 6 remains valid and eq 7 becomes' I+(t)) = (1
+ (i7/2h)V(t))-l
exp(-iKr/h) (1 - ( i 7 / 2 h ) W
studies, the fact that 7 is either imaginary or complex leads to a Gaussian dependence on (xf - xg) which radically limits the values of xf, xg for which the propagator makes a contribution. Then the propagator for the complete time t can be formed from the short-time propagators byI3J8 exp(-iHt/h) = [exp(-iHr/h)lN where 7
= t/N
(17)
is short enough to make use of eq 4 or 7. The resolution Jdxjlxj) (xjl of the identity in the coordinate representation is inserted between the N - 1 pairs of factors in eq 16. This yields
where xNequals xf and the index k denotes the coordinate basis used at time kr. This, of course, is the discretized (in time) path integral form of the propagator. From eq 4 we obtain the well-known expression =
(xf(exp(-iHt/h)lxg),RSSO
lim (mN/ 27rih t) N G / 2
N--
X
(16)
mNN j dxl...j dxN- exp1-ihi c (Xk 22 I
k=I\,
- 7))IW - 7 ) ) (8)
Note that eq 8 is correctly time ordered. We now point out that methods very closely related to the split operator approach have been used in the Feynman path integral approach to many-particle d ~ n a m i c s . I ~In- ~Feynman ~ path integrals, one makes use of the fact that the short-time-limit expression for the full propagator F(T),in the coordinate representation, is given approximately by the so-called Trotter formula:Iob F(xf,xg17) = ( m / 2 r i h ~ ) exp[ ~ / ~I (z ( x f - x g ) 2 h 27
- V(xf)r)]
(9) or by the more accurate symmetrized (or trapezoidal) expression13
and by eq 7 we obtain the new expression (xdexP(-ifft/ h) 1%)
KRMC
=
The above equation is readily expressed as (xdexp(-iHt/ h 11x0) KRMC
(2 i j (10) where we assume the system has G-independentdegrees of freedom (so each vector x has G components). We note that (xdexp(-iKr/h)lxg) = ( m / 2 ~ i h r ) exp[ ~ / ~&(xf 2h7
- x0)2]
(11) Now considering the xf, xg matrix element of eq 4 for the KRSSO method, we have (xdexp(-iH7/h)lxo)uRSS0 E F(XR%17)KRSSO (12) = exp(-i V(xf)T / 2 h ) (xf(exp(-iKr / h )1%) exp(-i V(xg) T / 2 h ) (13) which, as discussed above, is seen to be identical with the symmetrized propagator expression, eq 10. By similar arguments, one readily sees that the short-time KRMC propagator in the coordinate representation is (xdexp(-iHT/ h)l%)KRMC E F(Xf&117)KRMC (14) (xdexp(-iK.r/h)lxg) Expressions such as eq 10 have been used extensively in studies of equilibrium statistical mechanics using purely imaginary time, as well as of correlation functions using complex time.1632In such
where
Axk)= h arcsin (4h7V(Xk)/[4h2 + ~ ~ V ( x k )=] ) h arccos {[4h2- r2V(xk)]/[4h2+ T ~ V ( X ~ (22) )]] and H(N - 1 - k) is the Heavyside function: H(y)=O, y < o
(23)
=1, y 1 0
(24)
The above form is very interesting, becauseflxk) can be viewed as an effective potential.26 The expressions eqs 19 and 21 are the path integral forms of the KRSSO and KRMC algorithms, with eq 19 being identical with a previously used expression.2' Path integral expressions like eqs 19 and 21 are only useful for purely imaginary, or in certain cases,complex time dynamics. This is because for real time, the factors containing exp[(imN/2ht)(xk - x ~ - , ) oscillate ~] more and more rapidly as N increases (or T is made smaller), as the distance Ixk - xk-,l increases, and as the "particle" mass m increases. Indeed, a basic outstanding problem is how one may utilize path integrals combined with Monte Carlo importance sampling methods when one wishes to treat real-time dynamic^.^^-^^ In the Monte Carlo method of integration, one desires to evaluate an integral q ( x ) by random sampling of the integrand. The uncertainty is determined by E/v'L, where Z is the variance (defined as (1/L)[C1v,11q(xi)12 - (l/L) X ( x C , q ( ~ ~ of ) )the ~ ]numerically ) evaluated integral and L is the
9624 The Journal of Physical Chemistry, Vol. 96, No. 24, 1992
number of randomly sampled points. Clearly, the variance Z is larger for rapidly oscillating functions, and an important strategy is to factor the integrand into a positive definite, integrable function, say d(x), and another factor h(x). It is then possible to employ d(x) as a probability distribution which, combined with an appropriate acceptance criterion, biases the selection of random points xi at which to accumulate the factor h(x,). If h(x) is a smootherfunction of x than q(x) in the important region selected by the bias, then the variance is accordingly reduced and a much more efficient evaluation of the integral achieved. It is important to keep in mind that for real-time dynamics, h(x) will, in general, be oscillatory (obviously, since d ( x ) is positive definite and d(x) h(x) equals the oscillatory functionflx)). Thus, the Monte Carlo procedure may be used for oscillating integrands provided there is a sufficiently rapid decay of the integrand in any regions where the oscillations become extreme. In the case of the exact propagator in Feynman path integrals, there is no obvious decaying factor, and the oscillations become increasingly violent as Ix' X I increases, so that it is impossible to utilize the Monte Carlo integration technique in any obvious fashion. The principle difficulty to be addressed is, therefore, the factorization of the integrand into a product of a positive, normalizable function to serve as the Monte Carlo importance sampling function times a function that varies slowly in the region of importance. There has, however, been some progress made. Early contributors include and Filin~v.~" They introduced factors in the integrand to bias the sampling toward the stationary phase points of the integrand. Thirumalai and Berne19b,29 introduced a procedure involving numerical analytic continuation in order to obtain real-time dynamics from complex time dynamics, evaluated using Metropolis type Monte Carlo integration. Subsequently, Makri and Miller2' extended the ideas of Doll and Filinov so that the stationary phase or semiclassical approximation would be the worst or most extreme limit. This leads to an exponential factor of the integrand which decays with a rate proportional to the square of the derivative of the action, but it necessitates the calculation of the classical action. A clear summary and discussion of the strengths and limitations of the approach are given by Makri.26bThe effective propagator technique developed by Makri in ref 26a (seealso ref 5 ) has similar objectives to those which we pursue in this paper but makes use of a sharp momentum cutoff, P,,,, to essentially coarse grain or smooth out the rapid oscillations of the short time propagator. We subsequently and independently considered this type approach in our first paper introducing the DDAF and chose not to pursue it because of the very slow decay of the resulting effective propagator.6a The sharp momentum cutoff leads to the result that the integrand is proportional to sin [Pmx(xr- x o ) / h ] / ( x-r xo),which is normalizable (due to interference cancellations) but not positive definite. The slowly decaying factor l / ( x f- xo) is a direct consequence of the use of a Fourier transform, plus the introduction of a sharp, discontinuous cutoff, P,,, in momentum space. (This is exactly analogous to the well-known behavior of the Fourier transform of a square waveform.33) Because sin [P,,,(xf x o ) / h ] / ( x -r xo) is not positive definite, it cannot be used as an importance sampling function in a Monte Carlo integration scheme. Makri instead multiplied and divided the integrand by the Gaussian ( X / P ) ' / ~ exp[-X(x'- x ) ~and ] fixed X so that the Gaussian falloff resembles the l / ( x ' - x ) behavior in the region important for the integration.26a Despite the introduction of the Gaussian sampling function, the rate of convergence of the Monte Carlo integration of Makri's expression is still governed by the behavior of sin [Pmax(x'-x ) / h ] / ( x ' - x ) due to the effective propagator since the integrand divided by the Gaussian remains oscillatory, and as (x'- x ) becomes large, the integrand becomes infinite in exactly the same way that the Gaussian tends to zero. For the P,, = 12 atomic units used, e.g., by Makri in calculations, the effective propagator still displays oscillatory behavior, even in the region where the quantity ( x ' - x ) increases. By contrast, the DDAF procedure first presented in ref 6a focused not on approximating the free propagator but rather on the accurate fitting or approximating of a specific class of wave
Kouri et al. packets, which is denoted as the "DAF class". The DAF class refers to L2 wave packets which are accurately fitted by a polynomial of degree M,under the DAF envelope; they have a smooth momentum distribution which decays rapidly at high energies. The DDAFs were constructed by requiring first that any function expressible as a polynomial of degree M could be exactly fitted. Second, it was required that the DAFs be exactly and analytically propagated by the free particle evolution operator. The only functions we are aware of which can be so propagated are plane waves, coherent states,32and Gaussians. The plane waves suffer from the disadvantage of being totally delocalized, so that the description of compact, L2 functions requires strong interference effects. Coherent state path integrals also have recently been shown by marc hi or^^^ to lead to a Gaussian weight factor in the integrand. However, the potential is not easily dealt with when it is not a polynomial. (We shall be carrying out a comparative study of the present CDAF approach and the coherent state approach in collaboration with Marchioro.) The Gaussians are convenient because they decay exponentially both in coordinate and momentum space, as desired. To efficiently satisfy the conditions that the DAFs exactly fit all polynomials of degree M and be exactly and analytically freely propagated, Hermite functions were used to construct them. The fact that Hermite functions can be exactly and analytically freely propagated results from the special property that the Nth Hermite function can be written in terms of the Nth derivative of the Gaussian generating function, and the fact that the free-particle evolution operator commutes with these derivatives (for Cartesian Having an accurate approximation to wave packets which are accurately approximated as Mth-degree polynomials under the envelope of the DAF at time t = 0, it is analytically and exactly propagated to time t = T by action of the free-particle evolution operator. Since this evaluation is done analytically and exactly, the only error introduced comes from the deviation of the wavepacket from a degree-M polynomial, and this accumulates slowly. By appropriate choice of the initial DAF fit, one can ensure that the new wave packet remains in the DAF class and can be similarly evolved forward freely. Of course, to carry out a full-time evolution, expressions such as eqs 4 or 7 are used, and it is required that the action of exp(-iVr/2h) or (1 =F i~V/2h)*ldoes not remove the packet from the DAF class. Finally, one may also extract a DDAF class free propagator from the analysis, as is discussed in ref 6a, this being simply the propagated DDAF for a single time step T. We close our introductory remarks by emphasizing that the payoff for an efficient Monte Carlo integration scheme for path integrals will be high, because the ability to carry out efficient evaluation of Feynman path integrals for real-time dynamics would make the treatment of an extremely wide range of problems possible. These include the time-dependent SchrMinger equation and the timedependent von Neumann equation, solutions for both of which can be expressed in terms of short-time propagators. In addition, recent efforts have been initiated to carry out detailed (nonperturbative) calculations of systems described by quantum ~hromodynamics.~~ These studies also make use of analogous path integral methods and likewise are feasible currently only for purely imaginary time dynamics. We also note that time-dependent correlation functions at finite temperatures can be expressed in terms of such propagators, except with complex time. Thus, the successful development of efficient Monte Carlo importance sampling methods to calculate real-time path integrals can be hoped to impact the entire range of phenomena described by quantum mechanics. We believe the present approach holds promise for progress in the pursuit of this goal. This paper is organized as follows. In section 11, we extend the discretized distributed approximating functions (DAF) method6for obtaining the DAF class coordinate representation free propagator to the continuous case. In section I11 we indicate how real-time Feynman path integrals can be evaluated by Monte Carlo methods with a Gaussian importance sampling function. Section IV contains a computational example treated using the present CDAF approach for the Monte Carlo evaluation of real-time
Real Time Feynman Path Integrals quantum dynamics. Finally, we present our conclusions in section V.
11. Continuous Distributed Approximating Function-Class Free Propagator We begin by examining a system with one degree of freedom and then generalize to one with G degrees of freedom. Recently we have considered the problem of accurately evaluating the action of the short time free propagator in the coordinate representation! The work was motivated, in part, by the desire to take advantage of the fact that for short times T , the kinetic energy operator is almost local, while the potential is exactly local regardless of time. However, as is well-known, the exact coordinate representation short-time free propagator becomes increasingly highly oscillatory as T is reducedi3and in fact has constant modulus and thus is of infinite extent as a function of distance Ix - x’l. This is a consequence of the necessity to include infinitely high momenta with equal weight, in evaluating the exact (analytical) coordinate representation matrix element of e x p ( - i K ~ / h ) . ~ qHowever, ~~ as stated earlier, the action of the free propagator on Hermite functions (Le., the Hermite polynomials times the Gaussian generator of these polynomials) can be obtained exactly and analytically.6 The natural spreading of wave packets under propagation is conveniently included (exactly) through a complex, time-dependent width, u(T), for the generator and the corresponding Hermite polynomial^.^^ Furthermore, as will be important later, the spread of the Gaussian envelope is the same for all of the Hermite functions. It is well known that one obtains the minimum packet spreading for a Gaussian of which the envelope of the Hermite functions is a special case. The above factors strongly suggest that one use the Hermite functions to develop an efficient procedure for evaluating the action of the free propagator on particular wave packets which can be accurately represented in terms of distributed approximating functions constructed using these Hermite functions. However, standard spectral collocation techniques lead to a description in which the truncated basis representation of the wave packet is constrained to be exactly equal to the wave packet on a finite number of grid points, the value of the wave packet between grid points being obtained by interpolati~n.~~ This leads to a slow decay of the correlations between values of the function at any pair of grid points.6 That is, in the case, e.g., of Fourier series, it, combined with the sharp momentum cutoff, leads to the sin [P-(x’ - x ) / h ] / ( x ’ - x ) b e h a ~ i o r . ~As, ~discussed ~ above, in the DAF approach previously developed: we have deliberately avoided such an interpolation. Instead, we begin by positing approximating or fitting functions (DAFs) a.(x) such that a class of well-behaved functions, A x ) , including f i packets, can be accurately approximated by
where These functions are required to fit exactly any polynomial of degree M. This weak set of conditions does not uniquely define aj(x). Instead, it permits great flexibility in defining the “range of influence” of the grid valuef(xj) on the approximate value of the function at neighboring points. In general, one expects that if M is increased (and hence the class of functions being approximated becomes more inclusive) the approximation of eq 25 would become more nearly an equality. The strategy next is to use the inherent flexibility of the fitting functions to express aj(x) in terms of the Hermite functions so that the response of these fitting functions to the free propagator may be calculated readily. The effect of the exact free propagator on the discretized wave packet can be obtained then by reading off the coefficient of the packet value at x,, after evolving the packet by time step T . This will simply be the value of the fitting function at the new time T , evaluated at (xf - xi). The condition that all polynomials of degree M be expressible exactly in terms
The Journal of Physical Chemistry, Vol. 96, No. 24, 1992 9625 of the fitting functions serves to fix the coefficients of a finite Hermite function expansion of the u,(x) involving Hermite functions up to degree M. These coefficients are, in fact, not constant but are instead periodic functions on the grid. For each x value (0 Ix < Ax, where Ax is the grid spacing) they satisfy M linear equations. Each element in the M X M kernel matrix is given as a discrete sum of products of two Hermite polynomials and the Gaussian weight function. In principle (a) these linear equations could be solved for each x value and (b) the finite sum of the periodic coefficient functions thus obtained, multiplied by the appropriate Hermite functions, could then be reexpanded to yield an infinite sum of Hermite functions with constant coefficients. Instead, we approximate the matrix elements of the kernel, which are given as discrete sums over all the grid points, by integrals over the whole line. Because of the orthogonality of Hermite polynomials under the Gaussian weight, this obviates the need to solve any matrix equations numerically. The result is an asymptotic expansion in M. The optimum value for M is roughly fitted by the width of the Gaussian generating function relative to the grid spacing, but within a broad tolerance range the choice of M is not critical. The uj(x) which we obtain by this procedure we call discrete distributed approximating functions (DDAFs). As developed in our earlier papers, we have previously required only the DDAF class free propagator, which is given by the freely propagated DDAF evaluated at points on the grid.6 That is F’j(T)free
lF(T)freeajKx=x~)=
(Ax/gx,(O))
X
MI 2
exp(-(x’ - x , ) ~ / ~ u , ~ (C T )(-1/4)n[ux,(0)/u,,(~)12n+i ) x n-0
(2*(n!)2)-1/ZHzn[2-1/2(~,, - xj) /u,(T)] (27) It is seen that only Hermite functions of even index appear in this expression. The result immediately generalizes for multidimensional systems to
where j is a collective index and j(q) indicates the jth grid point for the Cartesian degree of freedom, xp. In eq 27, u,(O) is the original width of the Gaussian-generating function for the Hermite polynomials of Cartesian coordinate xq, and uq(r) is given by u:(T)
= u:(O)
+ ihr/m,
(29)
Thus, in the multidimensional case, uJ0) depends on the degree of freedom, q, and the mass, mq, is that of the particle whose coordinates are labeled by q. However, for the most common derivation of a Monte Carlo evaluation of the path integral form of the free propagator, it is desired to have the complete, continuously indexed matrix F(x,xlT)free.(In actual computations, it is convenient and much more efficient to revert back to the discretized or quadrature form of the spatial integrations. One may then randomly sample the integration variables just as if the continuum form of the integral were being evaluated, while interpolating between grid points to evaluate the contribution from that point in the random walk to the integral. This enormously reduces the computational effort and makes possible the efficient sampling of extremely large numbers of points. This is discussed further in section IV.) To obtain F(x,xlT)frcewe first note from eqs 25 and 27 that m
where F j t ( ~ I ~ ) fisr eobtained e by replacing xj is eq 27 by the continuous variable x. Thus Fy(xIT)free = ( h / u , ( O ) ) exp(-U’Ax - XI’/ MI2 2U?(T)) (-1 / 4 ) n [ U , ( o ) / U , ( T ) ] 2 n ” ( 2 * ( n ! ) 2 ) - i / 2 n=O
x
H 2 n [ 2 - ’ / 2 ( ~-j tx ~ ) / u , ( T ) ](31)
Replacement of the j’index by a continuous variable x’is more
9626 The Journul of Physicol Chemistry, Vol. 96, No.24, 1992
subtle. Since the DAFs provide a general fitting procedure for the class of functions to whichflxp) belongs they can be used to approximate tbe fuactianflx 4). Howeva, this is, of copr# equivalent to fittingflqo) on a now grid shifted by a distance e fram the OM. If we give all wch grids for 0 I e < Ax equal weight, w h m Ax is the grid spa-, tben we can write
-
0
Axld = ( l / N P-c dc Fy&IT)rloJ(xy(,)lo)
Sldf’
F ( x s e l r ) r d 7 x l O ) (32)
w h m j‘(t) is the j t h grid point on the grid shifted by t/Ax from the original. It is readily verified that j’(c) is given by j’(t) = j ’ + c / A x (33)
Kouri et al. potential parameters? D = 0,164 au a = 0.9374 au Gaussian pack& U, = 0.248421 33 au ~p = 4.188971 259 au DAF parameters for comparison to Fourier effective propagatof with , P = 12 au o(0) = 0.3, M = 14, m = m“/2 = 918.5 au, 7 = 10.3352 au for amparim to Fourier effective propapt& with ,P = 48 au 4 0 ) = 0.1. M = 10, m = 4 2 = 918.5 au, T = 10.3352 au for Monte Carlo c o l ~ ~ l a t i o ~ o(0) = 0.2S, d( = 10. m = m ~ / = 2 918.5 au. T = 10.3352 au interpolation grid Ax = 4 X l(TJ au -2 I x I 2 a
V(x) = D[exp(-2ax)
- 2 exp(-cu)].
- X , ~ ) ~ / ~ U :‘See ] . ref 26a.
‘W(XJ0)
= [*a:]-’/‘
exp[+x
DAF(M= 14,SigmaO=O 3),Makri(P=12),Trotter
400
I
I
and then denoting j’(c)Ax by x’, we have x ‘ = j’Ax
+t
(34)
It follows that the continuous DAF class free propagator, F(X,X’~T)~is given by F(X,XlT)rw
(1 / u x ( o ) ) exp(-(x‘-
X ) 2 / ( 2 u x 2 ( T ) ) )x
1112
x (- 1 /4)”[ ax(0)/
u x ( T ) ]2n+l(2*(n!)2)4/2H2J
n 4
2-1 /2( x ’ -
X ) / Q x ( d I (35)
The crucial feature of F ( x , x ~ T is) ~the ~ presence of the ]. has a bare G a w i a n factor exp[-(x’ - X ) ~ / ~ U : ( T ) which Gaussian modulus exp[-(x’- x ) ~ u , ~ ( O ) / ( ~ ( U+~ (h2?/m2))]. O) This, tben,will lead to thc desind Gaussian importance sampling function in the Feynman path integral expressions. It is also of interest to examine the oscillatory part of the complex Gaussian factor, exp[(i(x’ - ~)~h~)/(2m(u,’(O) h 2 ? / d ) ) ] ,to sa whether it creates problems of high-frequency owillations for small or large particle mass m. It is obvious that tbjs quantity, in fact, beoomes lets ascillatory as T is made smaller or m is made larga (provided, of cou~sc,that ux(0)is independent of T and m). For notational purpoees, it is convenient to d e f m the ‘shape polynomial” g d x ‘ X ~ Tby )
+
-400 1 -2 00
-1 00
0 00
100
2 00
x (bohr)
Figm 1. Comparison of the M = 14, o(0) = 0.3, CDAF c b propagator (solid w e ) ; the , P = 12, F 0 U r i ~ r - Wcffeaiv~PrOpegPtW of ref 26a (daabed curve); and the exact, symmetrized Trotter full propagator (dotted curve). The potential is that given in ref 26a.
as M becomes inftnte, thc &function representation is obtained.
For the DDAF one cannot take this limit because the large M amtribu~cannotbea&quatdyrepresffltalcmthegrid. Henct, the expansion in M is asymptotically convergent. (Although by increasing ux(0), arbitrary accuracy can be obtained.) This 1112 difficulty disappears in the continuous limit. The approximate B d X ’ - 47) = (1 /Ux(O)) c (-1 /4)n[~x(0)/gx(T)lMI x no0 expression for the fiee propagator of eq 35 is thus simply the (2r(n!)2)-1/2Hk[2-112(X’ - X)/Ux( T ) ] (36) propagated truncated expansion of the b function (which can be evaluated analytically as we have discussed). The exact free so that propagator approaches b(x - x’), as T goes to zero, in a highly oscillatory manner. In contrast, the approximate expression of F(XPIT)fm exp[-(x’- x ) 2 / 2 u x ( r ) ] g v ( x ’ (37) eq 35 approaches the truncated b function smoothly in this limit. For G degrees of freedom, this leads to These points will be explicated more fully in a later communication. F(WIT)rroc exp[-W - W 2 / 2 1 g d x ‘ - 47) (38) Finally in this section, we present a direct comparison of the where CDAFtlass propagator with Makri’s Fourier-based effective propagator.26 Parameter values for the comparison are given in G G Table I. (The comparison is for the short-time symmetric Trotter (t’- w2 = c [(Xql/Qq(T)) - ( x 4 / ~ 4 ( 7 ) ) 1 2= [(a,‘) - (ji,)12 4’ I 4-1 expression for the full propagator, and this is constructed, on one (39) had, in tgmsof theCDAFclass free propagator and 011 the& in tamr ofthe Fauier-basad d€&h fm propal~ator.)La I, the CDAF parametera are taken to b e M = 14 and 4 0 ) = 0.3, which products accurate agreement between the CDAF propagator and Makri’s , P = 12, Fourier-based propagator for small (2*(n(q)!)2)-’/2Hk(q, [ ( x i X J / f i U & T ) l I (4) (x’- x ) values. It is dear that the CDAF propagator damps to zcfo much more rapidly, which rhould be advantageous whm OM The details of choosing M,~ ~ ( 0and ) . T for the discrete case are attempts the Moa& Carloevaluath ofa path integral. In F i i dibavlsed in our earlier papers6 We simply note here that these 2, we give a similar comparison, but for the case where Makri hrs increased ,P to 48. This products much better agreement parameten are easily c h a m so as to produce an effective free propagator that accurately evolves wavepackeu. betwm Makri’~F d a - b a s e d C ~ ~ & V C -tor and the exBd Before closing this discussion, we take a retrospective look at Trotter form short-time fall propagator. Again, the CDAF-pacq 35, which contains the main result of this section. The exact ramctm, M = 10 and 4 0 ) = 0.1, were c b n to produce clo6c agreement between the CDAF-clw propagator and that of free propagator, F ( X , X ~ ~(i.e., ) , ~the propagator when T = 0), is, of cwrae, rigorously equal to b(x - x’). But, since cq 35 is Makri,- in the small (x’- x) region. H m , the dampingof the orcillations in the CDAF-cleae propagator, compared to the anlyappmximatc the right-hand side for T =Oh not exactlyequal Fourier-basbd dT& p r o ~ e ~ t oisre, y ~ mott l proaaurccd. This to the b function. However, it is a ttuncored Hermite function (not plynowtiof) repredentation of the b function. In the limit behavior is the direct a x t a c q u a ~of~the ~ use of Hamite functions
-
x
-
Real Time Feynman Path Integrals DAF(M- 10,Sigma0=0 l),Makri(P=AB).Trotter
400 r
'
-400 - 2 00
1
- 1 00
0 00
100
I 2 00
x (bohr)
Figure 2. Comparison of the M = 10, o(0) = 0.1, CDAF class propagator (solid curve); the P, = 48, Fourier-based effective propagator of ref 26a (dashed curve); and the exact, symmetrized Trotter full propagator (dotted curve). The potential is that given in ref 26a. TABLE Ik Sarird Prob8bWka for Gai" Wave heket Ropptioa in 8 M o m Pohthl, Computed by M8lrix MultipliatiorP propagator Fourier: DAF, a(0) = 0.3, Fourier, DAF, o(0) = 0.1, t P-.. = 12 M = 14 P-.. = 48 M 5 10 0.25 0.50 0.75 1.00 1.25 1-50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75 5.00
0.9714 0.9000 0.8100 0.7180 0.6321 0.5554 0.4884 0.4307 0.3816 0.3404 0.3063 0.2784 0.2563 0.2393 0.2269 0.2186 0.2142 0.2134 0.2163 0.2227
0.9709 0.8988 0.8078 0.7148 0.6283 0.5515 0.4850 0.4282 0.3801 0.3397 0.3062 0.2789 0.2570 0.2400 0.2276 0.2192 0.2147 0.2139 0.2167 0.2234
0.9714 0.9000 0.8102 0.7188 0.6334 0.5569 0.4897 0.4316 0.3821 0.3405 0.3061 0.2782 0.2561 0.2391 0.2267 0.2184 0.2141 0.2133 0.2161 0.2225
0.9715 0.9000 0.8102 0.7188 0.6334 o . 5 ~ 0.4897 0.43 16 0.3821 0.3405 0.3061 0.2782 0.2561 0.2391 0.2267 0.2184 0.2141 0.2133 0.2162 0.2225
ref 19 for this treatment of real-time path integration. bSee ref 26. The present results were obtained using our codes.
as a means of constructing localized polynomial-type approximations (as opposed, e.g., to interpolation) to the class of wave packets which are to be propagated.& The survival probabilities for a Gaussian wave packet in a Morse potential (see Table I) obtained by matrix multiplication-propagation for t h m four casea are given in Table 11. We see that sufficiently accurate results are obtained in all cases.
To compute, e.g., an S matrix element using eq 44,one transforms from the 1%) representation initial state to the desired initial ( r = 0) wave packet I$i(0)) appropriate to the process of interest, and from the (xd dual-space representation final state to the desired final quantum state (4d of interest. Thus
(i/h:&xk))(l k= I
+ (iT/2h)V(xN))-'(l - (iT/2h)v(%)) (50)
which again is simply proportional to a specific state-to-state, energy resolved S matrix element.
9628
Kouri et al.
The Journal of Physical Chemistry, Vol. 96, No. 24, 1992
IV. Computational Example of Gaussian Sampling Monte Carlo To illustrate the DAF approach to the evaluation of real-time quantum path integrals by Monte Carlo with Gaussian importance sampling, we have carried out calculations for a simple 1-D model problem. We consider a Gaussian wave packet moving in a Morse potential, with the center of the packet initially located in the repulsive region of the potential.26a-3sValues of the system and computational parameters have been given in Table I. We calculate the autocorrelation function for the wave packet, at N times, t = N r , N = 1 , 2,..., 16: A(t) = I m -m d x\k(xlO) \k(xlt)
(51)
TABLE III: Continuous Random-Walk Constant Interpolation Method with Purelv Gaussian SamDLind Monte Carlo time 1 2 3 4 5 6 7 8 9
IO 11
total sampling points 400 000 900 000 1 600 000 2 500 000 3 600 000 4 900 000 6 400 000 18 000 000 20000000 22 Do0 000 24000000 39 000 000 42000000 45000000 56000000 72000000 540000000 475 000000
variance 1.4E-3 1.5E-3 1.7E-3 2.OE-3 2.58-3 3.38-3 4.48-3 2.7E-3 3.88-3 5.5E-3 8.7E-3 9.88-3 1.5E-2 2.2E-2 3.3E-2 5.OE-2 2.58-2 4.28-2
IA(t)I2
0.9708 0.9075 0.8088 0.7 192 0.6330 0.5543 0.4919 0.4266 0.3874 0.3441 0.3051 0.2809 0.2516 0.2362 0.2357 0.2275 0.2254 0.2266
quadrature IA(t)lZ
0.9714 0.9000 0.8102 0.7187 0.6332 0.5567 0.4896 0.4316 0.3821 0.3406 0.3062 0.2782 0.2561 0.2391 0.2267 0.2184 0.2140 0.2162
where the initial wave packet is real and is the ground-state 12 harmonic oscillator eigenfunction of the quadratic fit to the bottom 13 14 of the Morse potential well. (However, because the center of the 15 packet is initially displaced from the minimum of the Morse 16 potential, this initial packet has nonzero overlap with excited states 17 of the Morse potential.) 18 To provide results against which to compare, we have also carried out a straightforward trapezoidal quadrature e v a l ~ a t i o n ~ ~ The survival probability, IA(t)I2,is computed for a Gaussian wave packet moving in a Morse potential. See Table I for the system paof A ( t ) for the same times (with the time step r taken to be 10 rameters. fsz6a). Equation 51 can be written as dimensional space). We use the simplest possible interpolation for getting the value of the integrand between two grid points, namely, a constant value, which is taken equal to the value determined by the ( N l)-dimensional random point of interest. This is analogous to using “rectangles” to approximate the area under a curve in a 1-D integral. More sophisticated interpolations can, of course, be used, if desired. The important point is that the result becomes more and more accurate as the grid mesh used is made finer. The results obtained using the constant interpolation combined with a random walk in the continuous ( N + 1)-dimensional space are shown in Table 111. In this approach, A ( t ) is written as
+
with ug being the initial variance of the Gaussian wave packet to be propagated, and xinbeing the average position of this initial wave packet (equal to -0.1 au). Once having written eq 52 for the autocorrelation of the wave packet, it is necessary to choose an appropriate sampling function. The objective is to choose a sampling function, d(xo,...,xN) so that
A ( t ) = JmdxW..JmdXo d(xo,...,xN) k ( ~ 0..., , xN)
(55)
d(X0,...,XN) = w ( x ~..., , xN)/K
(56)
-m
-m
where w(x0,...,XN) =
d(xO,...,xN) is positive definite, and the ratio of the original integrand to d(xO,..,,xN)vanes smoothly and is large (small) where d(
[email protected],) is large (small). In fact, the choice of the bias sampling function is important since it can dramatically affect the number of sampling points required to converge the integration. In this paper we examine the simplest possible choice for the bias sampling function, but there may be others that are better. Another important practical point is that one should arrange the calculation so as to minimize the amount of work required to evaluate the integrand at the sampling points. S ince these points are only generated as one does the random walk (we use the Metropolis procedure), it might appear that one must calculate the integrand separately at each random sampling point. However, much of the difficulty can be circumvented in the following way, thereby leading to a rather efficient sampling. One evaluates the DAF, exp[-(x - ~ ’ ) ~ / 2 u ~ ( r ) ] g- ~XI?) ( x on a finite grid of difference points, the wave packet, P(xlO), and potential V(x) on a finite grid of x points. The finite grid of difference points is obtained by introducing a grid of M discrete values, x,, for the single variable, x, and forming all possible differences between any two of the grid points. Clearly, the difference grid is twice as large as the basic grid of M points. The integrand (or any part of it) at any ( N + 1)-dimensional grid point can then be quickly and efficiently constructed by a few multiplications and table ‘lookups”. This reduces the required computational effort without requiring a large amount of storage. To do the integral, we use the grid values to interpolate to the particular random point generated in the random walk point (in the continuous N 1
+
the (real) width parameter,,:a
is given by
= 1/Re( 1/ u Z ( r ) )
(58)
K is the integral K = JmdXo
...I -dXN W(X0,...,XN) -m
(59)
ci = (i/2) Im (l/uz(r))
(61)
-m
and The Monte Carlo expression for the path integral using a total of L random points is then L
A ( t ) = (K/L) h(xo,p,...Jjv,p)
(62)
P-1
where the points ( x ~,..., , xNP), ~ p = 1, ..., L, are determined by the
The Journal of Physical Chemistry, Vol. 96, No. 24, 1992 9629
Real Time Feynman Path Integrals Metropolis procedure. In general, the points (xo,, ...,xN,) do not lie on the grid, but d(xO,p,...,xN,e)(and h(xo,p,...,xN.p). if [d(xq, ...,xN )/d(xOSI, ...,xNS!)] satisfies the random selection critena for dDetermining a contnbuting point) is evaluated at the closest grid point. As the mesh is made finer, this converges to the exact result, provided enough sampling points are taken. The variance, uA,of A(t) is calculated in the usual way from
I ,
h
l
so that uA decreases only as l / d L . It should be noted that
reasonable results may be obtained with substantially fewer points if one is content with larger variance. In addition, the computational effort required to carry out literally hundreds of millions of samples is rather small; the implementation of the approach on a massively parallel supercomputer would enable one to carry out billions of samples easily.
V. Discussion In this paper we have presented the theory and an example application of a new, fully quantal procedure for evaluating the multidimensional path integral form of the real time full propagator by means of Monte Carlo integration with Gaussian importance sampling. Two alternative expressions for the propagator are developed based respectively on the kinetic referenced symmetric split operator and the kinetic referenced modified Cayley forms of the short time full propagator. Both approaches require the short-time free propagator in order to construct the full propagator for both short and long times. Indeed, any treatment of the short-time full propagator which expresses it in terms of the short-time free propagator can be considered by our approach, since the key step involves replacing the exact, highly oscillatory coordinate representation matrix elements of the free propagator by the new CDAF-class free propagator? The CDAF-class free propagator is characterized by the natural Occurrence of a Gaussian factor times a “shape polynomial”, both of which are functions of the generalized distance N G
d
c ( 3 k . q - Zk-i,,)’
k=lq=I
(64)
This leads to an integrand for the multidimensional path integral which can be put in a form to display an exponential damping factor similar to that which occurs in equilibrium statistical mechanics, i.e., Gaussian, and suggests the possibility of treating a range of real-time quantum dynamical processes in a similar way. This includes ordinary scattering processes, relaxation phenomena described by finite temperature quantum time correlation functions, processes described by the von Neumann equation for the density matrix, and nonequilibrium processes described by reduced quantum equations. The major feature (and in fact, the blessing and bane) of the Monte Carlo approach is, of course, the fact that the uncertainty of the result decreases according to (L)-I/’, where L is the number of sampling points. To the extent that the integrand remains smooth with increasing dimensionality, this scaling with L is independent of the dimensionality of the integr~l.~’The uncertainty, however, also depends directly on the variance, 2. Methods to reduce Z generally involve choosing an importance sampling function such that the ratio of the complete integrand to the sampling function is smoothly varying in the region of dominant contribution (Le., the region to which the sampling procedure is biased). The sampling function is required to be positive definite and normalizable. The CDAF class free propagator leads to a much smoother integrand compared to the original integrand containing the exact coordinate representation free propagators. This was illustrated by a comparison of the CDAF class propagator with the exact short-time full propagator, and with a Fourier-based effective propagator.268It is, however, important to note that the CDAF class free propagator does exhibit increasingly rapid oscillations as Ix’- XI becomes large, but the modulus is driven to zero sufficiently rapidly by the real
-10
-5
o
5
10
Figure 3. Comparison of the modulus of the CDAF-class free propagator, the bare Gaussian envelope factor, and the extended Gaussian, as a function of ( x - x’). The bare Gaussian is the narrowest. The CDAF modulus has the indentation at ( x - x ? = 0 and extends out the second farthest. The extended Gaussian just “kisses” the CDAF envelope and dies to zero just outside it. The extended Gaussian has a variance equal to 1.155.
Gaussian envelope. This gives reason to hope that the present approach will be applicable to longer time propagations, by virtue of a smaller variance for the Monte Carlo evaluation of the DAF expression for the path integral? Next, in Figure 3, we compare the modulus of the CDAF class free propagator with the Gaussian envelope factor exp(-m(x ~ ) ~ / 4 h It~ is) .obvious that the complete CDAF class propagator extends somewhat beyond the Gaussian envelope. The reason for this is the behavior of the shape polynomial, g&x’- x~T),which becomes large in magnitude as (x’- x) increases. Of course, the Gaussian factor ultimately dominates to drive the complete CDAF class propagator to zero, but the shape polynomial could give rise to numerical difficulties in the wings. To avoid this, one could choose the modulus of the CDAF class propagator to be the sampling bias function, since this leads to an “effective shape function”,gM(x’- x~T),which is everywhere bounded by 1. One also can numerically accomplish something similar by introducing a second, broader Gaussian function to obtain F(x,x~T)~~,, = exp[-(x’- x)’A] e x p [ i ( h ~ / h )x (x - ~ ? ’ / ( 4 0 ) + ~ (h7/m)’)l&!&’-
~ 1 7 )(65)
where XlT) = exp[-(x - x?~[u(O)’/(U(O)~ + (hT/m)’) - A]]gM(x’- X ~ T ) (66)
~M(X’-
+
and [ u ( O ) ~ / ( U ( O ) ~ (hT/m)’) - A] > 0. By proper choice of A we can obviously always choose g&x’- xlr) to be bounded by unity. It then follows that a biased sampling function based on the “extended Gaussian”, exp[-(x’- x)’A], will avoid contributions of large magnitude arising from the very large values of g&’ - X ~ Tat) the extremities of the somewhat narrower original Gaussian. We now contrast the numerical procedure just outlined to an apparently similar procedure followed by Makri,’“ who also multiplied and divided by a Gaussian. The purpose, however, was to introduce a sampling function which narrows the distribution. To do this, an attempt was made to find a Gaussian which simulates the envelope of sin [P,,,(x’- x)/h]/(x’- x), which of course eventually decays very slowly compared to any Gaussian. In this procedure the oscillations are carried by exp[(x’- x)’A] sin[P,(x’x)/h]/(x’- x) (cf. gM of eq 53), which grows rapidly and without bound for any choice of A as 1x’- XI becomes large. Thus, the ad hoc introduction of a Gaussian sampling function cannot alter the basic fact that Ix’- XI-’ decays very slowly. By contrast, the above procedure introduces a Gaussian distribution which widens the distribution so that lgMl < 1 as is convenient in the Monte Carlo sampling.
9630 The Journal of Physical Chemistry, Vol. 96, No. 24, 1992 The CDAF class free propagator contains the effect of wave packet spreading through the width parameter uq(7),which is the width associated with coordinatexq after a time step 7. The fact that the CDAFs are naturally expressed in terms of Hermite functions, expressible in terms of the Gaussian generating function6 suggests that the CDAF class free propagator contains the minimum possible effects of spreading. This, in turn, controls the rate for decay of the Gaussian importance sampling function. Finally, we make some observations regarding some practical considerationsof a Monte Carlo approach to real-time path integral quantum dynamics. First, because the path integrals are most straightforward when Cartesian coordinates are used, it is probably m a t natural, when dealing with systems containing more than two particles, to avoid separating out the Eulerian angles to take advantage of conservation of total angular momentum. Thus,one probably would calculate the plane-wave state resolved differential scattering amplitude directly and from it state resolved differential and integral cross sections. To appreciate the challenge posed by an atom-diatom collision system, one would have G = 6 (after separating out the center of mass degrees of freedom). It would not be unusual to require several hundred time steps to complete the collision process, so one might have N = 200. The result would be a 1206-fold multiple integral to be done by Monte Carlo with Gaussian importance sampling. To calibrate things, the real-time path integral calculation carried out in this work using Monte Carlc-Gaussian sampling techniques, involved N = 1 to N = 18, but for a one-dimedional system (so G = l), leading at most to a 19-fold multiple integral. Thus, the leap, even to atom-diatom collisions, is enormous and presents a great challenge. At the same time, the Monte C a r l d a u s s i a n importance sampling approach has several extremely attractive features compared to more conventional real-time wave packet propagation or time-independent methods. For example, there are minimal storage requirements since the DAF class free propagator is highly banded and the potential is diagonal in the coordinate representation. One also only need compute the DAF class free propagator as a function of a single difference variable in order to generate it for the entire multidimensional integrand. The use of a grid combined with interpolation to obtain the integrand at any random multidimensional point saves a huge amount of computation effort. One can compute a large number of transition amplitudes at once since they differ only in the final state and can be accumulated simultaneously. Each calculation at a random point is independent of all the others so the procedure is perfectly suited for massively parallel computing. Thus, it would appear that the approach is suited ideally for implementation on the new, massively parallel processing computers, such as the CM5. The calculations are done completely in the coordinate representation, and there are no complicated Yexchangenmatrix elements of the potential associated with reactive scattering. For these reasons, continued efforts to develop sufficiently powerful, new Monte Carlo methods for computing path integrals are well justified. We believe the present approach, based on the DAF ideas, holds promise for such calculations. We are continuing to pursue both formal improvements in the DAF approach, as well as exploration computations. Of course, the ultimate goal is to begin exploring fully converged quantum dynamics calculations for systems involving more than the traditional two or three atoms to which theorists have been confined. Acknowledgment. The authors thank Professor N. Makri for providing us with the necessary parameter values to make the comparisons shown in Figures 1 and 2. Thanks are also expressed
Kouri et al. to the referees for helpful comments. D.J.K. and X.M. were supported in part under National Science Foundation Grant CHE89-07429, and R. A. Welch Foundation Grant E-608. B.M.P. was supported in part under R. A. Welch Foundation Grant E-1028.
References and Notes ( I ) Time Dependent Methods for Quantum Dynamics; Kulander, K. C., Ed.; Comput. Phys. Commun. 1991, 63, 1-584. The articles in this special volume contain an excellent set of references to work on this subject in nonrelativistic quantum mechanics. (2) Hoffman, D. K.; Sharafeddin, 0. A.; Judson, R. S.; Kouri, D. J. J . Chem. Phys. 1990, 92,4167. (3) Sharafeddin, 0. A.; Judson, R. S.; Kouri, D. J.; Hoffman, D. K. J . Chem. Phvs. 1990.93. 5580. (4) Judson, R. S.; McGarrah, D. B.; Sharafeddin, 0. A,; Kouri, D. J.; Hoffman, D. K. J . Chem. Phys. 1991, 94, 3577. ( 5 ) Sharafeddin, 0.A.; Kouri, D. J.; Nayar, N.; Hoffman, D. K. J . Chem. Phys. 1991, 95, 3224. (6) (a) Hoffman, D. K.; Nayar, N.; Sharafeddin, 0. A,; Kouri, D. J. J . Phys. Chem. 1991, 95, 8299. (b) Hoffman, D. K.; Kouri, D. J. Ibid. 1992, 96, 1179. (7) Kouri, D. J.; Hoffman, D. K. Chem. Phys. Lett. 1991, 186, 91. (8) (a) Tal-Em, H.;Kosloff, R. J . Chem. Phys. 1984, 81, 3967. (b) Kosloff, R.; Tal-Ezer, H. Chem. Phys. Lett. 1986, 127, 223. (9) See, e.g., the review by: Kosloff, R. J. Phys. Chem. 1988, 92, 2087 and many papers in ref 1. (10) (a) Feit, M. D.; Fleck, J. A. J. Chem. Phys. 1982, 78, 301; 1983, 79, 301; 1984.80, 2578. (b) Trotter, M. F. Proc. Am. Math. Soc. 1959, IO, 545. (c) Nelson, E. J . Math. Phys. 1964, 5, 332. (11) DeVries, P. Nato ASI Ser. 1988, B171, 113. (12) Kosloff, D.; Kosloff, R. J . Comput. Phys. 1983, 52, 35; J . Chem. Phys. 1983, 79, 1823. (13) (a) Feynman, R. P. Reu. Mod. Phys. 1948, 20, 367. (b) Feynman,
R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (14) Pechukas, P. Phys. Reu. 1969, 181, 166, 174. (15) Miller, W. H. Adu. Chem. Phys. 1974, 25.69 and references therein. (16) (a) Doll, J. D.; Meyers, L. E. J. Chem. Phys. 1979, 71, 2880. (b) Doll, J. D. J. Chem. Phys. 1984,81, 3536. (17) Chandler, D.; Wolynes, P. G. J . Chem. Phys. 1981, 74, 4078. (18) Schulman, L. S. Techniques and Applications of Path Integration; Wiley: New York, 1981. (19) (a) Thirumalai, D.; Bruskin, E. J.; Berne, B. J. J . Chem. Phys. 1983, 79, 5063. (b) Thirumalai, D.; Berne, B. J. Ibid. 1983, 79, 5029. (c) Thirumalai, D.; Berne, B. Ibid. 1984, 81, 2512. (20) (a) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. J . Chem. Phys. 1983, 79,4889. (b) Yamashita, K.; Miller, W. H. J . Chem. Phys. 1985,82, 5474. (c) Tromp, J. W.; Miller, W. H. Faraday Discuss. Chem. SOC.1987, 84,441; J . Phys. Chem. 1986, 90, 3482. (21) Coalson, R. D. J. Chem. Phys. 1985,83,688. (22) (a) Doll, J. D.; Freeman, D. L. Science 1986, 234, 1356. (b) Doll, J. D.; Coalson, R. D.; Freeman, D. L.J . Chem. Phys. 1987, 87, 1641. (23) (a) Behrman, E. C.; Wolynes, P. G. J . Chem. Phys. 1985,83,5863. (b) Cline, R. E.; Wolynes, P. G. J . Chem. Phys. 1988,88,4334. (24) Berne, B. J.; Thirumalai, D. Annu. Reu. Phys. Chem. 1986, 37, 401. (25) Hendrikson, N. E.; Heller, E. J. Chem. Phys. Left. 1988, 148, 277. (26) (a) Makri, N. Chem. Phys. Lett. 1989, 159, 489. (b) Makri, N. Computer Phys. Commun. 1991, 63, 389. (27) Makri, N.; Miller, W. H. Chem. Phys. Lett. 1987, 139, IO; J . Chem. Phys. 1989, 90, 904. (28) Sethia, A.; Sanyal, S.; Singh, Y. J . Chem. Phys. 1990, 93, 7268. (29) Thirumalai, D.; Berne, B. J . Comput. Phys. Commun. 1991,63,415. (30) Filinov, V. S. Nucl. Phys. 1986, 8271, 717. (31) Baym, G.; Mermin, N. D. J. Math. Phys. 1961, 2, 232. (32) Marchioro, T. J . Math. Phys. 1990, 31, 2935. Marchioro, T. L.; Beck, T. L. J. Chem. Phys. 1992, 96, 2966. (33) Powell, J. L.; Craseman, B. Quantum Mechanics; Addison-Wesley: Reading, MA, 1961. 134) DeTar. C. Bull. Am. Phvs. SOC.1991, 36. 1778. (35) See, e.& the very nice discussion in Berman; M.; Kosloff, R. Comput. Phys. Commun. 1991,63, 1. (36) See, e&: Mowrey, R. C.; Kouri, D. J. J . Chem. Phys. 1986,84,6466. (37) See, e.g.: Koonin, S . E. Computational Physics; Benjamin Cummings: Menlo Park, CA, 1986; pp 190-191. (38) We use the same potential and wave packet parameters as those used in ref 26a.
