Analytic banded approximation for the discretized ... - ACS Publications

8299

J. Phys. Chem. 1991, 95, 8299-8305 disagreement between the vibrational levels of it and those of the Bowman potential and for remaining uncertainties in the assignment of the observed levels. We have suggested in the above analysis a corrected well depth of 11075 cm-' for the Bowman potential. This is to be compared with the ab initio well depth of 1109 cm-'. There is clearly reasonable agreement among these two values and the 1062-cm-' well depth of the Bowman potential. There is, however, a large discrepancy between the zero-point levels of the two surfaces that arises primarily from what appears to be an overly large contribution (-440 cm-')from the bending vibration on the ab initio surface. This leads not only to a systematic shift in the computed vdW stretching levels compared to experiment, but to the prediction of fewer excited bending states. The fact that experimental estimates of the dissociation energy, which range from 718 to 742 cm-', are so closely in accord with the present calculation (724 cm-l) and estimation (-739 cm-') on the Bowman surface strongly suggests that the bending po-

Downloaded by UNIV OF NEBRASKA-LINCOLN on September 1, 2015 | http://pubs.acs.org Publication Date: October 1, 1991 | doi: 10.1021/j100174a052

-

tential is better described by that surface. This inference is independent of the manner in which the non-vdW stretching levels are assigned. The well depth of the ab initio surface would have to be in error by enough to allow a 126-cm-' increase in the dissociation energy. While Chakravarty and Clarylo are inclined to accept a 15% error in the ab initio well depth, one must expect such an error also to affect the bending potential.

Acknowledgment. We are grateful for several instances of communication via electronic mail with J. M.Bowman, who also kindly sent his computer subroutine for evaluating the potential energy function proposed by him and his co-workers. Registry No. Ar, 7440-37-1; OH, 3352-57-6.

Supplementary Material Available: Contour plots of the wave functions of the seven states designated as delocalized in Table I11 (7 pages). Ordering information is given on any current masthead page.

Analytic Banded Approximation for the Discretized Free Propagator David K. Hoffman, Naresb Nayar, Department of Chemistry and Ames Laboratory,' Iowa State University, Ames, Iowa SO01 1

Omar A. Sharafeddin,* and D. J. Kouri**S Department of Chemistry and Department of Physics, University of Houston, Houston, Texas 77204-5641 (Received: March 27, 1991)

A procedure for obtaining an approximate sparse (banded) matrix for the coordination representation of the free-particle propagator is presented. The technique takes advantage of the fact that the action of the free propagator on Hermite polynomials can be obtained analytically and represents the propagator in terms of its effect on the Hermite basis.

I. Introduction The wave packet that represents a physical molecular scattering system of interest evolves according to the time-dependent SchrGdinger equation (TDSE), which is a linear, first order (in time), differential equation. In recent years, there has been considerable interest in the development and application of numerical methods for solving dynamics problems using time-dependent methods.'* In part this is because the TDSE constitutes an initial value problem. This enables one to devise solution methods that scale slowly with problem size. Typically, solutions involve developing a matrix representation of the problem, which leads to an Nz scaling, where N determines the dimension of the matrices. This reflects the principal computational step, which is forming the product of a square matrix with a column vector. Great interest has centered on using Fourier methods"J2 for representing the action of the Hamiltonian on the wave packet (especially the kinetic energy operator), so that with a Fourier grid of N points, the scaling associated with the matrix/vector product is N log2 N . Alternatively, finite difference methods applied to direct solutions of the TDSE lead to sparse, or banded, matrices for the Hamiltonian. Then if the Hamiltonian matrix has dimension N X N , with a band width of 2n + 1, applying this matrix to the column vector wave packet involves (2n + l ) N multiplications. (It should be noted that the value of the dimension N, as well as the size of the time step possible in solving the TDSE,

'Ames Laboratory is operated for the US.Departmentof Energy by Iowa

State University under Contract No. 2-7405-ENG82.

tR. A. Welch Foundation Postdoctoral Fellow under Grant E-608. #Supported in part under National Science Foundation Grant CHE8907429.

0022-3654/91/209S-8299$02.50/0

varies from method to method.) The coordinate and momentum representations have played a particularly important role in the ~

~

~~

(1) Some recent reviews include: (a) For molecule-surface scattering: Gerber, R. B.; Kosloff, R.; Eerman, M.Comput. Phys. Rep. 1986,5,59. (b) For collinear gas-phase reactive scattering: Mohan, V.; Sathymurthy, N. Ibid.

1988, 7,213. (c) A general review: Kosloff, R. J . Phys. Chem. 1988, 92, 2087. (2) Mazur, J.; Rubin, R. J. Chem. Phys. 1959, 31, 1395. (3) McCullough, E. A.; Wyatt, R. E. J. Chem. Phys. 1971,54,3578,3592. (4) Zubert, C.; Kamal, T.; Zulike, L. Chem. Phys. Loti. 1975,36, 396. (5) Heller, E. J. J. Chem. Phys. 1975,62,1544; 197665,4979. Kulandcr, K. C.; Heller, E. J. Ibid. 1978,69,2439. Drolshagen, G.; Heller, E. J. Ibid. 1983, 79, 2072. (6) Askar, A.; Cakmak, A. S.1. Chem. Phys. 1978,68, 2794. (7) Kulander, K. C. J . Chem. Phys. 1978,69,5064. Gray, J. C.; Frascr, G. A.; Truhlar, D. G.; Kulander, K. C. Ibfd. 1980, 73, 5726. Orel, A. E.; Kulander, K. C. Chem. Phys. Left. 1988, f46,428. (8) LeForestier, C.; Bagerson, G.; Hibcrty, P. C. Chem. Phys. Lett. 1981, 84,385 (see also: LeForestier, C. Chrm. Phys. 1974,87,241). (9) Agrawal, P. M.;Raff, L. M. J . Chem. Phys. 1982, 77, 3946. (IO) Heller, E. J.; Sundberg, R. L.; Tannor, D. F. J . Phys. Chem. 1982, 86, 1822. (11) Feit, M.D.; Fleck, J. A. J . Chem. Phys. 19g2,78,301;lW, 79,301, 1984,80, 2578. DcVrics, P. In NATOASI Series, Vol. B 1988, 171, 113. (12) Kosloff, D.; Kosloff, R. J. Compur. Phys. 1983, 52. 35; J . Chem. Phys. 1983, 79, 1823. (13) LeForestier, C. Chem. Phys. 1984,87,241. In Theory of Chemical Reaction Dynamics, Clary, D. C., Ed.; Reidel: Dordrecht, 1986; p 235. IC Quere, F.; LeForstier, C. J. Chrm. Phys. 1990, 92, 247. (14) Drolshagen, G.; Heller, E. J. J . Chem. Phys. 1983, 79, 2072. (IS) Tal-Ezer, H.; Kodoff, R. J. Chrm. Phys. 1984,8f, 3967. (16) Kosloff, R.; Cerjan, C. J . Chem. Phys. 1984,81, 3722. (17) Yinnon,A. T.; Kosloff, R.; Gerber, R. B. Sur/. Sei. 1984,148, 148. Gerbcr, R. B.; Yinnon, A. T.; Kosloff, R. Chem. Phys. Lett. 1984,105,523. (18) Imre, D.; Kinsey, J.; Sinha, A,; Krenos, J. J. Phys. Chem. 1984,88, 3956.

