03964.com

文档资料库 文档搜索专家

文档资料库 文档搜索专家

View Online

Negative imaginary potentials in time-dependent quantum molecular scattering

Susanta Mahapatra and Narayanasami Sathyamurthy¤ Department of Chemistry, Indian Institute of T echnology, Kanpur 208 016, India

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

Re?ection or wrap around of the wavefunction from the grid edges is often avoided in time-dependent quantum mechanical calculations by using a negative imaginary potential (NIP) near the grid edges. The stability of the various (second-order di?erencing, split operator, Chebyshev polynomial and short iterative Lanczos) schemes used, in conjunction with the NIP, for time evolution is discussed using collinear (He, H `) collisions as a test case. It is shown that the difficulties encountered in obtaining 2 converged reaction probability [PR(E)] values at high energies for the system when NIPs are used, are avoided by using a properly chosen damping function externally.

1 Introduction

In time-dependent quantum mechanics the wavefunction W at time t is obtained by solving the time-dependent Schrodinger ? equation (TDSE). Numerically this is accomplished, in what is called the grid method, by representing the initial state of the system by W(0) on a discrete grid in coordinate space and following its evolution by slicing time into small intervals. Often one switches to the energy space to extract the desired dynamical attributes from the time-evolved wave packet (WP). Various aspects of this problem have been discussed.1h11 One serious problem in time-dependent quantum mechanical (TDQM) calculations is that, as time progresses, the WP reaches the grid edges and undergoes spurious re?ections or wraps around, depending on the choice of method used in evaluating the action of the kinetic energy operator on W. One way to avoid this problem is to use a very large grid which would delay the WP reaching the boundaries. However, this would mean a tremendous increase in computer memory and time requirements and that is not feasible in many practical applications. Leforestier and Wyatt12 used the Saxon?Woods potential13 W (R) \ [iV 0 1 ] exp[a(R [ R*)] (1)

in dampening the wavefunction before it reaches the grid edges. X and X de?ne the range in which the NIP is oper1I 2I ative. They also derived16 the criteria for selecting the optimum height (V ) and width (*X \ X [ X ) of the NIP I0 I 2I 1I as follows : *X J(8k)E3@2 ?E1@2 I t t OV O (3) I0 ? *X J(8k) I E in eqn. (3) represents the translational energy of the WP. t Similar NIPs have been considered by Child,18 Seideman and Miller19 and Vibok and Balint-Kurti.20 Monnerville et al.21 ? have made use of NIPs in their time-dependent reactive scattering studies. Zhang et al.22 have proposed a new scheme while obtaining the ?nal product-state distribution from a time-dependent wave packet calculation in the interaction representation. They have used a cut-o? function (F \ exp cut [[iV (x)]) derived from the NIP in eqn. (2), while propagating I the wavefunction. Such an approach has the advantage that it propagates the wave packet in the ?eld of a real potential and it preserves the Hermitian property of the Hamiltonian. More recently, they have applied the same method in their other time-dependent reactive scattering studies.23 Macias et al.24 have described a systematic inversion technique to optimize di?erent forms of NIP. In spite of their widespread use in recent years, the use of NIPs is known to cause problems in some of the numerical time-evolution schemes, since the Hamiltonian becomes nonHermitian in their presence. An alternative is to use real damping functions which retain the Hermitian property of the Hamiltonian. A comparison of the second-order di?erencing (SOD), split operator (SO), Chebyshev polynomial (CP) expansion and short iterative Lanczos (SIL) schemes for numerical time evolution in terms of their stability has been made by Leforestier et al.25 and also by Truong et al.26 For the sake of completeness we will review them brie?y here and then analyse their applicability when used with an NIP. In Section 2 we will describe the general set up of the spatial grid that will be utilised in the discussion in the rest of this paper. The time evolution of the WP is described in Section 3. Section 3A describes the properties of the timeevolution operator and Sections 3B?E discuss the stability of the four numerical time-evolution schemes (vide supra) that are commonly used, in the presence of an NIP. The reaction probability results are presented and discussed in Section 4 for the test case of collinear (He, H `) collisions and our conclu2 sions are presented in Section 5. J. Chem. Soc., Faraday T rans., 1997, 93(5), 773?779 773

in addition to the real potential of the system, while investigating multi-photon dissociation in diatomic molecules. Here R is the spatial variable, R* the point at which the negative imaginary potential is activated and V and a are parameters 0 that de?ne the maximum height and slope of the potential variation. Koslo? and Koslo?14 derived the absorbing boundary condition of the propagating waves for the TDSE and the acoustic wave equation and concluded that the WP is absorbed near the grid boundaries when NIPs are added to the Hamiltonian. Koslo? and Cerjan15 had made use of NIPs while investigating surface desorption phenomena. Neuhauser et al.16,17 demonstrated the usefulness of NIPs of the form : V (X) \ I

7

