- Email: [email protected]
Time-dependent wavepacket methods for the calculation of collinear atom-diatom exchange reaction probabilities
Computer Physics Communications 63 (1991) 209—215 North-Holland
209
Time-dependent wavepacket methods for the calculation of collinear atom—diatom exchange reaction probabilities N. Balakrishnan and N. Sathyamurthy Department of Chemistry, Indian Institute of Technology, Kanpur 208 016, India Received 7 february 1990
Different strategies available for extracting the translational energy dependence of state-selected reaction probabilities for a collinear atom—diatom exchange reaction are reviewed and the usefulness of some of them is examined using collinear He+ H~(v~ = 0) -~ HeH~+ H as a test case.
1. Introduction The time-dependent quantum mechanical (TDQM) wavepacket method has emerged as a powerful tool for the study of atomic and molecular collision processes in recent years. Its success emanates from its simplicity, ease of interpretation of the results and the detailed physical insight into the dynamics obtained. One can view snapshots of the wavepacket at every stage of time evolution and get a “feel” for the reaction, analogous to the classical description of the process but within the quantal framework. With the introduction of accurate algorithms for the propagation of wavepackets in space as well as time, the method could yield results of accuracy comparable with that of the time-independent coupled channel methods. The TDQM method has beeti applied successfully also to the study of eigenvalue problems, photodissociation, Raman scattering, modeselectivity in chemical reactions, gas—surface scattering etc. The method, its implementation and its application till 1987 have been reviewed elsewhere [1—3].For more recent applications, see refs. [4—28]. We refer to the TDQM method to mean that the time-dependent Schrodinger equation (TDSE) INSA research fellow (1989—1991). 0010-4655/91/$03.50 © 1991
—
is solved numerically for the problem under consideration; it is also referred to as the grid method. It consists of three parts: (i) choosing the initial wavefunction for the system at time t~,‘I’(t~),(ii) propagating it in space and time and (iii) analyzing the wavefunction, ~P(t), after time t to extract reaction attributes. Major developments in the methodology in the 1980’s concerned primarily part (ii): Kosloff and Kosloff [29] pointed out thai the Laplacian of the wavefunction could be evaluated using fast Fourier transform (FFT) aI~ gorithms; Askar and Cakmak [30] had earliet pointed out that the time-evolution could be car ried out by a second-order difference scheme pro. vided that the time increment met certain stabiliti criteria; its equivalence to the explicit metho (also called the Richardson scheme) [31] has beei discussed elsewhere [2,32,33]. Feit et al. [34] dem onstrated the utility of a split-operator time prop agation scheme. Tal-Ezer and Kosloff [35] showed that the time-evolution could be carried to mac hine accuracy by expanding the evolution operatc in terms of complex Chebyschev polynomiali Chattaraj et al. [36] proposed an unconditionall stable matrix inversion scheme but it turned out t be computer intensive. Recently Tal-Ezer et a [37] have introduced a low-order polynomial a~ proximation of propagators, particularly when tF Hamiltonian is time-dependent. As an outcome
Elsevier Science Publishers B.V. (North-Holland)
N. Balakrishnan, N. Saihyamurthy / Time-dependent wavepacket method
210
a CECAM workshop, Leforestier et al. [38] have carried out a critical analysis of the accuracy and stability of the different algorithms used for ternporal evolution, Regardless of the method used for the propagation of the wavefunction in space and time, there are some practical difficulties. For example, considering one of the simplest reactions possible, an atom—diatom exchange reaction constrained to a line, A
+
BC(v,)
—‘
AB
+
C,
(Ri)
the initial wavefunction would be, typically, a product of the vibrational wavefunction corresponding to the initial state v, of the diatom and a translational wavepacket, often a minimum uncertainty Gaussian. That is, there is a Gaussian distribution in the momentum space corresponding to the relative translational motion of the reactants. Choice of such a wavepacket becomes necessary as otherwise we have to start with a delta function in the momentum space in order to compute the reaction probability, pR, at a fixed Etrans. This would mean that we use an infinitely spread out wavepacket in the configuration space and thus an infinitely large grid. In practice, the dependence of pR Ofl Etrans is extracted from the wavefunction at the end of the time evolution by a deconvolution procedure [39]. Agrawal and Raff [40] had proposed the use of rectangular wavepackets to compute P’~at a given Etrans. Unfortunately, the method calls for time evol4ving a pair of wavepackets centered at Etrans at which the value of pR is sought and if the energy dependence of P’~is required, the calculations have to be repeated again and again at different Etrans. Therefore it is computationally not appealing. Yet another difficulty of these TDQM calculations is that, with increasing time, the wavepacket spreads over the reactant and product spaces and before the end of the calculation it gets reflected from the edges of the gnd causing artificial interference effects. One way to avoid such spurious effects is to increase the size of the grid which unfortunately increases the computational time. Some authors [41,42] have circumvented the prob1cm by using absorbing boundaries but this im-
plies loss of information on the product (reactant) states. Heather and Metiu [43] proposed a wavefunction splitting algorithm to keep the grid extended barely beyond the interaction region. Recently Zhang [44] has shown that by switching over to an interaction representation, finite boundary reflection is eliminated, with concomitant reduction in the number of grid points and increase in the time increment for the numerical integration. Neuhauser and Baer [24] have developed a novel projection operator formalism which eliminates the need for inclusion of the product channel. Their approach calls for a relatively small rectangular grid in the interaction region and a long one-dimensional grid in the approach coordinate thus reducing the computational time. It also provides for a choice of a wide wavepacket in the coordinate space enabling the computation of pR at a narrow range of Etrans. The authors demonstrated that they could reproduce the dip in PR(Etrans) curve for H + H 2(v, = 0). We have chosen to examine the usefulness of some of the available approaches using the collinear He + H~(v 0)—’ HeH~+H (R2) =
reaction as a test case since the recently reported [45,46] converged R-matrix calculations on an accurate ab-initio potential energy surface (PES) have revealed reactive scattering resonances in the energy range Etrans 0.82—1.35 eV. Also quasiclassical trajectory (QCT) calculations have been carried out [47,48] over a wide range of Etrans with pR varying from 0 to 1. In the following section we outline the different approaches we have used in computing PR(Etrans) by numerically solving the TDSE and also present the results obtained. A summary of our findings is given in section 3. =
2. Different strategies to compute pR ( Etrans) 2.1. Gaussian wavepacket (GWP) method
The basic philosophy of the grid method is to solve the TDSE ih
( \
~
=
~t
I
I-J’I’,
(1)
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepacket method
where H=
wherein the superscript refers to the number of 2
2
~
2ILA,BC) h
—
Bq~+ 0q
+ V(
-~—~
2
time steps. We have used Lit
qi q ‘
2)
(2)
in mass weighted coordinates for a collinear A BC system: 3rsc sin 9, q~ rAB + /
+
=
q2 $~‘BC cos with
(3)
=
mC(mA + mB) m~(m~+ mc)
1/2
(\ (m~+ mB)(mB mAmC + mc) /
=
In such a coordinate system ‘I’( t
=
(5)
0) is chosen as
0) ‘P(q1, q2) ~(q1)x(q2), (6) where ~(q 2) is vibrational the Morse oscillator wavefunction for the ground state of H~ and ~(q 1) is the minimum uncertainty Gaussian: =
=
=
=
=
2 ~q1,
~,
with 3 0.25 au. q10 refers to the centre of the wavepacket and k0 is related to its average translational energy,
—
~2 /4
—
‘~
xexp
~ )/
~
—
q’°~
—
2itCh/pAB~
46
/4
2 32 1A,BC
—1/4
I-
ik0q~
(10)
,
—
.
which contracts to a nMnimal Gaussian wavepacket after the dilation time t~ which was chosen to be (1000—2000) Lit, depending upon the grid size. The time evolution was followed by monitoring the reaction probability
=
f
—
~,
2\
h
2 34 !~A,BC 2
2
(q
1—q10) _iko~i) 46 (7)
—
=
=
=
~(qi)=(2~62)1/4exp(_
0.0269 fs. To ii-
=
(4)
and sin 9
=
tialize theusing temporal ‘I~was from ~I,0 eq. (1)evolution, by expanding the evaluated evolution operator in Taylor series in Lit and truncating after the second-order term. After preliminary investigations using a 64 X 64 grid in (q 1, q2) space with q~fl 1.0 and ~ 0.4 with Liq~ Liq2 0.1 au, we carried out the time evolution of ‘I’ at
=
211
product
I ~I’(qi~ q 2)
2
dq1 dq2
(11)
ILA,BC
where !~A,BC is the reduced mass of the reactants, The potential energy, V, for the system was taken to be the analytic fit of the McLaughlin—Thompson ab-initio PES [491 reported earlier from our laboratory [50]. The Laplacian of the wavefunction was evaluated using the Cooley—Tukey FFT algorithm [51]. The time evolution of ‘I’ was followed using the second-order finite-difference method [2]: 2, (9) [2i LitH/h]~”
as a function of time as illustrated in fig. 1 with the above integral being evaluated numerically using extended Simpson’s rule [52] at frequent time intervals. The reaction (R2) is substantially endothermic and therefore is expected to have the saddle point in the exit channel. It turns out that the ab-initio PES for the collinear reaction has no saddle point in the (q1, q2) space under investigalion. Therefore we chose a dividing line at q2 3.6 at a q au such to separate the products fromtothe 2 value the barrier thereactants reaction as is =
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepacket method
212
12
0.8
—
C 0
(ER) 0.0
1
2
3
4
Ui
04 ______________________________________ 00 0~6
NT/bOO
E
11009 Cay)
Fig. 1. Vanation of
Fig. 2. Translational energy distribution for the four GWPs used in the study.
where
dip in
f( ~
=
X exP[
I
(
j~
i/2~
ABC ‘IT
—46
6
~
~
)
j
1E1/2 1ran~
t2~ tran.r o 2~A.BC\(E1/2 _E1/2)2
(13) nearly overcome. That such a choice of dividing line was reasonable was apparent from contour plots of I ~Itj 2 in (q 1, q2) space at different time intervals. There is a clear separation in the I ‘P 2 plot in the neighborhood of q2 3.6 au. Shifting the dividing line to a larger value of q2 in the neighborhood had only a marginal effect on the (PR) computed. As was shown in an earlier pub-
The resulting (pR) 0.241 is in excellent accord with the (PR) obtained from the TDQM calcula=
tions using a 128
X
128 and a 256 X 256 grid:
=
lication [26], the sum of informatioq theoretic entropies (S~+ S~)in the coordinate and momenturn spaces reached a maximum value during the reaction and levelled off indicating that we had time evolved the wavefunction for a sufficiently long duration. We kept a check on the accuracy and stability of the time evolution by ensuring that the norm of the wavefunction was conserved. In addition, we computed (PR) by taking an average of the timeindependent quantal PR(Etrans) results [45,46]over the Eirans distribution for the GWP centered at 1.0 eV: (PR)
=
R
( Etrans )f( Eirans)
d Etrans,
(12)
-
He +H;(v,=o)_. OCT RWP KGwP) GWP
HeH + H
0 ..
0.8
-
x.•0
0.5 ,-
0.4-
‘•
•,
0.0 0.0
0.8 1
2
1.2 3
4
Etrans Cay) Fig. 3
P’~ values computed from TDQM calculations corn-
pared with the time-independent quantal results shown in the inset along with the QCT results over a wide range of Eirans.
values obtained directly from the GWPs are indicated as results obtained from eqs. (17)—(19)
(~)
are shown as (2<) GWP. Results from RWP (-) calculations include the error bars which indicate the L~Etrans range for which they are valid.
N. Balakrishnan, N. Sathyamurthy
0.258 and 0.251 respectively showing that the GWP method gives, on the average, an accurate description of the reaction dynamics. We did additional calculations using GWPs with (Etrans) 2.0, 3.0 and 4.0 eV, the distributions for which are given in fig. 2. There has been no other quantum calculation reported in this energy range. Therefore the (pR) values obtained from the TDQM calcula-. tions are compared with the available QCT results in fig. 3. Although they do not agree with each other quantitatively, they depict the same trend in PR(Etrans) and converge towards each other at high energy as would be expected. Different workers [39,53] have shown that state-to-state inelastic (I) and reaction probabilities could be obtained by projecting the wavefunction at the end of the time evolution onto the different reactant and product vibrational states. Since the details are given elsewhere [2 39 53] we only mention that in the absence of any dissociation =
PR(Etrans)
=
1
—
P’(Etrans),
(14)
where P’(Etrans)
Pon(Etrans)
(15)
~Pon(Etrans),
=
=
k \ k_)(2~IF(k)I2)
..~
/
Time-dependent wavepacket method
213
Alternatively, we could adopt any of the procedures (see for example, ref. [54]) used to obtain the excitation function from rate constant data. Since we have only limited (PR) information we adopted the following strategy. The threshold for the exchange reaction (R2) is 0.81 eV. Therefore we let x E 0 81 PR(E ) ax2 + bx =
—
=
trans
\
trans
‘
17 This is equivalent to expanding (PR) in terms of “moments” of Etrans upto and including the second: (pR)
=
J
(ax2
+
bx )f( Eirans) d Eirans.
(18)
By utilizing the additional information from the QCT calculations that pR 1 for Etrans 4.0 eV, =
2
a(4 —0.81) + b(4 —0.81) 1. (19) Using eqs. (18) and (19) we obtained the values of a and b and the resulting (PR) values in the neighborhood of 1 eV are shown in fig. 3. They fall within the range of the time-independent quantal results as shown in the inset. That is, although the method is not able to reproduce the oscillations due to resonance, it yields reliable results, on the average. =
2.2. Rectangular wavepacket (R WP) method
x
ff dq~dq2x~(q2)exp(ik~q1)’P(q1,q2)
.
(16) F(k) is the momentum distribution corresponding to the initial wavepacket and k0 and k~are the wavevectors corresponding to E~ rans for the initial and final states. Our attempts to use the above approach in determining PR(Etrans) were not successful at low (Etrans) as a substantial portion of the WP lingered around the interaction region for a long time making the deconvolution meaningless. With increase in (Etrans>, the problem became less severe. At the highest (Etr~s)the final wavefunc.tion was almost completely in the exit channel suggesting a pR 1, in excellent accord with the QCT results. —
In order to circumvent the problem of deconvolution, Agrawal and Raff [40] had proposed the use of the following wavepacket:
~I~’ ( q1)
ik0q1) sin Lik ~ q10) (q1 which has a rectangular distribution momentum space: =
exp(
—
—
(20).
—
.
F(k)
.
=
.
in the
exp(ikq10) V2 Lik k~ Lik < k k0 + Lik or k
—
~21 and \
‘Etrans!
‘~ =
2 ko
+
~~
\
2
/
h2 P’A.BC
.
22
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepackel method
214
In principle one could use a very small value of Lik to start with a delta-function distribution for the wavepacket in the momentum space and extract the reaction probability in a very narrow range of energy. But, as mentioned earlier, this will require an unrealistically large grid in the q1 direction. In the present set of calculations we have used Lik 2.0, for which the wavepacket was adequately represented in a 128 X 128 spatial grid. Other parameters used for time-evolving the RWPs were the same as those used for the GWPs. But we must reiterate that the use of the RWP method to compute pR at a narrow Eirans range calls for the use of a pair of wavepackets one centered at k0 k1 with Lik Lik1 and another at k0 k2 k1 + with Lik Lik2 Lik~+ e. If P1 and P2 are the corresponding average reaction probabilities, pR ai
(b) 8
(3> (2> (1>
-
(4)
6
-
=
4.
(a) 0.8
~
/ ~
-
(2>
—
=
=
=
=
=
=
-
~ “
2/2~ABC
0.4
-
0.0
1
Etrans=h2(kl+Likl+() NT /1000 LiEtrans
=
h2(k 1
+
Lik1
(23)
+ c)/~LA BC
is obtained from DR( I;’ ‘
~.‘-‘trans ±
A
77
Fig. 4. (a)
\ trans/
— —
P2 Lik~—
Lik2
—
P1 LikLik1 i
‘c” I
Time evolution of the RWPs was followed by monitoring the variation of (PR) with t. Similar to the behaviour for the GWPs, S + for the RWPs also increased with t, reached a maximum and levelled off as shown in fig. 4, indicating that the time evolution was complete. Values of pR obtained at Etran. 0.95, 1.0, 1.05, 1.33, 1.38, 2.67 and 3.56 eV are included in fig. 3. They fall within the range of the time-independent quantal results shown in the inset. They exhibit some structure around 1 eV but they do not reproduce the resonances presumably due to the choice of the Lik value. A further reduction in Lik would require the use of a larger grid. Clearly the RWP method is not computationally viable for computing pR5 when there are very narrow resonances as is the case for the system under investigation. Although the rectangular wavepacket results are nearly in quantitative agreement with the QCT results at higher Etrans, calculations using an RWF centered S
=
ingRWPs.
(IA\
‘
at
4 eV yielded a
4.64 eV
—
pR
greater than unity at
Etrans
=
clearly a non-physical result.
3. Summary and conclusion Although TDQM calculations using GWPs are not suitable to predict narrow resonances, our present calculations show that, on the average, they yield pR~5 in reasonable accord with the available time-independent quantal results at low energies and with the QCT results over the entire Etrans range under investigation. The RWP method yields pR values in better agreement with the QCT calculations. However, despite its success in reproducing resonances in collinear (H, H2) collisions, it proves to be too time-consuming and of no practical use in reproducing the narrow resonances in (He, H~) collisions. The only other method that could possibly do the job is the projection operator algorithm of Neuhauser and Baer which allows for a larger one-dimensional
N. Balakrishnan, N. Sathyamurthy
grid along the approach coordinate. We are in the process of testing its utility for the (He, H~) system. Acknowledgements We would like to thank Dr. Kulander for his kind invitation to contribute this article. We alsothank him and the anonymous referee for their comments on an earlier version of the manuscript. This work was supported in part by a grant from the INDO—US subcommission.
/
Time-dependent wavepacket method
[241D.
215
Neuhauser and M. Baer, J. Phys. Chem. 93 (1989)
2872; J. Chem. Phys. 91(1989) 4651.
[25] D. Neuhauser, Baer, R.S. Judson and D.J. Kouri, J. Chem. Phys. 90 M. (1989) 5882. [261 N. Balakrishnan and N. Sathyarnurthy, Chem. Phys. Lett. 164 (1989) 267. [27] S.E. Bradforth, A. Weaver, D.W. Arnold, RB. Metz and D.M. Neumark, J. Chem. Phys., in press. [28] A.E. Orel and K.C. Kulander, J. Chem. Phys. 91(1989) 6086. [29] D. Kosloff and R. Kosloff, J. Comput. Phys. 52 (1983) 35. [30] A. Askar and AS. Cakmak, J. Chem. Phys. 68 (1978) 2794. [31] Harmuth, J. Math, and70Phys. 36 4811. (1957) 269. [32] HF. R.J. Rubin, J. Chem. Phys. (1979)
[33] RE. Mickens, Phys. Rev. A 39 (1989) 5508.
References
[34] 47 M.D. Feit,412. J.A. Fleck, Jr. and A. Steiger, J. Comput. Phys. (1982)
[1] R.B. Gerber, R. Kosloff and M. Berman, Comput. Phys. Rep. 5 (1986) 59. [2] V. Mohan and N. Sathyamurthy, Comput. Phys. Rep. 7 (1988) 213. [3] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [4] D.J. Tannor and S.A. Rice, Adv. Chem. Phys. 70 (1988) 44. [5] T. Joseph, J. Manz, V. Mohan and H.-J. Schreier, Ber. Bunsenges. Phys. Chem. 92 (1988) 397. [6] N.E. Henriksen, J. Zhang and D.G. Imre, J. Chem. Phys. 89 (1988) 5607. [7] S.O. Williams and D.G. Imre, J. Phys. Chem. 92 (1988) 6636. [8] A.E. Orel and K.C. Kulander, Chem. Phys. Lett. 146 (1988) 428. [91J. Zhang and D.G. Imre, J. Chem. Phys. 90 (1989) 1666. [10] M. Messina and R.D. Coalson, J. Chem. Phys. 90 (1989) 4015. [11] R. Heather and H. Metiu, J. Chem. Phys. 90 (1989) 6903. [121 M. Founargiotakis, S.C. Farantos, G. Contopoulos and C. Polymilis, J. Chem. Phys. 91 (1989) 1389. [131 M. Jacon, 0. Atabek and C. Leforestier, J. Chem. Phys. 91 (1989) 1585. [14] V. Engel and H. Metiu, J. Chem. Phys. 91(1989)1596. [15] D.G. Imre and J. Zhang, Chem. Phys. 139 (1989) 89. [161 B. Hartke, J. Manz and J. Maths, Chem. Phys. 139 (1989) 123. [171 R. Kosloff, S.A. Rice, P. Gaspard, S. Tersigni and D.J. Tannor, Chem. Phys. 139 (1989) 201. [18] S.-Y. Lee, W.T. Pollard and R.A. Mathies, Chem. Phys. Lett. 160 (1989) 531. [191 B. Hartke, Chem. Phys. Lett. 160 (1989) 538. [20] J. Zhang, D.G. Imre and J.H. Frederick, J. Phys. Chem. 93 (1989) 1840. [21] Y. Murayama, Phys. Lett. A 140 (1989) 469. [221 S. Thareja and N. Sathyamurthy, J. Indian Chem. Soc. ~ (1989) 596. [23] P. Pernot, R.M. Grimes, W.A. Lester Jr. and C. Cedan, Chem. Phys. Lett. 163 (1989) 297.
[351H.
Tal-Ezer and R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [36] P.K. Chattaraj, S. Rao Koneru and B.M. Deb, J. Comput. Phys. 72 (1987) 504. [37] H. Tal-Ezer, R. Kosloff and C. Cerjan, J. Comput. Phys., to be published. [38] C. Leforestier, R. Bisseling, C. Cerjan, M.D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, 0. Roncero and R. Kosloff, J. Comput. Phys., in press. [39] C. Leforestier, Chem. Phys. 87 (1984) 241. [40] P.M. Agrawal and L.M. Raff, J. Chem. Phys. 74 (1981) 5076. [41] R. Kosloff and D. Kosloff, J. Comput. Phys. 63 (1986) 363. [42] D. Neuhauser and M. Baer, J. Chem. Phys. 90 (1989) 4351. [43] R. Heather and H. Metiu, J. Chem. Phys. 86 (1987) 5009. [44] J.Z.H. Zhang, Chem. Phys. Lett. 160 (1989) 417. [451T. Joseph and N. Sathyamurthy, J. Indian Chem. Soc. 62 (1985) 874. [46] N. Sathyamurthy, M. Baer and T. Joseph, Chem. Phys. 114 (1987) 73. [471J.E. Dove, ME. Mandy, V. Mohan and N. Sathyamurthy, J. Chem. Phys. 92 (1990) 7373. [48] V. Balasubramaman, B.K. Mishra and N. Sathyamurthy, J. Chem. Phys., submitted. [49] DR. McLaughlin and D.L. Thompson, J. Chem. Phys. 70 (1979) 2748. [501 T. Joseph and N. Sathyamurthy, J. Chem. Phys. 86 (1987) 704. [511 J.W. Cooley and J.W. Tukey, Math. Comput. 19 (1965) 297. [521 M. Abramowitz and l.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1968) p. 886. [53] (a) E.J. Heller, J. Chem. Phys. 62 (1975) 1544. (b) K.C. Kulander, J. Chem. Phys. 69 (1978) 5064; Nucl. Phys. A 353 (1981) 341c. [541A. Menon and N. Sathyamurthy, J. Phys. Chem. 85 (1981) 1021.
209
Time-dependent wavepacket methods for the calculation of collinear atom—diatom exchange reaction probabilities N. Balakrishnan and N. Sathyamurthy Department of Chemistry, Indian Institute of Technology, Kanpur 208 016, India Received 7 february 1990
Different strategies available for extracting the translational energy dependence of state-selected reaction probabilities for a collinear atom—diatom exchange reaction are reviewed and the usefulness of some of them is examined using collinear He+ H~(v~ = 0) -~ HeH~+ H as a test case.
1. Introduction The time-dependent quantum mechanical (TDQM) wavepacket method has emerged as a powerful tool for the study of atomic and molecular collision processes in recent years. Its success emanates from its simplicity, ease of interpretation of the results and the detailed physical insight into the dynamics obtained. One can view snapshots of the wavepacket at every stage of time evolution and get a “feel” for the reaction, analogous to the classical description of the process but within the quantal framework. With the introduction of accurate algorithms for the propagation of wavepackets in space as well as time, the method could yield results of accuracy comparable with that of the time-independent coupled channel methods. The TDQM method has beeti applied successfully also to the study of eigenvalue problems, photodissociation, Raman scattering, modeselectivity in chemical reactions, gas—surface scattering etc. The method, its implementation and its application till 1987 have been reviewed elsewhere [1—3].For more recent applications, see refs. [4—28]. We refer to the TDQM method to mean that the time-dependent Schrodinger equation (TDSE) INSA research fellow (1989—1991). 0010-4655/91/$03.50 © 1991
—
is solved numerically for the problem under consideration; it is also referred to as the grid method. It consists of three parts: (i) choosing the initial wavefunction for the system at time t~,‘I’(t~),(ii) propagating it in space and time and (iii) analyzing the wavefunction, ~P(t), after time t to extract reaction attributes. Major developments in the methodology in the 1980’s concerned primarily part (ii): Kosloff and Kosloff [29] pointed out thai the Laplacian of the wavefunction could be evaluated using fast Fourier transform (FFT) aI~ gorithms; Askar and Cakmak [30] had earliet pointed out that the time-evolution could be car ried out by a second-order difference scheme pro. vided that the time increment met certain stabiliti criteria; its equivalence to the explicit metho (also called the Richardson scheme) [31] has beei discussed elsewhere [2,32,33]. Feit et al. [34] dem onstrated the utility of a split-operator time prop agation scheme. Tal-Ezer and Kosloff [35] showed that the time-evolution could be carried to mac hine accuracy by expanding the evolution operatc in terms of complex Chebyschev polynomiali Chattaraj et al. [36] proposed an unconditionall stable matrix inversion scheme but it turned out t be computer intensive. Recently Tal-Ezer et a [37] have introduced a low-order polynomial a~ proximation of propagators, particularly when tF Hamiltonian is time-dependent. As an outcome
Elsevier Science Publishers B.V. (North-Holland)
N. Balakrishnan, N. Saihyamurthy / Time-dependent wavepacket method
210
a CECAM workshop, Leforestier et al. [38] have carried out a critical analysis of the accuracy and stability of the different algorithms used for ternporal evolution, Regardless of the method used for the propagation of the wavefunction in space and time, there are some practical difficulties. For example, considering one of the simplest reactions possible, an atom—diatom exchange reaction constrained to a line, A
+
BC(v,)
—‘
AB
+
C,
(Ri)
the initial wavefunction would be, typically, a product of the vibrational wavefunction corresponding to the initial state v, of the diatom and a translational wavepacket, often a minimum uncertainty Gaussian. That is, there is a Gaussian distribution in the momentum space corresponding to the relative translational motion of the reactants. Choice of such a wavepacket becomes necessary as otherwise we have to start with a delta function in the momentum space in order to compute the reaction probability, pR, at a fixed Etrans. This would mean that we use an infinitely spread out wavepacket in the configuration space and thus an infinitely large grid. In practice, the dependence of pR Ofl Etrans is extracted from the wavefunction at the end of the time evolution by a deconvolution procedure [39]. Agrawal and Raff [40] had proposed the use of rectangular wavepackets to compute P’~at a given Etrans. Unfortunately, the method calls for time evol4ving a pair of wavepackets centered at Etrans at which the value of pR is sought and if the energy dependence of P’~is required, the calculations have to be repeated again and again at different Etrans. Therefore it is computationally not appealing. Yet another difficulty of these TDQM calculations is that, with increasing time, the wavepacket spreads over the reactant and product spaces and before the end of the calculation it gets reflected from the edges of the gnd causing artificial interference effects. One way to avoid such spurious effects is to increase the size of the grid which unfortunately increases the computational time. Some authors [41,42] have circumvented the prob1cm by using absorbing boundaries but this im-
plies loss of information on the product (reactant) states. Heather and Metiu [43] proposed a wavefunction splitting algorithm to keep the grid extended barely beyond the interaction region. Recently Zhang [44] has shown that by switching over to an interaction representation, finite boundary reflection is eliminated, with concomitant reduction in the number of grid points and increase in the time increment for the numerical integration. Neuhauser and Baer [24] have developed a novel projection operator formalism which eliminates the need for inclusion of the product channel. Their approach calls for a relatively small rectangular grid in the interaction region and a long one-dimensional grid in the approach coordinate thus reducing the computational time. It also provides for a choice of a wide wavepacket in the coordinate space enabling the computation of pR at a narrow range of Etrans. The authors demonstrated that they could reproduce the dip in PR(Etrans) curve for H + H 2(v, = 0). We have chosen to examine the usefulness of some of the available approaches using the collinear He + H~(v 0)—’ HeH~+H (R2) =
reaction as a test case since the recently reported [45,46] converged R-matrix calculations on an accurate ab-initio potential energy surface (PES) have revealed reactive scattering resonances in the energy range Etrans 0.82—1.35 eV. Also quasiclassical trajectory (QCT) calculations have been carried out [47,48] over a wide range of Etrans with pR varying from 0 to 1. In the following section we outline the different approaches we have used in computing PR(Etrans) by numerically solving the TDSE and also present the results obtained. A summary of our findings is given in section 3. =
2. Different strategies to compute pR ( Etrans) 2.1. Gaussian wavepacket (GWP) method
The basic philosophy of the grid method is to solve the TDSE ih
( \
~
=
~t
I
I-J’I’,
(1)
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepacket method
where H=
wherein the superscript refers to the number of 2
2
~
2ILA,BC) h
—
Bq~+ 0q
+ V(
-~—~
2
time steps. We have used Lit
qi q ‘
2)
(2)
in mass weighted coordinates for a collinear A BC system: 3rsc sin 9, q~ rAB + /
+
=
q2 $~‘BC cos with
(3)
=
mC(mA + mB) m~(m~+ mc)
1/2
(\ (m~+ mB)(mB mAmC + mc) /
=
In such a coordinate system ‘I’( t
=
(5)
0) is chosen as
0) ‘P(q1, q2) ~(q1)x(q2), (6) where ~(q 2) is vibrational the Morse oscillator wavefunction for the ground state of H~ and ~(q 1) is the minimum uncertainty Gaussian: =
=
=
=
=
2 ~q1,
~,
with 3 0.25 au. q10 refers to the centre of the wavepacket and k0 is related to its average translational energy,
—
~2 /4
—
‘~
xexp
~ )/
~
—
q’°~
—
2itCh/pAB~
46
/4
2 32 1A,BC
—1/4
I-
ik0q~
(10)
,
—
.
which contracts to a nMnimal Gaussian wavepacket after the dilation time t~ which was chosen to be (1000—2000) Lit, depending upon the grid size. The time evolution was followed by monitoring the reaction probability
=
f
—
~,
2\
h
2 34 !~A,BC 2
2
(q
1—q10) _iko~i) 46 (7)
—
=
=
=
~(qi)=(2~62)1/4exp(_
0.0269 fs. To ii-
=
(4)
and sin 9
=
tialize theusing temporal ‘I~was from ~I,0 eq. (1)evolution, by expanding the evaluated evolution operator in Taylor series in Lit and truncating after the second-order term. After preliminary investigations using a 64 X 64 grid in (q 1, q2) space with q~fl 1.0 and ~ 0.4 with Liq~ Liq2 0.1 au, we carried out the time evolution of ‘I’ at
=
211
product
I ~I’(qi~ q 2)
2
dq1 dq2
(11)
ILA,BC
where !~A,BC is the reduced mass of the reactants, The potential energy, V, for the system was taken to be the analytic fit of the McLaughlin—Thompson ab-initio PES [491 reported earlier from our laboratory [50]. The Laplacian of the wavefunction was evaluated using the Cooley—Tukey FFT algorithm [51]. The time evolution of ‘I’ was followed using the second-order finite-difference method [2]: 2, (9) [2i LitH/h]~”
as a function of time as illustrated in fig. 1 with the above integral being evaluated numerically using extended Simpson’s rule [52] at frequent time intervals. The reaction (R2) is substantially endothermic and therefore is expected to have the saddle point in the exit channel. It turns out that the ab-initio PES for the collinear reaction has no saddle point in the (q1, q2) space under investigalion. Therefore we chose a dividing line at q2 3.6 at a q au such to separate the products fromtothe 2 value the barrier thereactants reaction as is =
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepacket method
212
12
0.8
—
C 0
(ER) 0.0
1
2
3
4
Ui
04 ______________________________________ 00 0~6
NT/bOO
E
11009 Cay)
Fig. 1. Vanation of
Fig. 2. Translational energy distribution for the four GWPs used in the study.
where
dip in
f( ~
=
X exP[
I
(
j~
i/2~
ABC ‘IT
—46
6
~
~
)
j
1E1/2 1ran~
t2~ tran.r o 2~A.BC\(E1/2 _E1/2)2
(13) nearly overcome. That such a choice of dividing line was reasonable was apparent from contour plots of I ~Itj 2 in (q 1, q2) space at different time intervals. There is a clear separation in the I ‘P 2 plot in the neighborhood of q2 3.6 au. Shifting the dividing line to a larger value of q2 in the neighborhood had only a marginal effect on the (PR) computed. As was shown in an earlier pub-
The resulting (pR) 0.241 is in excellent accord with the (PR) obtained from the TDQM calcula=
tions using a 128
X
128 and a 256 X 256 grid:
=
lication [26], the sum of informatioq theoretic entropies (S~+ S~)in the coordinate and momenturn spaces reached a maximum value during the reaction and levelled off indicating that we had time evolved the wavefunction for a sufficiently long duration. We kept a check on the accuracy and stability of the time evolution by ensuring that the norm of the wavefunction was conserved. In addition, we computed (PR) by taking an average of the timeindependent quantal PR(Etrans) results [45,46]over the Eirans distribution for the GWP centered at 1.0 eV: (PR)
=
R
( Etrans )f( Eirans)
d Etrans,
(12)
-
He +H;(v,=o)_. OCT RWP KGwP) GWP
HeH + H
0 ..
0.8
-
x.•0
0.5 ,-
0.4-
‘•
•,
0.0 0.0
0.8 1
2
1.2 3
4
Etrans Cay) Fig. 3
P’~ values computed from TDQM calculations corn-
pared with the time-independent quantal results shown in the inset along with the QCT results over a wide range of Eirans.
values obtained directly from the GWPs are indicated as results obtained from eqs. (17)—(19)
(~)
are shown as (2<) GWP. Results from RWP (-) calculations include the error bars which indicate the L~Etrans range for which they are valid.
N. Balakrishnan, N. Sathyamurthy
0.258 and 0.251 respectively showing that the GWP method gives, on the average, an accurate description of the reaction dynamics. We did additional calculations using GWPs with (Etrans) 2.0, 3.0 and 4.0 eV, the distributions for which are given in fig. 2. There has been no other quantum calculation reported in this energy range. Therefore the (pR) values obtained from the TDQM calcula-. tions are compared with the available QCT results in fig. 3. Although they do not agree with each other quantitatively, they depict the same trend in PR(Etrans) and converge towards each other at high energy as would be expected. Different workers [39,53] have shown that state-to-state inelastic (I) and reaction probabilities could be obtained by projecting the wavefunction at the end of the time evolution onto the different reactant and product vibrational states. Since the details are given elsewhere [2 39 53] we only mention that in the absence of any dissociation =
PR(Etrans)
=
1
—
P’(Etrans),
(14)
where P’(Etrans)
Pon(Etrans)
(15)
~Pon(Etrans),
=
=
k \ k_)(2~IF(k)I2)
..~
/
Time-dependent wavepacket method
213
Alternatively, we could adopt any of the procedures (see for example, ref. [54]) used to obtain the excitation function from rate constant data. Since we have only limited (PR) information we adopted the following strategy. The threshold for the exchange reaction (R2) is 0.81 eV. Therefore we let x E 0 81 PR(E ) ax2 + bx =
—
=
trans
\
trans
‘
17 This is equivalent to expanding (PR) in terms of “moments” of Etrans upto and including the second: (pR)
=
J
(ax2
+
bx )f( Eirans) d Eirans.
(18)
By utilizing the additional information from the QCT calculations that pR 1 for Etrans 4.0 eV, =
2
a(4 —0.81) + b(4 —0.81) 1. (19) Using eqs. (18) and (19) we obtained the values of a and b and the resulting (PR) values in the neighborhood of 1 eV are shown in fig. 3. They fall within the range of the time-independent quantal results as shown in the inset. That is, although the method is not able to reproduce the oscillations due to resonance, it yields reliable results, on the average. =
2.2. Rectangular wavepacket (R WP) method
x
ff dq~dq2x~(q2)exp(ik~q1)’P(q1,q2)
.
(16) F(k) is the momentum distribution corresponding to the initial wavepacket and k0 and k~are the wavevectors corresponding to E~ rans for the initial and final states. Our attempts to use the above approach in determining PR(Etrans) were not successful at low (Etrans) as a substantial portion of the WP lingered around the interaction region for a long time making the deconvolution meaningless. With increase in (Etrans>, the problem became less severe. At the highest (Etr~s)the final wavefunc.tion was almost completely in the exit channel suggesting a pR 1, in excellent accord with the QCT results. —
In order to circumvent the problem of deconvolution, Agrawal and Raff [40] had proposed the use of the following wavepacket:
~I~’ ( q1)
ik0q1) sin Lik ~ q10) (q1 which has a rectangular distribution momentum space: =
exp(
—
—
(20).
—
.
F(k)
.
=
.
in the
exp(ikq10) V2 Lik k~ Lik < k
—
~21 and \
‘Etrans!
‘~ =
2 ko
+
~~
\
2
/
h2 P’A.BC
.
22
N. Balakrishnan, N. Sathyamurthy / Time-dependent wavepackel method
214
In principle one could use a very small value of Lik to start with a delta-function distribution for the wavepacket in the momentum space and extract the reaction probability in a very narrow range of energy. But, as mentioned earlier, this will require an unrealistically large grid in the q1 direction. In the present set of calculations we have used Lik 2.0, for which the wavepacket was adequately represented in a 128 X 128 spatial grid. Other parameters used for time-evolving the RWPs were the same as those used for the GWPs. But we must reiterate that the use of the RWP method to compute pR at a narrow Eirans range calls for the use of a pair of wavepackets one centered at k0 k1 with Lik Lik1 and another at k0 k2 k1 + with Lik Lik2 Lik~+ e. If P1 and P2 are the corresponding average reaction probabilities, pR ai
(b) 8
(3> (2> (1>
-
(4)
6
-
=
4.
(a) 0.8
~
/ ~
-
(2>
—
=
=
=
=
=
=
-
~ “
2/2~ABC
0.4
-
0.0
1
Etrans=h2(kl+Likl+() NT /1000 LiEtrans
=
h2(k 1
+
Lik1
(23)
+ c)/~LA BC
is obtained from DR( I;’ ‘
~.‘-‘trans ±
A
77
Fig. 4. (a)
\ trans/
— —
P2 Lik~—
Lik2
—
P1 LikLik1 i
‘c” I
Time evolution of the RWPs was followed by monitoring the variation of (PR) with t. Similar to the behaviour for the GWPs, S + for the RWPs also increased with t, reached a maximum and levelled off as shown in fig. 4, indicating that the time evolution was complete. Values of pR obtained at Etran. 0.95, 1.0, 1.05, 1.33, 1.38, 2.67 and 3.56 eV are included in fig. 3. They fall within the range of the time-independent quantal results shown in the inset. They exhibit some structure around 1 eV but they do not reproduce the resonances presumably due to the choice of the Lik value. A further reduction in Lik would require the use of a larger grid. Clearly the RWP method is not computationally viable for computing pR5 when there are very narrow resonances as is the case for the system under investigation. Although the rectangular wavepacket results are nearly in quantitative agreement with the QCT results at higher Etrans, calculations using an RWF centered S
=
ingRWPs.
(IA\
‘
at
4 eV yielded a
4.64 eV
—
pR
greater than unity at
Etrans
=
clearly a non-physical result.
3. Summary and conclusion Although TDQM calculations using GWPs are not suitable to predict narrow resonances, our present calculations show that, on the average, they yield pR~5 in reasonable accord with the available time-independent quantal results at low energies and with the QCT results over the entire Etrans range under investigation. The RWP method yields pR values in better agreement with the QCT calculations. However, despite its success in reproducing resonances in collinear (H, H2) collisions, it proves to be too time-consuming and of no practical use in reproducing the narrow resonances in (He, H~) collisions. The only other method that could possibly do the job is the projection operator algorithm of Neuhauser and Baer which allows for a larger one-dimensional
N. Balakrishnan, N. Sathyamurthy
grid along the approach coordinate. We are in the process of testing its utility for the (He, H~) system. Acknowledgements We would like to thank Dr. Kulander for his kind invitation to contribute this article. We alsothank him and the anonymous referee for their comments on an earlier version of the manuscript. This work was supported in part by a grant from the INDO—US subcommission.
/
Time-dependent wavepacket method
[241D.
215
Neuhauser and M. Baer, J. Phys. Chem. 93 (1989)
2872; J. Chem. Phys. 91(1989) 4651.
[25] D. Neuhauser, Baer, R.S. Judson and D.J. Kouri, J. Chem. Phys. 90 M. (1989) 5882. [261 N. Balakrishnan and N. Sathyarnurthy, Chem. Phys. Lett. 164 (1989) 267. [27] S.E. Bradforth, A. Weaver, D.W. Arnold, RB. Metz and D.M. Neumark, J. Chem. Phys., in press. [28] A.E. Orel and K.C. Kulander, J. Chem. Phys. 91(1989) 6086. [29] D. Kosloff and R. Kosloff, J. Comput. Phys. 52 (1983) 35. [30] A. Askar and AS. Cakmak, J. Chem. Phys. 68 (1978) 2794. [31] Harmuth, J. Math, and70Phys. 36 4811. (1957) 269. [32] HF. R.J. Rubin, J. Chem. Phys. (1979)
[33] RE. Mickens, Phys. Rev. A 39 (1989) 5508.
References
[34] 47 M.D. Feit,412. J.A. Fleck, Jr. and A. Steiger, J. Comput. Phys. (1982)
[1] R.B. Gerber, R. Kosloff and M. Berman, Comput. Phys. Rep. 5 (1986) 59. [2] V. Mohan and N. Sathyamurthy, Comput. Phys. Rep. 7 (1988) 213. [3] R. Kosloff, J. Phys. Chem. 92 (1988) 2087. [4] D.J. Tannor and S.A. Rice, Adv. Chem. Phys. 70 (1988) 44. [5] T. Joseph, J. Manz, V. Mohan and H.-J. Schreier, Ber. Bunsenges. Phys. Chem. 92 (1988) 397. [6] N.E. Henriksen, J. Zhang and D.G. Imre, J. Chem. Phys. 89 (1988) 5607. [7] S.O. Williams and D.G. Imre, J. Phys. Chem. 92 (1988) 6636. [8] A.E. Orel and K.C. Kulander, Chem. Phys. Lett. 146 (1988) 428. [91J. Zhang and D.G. Imre, J. Chem. Phys. 90 (1989) 1666. [10] M. Messina and R.D. Coalson, J. Chem. Phys. 90 (1989) 4015. [11] R. Heather and H. Metiu, J. Chem. Phys. 90 (1989) 6903. [121 M. Founargiotakis, S.C. Farantos, G. Contopoulos and C. Polymilis, J. Chem. Phys. 91 (1989) 1389. [131 M. Jacon, 0. Atabek and C. Leforestier, J. Chem. Phys. 91 (1989) 1585. [14] V. Engel and H. Metiu, J. Chem. Phys. 91(1989)1596. [15] D.G. Imre and J. Zhang, Chem. Phys. 139 (1989) 89. [161 B. Hartke, J. Manz and J. Maths, Chem. Phys. 139 (1989) 123. [171 R. Kosloff, S.A. Rice, P. Gaspard, S. Tersigni and D.J. Tannor, Chem. Phys. 139 (1989) 201. [18] S.-Y. Lee, W.T. Pollard and R.A. Mathies, Chem. Phys. Lett. 160 (1989) 531. [191 B. Hartke, Chem. Phys. Lett. 160 (1989) 538. [20] J. Zhang, D.G. Imre and J.H. Frederick, J. Phys. Chem. 93 (1989) 1840. [21] Y. Murayama, Phys. Lett. A 140 (1989) 469. [221 S. Thareja and N. Sathyamurthy, J. Indian Chem. Soc. ~ (1989) 596. [23] P. Pernot, R.M. Grimes, W.A. Lester Jr. and C. Cedan, Chem. Phys. Lett. 163 (1989) 297.
[351H.
Tal-Ezer and R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [36] P.K. Chattaraj, S. Rao Koneru and B.M. Deb, J. Comput. Phys. 72 (1987) 504. [37] H. Tal-Ezer, R. Kosloff and C. Cerjan, J. Comput. Phys., to be published. [38] C. Leforestier, R. Bisseling, C. Cerjan, M.D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, 0. Roncero and R. Kosloff, J. Comput. Phys., in press. [39] C. Leforestier, Chem. Phys. 87 (1984) 241. [40] P.M. Agrawal and L.M. Raff, J. Chem. Phys. 74 (1981) 5076. [41] R. Kosloff and D. Kosloff, J. Comput. Phys. 63 (1986) 363. [42] D. Neuhauser and M. Baer, J. Chem. Phys. 90 (1989) 4351. [43] R. Heather and H. Metiu, J. Chem. Phys. 86 (1987) 5009. [44] J.Z.H. Zhang, Chem. Phys. Lett. 160 (1989) 417. [451T. Joseph and N. Sathyamurthy, J. Indian Chem. Soc. 62 (1985) 874. [46] N. Sathyamurthy, M. Baer and T. Joseph, Chem. Phys. 114 (1987) 73. [471J.E. Dove, ME. Mandy, V. Mohan and N. Sathyamurthy, J. Chem. Phys. 92 (1990) 7373. [48] V. Balasubramaman, B.K. Mishra and N. Sathyamurthy, J. Chem. Phys., submitted. [49] DR. McLaughlin and D.L. Thompson, J. Chem. Phys. 70 (1979) 2748. [501 T. Joseph and N. Sathyamurthy, J. Chem. Phys. 86 (1987) 704. [511 J.W. Cooley and J.W. Tukey, Math. Comput. 19 (1965) 297. [521 M. Abramowitz and l.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1968) p. 886. [53] (a) E.J. Heller, J. Chem. Phys. 62 (1975) 1544. (b) K.C. Kulander, J. Chem. Phys. 69 (1978) 5064; Nucl. Phys. A 353 (1981) 341c. [541A. Menon and N. Sathyamurthy, J. Phys. Chem. 85 (1981) 1021.