0 1991 American Chemical Society

8300 The Journal of Physical Chemistry, Vol. 95, No. 21, 1991


methods studied to date. For Cartesian variables, the kinetic energy is evaluated in the momentum representation, in which

(19) Skodje, R. T. Chem. Phys. Lett. 1984,109, 227. (20) Mowrey, R. C.; Kouri, D. J. Chem. Phys. Lett. 1985, 119, 285. (21) Jackson, B.; Metiu, H. J . Chem. Phys. 1985,83, 1952. Sawada, S.; Heather, R.;Jackson, B.; Metiu, H. Ibid. 1985, 83, 3009. (22) Zhang, Z. H.; Kouri, D. J. Phys. Rev. 1986, A34, 2687. (23) Jackson, B.; Metiu, H. J. Chem. Phys. 1986, 84, 3535; Ibid. 1986, 85, 4192; Ibid. 1987,86, 1026. (24) Mowrey, R. C.; Kouri, D. J. J . Chem. Phys. 1986.84.6466, (25) Park, T. J.; Light, J. C. J . Chem. Phys. 1986, 85, 5870. (26) Mohan, V.; Sathymurthy, N. Curr. Sci. 1986,55, 115. Thareja, S.; Sathymurthy, N. J . Phys. Cham. 1987, 91, 1970. (27) Mowrey, R. C.; Kouri, D. J. J . Chem. Phys. 1987,86.6140. (28) Heather, R.; Metiu, H. J . Chem. Phys. 1987, 86, 5009; 1989, 90, 6116; 1989, 91, 1596. (29) Kouri, D. J.; Mowrey, R. C. J . Chem. Phys. 1987,87,2087. (30) Sun, Y.; Mowrey, R. C.; Kouri, D. J. J . Chem. Phys. 1987,87,339. (31) Mowrey, R.C.; Bowen, H. F.; Kouri, D. J. J . Chem. Phys. 1987,86, 2441. Mowrey, R. C.; Bowen, H. F.; Yinnon, T.; Gerber, R. B. Ibid. 1988, 89, 3925; and manuscript in preparation. (32) Huber, D.; Heller, E. J. J . Chem. Phys. 1987, 87, 5302; 1988, 89, 4752. Huber, D.; Heller, E. J.; Littlejohn, R. G. J . Chem. Phys. 1987, 87, 5302. Huber, D.; Ling, S.;Imre, D. G.; Heller, E. J. Ibid. 1989, 90, 7317. (33) Sun,Y.; Kouri, D. J. J . Chem. Phys. 1988, 89, 2958. (34) Jackson, B. J . Chem. Phys. 1988,88, 1383; 1988,89, 2473; and in

press.

(35) 6667. (36) (37) (38) (39)

Chaseman, D.; Tannor, D. J.; Imre, D. G. J . Chem. Phys. 1988,89,

Williams, S.0.;Imre, D. J. J. Phys. Chem. 1988, 92, 6648. Tannor, D. J.; Rice, S.A. Adu. Chem. Phys. 1988, 70,441. Sun, Y.; Judson, R. S.;Kouri, D. J. J . Chem. Phys. 1989, 90, 241. Neuhauser, D.; Baer, M. J . Chem. Phys. 1989, 90, 4351; J . Phys. Chem. 1989, 93, 2872; J . Chem. Phys. 1989, 91,4651. (40) Neuhauser, D.; Baer, M.; Judson, R.S.;Kouri, D. J. J . Chem. Phys. 1989, 90,5882; 1990, 93, 312; Comput. Phys. Commun., in press. (41) Jiang, X.-P.; Heather, R.; Metiu, H. J . Chem. Phys. 1989,90, 2555. (42) J a m , M.;Atabek, 0.;LeForestier, C. J. Chem. Phys. 1989,91,1585. (43) Zhang, J.; Imre, D. G. J . Chem. Phys. 1989, 90,1666. (44) Dixon, R. N. Mol. Phys. 1989, 68, 263. (45) Judson, R. S.;Kouri, D. J.; Ncuhauscr, D.; Baer, M. Phys. Reo. 1990,

A42, 351. (46) Zhang, J. Z. H. Chem. Phys. Lett. 1989, 160,417; J. Chem. Phys. 1990, 92, 324. (47) Hoffman, D. K.; Sharafeddin, 0.;Judson, R.S.;Kouri, D. J. J. Chem. Phys. 1990, 92, 4167. (48) Das, S.; Tannor, D. J. J . Chem. Phys. 1990, 92, 3403. (49) Sharafeddin, 0.;Judson, R.; Kouri, D. J.; Hoffman, D. K. J . Chem. Phys. 1990, 93, 5580. (50) Judson, R. S.;McGarrah, D. B.; Sharafeddin, 0.; Kouri, D. J.; Hoffman, D. K. J . Chem. Phys. 1990,94, 5580. (51) Kouri, D. J.; Hoffman, D. K. Chem. Phys. Lett., submitted. (52) Neuhauser, D.; Judson, R. S.J. Chem. Phys., submitted. (53) Neuhauser, D.; Judson, R.S.;Baer, M.; Kouri, D. J. Chem. Phys. Left. 1990, 169, 372. (54) Sharafcddin, 0. A,; Bowen, H. F.; Kouri, D. J.; Hoffman, D. K. J .