[iV I0 0

C

X[X 1I X [X 2I 1I

D

if X O X O X 1I 2I otherwise

(2)

¤ Honorary Professor, S. N. Bose National Center for Basic Sciences, Calcutta, India.

View Online

2 Spatial grid

Let us assume that the coordinate space is discretized into a set of N points. If the spacing between two successive points is ü *x, then the eigenvalue of the position operator x at each point on the grid is given by x \ (i [ 1)*x ; i \ 1, . . . , N (4) i The corresponding eigenvectors are denoted by o x T. The i orthogonality and the completeness relations on this discrete grid are given by : ; *xSx o x T \ d i j ij i and N I? \ ; o x [ *x \ x o (6) x i i i/1 where I? is the identity operator. The functions represented at x the grid points are given by Sx o /T \ /(x ) (7) i i The continuous normalization integral, /= /*(x)/(x) dx \ 1, ~= on this grid transforms to a discrete sum : N ; /*(x )/(x )*x \ 1 (8) i i i/1 The maximum length L of this grid along the spatial coordinate is, N*x. This length determines the spacing between two successive points in the reciprocal k space : *k \ 2n N*x (9) (5)

3 Time evolution of a wave packet

A Time-evolution operator

The solution of the TDSE is given by27 W(t) \ P exp [i/? ?

C P

T

H(t@ ) dt@ W(0) ?

D

(13)

0 where W(t) is the wavefunction at time t, P is the time? ordering operator and H is the Hamiltonian operator. For an ? explicitly time-independent Hamiltonian the above equation simpli?es to : W(t) \ exp[[iHt/?]W(0) ? (14)

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

The exponential operator on the right-hand side of the above equation forms a continuous dynamical group where time t is a parameter, and is known as the (time) evolution operator denoted by U(t, t ). For t \ 0, ? 0 0 U(t, t ) \ exp([iHt/?) ? ? (15) 0 In actual computation, t is sliced in smaller steps of length, *t, and the entire time evolution is accomplished through : Nt~1 ? (16) U(t) \ < U[(n ] 1)*t, n*t] ? n/0 where N is the total number of time steps and *t \ t/N . t t U(t, t ) is unitary : ? 0 UUs \ UsU \ 1 ? ? ? ? (17)

B

Second-order di?erencing scheme

In k space, the grid is centred at zero and all other points are distributed symmetrically on either side. If p (\?k ) repmax max resents the maximum momentum in the k space then the N points are distributed in the interval M[p , . . . , 0, . . . , p N. max max Hence the total length of the grid in this space is 2o p o. Now max *x can be written as : (10) ok o max Therefore, the total volume v of phase space concerned is given by \ Nh (11) max is decided by the number of grid points. The maximum energy represented on the grid through this discretization is given by E max \T ]V max max p2 \ max ] V max 2m ?2o k o2 max ] V \ max 2m n2?2 \ ; ]V (12) max 2m (*x )2 i i i/1 where m , *x are the mass and grid spacing along the channel i i coordinate i, respectively, and V is the maximum value of max the potential energy represented on the grid. The minimum energy E (\V ) is equal to the minimum value of the min min potential energy. Our subsequent discussion will follow upon this general set up of the discrete grid. 774 J. Chem. Soc., Faraday T rans., 1997, V ol. 93 v \ 2L p *x \ n

The time-evolution operator and its adjoint for a time step of length *t can be written as : U(*t) \ exp([iH*t/?) ? ? Us(*t) \ exp(iH*t/?) ? ? Therefore, W(t ] *t) \ U(*t)W(t) ? W(t [ *t) \ Us(*t)W(t) ? Using the above equations one can write W(t ] *t) [ W(t [ *t) \ [exp([iH*t/?) [ exp(iH*t/?)]W(t) ? ? (21) Expanding the exponential terms in the above equation in Taylor series and keeping only the terms up to second order in *t, the equation becomes W(t ] *t) \ W(t [ *t) [ 2iH*t ? W(t) ? (22) (19) (20) (18)

Thus the SOD scheme28,29 evaluates W(t ] *t) from its values at time t and t [ *t. Computationally, to initialize the iteration process, the value of W at the ?rst step is computed either by a fourth-order Runge?Kutta scheme or by using a higher-order Taylor series expansion of U(*t). The above ? equation can also be written as : W(t ] *t) \ exp(iH*t/?)W(t) [ ? \ QW(t) 2iH*t ? W(t) ? (23)

View Online

The scheme will be unitary if QQs \ QsQ \ 1. From the above it follows that ? QQs \ exp(iH*t/?) [ \1[4

C

DC A B A B A B

2iH*t ? ?

exp([iH*t/?) ] ?

2iH*t ? ?

D

that containing the kinetic-energy operator is sandwiched in between. The time evolution of the wavefunction in the potential-referenced SO scheme is given by W(t ] *t) \ exp([iT? *t/2?)exp([iV? *t/?)exp([iT? *t/2?)W(t) \ QW(t) Clearly (30)

*tH ? *tH ? *tH 2 ? sin ]4 ? ? ? (24)

*t4H4 \ 1 ] (4/6) ?4

QQs \ exp([iT? *t/2?)exp([iV? *t/?)exp([iT? *t/2?) ] exp(iT? *t/2?)exp(iV? *t/?)exp(iT? *t/2?) \ 1 (31)

neglecting higher-order terms in the expansion of the sine function. The conservation of the norm of the wavefunction can be checked through : SW(t ] *t) o W(t ] *t)T \ SQW(t) o QW(t)T \ SW(t) o QsQ o W(t)T

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

\ SW(t) o W(t)T ] (2/3)(*t/?)4 ] SH2W(t) o H2W(t)T ? ? (25)

and the scheme is unitary, and the norm of the wavefunction is conserved. However, because of the non-commutability of the kinetic- and the potential-energy operators the scheme does not conserve energy. The scheme is unconditionally stable and it does not depend on the kinetic-energy spectrum of the grid. However, the size of the time step is selected based on the maximum value of the potential energy on the grid. Good results are obtained when n? (32) 3*V max where *V (\V [ V ) is the maximum excursion of the max max min potential. This step size is so chosen as to accommodate the entire energy level spectrum of the system under investigation.31 In the presence of an NIP, Q@ becomes : *t \ Q@ \ exp([iT? *t/2?)exp([iV? *t/?) ] exp([V *t/?)exp( [ iT? *t/2?) 0 The unitary relation now becomes (33)

Therefore, it becomes clear that, by choosing *t sufficiently small, the unitarity can be maintained throughout the time evolution. It can be shown that29 the algorithm is stable if *t @ ?/E . In practice, a choice of *t \ 0.2?/E is found max max to yield good results.29 Since the Hamiltonian for the present case commutes with itself, [Q, H]W \ 0 and energy is conserved throughout the ? time evolution. If a negative imaginary potential (C \ [iV ) is added to the 0 real potential of the system, the Hamiltonian becomes complex (H@ \ H ] C) and Q@ can be written as ? ? Q@ \ exp[i*t(H [ iV )/?] [ ? 0 2i*t(H [ iV ) ? 0 ? (26)

Right-multiplying the above by Q@s and expanding the exponential terms in Taylor series and ?nally rearranging the resulting expression, one obtains : Q@Q@s \ exp(2V *t/?) [ 4 0 Hence, SW(t ] *t) o W(t ] *t)T

A B

V *t 0 ? *tV 0 ?

Q@Q@s \ exp([iT? *t/2?)exp([2V *t/?)exp(iT? *t/2?) (34) 0 In the limit of the commutation error becoming zero, the last two terms can be interchanged and the norm of the wavefunction at time t ] *t becomes : SW(t ] *t) o W(t ] *t)T \ SW(t) o Q@sQ@ o W(t)T \ SW(t) o exp([2V *t/?) o W(t)T (35) 0 It is clear that the norm decreases exponentially with increase in time.

(27)

\ SW(t) o exp(2V *t/?) [ 4 0

C

A BD

o W(t)T

(28)

D Chebyshev polynomial expansion scheme CPs are found to be superior to many other polynomials and are optimal for a scalar function F(x) bounded in the interval [[1, 1]. So, a scalar function such as exp(ax) can be expressed in terms of these polynomials in the interval [1 O x O 1 as = exp(ax) \ ; (2 [ d )J (a)T (x) (36) n0 n n n/0 where d is the Kronecker delta and a \ *E*t/2?. J (a) are n0 n the modi?ed Bessel functions of order n. T (x) are the CPs of n order n, calculated using the recursion relation32 T (x) \ 2xT (x) [ T (x) (37) n`1 n n~1 with T (x) \ 1 and T (x) \ x. 0 1 The evolution operator is a function of an operator. It has been shown7 that a function of an operator can be expressed as a function of a scalar in the complete basis of the operator. Hence, the function of the operator can be approximated in the Chebyshev series, provided the domain of the operator is con?ned to the interval [[1, 1] in which the CPs are optimal. In case of a Hamiltonian that is self-adjoint, the eigenvalues lie on a real axis, and they can be positioned from [1 to 1 by J. Chem. Soc., Faraday T rans., 1997, V ol. 93 775

It is clear that the error in the norm would grow exponentially. It is worth mentioning here that the use of an NIP shifts the eigenvalue spectrum to the complex plane and methods, which are conditionally stable, and which depend on the eigenvalue spectrum of the Hamiltonian become unstable under such conditions. C Split operator method

In the application of the evolution operator, there is an error arising from the fact that the kinetic- and potential-energy operators do not commute. However, by splitting the time evolution operator symmetrically the commutator error is reduced to third order. Such an approach is known as the second-order split-operator (SO) scheme.30 It can be either potential referenced or kinetic referenced. The potentialreferenced SO scheme is given by : exp([iH*t/?) \ exp([iT? *t/2?)exp([iV? *t/?) ? ] exp([iT? *t/2?) ] O(*t3) (29)

In the kinetic-referenced SO scheme the exponential containing the potential-energy operator is symmetrically split and

View Online

renormalizing the Hamiltonian as follows : H ? \2 norm H [ I? H ? 1 *E (38)

The evolution operator becomes N exp([iHt/?) \ exp([iHt/?) ; ([i)n(2 [ d )J (*Et/?)T (H ? 1 ? ) n0 n n norm n/0 (45) Mandelshtam et al.37 have made a simple analytic continuation of the CPs while keeping their properties unaltered. They have used an exponential damping factor in the de?nition of the CPs without disturbing the Hamiltonian. This preserves the renormalization step, as in the case of the Hermitian Hamiltonian in the presence of a NIP. They have obtained an expression for the time-evolution operator as : N U(t) \ ; a (t)Q (H ? ? ; cü ) (46) n n norm n/0 where a (t) \ [(2 [ d )exp([iHt/?)([1)nJ (*Et/?)] 1 and n n0 n Q (H ? , cü ) are the analytic continuation of the CPs satisfying n norm the following recursion relation : exp([cü )Q (H ? ; cü ) ] exp(cü )Q (H ? ; cü ) n~1 norm n`1 norm \ 2H ? Q (H ? ; cü ) (47) norm n norm with Q (H ? ; cü ) \ I? and Q (H ? ; cü ) \ exp([cü )H ? . The 0 norm 1 norm norm operator cü is dimensionless and it de?nes the damping factor. c is set to zero in the strong-interaction region and rises slowly as the asymptotic region is approached. In this prescription, the NIPs are used externally as damping functions. Therefore, their complex nature does not interfere with the eigenvalue spectrum of the Hamiltonian. Another important point to be mentioned here is that the CP expansion method is conditionally stable because the scheme will be unstable if the energy range of the Hamiltonian *E is underestimated. The scheme does not conserve norm and energy and the resulting deviation is used as a measure of the error.1,7 Very recently Neuhauser38 has investigated the applicability of the NIPs in TIQM scattering. He proposes anomalyfree very short-range imaginary potentials that cover only 1 or 2 grid points. His method involves the computation of the full S-matrix by diagonalization/inversion of the complex Hamiltonian. E Short iterative Lanczos method

where H \ (E ] E )/2 and *E \ (E [ E ). In terms 1 max min max min of this renormalized Hamiltonian, H ? , the evolution opernorm ator can be written as follows : exp([iH*t/?) \ exp([iH*t/?)exp([iaH ? ? ? ) (39) norm The ?rst term in the above expression is the phase shift due to the shift of the energy scale. The second term is approximated by the Chebyshev series as33 = ) \ ; (2 [ d ) J (a)U ([iH ? ) (40) norm n0 n n norm n/0 where U ([iH ? ) are the complex CPs of order n satisfying n norm the recursion relation : exp([iaH ?

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

\ [2iH ? U ]U (41) n`1 norm n n~1 where U \ 1 and U \ [iH . Therefore, the evolution of 0 1 norm W(t) in this scheme on a discrete grid is given by : U N W(t ] *t) \ exp([iH*t/?) ; (2 [ d ) J (a)U ([iH 1 ? )W(t) n0 n n norm n/0 (42) The number of terms to be used in the above expansion is estimated from the time?energy space volume a. In practice, the number of terms used is slightly larger than this estimate for a good convergence. Since the evolution operator is expanded in a series of polynomials in the Chebyshev method, the scheme is strictly not unitarity. The deviation from the unitary corresponds to the remainder term in the expansion. This deviation is used as an accuracy check of the scheme. The errors are uniformly distributed in the bounded interval.1,7 Since Bessel functions show exponential convergence for n [ a, the error is usually very small. The instability of the CP scheme in the presence of an NIP was ?rst noticed by Mowrey.34 It has been realised subsequently7,35h38 that the instability is caused by the NIP. The main point to note here is that the Hamiltonian ceases to be Hermitian once the NIP is added. The eigenvalue spectrum of the Hamiltonian shifts to the complex energy plane and the renormalization of the Hamiltonian introduced above becomes invalid. Koslo? and co-workers7,39 have proposed the use of Newton interpolating polynomials in which the interpolation points are predetermined on a boundary curve, which encloses the spectrum of the Hamiltonian. This is done by conformal mapping. For a Hermitian Hamiltonian the interpolation points become the zeros of the CP, which is not true for a non-Hermitian Hamiltonian. Huang et al.36 proposed the use of generalized Faber polynomials of which CPs constitute a special case. The latter are generated by a one-to-one conformal mapping of the exterior of a disk to the exterior of a simple closed curve, L . In the case of CPs, L becomes an c c ellipse. The renormalization of the Hamiltonian is done in such a way as to account for the shift of its eigenvalue spectrum to the complex plane in the non-Hermitian case. These authors used the renormalization H ? \2 norm H [ I? H ? 1 2c (43)

The SIL scheme, adapted to a numerical solution of the TDSE by Park and Light40 computes the action of the time evolution operator on the wavefunction by forming a reduced subspace (Krylov space) of the Hamiltonian matrix. The orthonormal basis set, q ( j \ 0, . . . , N [ 1) (Krylov basis set), j which is spatially and temporally tailored, is generated from the initial vector q \ t(0) as follows40h42 0 Hq \ a q ] b q ? (48) 0 0 0 0 1 Hq \ b ? q ]a q ]b q ; jP1 (49) j j~1 j~1 j j j j`1 with a \ Sq o H o q T ? j j j and \ Sq oHoq T ? (51) j~1 j~1 j The Hamiltonian matrix becomes tridiagonal43 in this reduced subspace : b (50)

where 2c P *E. The eigenvalues of H ? are no longer norm bounded by one but are mapped conformally from a unit disk. 2c equals the sum of the major and the minor axes. H in eqn. 1 (43) is given by H \ 1/2[E 1 776 max ]E ] [ iV min 0 (44)

H \ ? NL

a

a 0 b 0 0 < 0 0

b 0 a 1 b 1 < 0 0

0 b 1 a 2 < 0 0

... ... 0 ... b ... 2 } < ... a NL~2 ... b NL~2

0 0 0 <

b NL~2 a NL~1

b

(52)

J. Chem. Soc., Faraday T rans., 1997, V ol. 93

View Online

It is of size N ] N compared to N ] N (where N is the total L L number of grid points) of the original Hamiltonian. The Her? miticity of H is retained in H . H is diagonalized by the ? ? NL NL unitary transformation matrix Z, the column of which contains its eigenvectors : (53) H \ ZsD Z ? NL NL where D is the diagonal matrix of eigenvalues and Zs is the NL conjugate transpose of Z. The evolution operator can now be written as25 U(*t) \ exp([i*tH /?) ? ? NL \ Zs exp([i*tD /?)Z NL The time-evolved wavefunction then becomes

where the mass-scaled Jacobi coordinates R (corresponding to He, H ` translation) and r (corresponding to H ` vibration) 2 2 are used. F(R) is a minimum uncertainty Gaussian wave packet (GWP) : F(R) \ (2nd2)~1@4 exp [

C

(R [ R )2 0 [ ik R 0 4d2

D

(58)

(54)

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

W(t ] *t) \ Zs exp([i*tD /?)ZW(t) NL N~1 \ ; U(*t)q ? (55) n n/0 The summation here arises because the ?rst basis function in the Krylov space, q is identical to W(0). The above pro0 cedure is repeated successively N times. For small *t (ca. 0.1 fs), typically 5?10 basis functions are sufficient to construct the tridiagonal matrix. This method is particularly useful for small *t and since the scheme involves the exponential operation in the reduced subspace, it conserves norm, as in the case of the SO scheme. It conserves energy also as the exponential operator commutes with the Hamiltonian, unlike the SO scheme, which involves splitting the Hamiltonian into kinetic- and potential-energy parts. The scheme is unconditionally stable. In the presence of an NIP, the SIL scheme is stable. However, the basis vectors often become non-orthogonal and one needs to resort to (Gram?Schmidt) orthogonalization at regular intervals. The loss in orthogonality arises from the truncation error introduced in each time step. Error starts building up when the recursion vectors span outside the reduced space. Park and Light40 estimated this error by multiplying together the o?-diagonal elements of the tridiagonal matrix : ([i*t)NL~1 o q (*t) o B NL (N [ 1) ! L where o q (*t) o is the magnitude of NL outside the Krylov space. NL~1 <b (56) n n/0 the ?rst vector lying

Here d is the width parameter of the GWP, R and k denote 0 0 the location of its maximum in the coordinate and the momentum space respectively. / (r) are the vibrational eigenv functions of H ` molecule corresponding to its vibrational 2 state v of the extended Rydberg potential-energy curve,46 computed by the Fourier-grid Hamiltonian method.47 We have used a 256 ] 256 grid in (R, r) space originating at (2.2, 0.4) a and with a spacing *R \ *r \ 0.05 a . The 0 0 spatial evolution of the WP is carried out by the fast Fourier transform (FFT) method and the temporal evolution is carried out by the SO method. Since the SO method revealed an exponential damping of the WP norm in the presence of a NIP (Section IIIC) we use it in the present calculation in order to demonstrate the e?ect of NIP on the PR(E) of reaction (I). v The length of the time step *t is chosen to be 0.1616 fs. The PR(E) values are computed through : v LW(R, r , E) ? I (59) PR(E) \ Im W(R, r , E) I v Lr k

CT C

K

UD

where k \ [(m m m )/(m ] m ] m )]1@2 is the three-body He H H He H H reduced mass. In Fig. 1(a) we show the PR(E) values for reaction (I) (solid 0 curve) computed in presence of an NIP of the type48 V (X) \ I

7

[iV I0 0

X[X 2 1I X [X 2I 1I

D

if X O X O X 1I 2I (60) otherwise

4 Test case

In Section 3 we discussed the applicability of NIPs in di?erent time-evolution schemes. Here, we report the vibrational (v) state-selected reaction probability [PR(E)] values for the colv linear He ] H `(v \ 0) ] HeH` ] H (I) 2 reaction computed in the presence and in the absence of an NIP. The variation of the reaction probability with energy for reaction (I) has been shown44 to be highly oscillatory and the oscillations were identi?ed in terms of a large number of closely spaced narrow resonances arising from the quasibound states supported by the vibrationally adiabatic potentials for the system [for example see ref. 44(d)]. More recently, we have identi?ed a large number of these resonances arising from the quantized transition states of the system.45 Since resonances are identi?ed from the oscillatory variation of the PR(E) and v reaction (I) is dominated by them, it would be an ideal testing ground for examining the e?ect of NIP on the dynamics. While carrying out the dynamical study, we have used the following initial WP : W(R, r, t \ 0) \ F(R)/ (r) v (57)

The height and the width of the NIP are chosen depending on the translational energy of the WP using eqn. (3). The timeindependent quantal PR(E) values reported by Sakimoto and 0 Onda49 are superimposed on Fig. 1(a) in the form of a dashed curve. We have performed several WP calculations by varying the height and width of the NIP in order to check the convergence of the results. The results of di?erent calculations di?er slightly in magnitude and the probability decreases when the height of the NIP exceeds some critical value. In each calculation, the PR(E) values show arti?cial oscillations 0 at higher energies, near the onset of collision-induced dissociation.

Fig. 1 State-selected reaction probabilities for the reaction He ] H `(v \ 0) ] HeH` ] H in the energy range 2.0?2.78 eV com2 puted by (a) adding an NIP to the real potential (solid curve) and (b) using a damping function (solid curve). The TIQM results (ref. 49) obtained without the use of the NIP or the damping function are plotted as dashed curves in both panels. The results at lower energies obtained using the NIP or the damping function are indistinguishable from each other and also from the TIQM results and, hence, are not included.

J. Chem. Soc., Faraday T rans., 1997, V ol. 93

777

View Online

In the second set of our calculation we have used a damping function in place of the NIP (but without actually adding it to the potential), to remove the WP that reaches the grid edge. The use of these functions is well known in the literature. We have used a damping function of the following type : n (X ] *X [X) mask mask i ; X PX (61) f (X ) \ sin i mask i 2 *X mask activated outside the dividing line in the product channel and also in the asymptotic reactant channel. X is the point at mask which the masking function is initiated and *X (\X mask max [X ) is the width over which the function decays from 1 mask to 0 with X being the maximum length of the grid along a max particular channel. The time evolved WP is multiplied by f (x ) i in each channel resulting in a sin2 masking. The PR(E) values computed using the above function do not 0 reveal any spurious oscillations at the higher energies and are presented in Fig. 1(b) (solid curve). These PR(E) values are in 0 excellent agreement with those obtained by Sakimoto and Onda (dashed curve). The PR(E) values obtained by Balakrishnan and 0 Sathyamurthy50 using an FFT-SIL propagator also do not reveal any oscillation at higher energies and are in excellent agreement with our results obtained using the damping function. Balakrishnan and Sathyamurthy used a damping function derived from the Neuhauser and Baer?s linear absorbing potential16 and multiplied it by the time-evolved WP at the end of each time step without additing it directly to the real potentials at the grid points. The damping function derived in this way is real and exponential and it preserves the Hermitian property of the Hamiltonian. The implication of the above results is that the artefacts in the computed result, particularly at higher energies, arise when NIPs are added to the real potentials at the grid points. The artefacts disappear when an exponential damping function resulting from NIPs [exp([iC*t/?), C \ [iV ] or any 0 other damping function [e.g. the one in eqn. (61) or the one used by Heather and Metiu51] is used. This method of damping can be achieved more easily than by the direct use of NIPs. Particularly for systems with many narrow resonances, use of the damping function seems more appropriate, as it yields well converged results. The prescription outlined by Taylor and co-workers37 on the use of NIPs in conjunction with the CP scheme follows the latter idea and is free from anomaly since the Hamiltonian remains Hermitian. A nonHermitian Hamiltonian also would not preserve the timereversal symmetry of the TDSE.

C

D

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

ence of an NIP. The SIL scheme also involves an exponential operation and is unconditionally stable, like the SO scheme. It conserves norm as well as energy. This also leads to the same damping factor as the SO, in the presence of an NIP. The accuracy of the scheme depends on the length of the time step. The additional complication of the non-orthogonality of the basis vectors in the presence of an NIP can be taken care of by using an orthogonalization scheme. An additional point to be mentioned here is that the Hamiltonian that includes an NIP does not preserve the time-reversal symmetry of the TDSE. This raises a fundamental question on the wisdom of the use of NIPs in solving the TDSE. For the test case of collinear He, H ` collisions we have 2 shown that arti?cial oscillations in PR(E) arise at higher ener0 gies, if an NIP is included, and that the magnitude of the oscillations depends on the height and width of the NIP. The spurious oscillations are avoided by using a properly chosen damping function externally. Such a procedure works with ease and it also preserves the Hermitian property of the Hamiltonian and the associated basic conservation rules of quantum mechanics. This study was supported in part by a grant from the Commission of European Communities. References

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 R. Koslo?, J. Phys. Chem., 1988, 92, 2087. V. Mohan and N. Sathyamurthy, Comput. Phys. Rep., 1988, 7, 213. Comput. Phys. Commun., 1991, 63, 1?582. J. Broeckhove and L. Lathouwers, T ime Dependent Methods for Quantum Dynamics, NATO ASI Ser. B229, Plenum Press, New York, 1992. G. G. Balint-Kurti, R. N. Dixon and C. C. Marston, Int. Rev. Phys. Chem., 1992, 11, 317. Numerical Grid Methods and T heir Applications to Schro? dinger?s Equation, ed. C. Cerjan, NATO ASI Ser. C412, Kluwer, Dordrecht, 1993. R. Koslo?, Annu. Rev. Phys. Chem., 1994, 45, 145. D. G. Truhlar, Comput. Phys. Commun., 1994, 84, 78. R. Koslo?, in Dynamics of Molecules and Chemical Reactions, ed. R. E. Wyatt and J. Z. H. Zhang, Marcel Dekker, New York, 1996. B. M. Garraway and K-A. Suominen, Rep. Prog. Phys., 1995, 58, 365. N. Balakrishnan, C. Kalyanaraman and N. Sathyamurthy, Phys. Rep., in the press. C. Leforestier and R. E. Wyatt, J. Chem. Phys., 1983, 78, 2334. E. Segre, Nuclei and Particles, Benjamin, New York, 1965, p. 467. R. Koslo? and D. Koslo?, J. Comput. Phys., 1986, 63, 363. R. Koslo? and C. Cerjan, J. Chem. Phys., 1984, 81, 3722. D. Neuhauser and M. Baer, J. Chem. Phys., 1989, 90, 4351. D. Neuhauser, M. Baer, R. S. Judson and D. J. Kouri, J. Chem. Phys., 1989, 90, 5882 ; D. Neuhauser, M. Baer and D. J. Kouri, J. Chem. Phys., 1990, 93, 2499. M. S. Child, Mol. Phys., 1991, 72, 89. T. Seideman and W. H. Miller, J. Chem. Phys., 1992, 96, 4412. A. Vibok and G. G. Balint-Kurti, J. Phys. Chem., 1992, 96, 8712. ? ? M. Monnerville, P. Halvick and J. C. Rayez, Chem. Phys., 1992, 159, 227. D. H. Zhang, O. A. Sharafeddin and J. Z. H. Zhang, Chem. Phys., 1992, 167, 137. D. H. Zhang and J. Z. H. Zhang, J. Chem. Phys., 1993, 99, 5615 ; 1994, 100, 2697 ; 1994, 101, 1146. D. Macias, S. Brouard and J. G. Muga, Chem. Phys. L ett., 1994, 228, 672. C. Leforestier, R. Bisseling, C. Cerjan, M. D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero and R. Koslo?, J. Comput. Phys., 1991, 94, 59. T. N. Truong, J. J. Tanner, P. Bala, J. A. McCammon, D. J. Kouri, B. Lesyng and D. K. Ho?man, J. Chem. Phys., 1992, 96, 2077. E. Merzbacher, Quantum Mechanics, Wiley, New York, 2nd edn., 1970. A. Askar and A. S. Cakmak, J. Chem. Phys., 1978, 68, 2794. D. Koslo? and R. Koslo?, Comput. Phys. Commun., 1983, 30, 333 ; J. Comput. Phys., 1983, 52, 35.

5 Summary and Conclusion

We have presented an overview of the use of NIPs in timedependent wave packet calculations. We have highlighted the usual difficulties one encounters while using NIPs. The SOD scheme is unstable in the presence of an NIP. The SO scheme is unconditionally stable. It conserves norm but the noncommutability of T? and V? results in non-conservation of energy. The exponential operation gives rise to a damping factor in the presence of an NIP and it can be used successfully for damping the WP near the grid edges. Since the method does not depend upon the spectral range of the Hamiltonian, it is stable in the presence of an NIP also. The CP scheme is conditionally stable, and it does not, in principle, conserve either norm or energy. The deviations are used as a check on the accuracy of the method. The stability of the scheme depends strongly on the energy range of the Hamiltonian. The scheme becomes unstable in the presence of an NIP. This can be made stable by adjusting for the shift of the eigenvalue spectrum of the Hamiltonian to the complex plane in the pres778 J. Chem. Soc., Faraday T rans., 1997, V ol. 93

26 27 28 29

View Online

30 31 32 33 34 35 36 37 38 39 40 41 42 43

Downloaded on 14 June 2011 Published on 01 January 1997 on http://pubs.rsc.org | doi:10.1039/A605778K

J. A. Fleck Jr., J. R. Morris and M. D. Feit, Appl. Phys., 1976, 10, 129. M. D. Feit, J. A. Fleck Jr. and A. Steiger, J. Comput. Phys., 1982, 47, 412. G. Arfken, Mathematical Methods for Physicists, Academic Press, Bangalore, 2nd edn., 1994. H. Tal-Ezer and R. Koslo?, J. Chem. Phys., 1984, 81, 3967. R. C. Mowrey, J. Chem. Phys., 1993, 99, 7049. G. Wiesenekker, G. J. Kroes, E. J. Baerends and R. C. Mowrey, J. Chem. Phys., 1995, 102, 3873. Y. Huang, D. J. Kouri and D. K. Ho?man, Chem. Phys. L ett., 1994, 225, 37 ; J. Chem. Phys., 1994, 101, 10493. V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys., 1995, 103, 2903. D. Neuhauser, J. Chem. Phys., 1995, 103, 8513. M. Berman, R. Koslo? and H. Tal-Ezer, J. Phys. A, 1992, 25, 1283. T. J. Park and J. C. Light, J. Chem. Phys., 1986, 85, 5870. D. J. Tannor, A. Besprozvannaya and C. J. Williams, J. Chem. Phys., 1992, 96, 2998. D. Kohen and D. J. Tannor, J. Chem. Phys., 1993, 98, 3168. C. Lanczos, J. Math. Phys., 1938, 17, 123.

44

45 46 47 48 49 50 51

(a) D. J. Kouri and M. Baer, Chem. Phys. L ett., 1974, 24, 37 ; (b) J. T. Adams, Chem. Phys. L ett., 1975, 33, 275 ; (c) T. Joseph and N. Sathyamurthy, J. Indian Chem. Soc., 1985, 62, 874 ; (d) N. Sathyamurthy, M. Baer and T. Joseph, Chem Phys., 1987, 114, 73 ; (e) J. D. Kress, R. B. Walker and E. F. Hayes, J. Chem. Phys., 1990, 93, 8085 ; ( f ) J. Z. H. Zhang, D. L. Yeager and W. H. Miller, Chem. Phys. L ett., 1990, 173, 489 ; (g) B. Lepetit and J. M. Launay, J. Chem. Phys., 1991, 95, 5159 ; (h) N. Balakrishnan and N. Sathyamurthy, Chem. Phys. L ett., 1993, 201, 294. S. Mahapatra and N. Sathyamurthy, J. Chem. Phys., 1995, 102, 6057. T. Joseph and N. Sathyamurthy, J. Chem. Phys., 1987, 86, 704. C. C. Marston and G. G. Balint-Kurti, J. Chem. Phys., 1989, 91, 3571. T. Seideman and W. H. Miller, J. Chem. Phys., 1991, 95, 4412. K. Sakimoto and K. Onda, Chem. Phys. L ett., 1994, 226, 227. N. Balakrishnan and N. Sathyamurthy, Chem. Phys. L ett., 1995, 240, 119. R. Heather and H. Metiu, J. Chem. Phys., 1987, 86, 5009.

Paper 6/05778K ; Received 19th August, 1996

J. Chem. Soc., Faraday T rans., 1997, V ol. 93

779

赞助商链接

相关文章:

- 错误提示及解决方法
- Often, a becomes only slightly
*negative*(due to...7000—7999 Solvers and Preconditioners TABLE 2-5...7147 Out of memory*in**time-dependent*solver The...

- Imaginary Number
- "positive" direction going up; "positive"
*imaginary*numbers then increase*in*magnitude upwards, and "*negative*"*imaginary*numbers increase*in*magnitude ...

更多相关标签:

- Time Dependent Supersymmetry in Quantum Mechanics
- Real-Time Dynamics from Imaginary-Time Quantum Monte Carlo Simulations Tests on Oscillator
- Time-dependent scattering theory of N-body quantum systems
- Differential Harnack Estimates for Time-dependent Heat Equations with Potentials
- Complex Envelope Soliton in Bose-Einstein Condensate with Time Dependent Scattering Length
- Close-Coupling Time-Dependent Quantum Dynamics Study of the H + HCl Reaction
- Local Scattering Operators for P(phi)_2 Models and the Time-dependent Schrodinger Equation
- Resonance Raman scattering in semiconductor quantum dots Adiabatic vs. time-dependent pertu
- The BCS Critical Temperature for Potentials with Negative Scattering Length
- Time delay for one-dimensional quantum systems with steplike potentials
- Noether's Theorem and time-dependent quantum invariants
- Scattering matrix ensemble for time-dependent transport through a chaotic quantum dot
- Polymer dynamics in time-dependent periodic potentials
- Differential Harnack Estimates for Time-dependent Heat Equations with Potentials
- Negative energies and time reversal in Quantum Field Theory