Hoffman et al. it is diagonal. Then the action of the N X N matrix that transforms from the momentum to the coordinate representation is made efficient by using fast Fourier transforms (FFTs). From the point of view of formulating a solution method based on making the Hamiltonian matrix sparse, it it natural to think in terms of using a single representation throughout the calculation. The most obvious candidates are the coordinate representation and the momentum representation. In the coordinate representation, the change of the packet in time at a point in configuration space is governed by the potential energy of the system at that point, as well as the instantaneous kinetic energy of the packet at the configuration point. The first is a strictly local (diagonal) operator and the second is ‘nearly local” (almost diagonal) in the sense that it is defined as the limit as the neighborhood goes to zero of a second difference function defined on the neighborhood around the configuration point of interest. That is, the time evolution of the packet at a point in configuration space is governed by the value of the function at neighboring points. By contrast, in the momentum representation the kinetic energy is exactly local, but the potential energy is nonlocal. Furthermore, the latter is an integral operator and the extent of its nonlocal character is typically much greater than is the case for the coordinate representation of the kinetic energy. For this reason the coordinate representation is likely to be the more profitable representation to use for developing “banded” propagator matrix structures. There is a second compelling reason for seeking to develop sparse, or banded, matrix approaches to solving the TDSE. It is highly likely that most of the future major improvements in computational power will involve massively parallel systems. To take advantage of this, it should be possible in the coordinate representation, at least in principle, to assign a configuration space point to each processor and thus to form a grid on which the scattering calculation could be made to mimic the physical process (in that the calculation would proceed by each grid point “talking to” its close neighbors). Of course, at the present time the number of configuration grid points needed for quantitative work, for all but the simplest systems, far exceeds the number of processors available on any parallel machine, and one is forced to map many grid points onto a single processor. For certain parallel architectures, even if this were not required (say for some simple I-D systems), it would still be desirable because communication on a single processor is much faster than communication between processors. Thus, the guiding philosophy of making the numerical algorithm as “local” as possible in the coordinate representation so that the processor grid can be made to map closely onto the configuration grid seems to be a good one for making most efficient use of parallel machines. 11. Formal Development

Coordinate Space Calculations. Consider a 1 D scattering problem on the infinite line. It is well-known that an analytic expression exists for the free propagator, F(T), in the coordinate repre~entation.~~ Here

Comput. Phys., in press. (55) Viswanathan, R.;Shi, S.;Villalonga, E.; Rabitz, H. J . Chem. Phys.

1989, 91, 2333. (56) Neuhauser, D. J . Chem. Phys., in press. (57) Sharafeddin, 0. A.; Kouri. D. J.; Hoffman, D. K. Can. J . Chem.,

submitted. (58) Hoffman, D. K.; Sharafeddin, 0.A.; Kouri, D. J.; Carter, M.; Nayar, N.; Gustafson. J. Theor. Chim. Acta 1991. 79, 297. (59) Sharafeddin, 0. A.; Kouri, D. J.; Nayar, N.; Hoffman, D. K. J . Chem. Phys., in press. (60) Heller, E. J. J . Chem. Phys. 1991, 94, 2723. (61) Bandrauk, A. D.; Shen, H. Chem. Phys. Lett. 1991, 176, 428. (62) Wiliams, C. J.; Qian, J.; Tannor, D. J. J. Chem. Phys., submitted. (63) Dixon, R. N.; Marston, C. C.; Balint-Kurti, G . G. J . Chem. Phys. 1990, 93, 6520. (64) Founargiotakis, M.; Light, J. C. J. Chem. Phys. 1990, 93, 633. (65) Carter, M.; Nayar, N.; Gustafson, J.; Hoffman, D. K.; Sharafeddin, 0. A.; Kouri, D. J., to be published.

Although the kinetic energy operator itself has a “nearly local” character in the coordinate representation (which is easily seen on the basis of finite difference schemes), the free-particle propagator is not obviously nearly local. Indeed, for a long propagation time, it is far from local since it must translate a wave (66) Hoffman, D. K.; Ma, X.;Kouri, D. J. J . Chem. Phys., manuscript in preparation. (67) Truong, T. N.; Tanner, J. J.; Bala, P.; McCammon, J. A.; Kouri, D. J.; Lesyng, B.; Hoffman, D. K. J . Chem. Phys., submitted. (68) Parker, G. A.; Hoffman, D. K.; Kouri, D. J., manuscript in preparation. (69) Feynmann, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965.

The Journal of Physical Chemistry, Vol. 95, No. 21, 1991 8301

Discretized Free Propagator


packet in accordance with Ehrenfest’s theorem.70 That is d/dt ( x ) = ( p ) / m

(3)

d/dt ( p ) = -(dV/dx)

(4)

where V = V ( x ) is the potential. Thus, as time progresses, the wave packet must translate such that the expectation values of x and p change in accordance with the quantum average of Newton’s equations. However, for short enough time steps, one expects to recover “nearly local” behavior just on the basis of the Taylor series representation of F(T).Although for T 0, F(T) is proportional to the Dirac 6 function 6(x - x?, it is a well-known fact that for T # 0 (no matter how small) F(T)does not manifest the ‘nearly local” character leading to a narrow “bandwidth” (Le,, the exact matrix element F ( x , x ~ Tdoes ) not tend to zero rapidly as 1% - X I increases). Instead, F ( x , x ~ Tbecomes ) increasingly highly oscillatory as the magnitude of the coordinate difference increases and the time step decreases (but still is nonzero). The difficulties for real-time propagation caused by this behavior have been intensively documented and ~ t u d i e d . ~ ’ “ ~ Recently we have shown that the exact propagator can be “coarse grained” by eliminating its high-momentum components to obtain a banded operator.sg This was achieved by an approximate numerical quadrature evaluation of F ( x , x ~ Tso ) that (a) a discrete grid in k space was employed and (b) a finite cutoff k,,, was introduced. The results might have been thought to depend critically on the choice of the momentum cutoff; however, although the coarse-grained propagator operators created in this way bear very little resemblance to the exact propagator as a function of the coordinates, they can be used to propagate the wave packet without an unacceptable loss of accuracy provided that the momentum cutoff is not too drastic. In the present paper, we explore an alternative approximate evaluation of the coordinate representation freepropagator action on an arbitrary wave packet. Ultimately, our results can be viewed as a different sort of coarse-graining approximation. Before the explicit details are given, it is useful to summarize the strategy to be employed. Essentially, we propose to take advantage of the fact that the action of the free-particle propagator on a Gaussian wave packet is known analytically. Thus, we do not rely on an approximate quadrature treatment of the x,x’ matrix element of the propagator as in our earlier work.s9 We also do not directly use the exact F ( x , x ~ T(which ) is highly oscillatory, highly nonlocal, and computationally inconvenient). Instead we work directly with the exact analytic result of F(T)acting on a Gaussian and related Hermite polynomials, which are orthogonal under a Gaussian weight function and can be obtained via a Gaussian generating function. As a consequence, the action of the free propagator F(T) on the Hermite polynomials can also be obtained analytically, However, we do not simply expand the wavepacket in terms of (70) Messiah, A. Quanrum Mechanics; Wiley: New York, 1966; pp 216-242. Powell, J. L.; Craseman, B. Quantum Mechanics; Addison-Wesley: Reading, MA, 1961; pp 69-84, 98-100. (71) Pechukas, P. fhys. Rev. 1969, 18I, 166, 174. (72) Miller, W. H. Adu. Chem. fhys. 1974,25,69, and references therein. (73) Marcus, R. A. Faraday Discuss. Chem. Soc. 1973, 55, 34 and references therein. (74) Doll, J. D.; Myers, L. E. J . Chem. fhys. 1979, 71, 2880. Doll, J. D. Ibid. 1984, 81. 3536. (75) Chandler, D.; Wolynes, P. G. J . Chem. fhys. 1981, 74, 4078. (76) Thirumalai, D.;Bruskin, E. J.; Berne, B. J. J. Chem. fhys. 1983,79, 5063. Thirumalai. D.;Berne, B. J. Ibid. 1984,81, 2512. (77) Miller, W. H.; Schwartz, S. D.;Tromp, J. W. J. Chem. fhys. 1983, 79, 4889. Yamshita, K.;Miller, W. H. Ibid. 1985,82, 5474. Tromp, J. W.; Miller, W. H. Faraday Discuss. Chem. Soc. 1987,84,441; J . fhys. Chem. 1986, 90, 3482. Chang, J.; Miller, W. H. J. Chem. fhys. 1987,87, 1648. (78) Coalson, R. D. J. Chem. Phys. 1985,83, 688. (79) Doll, J. D.; Freeman, D.L. Science 1986, 234, 1356. Doll, J. D.; Coalson, R. D.; Freeman, D. L. J . Chem. fhys. 1987,87, 1641. (80) Behrman, E. C.; Wolynes, P. G.J . Chem. fhys. 1985, 83, 5863. Cline, R. E.; Wolynes, P. G. Ibid. 1988, 88, 4334. (81) Hendriksen, N. E.; Heller, E. J. Chem. Phys. Len. 1988, 148, 277. (82) Makri, N. Chem. Phys. Lctr. 1989, I59,489. Makri, N.; Miller, W. H. J . Chem. fhys. 1989, 90,904. (83) Sethia, A.; Sanyal, S.;Singh. Y . J . Chem. fhys. 1990, 93, 7268.

a Gauss-Hermite basis. Instead, an approximating function is introduced for representing an arbitrary wave packet in terms of the values of the packet on a discrete (infinite) grid, times fitting functions of x associated with each grid point, x,, These fitting functions themselves can be represented as an expansion in terms of Hermite polynomials, so that propagation of the wave packet is expressed in terms of the propagation of the fitting functions, evaluated in terms of the propagated Hermite polynomials. Then the effect of the propagator on the wave packet, which determines the amplitude at the grid point xy due to the particle starting at point xj, can be obtained by simply reading off the coefficients of the packet at grid point x . when the packet is evolved forward in time by an amount T . ?chis turns out to be the value of the fitting function at time T , evaluated at the point ( x r x j ) . The fact that our fitting function yields only an approximation to the wave packet (even at the grid points) implies a great deal of latitude in its functional form. The desire ultimately is to choose an approximating function such that the matrix element of the propagator for the contribution of the packet at xi and time t T to its value at xf and time t becomes very small as Ixy-x,l becomes larger. We speak of this as limiting the “range of influencen of the grid point x . during the time interval T . As the matrix representation of F ( Tbecomes ~ sparser, or more banded, the range of influence decreases. Finally we wish to develop a banded approximation to F(T)that will be accurate for propagating a general wavepacket. To achieve this, the fitting function will be required to fit exactly all polynomials up to a given order. This is combined with the Hermite expansion of the fitting function to derive a linear set of algebraic equations for the expansion coefficients, which are made approximate by truncation of the expansion. The coefficient matrix of these algebraic equations has elements that are expressed as an infinite, discrete sum of products of two Hermite polynomials and a Gaussian weight function. If this discrete sum is approximated by an integral, then we find that the expansion coefficients for the hermite representation of the fitting function can be evaluated analytically. Finally, having these expansion coefficients,the approximate fitting function at time t + T is readily evaluated, and its value at the grid points (xy-xj) yields the discretized approximation to the propagator, FyJ(7). We now turn our attention to discuss the explicit details of the procedure. Coordinate Space Approximation Formula for an Arbitrary Wave Packet. Let us assume we know the wave packet on an infinite grid of equally spaced points. We first wish to develop an approximation formula of the following form: 0

AX) = IC a j ( x ) fix,) .--

(5)

where aj(x) is termed a “fitting function”. Thus associated with each grid point xj.is a fitting function u,(x), which will be chosen such that its dominant contribution tofix) comes from the locus of points close to x,, It is important to keep in mind that eq 5 does not represent an interpolation formula. The importance of this distinction lies in the fact that aj(xk)is not rigorously equal to 6. This has important consequences on the banded structure of ti: free propagator as we later discuss. We assume aj(x) to have a translational property defined by

.

a,(x) = a&

- Xj)

(6)

From elementary considerations the fitting function ao(x) must be symmetric. Obviously, however, these requirements do not uniquely determine its functional form, even if one also specifies the ‘range of influencen of the grid point value, Ax,), in the approximation formula (Le., the extent to which aj(x) is nonzero away from the grid region around x,). The explicit form of a&) depends on the approximation scheme one wishes to employ, which is itself governed by the gross features of the packet one wishes to propagate. For example, with only very mild restrictions, any wave packet with momentum components less than r/Ax can be expressed exactly as a Fourier series, in which case u,(x) = sin

(A(X

- xj))/(n(x- xj))

(7)

8302 The Journal of Physical Chemistry, Vol. 95, No.21, 1991 This, however, is an interpolation formula for which the range of influence of Ax,) falls off very slowly (i.e., as the inverse of the separation distance from the grid point x,). As a practical matter, to approximate typical wave packets to sufficient accuracy, fitting functions with much narrower band widths will suffice. As we shall see,this can be exploited to develop a sparse or banded approximation to the free propagator. A. Qlipmid Basis. The free propagator, F(T),is defined by F(T)= exp(-iTK/h)

(8)

where K is the Cartesian kinetic energy operator. A Gaussian function with standard deviation u(0) evolves under this propagator according to the expression70

F(r)exp(-(x

P = X/(W (16) where Ax is the x-grid spacing. This new variable is convenient because the x-grid points are mapped onto the pgrid points, which fall on the integers. Then to specify the fitting function %(x) (and the concomitant translated fitting functions a,(x)), we require our fitting function to exactly fit all polynomials of order 0 In I m. According to eqs 5 and 6 we must have m

- x0)2/2.2(7))

p" =

(9)

C jx-m

ao(p-j)j"

0 I n Im

(17)

or equivalently, by making use of the symmetry of a.

where $(T)


for the coefficients, bznrin eq 14. The results in the next section suggest that for an analytic function, this approximation can be made to hold to arbitrary accuracy by properly controlling the spread of the Gaussian envelope of the fitting function. We begin by making a variable change

- X ~ ) ~ / ~ U ~=( O ) ) [ 4 0 ) / 4 7 ) 1 exp(-(x

= $(O)

+ ihr/m

m

exp(-x2)H,,(x) = (-l)"(d"/dx") exp(-x2)

0 I n < = (11)

It immediately follows that these functions can also be propagated analytically in lD, since the nth derivative with respect to x commutes with F(T),and exp(-$) is a special case for eq 9. Now if we define a scaled variable, ~ ( x , T )by , y(X,T) = (2)-'/2X/U(T)

(12)

then application of F(T) to eq 11 above leads to F(7) exp(-y2(x,o))HnCv(x.o)) = (-1 )"(d"/dy(x,O)n)F(r) X ~ x P ( - Y ~ ( x= ~ )[)d o ) / u ( ~ ) I & ' exp(-YZ(x,r))HnCv(x,T)) (13)

Thus, if we write the function ao(x) in terms of the Hermite basis set of eq 11 (see below), we find that the fitting function can be propagated according to

Cb2n[fl(o)/dT)I2"+'

n=O

(-p)" = j x - m ao(j + p)jn

(10)

and m is the particle mass. This is a special case of the well-known analytic result for a freely propagated Gaussian wave packet.70 (The parameter xois not a function of time because the momentum of the Gaussian function, considered as a Gaussian wave packet, is zero.) We wish to take advantage of this by introducing a basis that can be evolved analytically by use of eq 9. The Hermite polynomials form a complete set of orthogonal polynomials under a Gaussian weight function and are given in terms of their generating function by

Q(X,T) =

Hoffman et al.

exp(-y2(x,T))HznCv(x,T))

(14) where ao(x,O) ao(x) and the coefficients, bzn, are constants (albeit, unknown since we do not yet know the functional form of ao(x)). Only even Hermite polynomials appear because ao(x) is symmetric. Equations 5,6, and 14 together yield an expression for the propagated function. Further, because of the translational property of eq 6, the translation of any of the a,(x) follows from eq 14. An analytic expression for the discretized propagator is then given by FyJ'I) = ao(xpx,J) (15) since multiplicationof ao(x,-x,,~) by theAx ), followed by summing over 1 from -= to = yields exactly Axf,.), which when combined with eq 5 clearly yields Ax,T). We now develop an analytic approximation for this quantity. To do this, we impose restrictions on the fitting function, %(x), by requiring that it exactly fit a certain class of functions. Derivation of a O ( x ) . We wish to derive an explicit form for ao(x) that ultimately we shall confine, spatially, to the greatest extent possible, so as to obtain a minimum bandwidth for the propagator. Since we are deriving a fitting function assuming knowledge only of values on a grid of the function to be propagated, there is clearly no unique expression for ~ ( x ) However, . we now outline a general approach that yields explicit expressions

0I n I m

(18)

+

Writingj in the form j = (j p) - p and substituting into (18), we find that this set of equations can be reduced to the set m

+

= j.-m ao(j p ) ( j

+ p)"

0I n I m

(19)

We now convert these equations to an equivalent set written in terms of Hermite polynomials (in the scaled variable of eq 12 rather than in powers of (j + p ) ) . These equations have the form m

hJn) = j aC ao(j + p)H,,Cv(j + p,O)) -m

0 In Im

(20)

where

ho(") = (-1)"h!/(n/2)!, n even = 0 , nodd (21) is the coefficient of the constant term in the nth Hermite polynomial. Proceeding further, it is convenient to write ao(p)in the form sob) =

?bib) exp(-y*(p,o))HnOt,(p,O))

n=0

(22)

where the sum is over afinite number of Hermite functions and b',,(p) is a periodic function of p, which is symmetric for even values of n and antisymmetric for odd n (c.f., eq 14; the product of b',,(p)with H,,Q(p,O)) is even in p when n is odd). Substituting eq 22 into eq 20 and making use of the periodicity of b',,(p), we obtain the set of linear equations

where the (symmetric) matrix element, C,,&), is given by

Cn&)

=

2 exp(-y2~+ p,o))HnQ(j + P,o)) H ~ ~ C V+OP,o)) '

ja-m

(24)

Equation 23 uniquely determines the functions, b',,@). To complete the evaluation of the b,, coefficients in eq 14, we could proceed at this point by reexpanding the product function, b ' , , @ ) H n ~ ( i -p,O)),in Hermite polynomials of even index; however, we choose a slightly different approach. If u(0) is sufficiently large and the order of the largest Hermite polynomial required is not too high, the sum over j in eq 24 can be replaced by an integral. (These qualitative statements will be quantified in the next section.) In this case, one finds that the C,,d are independent of p , so that the functions, b',(p), become constants that are trivial to evaluate because of the orthogonality of the Hermite polynomials. For odd n the constants are zero; for even n they will be used as the constants of eq 14, thereby yielding an approximation to the fitting function that fits all


Discretized Free Propagator polynomials of degree m or smaller. Thus, we finally obtain the explicit result

I4

l2


(25) As will be illustrated in the next section, this expression displays a kind of asymptotic convergence. For a fixed a(0) the fit gets better with increasing m up to a point and and then steadily grows worse. The conditions required for eq 25 to give an accurate approximation to ao@)can always be met by broadening the envelope of ao@)as we presently discuss. Alternatively, one can, of course, choose a smaller grid spacing, Ax, thus reducing the number of Hermite polynomials needed to approximate the function adequately. The pertinent consideration here is the *smoothness” of the function being fitted which is ultimately controlled by the grid spacing. Combining eq 15 with eq 25, we finally obtain the analytic approximation m/Z

F/,/(r) =

CKZnO”-AT) n-0

(26)

for the discretized propagator, where

Using the standard recursion relation for the Hermite polynomials

= 2 x q w - 2m..I(X)

e 0

Figure 1. Measure of the number of significant figures, e = -log V; @) -f@)l, plotted versus the number (m 1) of Hermite polynomials. %e highest Hermite order included is 2m. The diamonds denote o(0) = 3.54 points, the plus signs are for o(0) = 2.36, the squares are for a(0) = 1.77,

+

and the crosses are for a(0) = 1.41.

j,T)Z)HznCvO” - j , T ) ) / [(2*)1/2n!I (27) H,+l(X)

t

(28)

it is a straightforward matter to obtain also a recursion relation for K&)T) (although to do so it is necessary to extend the definition, eq 27, for the K functions to odd index in the obvious way). After some algebra we find that

~n+~@b) = + ~ @ I ~ ) g j [a(O)/d~)l~n@I7) +l + Yz(n/(n + l))[a(o)/a(T)l2K,l@IT) (29)

given by eq 25, is able to fit a given function&), where we shall use the variablep (seeeq 16). We take as the function to be fitted &)

= exp(-j9/32) cos ( 0 . 7 ~ )

cos (7%)in This corresponds to the function exp(-~~/[2(0.4)~]) terms of the Cartesian variable x. We study first how accurately it is described by eqs 5 and 25, as a function of the number of Hermite polynomials included for ao(x), and the width of the Gaussian weight function a(0). Sample results for an arbitrarily chosen point, p = 0.08, and a(0) equal to 3.54, 2.36, 1.77, and 1.41 are shown in Figure 1. Plotted is roughly the number of significant figures, e, obtained in the fit, versus the number of Hermite functions used for the 4 values of a(0). Here

where

= -log K - , @ ~ T )=:

0

Ax

KO@lT)

= 75[0(0)/U(T)1 exP(-Y%w))/(2*)”2 g] = (r/2)1/Z

(31) (32)

2/[(n + 1)gnl (33) Thus, our analysis provides an analytic, closed-form expression for the discretized free propagator expressed as a finite sum of terms that can be generated from a three-term recursion relation. Application of the banded propagator matrix, F/J(T),to the column vector of discrete values of the wavepacket at the xrgrid points generates the j’component of the wave packet at a time interval T later. The evolution of the packet then consists of repeated multiplications of the matrix F’JT) to the wavepacket. Extension of these results to nodartesian variables and to higher dimension is currently under investigation and will be the subject of subsequent communications (see section IV). However, it is worth remarking that the propagator for an s-wave in polar coordinates is given by gn+l

o’/j?[F’,, - F’.-,I, j’, j 1 0 (34) where again we have assumed that the coordinates have been scaled so that the radial grid falls on the integers. (Although this expression is indeterminant for j ’ = 0,the matrix element vanishes when either j ’ or j is equal to zero). 111. Numerical Results

We now illustrate some features of the present approach to discretizing and banding the free propagator. To begin, we shall examine how accurately the approximate fitting function, Q,@)

(35)

If@) -f*,@)l

(36)

Since only even Hermite polynomials are involved, the number of Hermites is times the order of the highest polynomial, +1. The Ax spacing of the grid is Ax = 0.01, and x = p(Ax). We note first that the quality of the fit varies smoothly with the number of Hermites. It is clear that for each u(O), the quality of the fit fmt rises as the number of Hermite functions is increased, but past a certain point, the fit becomes progressively worse (for u(0) = 3.54, this point occurs off the plot at 78 Hermites). However, no stability problems are encountered. The apparently asymptotic nature of the fit presumably results from a breakdown of the approximation of replacing the sum over j in eq 24 by an integral. As a(0) increases, more Hermite polynomials can be included and higher accuracy achieved. As the width of a(0) decreases, the accuracy of the fit deteriorates, but the number of Hermites required to obtain maximum accuracy also decreases. In Figure 2 we consider the effect of making the function&) to be fitted more rapidly oscillatory. Thus, we consider fits of the class of functions &)

= exp(-p2/32) cos (k.4

(37)

where k (k is the momentum conjugate to the coordinate p ) is taken to be 0.4,0.7, 1.0, 1.3, and 1.6. The width of the Gaussian weight for the fitting function is taken to be u(0) = 2.36. Again, e is plotted versus the number of Hermite functions used in approximating the fitting function ao@) via eq 25. It is rather encouraging that the number of Hermite functions yielding optimum results is almost independent of the function being fitted (although all of the functionsm) are of the same generic form). It is also to be noted that the less oscillatoryfi) are more accurately fitted by orders of magnitude with 25 Hermite functions included in ao(x) than are the more oscillatoryflp). This loss

Hoffman et al.


8304 The Journal of Physical Chemistry, Vol. 95, No. 21, 1991

TABLE I: Fitting Fuoction a,@) Evaluated at 20 Grid P o h b with ~ ( 0 =) 3.54 Using ”e Different Numbers of Hermite PolyaomiaW P m=7 m = 14 m = 21 5.938 I899968361 61D-01 0 3.545468720829894D-01 4.890320208469076D-01 2.985996427495615D-01 1 2.8005240018761 80D-01 3.1 18726962120446D-01 2 1.166662322805951D-01 1.032720775207557D-02 -8.1 574321 31 881 334D-02 -8.8242827 19922326D-02 -0.565012449348001D-01 3 -1.708915123913788D-02 4 -5.594533355259128D-02 -8.598632392432148D-03 5.329960270324258D-02 5 -2.636458957103059D-02 3.805056693604263D-02 4.442974231 376523D-03 6 8.908390103656035D-03 6.285991575002825D-03 -0.2557 14672104994D-01 7 1.717733615474377D-02 -1.643730052365748D-02 7.360547089491 254D-03 8.4259531402621 15D-03 8 6.731429486037197D-03 -3.9935241 90909943D-03 -5.856816350328715D-03 9 -2.95 1952260409726D-03 6.429668495366497D-03 IO -4.436094467450371 D-03 2.177809876667323D-03 -1.4457725503941 5 5 W 3 I1 -1.60298491 5 l67141D-03 -2.150356794584995D-03 2.599357404664657D-03 -1.004826503308220W3 -2.2013 12453450756D-04 12 5.735 l989202097lOD-04 -7.924074646455041 D-04 13 8.685664792910462D-04 5.76238591 161 1586D-04 14 3.469171579956059D-04 3.85 I78 1708818005D-04 2.558888852975339D-04 15 -3.950936235 15 14991)-05 -1.078191 154504618D-04 1.706112882279607D-04 16 -1.158794342765743D-04 -1.195077493675571D-04 -1.01 74796721271 10D-04 17 -5.96236401 1275 l02D-05 6.451421609929189D-06 -2.436000727445084D-05 2.87205 19 18579976D-05 2.72891 5896957396D-05 18 -9.039952658843801I6 19 7.342911218249110D-06 4.1 566393838715 IOD-06 1S16085265855748D-06 20 6.2 14044608340459D-06 -4.868565 122035786D-06 -5.61 71 17535089199D-06 sum 9.99995252688 1 8 4 8 W 1 1.000002067157639D-00 9.99997992822638 1D-01

“The values are symmetric about zero. The sum is from -20 to 20. bD-Ol means times IO-’. TABLE II: Fitting Function a ( @ ) Evaluated at 20 Grid Points with u(0) = 2.36 Using Wee Different Numbers of Hermite PolyoomiaV 4 m=7 m = 14 m = 21 0 5.3 1820308l244842D-0l 7.335480312703613D-01 8.907284995254242D-01 1 3.0302829 19679426D-01 2.26209702451 1312D-01 1.025883301145940D-01 2 -2.563372685870682D-02 -1.323642407988349D-01 -8.4751 86740220014D-02 -6.850146844348796D-02 0.409427383647223D-01 6.127774087744189D-02 3 4 1.336258515548405D-02 9.428987362504238D-03 -0.3835720081 57491D-01 5 1.941468361736353W2 -1.876151436223907D-02 2.036325986820049D-02 6 -4.4279283906 14589D-03 -8.785224525493073D-03 9.644502743049746D-03 7 -1.117833516903085D-03 2.748791873010633D-03 -4.762510810541 567D-03 8 8.6027983803 14565D-04 -1 507239754962329D-03 -3.3019686801 761 35D-04 9 9.485078449949100D-04 -2.827944669938354D-04 9.465648664652409D-04 IO -5.926404352727259D-05 -1.61 7286731756929D-04 2.559 1693234194IOD-04 11 -1.390438375972854D-04 -7.265633893225665D-05 -1.241017341765646D-04 12 -1.355992898826570D-05 4.308077877869964D-05 4.093373845436095D-05 I . 165224007758289D-05 13 -4.42001 3523298274D-06 -8.542055347522355D-06 -2.73368443758 17 17D-06 3.912698806492329D-07 3.5861 2956183 1445D-06 14 8.235199453660597D-07 4.794362662691727D-07 -3.393530468472720D-11 15 5.4755018 19076569D-08 -2.03861 2660060736D-07 -2.202383054568 160D-07 16 -6.367536723 133574D-08 17 -4.86091 3999751044D-08 4.l012l0900972738D-08 -1.03625 1407321909D-08 5.529771 117767747D-10 -2.443868753959866D-09 I8 19 -1 . I 533541981 32456D-09 1.795358066172557D-09 -9.6508495241 20983D-10 -9.4107351 153433761)-11 2.76964622436902OD-10 20 -1.515341694122139D-12 1.000000000012242 1.000000000099994 1.000000000047557 sum

‘The values are symmetric about zero. The sum is from -20 to 20. of accuracy in fitting the more oscillatory functions can be addressed by increasing the width of u(0) of the Gaussian weight function, or by decreasing the Ax of the grid in x. Finally, it appears that a great deal of fine tuning of the fitting function u&) can be done before actually propagating the full solution of the time-dependent Schradinger equation. Of even greater interest is the question of how strongly peaked the approximate fitting function ao(p) is about p = 0, since by eq 15, this will determine the degree to which the discretized propagator matrix, FiJ,is banded. We have carried out calculations of a,@) for several choices of the parameter u(0) (fixing the width of the Gaussian weight function at time t = 0) and the number of Hermite functions included in approximating uo@). The results are tabulated in Tables 1-111. In Table I, u(0) = 3.54, and the number of Hermite functions are 7, 14, and 21, respectively. The results are given for &), with integer j ranging from 0 to 20. By the translation property, these are also the results ) - j’) with 0’”- j ’ ) equal to 0 I j ‘ l - j ‘ I 20. for ~ ~ 0 ”=’aoC” We see that the degree to which a0@) is peaked about p = 0 increases with the number of Hermite functions. Further, in Table 11, u(0) is decreased to 2.36, and that it is found that generally a&) becomes substantially more banded compared to the results

with u(0) = 3.54. However, the bandedness is not a simple monotonic function either of the number of Hermite functions used for a fixed value of u(O), or of the u(0) for a particular number of Hermite functions. This is the consequence of a grid point falling near a node of the fitting function, which gives rise to an apparent somewhat eratic behavior. In Table 111, ao@)is given for u(0) = 1.41 and 4 and 8 Hermite functions. We note that the degree of bandedness is primarily dependent on U ( T ) as opposed to the number of Hermite polynomials used in the approximate fitting function. This reflects two effects. First, at T = 0, it is the Gaussian weight function in a&) that exerts the dominant effect causing ao@) to decay as a function of p . For time evolution, U ( T ) depends on u(0) and T according to eq 10. This implies spreading of the Gaussian weight function with time that affects the bandedness. Clearly, the smaller is u(O), the more sharply peaked is the Gaussian (see eqs 9 and 25) as a function of p. Comparing results in Tables 1-111, we find that those with u(0) = 1.41 are by far the most strongly banded. Finally, we note that it has been found possible to freely propagate wave packets with the present method for arbitrarily long times without significant loss of accuracy. Further, such a free propagation can be achieved in a single time step or by a


Discretized Free Propagator


propagating wave packets solely in the coordinate representation. The critical requirement for this to prove efficient is that a banded (or sparse) matrix representation of the free-particle propagator can be developed so that the application of this propagator will involve as few multiplication operations as possible. There remains a number of problems that must be addressed in order to complete the development of a general wave packet propagation procedure. We discuss briefly below some of these that are under current study. Dynamical Basis Set in Higber Dimensions. The extension of the analysis leading from eq 5 to eq 25 for a Cartesian grid in more than 1D is trivial and requires no elaboration here. For a nodartesian grid (e.g., spherical polar) the situation is less obvious and requires exploration. We want to develop an a p proximation formula like that of eq 5, with the sum extending over the whole grid, and to introduce an invariant fitting function as in eq 6 that possesses the translational property with respect 0 5 10 15 20 25 30 35 40 45 50 to every grid point. With use of a complete set of functions that N e d - is a direct product of those used in the 1D case, a propagation formula analogous to eq 14 can be derived. (Drawing on the Figure 2. Measure of the number of significant figures, e, plotted versus familiar result that, for example, the wave functions for the isothe number ( m + 1) of Hermites for a family of oscillator functions of the form, exp(-d/32) cos (kp). The cro13ses are for k = 0.4, the triangles tropic oscillator in three dimensions can be obtained in both are for k = 0.7, the diamonds are for k = 1 .O, the plus signs are for k Cartesian and spherical coordinate systems with a simple coor= 1.3, and the squares are for k = 1.6. dinate transformation between them, it seems intuitively clear that the higher dimensional analogy to eq 14 need not necessarily be TABLE III: Fitting Function a,@) Evaluated at 20 Grid Points with expressed in Cartesian coordinates.) Following through the aru(0) = 1.41 Using Two Different Numbers of Hennite Polynomialsa gument, however, it would seem that there is no obvious analogue P m=4 m=8 to eq 23 unless one has the periodicity provided by a Cartesian grid; however, this still might not ultimately defeat the method 0 6.942176516310283D-01 9.417651 289704407D-01 1 2.3 1746OO80599210D-01 5.217096057803501 D-02 if one is able to extend the approximation of replacing the sum 2 -9.053451278378856D-02 -3.733952448461333D-02 over the mid bv an intearal. This is currently under investigation. 3 4.5731 251 34941802D-03 2.104712395952302D-02 hop&& shies f0; Coordite Space-Computations: NU4 9.701 104453320554D-03 -9.0985251 13136650D-03 merical experimentation for some example scattering and time5 -2.327 101744299739D-03 2.8767 18 162524704D-03 dependent dynamics problems will be carried out to test the ac6 -3.843057859389741D-04 -6.028029345129035D-04 curacy and efticiency of the propagating procedure outlined above 7 9.789646730961092D-05 6.099350358227510D-05 and to explore the relationship between bandwidth, grid spacing, 8 1.790856175646545D-05 4.126186269366692D-06 and time step. Also application of generalized methods to colliisions 9 1.024596621390184D-06 -1.71 1779195931949D-06 involving internal degrees of freedom will be carried out and 10 2.686232400944625D-08 6.638 118597644133D-08 compared both to alternative wave packet methods and to time11 3.605032724698174D-10 1.200224499529556D-08 independent approaches. 12 2.60821 991 961 7377D-I2 3.5251 66260291 3 18D-IO 13 1.047884667268035D-14 4.096653936512238D-12 Roposed Studies for pualkl Rocxdng. For reasons previously 14 2.382626388774752D-17 2.244938676023936D-14 described, coordinate space methods have the potential of being 15 3.1 0633280431 261 8D-20 6.261 59639213675OD-I7 particularly convenient and efficient for performing time-de16 2.344263092093616D-23 9.278109002105353Q-20 pendent scattering calculations on parallel processing machines. 17 1.0314131 136898521)-26 7.505654295636144D-23 Furthermore, the propagation procedure outlined above has the 18 2.660318137558049D-30 3.378089603512622D-26 design feature that the computations for each time step require 19 0.0 0.0 a very small number of operations. Indeed, one of the primary, 20 0.0 0.0 initial motivations for pursuing such an approach based entirely sum 1.000000000001604 1.000000002607525 on using the coordinate representation for all operators was the plan to implement the method on massively parallel computers. OThe values are symmetric about zero. The sum is from -20 to 20. The applications mentioned above will be coded for and tested sequence of them. The present formulation can be used in any on the two parallel processing machines (a MasPar with 8192 of the existing algorithms that involve the free propagator, inprocessors and an NCUBE 6400 with 64 processors) currently cluding the split operator,"S6' the interaction p i ~ t u r e ~and ~ ~ * * ~available ~ at the ISU Scalable Computing Facility. The results the kinetic energy referenced integral equation methods.47*49~50J9 will be compared with those obtained by Fourier transform based We have tested this method for some standard problems and algorithms on the same computers. Once this is completed, several obtained excellent results. calculations of increasing complexity will be undertaken: first, These results demonstrate that it is indeed possible to develop a 1D electron tunneling problem in the solid state, with an applied highly banded approximations to the action of the free propagator time-varying electric field, will be studied; second, a 3D electron on a general wavepacket. Such banded matrices are of interest tunneling problem in the solid state relevant to nanostructures because they greatly reduce the work involved in applying the free will be treated. In all cases the results and concomitant efficiencies propagator to an arbitrary packet. Since this is a key computawill be compared with calculations carried out on a serial machine. tional component in a number of numerical methods for solving Particular attention will be paid to finding the most efficient way time-dependent uantum mechanical problems (e.g., the split to assign grid points to processors for realistic problems where operator method! and the recently developed time-dependent the number of grid points greatly exceeds the number of processors. integral equation method^".^^^^.^^^^^), we are optimistic that the Acknowledgment. D.K.H. and D.J.K. express their great adapproach will prove useful. miration of and gratitude for Dick Bernstein. He epitomized all IV. Future Work that is best in science and was an inspiration to us. Lastly, he was a valued friend and encourager. We shall miss him sorely. The ideas presented in this paper represent an important step in our overall strategy to develop a new, efficient method for Registry No. H20, 7732-18-5.

Analytic banded approximation for the discretized ... - ACS Publications

Recommend Documents