The Chebyshev propagator for quantum systems

Computer Physics Communications ELSEYIER

Computer Physics Communications 119 (1999) 19-31

The Chebyshev propagator for quantum systems Rongqing Chen, Hua Guo ’ Department

of Chemistry,

University

of New Mexico,

Albuquerque,

NM 87131.

USA

Received 3 August 1998

Abstract This paper provides an overview of some recent developments in quantu_m dynamic and spectroscopic calculations using the Chebyshev propagator. It is shown that the Cheby?hev opera?: (Tk( H) ) can be considered as a discrete cosine type pr:Ipagaror ‘(cos~k?Y~j, ‘m titirdn ‘ine an&e oper%ror 1% = arccos ~1 ‘IS a drngm-va’rueir mappmg rr’r ‘tne sc&etr Rarr&oriran (H) and the order (k) is an effective time. Properties in the angle domain (thus in energy domain) can be obtained from those in the order domain via a cosine Fourier transform. The equivalence and parallelism between the order-angle formulation and the time-energy formulation allow one to transplant existing time-dependent quantum mechanical methods. Compared with the time propagator, the Chebyshev propagator has a number of numerical advantages. Because of its polynomial form, the action of the Chebyshev propagator can be evaluated exactly except for roundoff errors. The threeterm recursion formula for Chebyshev polynomials is stable and requires minimal memory. The discretization error in time propa_paium can k avdxku &cause. ‘Ihe
Quantum propagator; Eigenproblem; Dynamics; Fourier transform; Chebysbev potynomiais

where fi is the ro-vibrational Hamiltonian of the system, which is typically discretized using a basis

1. Introduction

The quantum mechanical characterizationof a

expansion. (The spatial discretization of molecu-

timeindependently or time-dependently [ 11. Although the two formulations are formally equivalent for timeindependent observables, their numerical efficiency may differ The calculation of mukcuktr ru-vibrmiunal energy levels serves as a good example to illustrate the differences between the two approaches.The traditional time-independent approach is based on the time-independent Schrijdinger equation,

lar Hamiltonian has been reviewed by a number of authors [2,3] and is not discussedhere.) The eigenvalues and eigenfunctions of H can be obtained using the Rayleigh variational principle which leads to the

moIecuIar

system

EpP) = EpP) , ’ CCorresponzIinp

can be formulated

either

hi&7

dkgmi~~2~&l

auihor.

OOlO-4655/99/$ - see front matter @ 1999 Elsevier Science B.V. All rights reserved. PIISOOlO-4655(98)00179-9

of the Fkani~tcmirm

matiix

141.

The resultant eigenvalues provide information on the ro-vibrational spectrum of the system. Because the Hamiltonian matrix is modified in the diagonalization process using, for instance, the method of Householder [5], this approach requires O( N2> memory and scales as O( N3) in CPU time, where N is the total number of basis functions used in the

20

R. Chen. H. &o/Computer

Physics

calculation. This approach, though very stable and accurate, is unsuitable for systems with large dimensions, even on today’s state-of-the-art computers. In the past decades, various truncation and contraction schemes have been proposed to reduce the size of the basis [6], but significant improvements beyond the current level of sophistication appear to be difficult. The time-dependent approach, on the other hand, takes advantage of the time-to-energy Fourier transformation. In the so-called spectral method [ 7-91, for example, an initial wave packet /@e) is propagated in time by solving the time-dependent Schrodinger equation (atomic units are used in this review unless stated otherwise), (2) where fi is assumed to be independent of time. For systems with time-dependent Hamiltonians, one can employ the (t, t’) formulation [ 10,11], which treats time as an extended spatial dimension. The Fourier transform of the auto-correlation function C(t) = .(@aJe-‘HrI@c) gives the spectrum of the system, G(E)

1 +CC P = k / eiEtC(t)

dt.

(3)

-m

The wave function at a given energy is then assembled from the time-dependent wave packets I@(t)) =

(4) -cc

In practice, the propagation is carried out on a discrete time grid and truncated at a finite time. Under such circumstances, discretization errors may occur in the above integrals. The time-dependent method is an energetically global method: the propagation of a single wave packet yields information at all energies. Of course, the resolution of the spectrum is restricted by the total time of propagation due to the uncertainty principle. In practical calculations, the propagation is truncated at a finite time and window functions can be incorporated into the transformation to reduce noise [ 8,9,12]. An added advantage of the time-dependent method

Communications

119 (1999)

19-31

is its relevance to time-dependent processes [ 13,141. We, however, focus in this review on the calculation on energy domain properties. A major drawback of the time-dependent method is the inability to directly calculate the action of the time propagator on a wave function. Approximate methods for propagating a quantum wave packet have been developed in several directions [ 13,151 and can be classified into two categories. Short-time schemes, such as the split-operator method (SO) [ 8,9], the second-order difference (SOD) method [ 161, and the short-iterative Lanczos method (SIL) [ 171, provide local approximations to the time propagator, U(t, E) = f+? As a result, their accuracy and efficiency depend sensitively on the time interval used in the calculation. Long-time approximations are based on global polynomial interpolations of the time propagator [ 13,1SA19]. Among them, Chebyshev polynomials (Tk( H) ) are one of the most accurate and efficient for Hermitian Hamiltonians [ 181, e -at _-2

(2

- ~kO)(-i)kJk(t>Tk(@)

9

(5)

k=O

where Jk are the Bessel functions of the first kind and i? is scaled to [ -1, 11. Eq. (5) has found a wide range of applications, thanks to the accuracy and stability of the Chebyshev recursion [ 191. It is interesting to note that two transformations are employed in calculating energy properties. Chebyshev polynomials are first used to provide an approximation to the time propagator. Then, information in time domain is transformed to energy domain. The time-dependent method for solving eigenproblems belongs to the class of recursive methods that generate a set of states, or wave packets, using specific recursive formulae. These states eventually span a Hilbert space in which the system Hamiltonian becomes invariant. The eigenproblem can then be solved via transformation. A canonical example of the recursive methods is the Lanczos algorithm [ 201, which generates orthonormal states from an arbitrary initial state using a three-term recursion. The Hamiltonian is reduced to a tridiagonal form in the Lanczos basis and approximate eigenvalues can be readily obtained by diagonalization. There are several distinct advantages associated with the recursive approach. First, approximate solutions can be obtained at any step during the

R. Chen, H. Guo/Computer

Physics

recursion. They may not be sufficiently accurate, but one does not have to wait until the end of the calculation to learn some useful properties of the system. The approximation can be systematically improved by running more recursive steps. In the Lanczos approach, for example, one can expect convergence of extremal eigenvalues before the entire spectrum is resolved. Secondly, unlike conventional diagonalization where direct modifications of the Hamiltonian matrix are often unavoidable, a recursive process typically involves a simple multiplication of the Hamiltonian matrix onto a vector. Because of its sparsity, the Hamiltonian matrix is often not stored and can be generated on-the-fly. Finally, the recursion usually needs only a small number of previous states. Thus, the storage requirement, which is proportional to O(N), is much less demanding. In this paper, we concentrate on a recursive method based on Chebyshev polynomials for calculating energy domain properties. Kosloff and coworkers were the first to use Chebyshev polynomials as building blocks in approximating the time propagator and its variants [ 13,18,21-231. The pioneering work by Kouri, Hoffman and their coworkers demonstrated the efficiency and accuracy in approximating timeindependent operator functions, such as the Green and delta functions [ 24-291. Extensions of this approach to include absorbing boundary conditions were advanced by Mandelshtam and Taylor [ 30-321. Other research groups [ 33-371 including us [ 38-441 have also contributed to the understanding of the Chebyshev approach. In this review, we show that the Chebyshev operator possesses many properties similar to those of the time propagator. In particular, a cosine Fourier transformation exists between the angle and order domains of the propagator, similar to the Fourier transformation between time and energy. The transformation can be used to extract information in energy domain in an efficient manner. The angle-order conjugacy is not only important in understanding the efficiency and accuracy of Chebyshev polynomial based approaches, but also instrumental in developing new and more general methods. In our discussion presented here, we do not follow strictly the chronological order of events leading to the method. Rather, we attempt to provide a coherent, self-contained, and uniform prospective of the approach. This paper is organized as follows. In the next section (Section 2))

Cornrnunications

119 (1999)

19-31

21

we first review the properties of Chebyshev polynomials, followed by the discussion of the Chebyshev propagator. The discretization of both the order and angle domains is considered next. Finally, the relationship with other methods is explored. In Section 3, several applications of the Chebyshev propagator are surveyed. A brief conclusion is given in the final section (Section 4).

2. Chebyshev

propagator

2.1. Chebyshev polynomials

The Chebyshev polynomials of the first kind are a unique class of orthogonal polynomials based on the cosine function [ 451, Tk(w) c cos(k arccosw) = cos( k4) .

(6)

They are polynomials of o, and at the sametime cosine functions of 4. Chebyshev polynomials satisfy the following recursion formula: Tk+l t@> = 2wT,(w> - Tk--l tw) 7 for k > 1 , (7a) with To(o) = 1,

TI(w) =m,

t7b)

where the recursion (Eq. (7a)) can be consideredas a variant of the trigonometric identity, cos[(k+

l>qbl +cos[(k-

= 2cos(+)

1)41

cos(k4).

(8)

It can be readily seen that Chebyshev polynomials are actually cosine functions in disguise. Computationally, fast Fourier transform can thus be implemented [ 51. The trigonometry and polynomial viewpoints are equivalent and complementary. The duality asboth cosinefunctions and polynomials plays an important role in our discussionbelow. 2.2. Propagation

in Chebyshev

order domain

Now, we replacethe scalar variables in the last section by operators.In particular, we are primarily interestedin Hermitian Hamiltonian which have real eigenvalues. Accordingly, the Chebyshev operator, which

22

R. Chen, H. Guo/Conzputer

is a polynomial function of the Hamiltonian, expressed in the cosine form, Fk(E)

z cos(kG)

Physics

can be

,

(9)

where the angle operator o^ is a single-valued, nonlinear mapping of the Hamiltonian, G = arccos a.

but (10)

Apparently, the Hamiltonian has to be scaled to [ -1, 11, which is achieved as below [ 181: fi=

(ii0 -fi)/AH,

(11)

where fia is the original Hamiltonian, whose spectral middle point and half-span are given by z = (H,,, + Hmin) /2 and AH = (H,,, - Hkn) /2. From here on, all the Hamiltonians are assumed to be properly scaled. An eige_nstate of the Hamiltonian is also an eigenstate of 0. The corresponding eigenvalue is B = arccos E and 8 E [ 0, z-1. The Chebyshev angle operator can thus be considered as an effective Hamiltonian [ 391. The cosine form of the Chebyshev operator implies a subtle connection to two exponential operators, cos( kO^) = ( e-ikZ + dkO^) /2,

(12)

where the terms on the right-hand side resemble the forward and backward time evolution operators. Note that k 2 0. The former propagator e@@, for example, leads to a differential equation similar to the timedependent Schrodinger equation [ 37,391, i$D(k))

=&P(k)),

(13)

if k can be treated as a continuous variable. It is thus reasonable to regard k as the effective time [ 391. Accordingly, the Chebyshev states are the analogs of the time-dependent wave packets, Chebyshev propagation

Discrete time propagation

PO) 7 I@,) = cos(Gi> l@o) , I@*) = cos(2G)~@lJ),

PO) 3 I#) = e-iq@o)

Communications

119 (1999)

19-31

where r is the time interval and ]@a) is the initial state. Consequently, the Chebyshev operator can be conveniently considered as a discrete cosine type evolution operator, or a propagator in the Chebyshev order domain [ 391. The Chebyshev states are the orderdependent wave packets. The analogous time-energy version of the cosine propagator has been extensively discussed by Gray [ 461. In contrast to the exponential time propagator which cannot be calculated directly, the Chebyshev propagator can be computed easily and exactly since it is a polynomial of fi. By this we mean that no approximation (such as interpolation) to the propagator is necessary. In practice, the recursion of the Chebyshev operator is not explicitly ca@ed out. Instead, the Chebyshev states ]@k) = Tk( H) /@a) are calculated recursively [ 181, I&+1) = %@b)

- I&-l)

,

for k 2 1 ,

(144

and PO) = @Jo),

I@)

=

~l@O)~

(14b)

where Fis the identity operator. The above recursion is accurate and stable for H E [ - 1, l] , which allows the propagation in the Chebyshev order domain for tens of thousands of steps without accumulating significant errors. Furthermore, since only two previous vectors need to be stored and Hamiltonian is usually sparse, the propagation is very memory friendly. The fact that the Chebyshev operator is an efficient propagator is useless unless an efficient and accurate connection to energy domain can be established. Previously, such a link comes from the Chebyshev approximation of the time propagator (Eq. (5)) and the subsequent time-to-energy transformation. Actually, a more direct relationship can be found between the Chebyshev order domain and energy domain. Indeed, the cosine Fourier transform of the Chebyshev states yields a wave function in the Chebyshev angle domain, (2 - &o) cos(Wl@k)

,

I&) = e--iE271@o) ,

k=O

=-

(2 - &,j) cos( k0) cos( k6) /Go) :g

= 6(8 - o^)J@o) = p(e))

t

(15)

R. Chen, H. Guo/Computer

Physics

where the completeness of the cosine functions is used. Unlike the time-energy Fourier transform, the transformation here is naturally discrete. The inverse transform is x

I@k)= / decos(ke)p(e)).

(16)

0

Finally, the wave-function in the Chebyshev angle domain can be readily transformed to the energy domain via a change of variables, E = costi. The interesting duality of the Chebyshev operator as a cosine type propagator and a recursive polynomial function distinguishes itself from other generalized propagators. Lanczos [47] has realized long ago that the cosine transform between the two domains can be used for spectral analysis. The cosine transform between the angle and order domains indicates that k and 13are a pair of conjugate variables, analogous to time and energy. In fact, the k-0 formulation is equivalent to the better known t-E formulation. This equivalence has two important implications. First, methods developed for the t-E formulation can be easily adopted for the k-8 formulation, as explained below. Secondly, one can choose the most appropriate formulation for the problem under investigation. While the t-E formulation may be the method of choice for time-dependent phenomena, the k-0 formulation is more efficient for time-independent properties. It should also be pointed out that the two formulations are not independent. Properties in one formulation can be readily expressed in the other formulation, e.g. Eq. (5). We note in passing that the propagation scheme can be adapted to take advantage of the molecular symmetry. For example, propagation states and/or correlation functions belonging to different symmetry species can be obtained with the propagation of a single wave packet [ 48,491. Information obtained in this approach is identical to that from multiple wave packet propagation. The trick here is to extract different symmetry species at a given propagation step using symmetry projection operators. This is possible because the system Hamiltonian and symmetry operators commute. The Chebyshev propagation retains many favorable properties of time propagation. It remains an initialvalue problem and requires a relatively small memory. Like time propagation, it is an energetically global

Cotnmunications

119 (1999)

19-31

23

method since single propagation is capable of yielding information at all energies. However, its numerical performance is superior to methods based on the time-to-energy transformation for a number of reasons. First, its polynomial form makes possible for direct computation with no interpolation. Second, it does away with the t-E transformation and avoids discretization errors such as those present in the time integration. 2.3. Discrete energy representation An immediate observation from the two equations above reveals the inequality between the Chebyshev order and angle: the former is discrete and the latter is continuous. In practical applications, it is convenient to work with discrete representations. To this end, we can select a set of equidistant points in the angle domain 8,, = (m+ l/2) m/K and approximate an inner product as a discrete sum over these points,

The above approximation is the same as the GaussChebyshev quadrature if one changes the angle to energy (E,,, = cos e,,). The points chosen in the above equation are the same as the Gauss-Chebyshev quadrature points in energy domain, which are the zeros of the Kth order Chebyshev polynomial. This set of points plays an important part in defining the discrete energy representation (DER) [ 421 for the Chebyshev propagator, which is discussed in more detail below. Note that the DER points are not equally spaced in energy domain due to the nonlinear nature of the mapping, but they are distributed globally and interlace. The above quadrature is exact if F; (0) Fz( 0) is a polynomial of E = cos 0 with an order of 2K - 1 or less. We note that Chebyshev polynomials are both orthogonal and complete on the DER grid, cos( kf?,) m=O

COS(

k’&,)

= 8kk’ , (18a>

24

R. Chen, H. Guo/Computer

Physics

Communications

k=O (18b)

The above equations can be deduced from either the properties of Chebyshev polynomials or those of trigonometric functions. Hence, the angle and order domains are completely equivalent and symmetric in their discrete representations. As a result, an analytical function of the Hamiltonian can be expanded in terms of the Chebyshev propagator,

K-l fk cos( ki+) .

(19)

k=O

From the orthogonality of the cosine functions or the Chebyshev polynomials, one can obtain the expansion coefficients, K-l

F(cos e,,,) cos( ken,) .

fk2+C

(20)

ni=o

This equation can be efficiently implemented using the fast Fourier transform [ 51. Furthermore, we can show that DER and the so-called generalized time representation (GTR) are isomorphic [42], where the latter is identified with the effective time k. If Eq. (20) is recast in the DER/GTR formulation, we have the matrix form [fkl

=

[Tknzl[Fnrl 9

(21)

where fk = F,,, = F(G)

1

) kk’

+2-

- SkO)

c

cos( ke,) cos( k/e”,>

ndl = akk’

,

(T+T)n,,n/ =

(234 (2 - 8kO) cos( ken,) cos( ken,)

k=O

[fkl

=

[TnS;kl

[fkl

.

(24)

Hence, [F,,l and [fk] are the conjugate representations of F(H) in DER and GTR, respectively. The concept of DER and GTR bears close resemblance to the well-known discrete variable representation (DVR) [ 2,505 1] and the Fourier method [ 31 in spatial discretization. In fact, all thesemethodsare examplesof the orthogonal collocation scheme [52] that uses a set of orthogonal functions to establish a discrete approximation. The collocation allows the calculation of a given operation in the most appropriate representation. The collocation based on Chebyshev polynomials is particularly important from the viewpoint of interpolation. In approximating an arbitrary function, the Chebyshev interpolation is very close to the best possible polynomial approximation to the function. Indeed, Chebyshev polynomials are the polynomials of choice for the global fitting of an arbitrary non-periodic function [ 451. The usefulnessof Eq. (19) was first recognized by Tal-Ezer and Kosloff [ 181 in time propagation. Their treatmentof Chebyshevpolynomials asoperatorshasa profound influence on later research.The adoption for time-independentoperatorfunctions was pioneeredby Kouri, Hoffman and coworkers [ 24-271. We give below the explicit expressionsfor two important operator functions, the Green function [24] and the Dirac delta function [ 271, -i

cc

(2 - &+)CikBcos(k6), sin e ck=O

S(E-$=--&g (2 -

(25)

a,&,) COS( ke) cos( k&j).

k=O

(26)

K-l (TT+

[Fml = [T,:l

1

The transformation matrix is unitary (+ denotes the Hermitian conjugate),

(23b)

which implies that DER and GTR contain the same amount of information and are completely equivalent. The inverse transformation of Eq. (2 1) also exists,

y=E - id

2 - s/& ~ COS( ke,,,) . K

Tkm=

19-31

= ~“wd,

f y(2 - SK))cos(k&,)cos(k&t) = a,,,/ .

F( ii) ?Z c

119 (1999)

It is interesting to note that the expansioncoefficients in both expressionsdo not decay as k increases.A truncation at a finite k may result in sine-like oscillations in the sidebands.Another interesting observation is that the energy information is all contained in the expansioncoefficients, which can be readily calculated

R. Chen, H. Guo/Cmpoter

Physics

as the energy parameter sweeps through the spectrum. Provided the initial state contains all the energy components, the same set of Chebyshev states can be used for all the energies of interest since the propagator is independent of E. This property is true for the Chebyshev expansion of arbitrary operator functions. Since Chebyshev polynomial based expansions are defined on the real axis, they work very well for Hermitian Hamiltonians which have only real eigenvalues. However, many scattering problems require the imposition of boundary conditions that renders the Hamiltonian non-Hermitian. Consequently, the Chebyshev propagation may become unstable. Although other propagators based on Legendre [ 531, Faber [ 25,261 and Newton 154-561 polynomials have been proposed, Mandelshtam and Taylor [ 30,311 showed that a simple modification of the recursion formula (Eqs. ( 14) ) is sufficient for this class of problems, c+,(H)

=eeY[2iic(fi)

-epYQe’_,(fi)],

(27a)

with

(27b) The phenomenological damping term can be identified with the absorbing boundary conditions for scattering problems and can also be thought as a window function in propagation which converges the expansions such as Eqs. (25) and (26). This modification greatly increased the scope of problems amenable to the Chebyshev approach. Similar ideas have also been proposed by others [ 29,36,40]. 2.4. Other related methods The use of polynomial expansions to approximate operator functions is not restricted to Chebyshev polynomials. Legendre [53], Faber [ 25,261 and Newton polynomials [ 54-561 have been successfully applied to achieve the same objectives. These types of polynomials are particularly important in dealing with non-Hermitian operators. One can also understand the polynomial approaches from the Krylov subspace viewpoint [ 57,581. The K-dimensional Krylov subspace is spanned by the Krylov states {skl@c): , K - l}, which is the same as the space k=O,l,... spanned by the Chebyshev states {I&) = Tk( g) I@n): k=O,l,... , K - I}. This is true for all polynomial

Comnunications

119 (1999)

19-31

25

propagators as long as the corresponding polynomials are linearly independent and up to an order of no more than K - 1. An interesting member of the Krylov subspace based methods is that of Lanczos [20], Using a three-term recursion, the Lanczos algorithm not only generates the Krylov vectors, but also orthogonalizes them at the same time. This results in a tridiagonal Hamiltonian matrix which can be diagonalized with relative ease. The final eigenstates can be obtained by a unitary transform of the Lanczos vectors. A distinct difference from the Chebyshev and other polynomial based methods is that the Lanczos propagation depends sensitively on the initial vector. As a result, the transformation matrix is not known prior to the propagation. The Lanczos method has been used to extract energy localized states using the Green function as a filter [ 59-631. In exact arithmetic, the Lanczos method should be equivalent to the polynomial based methods. However, the two behave quite differently in finite arithmetic. It is well known that the Lanczos vectors lose global orthogonality as K increases because of roundoff errors. The loss of global orthogonality may lead to the emergence of spurious eigenvalues. Although these spurious eigenvalues can be identified, for example, by the Cullum-Willoughby test [ 641, the effectiveness of generating Krylov vectors is to certain extent compromised. On the other hand, the Chebyshev states are not orthogonal in coordinate space (although they are in principle linearly independent). Hence, a Hermitian Hamiltonian in Chebyshev basis is usually full, in contrast to the extremely sparse and structured tridiagonal form in a Lanczos basis. Recent calculations of large systems [ 65-701 indicated that the performance of the Lanczos method is comparable to that of the Chebyshev method. It is, however, not clear at present which one, the Lanczos algorithm or the Chebyshev propagation, is more efficient in expanding the Krylov space in finite arithmetic. Our preliminary tests indicated that the answer would depend on the accuracy of the calculation.

26

R. Chen, H. Guo/Cotnputer

Physics

3. Applications 3. I. Spectral method In the spectral method, eigenvalues and eigenstates are obtained directly from a transformation between conjugate domains. For Chebyshev order-angle formulation, the spectrum in the angle domain is expressed as a cosine Fourier transform of the autocorrelation function in the order domain, G(B) = (cDo18(0 - o^+Do) =-:E

(2 - S@) COS(k~)C~)

where the order-dependent is defined as

auto-correlation

Ck = (@()I cos( k6) I@(J) .

= &j

function (29)

G(B) 1

IL@(S)) =S(O - G>l@o) (2 - ho) COS(Wl@k)~

(31)

The energy wave function can be obtained using a formula similar to Eq. (30). The resolution in the angle domain is inversely proportional to the total number of propagation steps in a truncated expansion. To improve the resolution, one can introduce window functions in the above transformation. The similarities to the time-energy formulation are apparent. In calculating the correlation function in Eq. (29)) one can take advantage of the symmetry of the propagator to double the propagation steps, C*k 5z (@ol cos( 2.M) po) = (@a\(2Cos*(ke) = 2(@&&)

- l)l@e)

- co.

C*k+l = q@k+l(@k)

- Cl .

19-31

Compared with the spectra1 method in t-E formulation, the correlation function here is directly obtained from the exact Chebyshev states while the time correlation function is obtained using approximate timedependent wave packets. If the time propagator is approximated by Chebyshev polynomials (Eq. (5) ) , many Chebyshev terms near the end of the expansion have very small expansion coefficients and are wasted in order to achieve convergence. In contrast, even the last state in the Chebyshev propagation contributes equally as the earlier states. 3.2. Filter-diagonalization

where E = cos 8. Similarly, an eigenstate of the angle operator and thus of the Hamiltonian can be assembled from the Chebyshev order-dependent wave-packets,

=-:g

I19 (1999)

method

(28)

The energy spectrum, normalized on the energy scale, is thus c(E)

Communications

The filter-diagonalization (FD) method was first introduced by Neuhauser [71,72] to circumvent the uncertainty limit in the spectra1 method. The idea is to construct a set of primitive filtered states in a prespecified energy range [El,,, EuP] using the spectral method. Since these states do not need to be completely resolved in energy, relatively short propagation suffices. The final resolution of the eigenvalues/eigenstates is achieved variationally in a subspace spanned by the filtered states. To this end, a generalized eigenequation HB = ESB

(33)

is set up and the elements of the Hamiltonian and overlap matrices are computed from the filtered states IWE/)) = F@; El) PO),

H/P= (Wdl&Wd) SW= (AIL)

9 >

(34a) C-b)

where El E [ EIOw, Eup] . If the Hamiltonian becomes invariant in the subspace, the eigenvalues E,, can be obtained from the solution of Eq. (33) and the eigenstates are approximated by the linear combination of the filtered states, L

IEn)= c &nl’f’(E~)). /=I In the original implementation, Neuhauser [ 7173] and others [ 68,741 used the time-energy spectra1 method to construct the filtered states. Subsequently, it was shown that the Chebyshev order-angle formulation provides a more efficient and accurate platform

R. Chen, H. Guo/Comnputer

Physics

for filtering [ 27,28,38,75,76]. To this end, a function of the Hamiltonian, or a filter, is first designed to amplify amplitudes of eigenstates near El = cos Br and suppress those elsewhere. The filter is then directly expanded in terms of the Chebyshev propagator, K-l

F(ii;E,)

= c fk(O/)cos(kO^) k=O

.

(36)

The filter can be quite general and tailored to suit the problem [ 381. Examples of the filter operator include zero width functions such as the Green [24,30,59] and delta functions [ 271, as well as finite width functions [ 21,22,38,77]. In practice, a number of filtered states at different El in [El,, Eup] are assembled during the Chebyshev propagation. Typically, the number of filtered states (L) is much smaller than the number of propagation states (K). Therefore, filtering according to Eq. (36) is essentially a projection of the entire Krylov subspace spanned by the Chebyshev states to a much smaller subspace spanned by the filtered states. Although completeness in [EI,,, EUP] can be expected for the filtered states with a sufficiently large L, linear independence is not guaranteed. In fact, the filtered states are often highly dependent. As a result, effort has to be directed to the removal of the linear dependency, usually via orthogonalization. In studying resonances, FD is still applicable, but requires several modifications. The inner product in Eq. (17) should be replaced by the complex symmetric product [78] to accommodate the non-Hermitian Hamiltonian. The Chebyshev propagator can be modified to satisfy the outgoing boundary conditions. This is typically done by removing the wave packet in the asymptotic region using, for example, the method of Mandelshtam and Taylor [ 30,3 11, cf. Eqs. (27). The variational solution of Eq. (33) yields complex eigenvalues, which determine the energies and lifetimes of the resonances. If one is only interested in the spectrum of the system, the explicit calculation and storage of the filtered states are unnecessary. The situation is similar to the Lanczos algorithm, where all scalar quantities are determined by the tridiagonal Lanczos matrix [79]. As was first realized by Wall and Neuhauser [ 801, the only information required in calculating the elements of the Hamiltonian and overlap matrices is the cor-

Communications

119 (1999)

19-3I

27

relation function. They also showed that the eigenvalues of the Hamiltonian, although represented by a full matrix in the propagation states, can be obtained efficiently using the so-called low-storage filter diagonalization (LSFD) . Later, Mandelshtam and Taylor [ 81-831 further improved the efficiency by using discrete versions of the Chebyshev/time propagators and truncated delta filters. Compared to other filters, the truncated delta filter has simple and analytically known expansion coefficients, and their recurrent relations afford a dramatic simplification of the expression. LSFD is a powerful method to extract frequency information from signals in the “effective” time domain. Since the signal is not restricted to quantum correlation functions, this method is applicable in diverse areas, including semiclassical eigenvalue determination [ 84,851, NMR [ 86,871, and molecular dynamics [ 881. One can also take advantage of the transformation between the angle and order domains to calculate the elements of the Hamiltonian and overlap matrices [ 411. Specifically, both matrices can be expressed in DER where the filters are diagonal,

2K-

1 F(E,,,;

=c

G)E,X(&;

EP)G(E,,,)

,

(37a)

m=O SW

= (PC&)

t~v(Er))

ZK-I F(E,,;E/)F(E,,,;E/l)G(En,).

=c

(37b)

m=o

The spectral function G can be directly calculated from the cosine Fourier transform of the order-dependent correlation function according to Eq. (28). The overlap between an eigenstate and the initial wave packet can also be calculated from the correlation function, L

K-l

(Enl@o) =c B/n c fk(E1)Ck. I=1

(38)

k=O

The summation in Eqs. (37) can usually be truncated to include only the DER points inside the energy range of interest and some in its immediate surroundings. Typically, the number of DER points is less than 1OL. It can be shown that our scheme of calculating the matrix elements is equivalent to the method of Mandelsh-

28

R. Chen,

H. Guo/Computer

Physics

tam and Taylor [ 81-831 for truncated delta filters. Because the simplicity of the delta filters, the latter approach is more efficient. However, our method outlined above is applicable to arbitrary filters. Further, it provides a different and complementary viewpoint. 3.3. Photodissociation spectrum

I19

(1999)

19-31

This is possible because the order-dependent wave packet, like the time-dependent wave packet, is a coherent superposition of the scattering states. The resonance Raman amplitude, in the formulation of Kramer, Heisenberg and Dirac (KHD), is writ-

ten as a matrix element of the Green function operator [ 94,951,

and resonance Raman

The time-dependent and resonance Raman Heller and coworkers Fermi’s golden rule, the section (T(W) c( wZ( w) nential Fourier transform function [ 71,

Communications

theory of photodissociation spectrum was pioneered by [7,89-921. Based on the total photodissociation cross can be expressed as an expoof the time auto-correlation

1 EiftLWl-@+ir

(43)

where g is the excited state Hamiltonian and r represents the damping factor. Using Eq (25)) we can express the amplitude as an exponential Fourier transform of the cross-correlation function C,f’i = (Gfl cos( kO^) I@;) [ 13,441, .

af;(w,)

= 3

M

x(2

- 8kO)e-ikHckfi,

(44)

k=O

= ( @OI 7 eicE-‘)’

dti@o)

--co 03 =

J

eiEtC( t) dt ,

(39)

--DC)

where E = Ei + tiw. In the order-angle formulation, the cross section can be directly expressed in terms of the auto-correlation function in the Chebyshev order domain [ 431, (2

- Sk01 COS(WCk,

3.4. Collisional dynamics (40)

where we have used Eq. (26). This equation is almost

identical to Eq. (30) except that the initial wave packet in this case is well defined, I@O)

= &lpi)

3

where the damping factor is assumed to be small. This formula also resembles the time-dependent version derived by Heller and coworkers [90,92]. The investigation of photodissociation dynamics and resonance Raman spectrum can be bundled into a single calculation. The major task of the calculation is the propagation of the excited state wave packet in the Chebyshev order domain. Since the initial wave packet is always real, the entire propagation can be carried out in real space.

(41)

where IPi) is the ro-vibrational state of the parent molecule in its ground electronic state and pu,, is the transition dipole moment. The partial cross section of a photodissociation process can be determined by projecting the orderdependent wave packet onto the asymptotic basis [ 43,931, (42)

The first application of the Chebyshev propagator in time-dependent scattering calculations was carried out by Huang et al. [24]. Subsequently, these au-

thors have explored the approach in many different contexts [ 2.5,26,29,96,97]. It was pointed out that the scattering problem can be formulated in terms of timeindependent wave packets, in which the major numerical task is to evaluate the action of the Green operator on the initial wave packet, subjected to appropriate boundary conditions. These steps can be efficiently accomplished by the Chebyshev expansion of the Green operator [98]. In practice, an initial wave packet near the edge of the interaction region in the reactant channel is propagated using the Chebyshev propagator. These propagation states can then be assembled to yield filtered states at various energies of

R. Chen, H. Guo/Cotnputer

Physics

interest, subjected to appropriate boundary conditions. Scattering information is extracted from the propagation states. It is worth noting that Gray and BalintKurti [ 371 have recently shown that the S-matrix can be calculated from an exponential Fourier transform of the product flux in the Chebyshev order domain. Very recently, Mandelshtam [ 991 pointed out that it is possible to extract all the S-matrix elements from crosscorrelation functions obtained in a single Chebyshev propagation using the idea of filter-diagonalization. For problems involving continua, negative imaginary potentials are usually included in the Hamiltonian to impose the boundary conditions [ 1001. However, they destroy the Hermiticity of the Hamiltonian and may result in instability in the Chebyshev propagation. This problem can be alleviated by removing outgoing wave packets near the edges of the grid using damping potentials [ 30,3 1,36,40]. By doing so, appropriate boundary conditions are enforced without sacrificing the Hermiticity of the Hamiltonian. Its time-dependent version has long been realized [ 1011. Mandelshtam and Taylor [ 30,3 1 ] pointed out that the damping potential can be considered as a convergence factor in the expansion, which is operative only in the asymptotic region, much like the treatment of the E factor in the Green function by Seideman and Miller [ 1021. As mentioned several times above, a distinct feature of the Chebyshev propagation is the possibility of working with real wave packets. This aspect has been explored by Kroes, Neuhauser and coworkers [ 103,104] and by Gray and Balint-Kurti [ 371. In cases where resonances overlap with direct scattering wave functions, a hybrid method has been suggested by Neuhauser [ 721, which combines the explicit propagation for the direct scattering and filterdiagonalization for extracting the resonances.

4. Conclusions The Chebyshev propagation represents an efficient and accurate method in calculating time-independent properties of a quantum mechanical system. This approach is fundamentally different from traditional time-independent methods. In fact, it is more closely related to time-dependent methods despite the absence of the time variable. A key characteristics of this approach is the fact that the angle and order of the

Cornrnunications

119 (1999)

19-31

29

Chebyshev propagator form a pair of conjugate representations that are related to each other via a cosine Fourier transformation. The angle-order conjugacy resembles in many ways the well-known energytime conjugacy, which is based on an exponential Fourier transformation. One can choose to generate information in one domain and transform it to the other domain. Compared with the time propagator, the Chebyshev propagator has a number of advantages. It is naturally discrete and can be computed free of interpolation errors because of its polynomial form. For a Hermitian Hamiltonian and a real initial state, the entire propagation can be carried out in real space. The propagation in the Chebyshev order domain requires relatively small memory and is very stable. Single propagation yields information at all energies and the global accuracy of the Chebyshev propagation is guaranteed by the minimax properties of Chebyshev polynomials. The Chebyshev propagator is destined to become a versatile and effective tool in quantum mechanical studies of molecular systems. Acknowledgements This work was supported by the National Science Foundation and by the Petroleum Research Fund administrated by the American Chemical Society. References [l]

G.C. Schatz, M.A. Ratner, Quantum Mechanics in Chemistry (Prentice Hall, Englewood, NJ, 1993). [Z] J.C. Light, 1.P. Hamilton, J.V. Lill, J. Chem. Phys. 82 (1985) 1400. [3 1 R. Kosloff, in: Numerical Grid Methods and their Applications to Schrodinger’s Equation, C. Cerjan, ed. (Kluwer, Dordtecht, 1993). [4] G.D. Carney, L.L. Sprandel, C.W. Kern, Adv. Chem. Phys. 37 (1978) 305. [S] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes, 2nd ed. (Cambridge University, Cambridge, 1992 ) [6J Z. Bacic, J.C. Light, Ann. Rev. Phys. Chem. 40 (1989) 469. 171 E.J. Heller, J. Chem. Phys. 68 (1978) 2066. [8] M.D. Feit, J.A. Fleck, A. Steger, J. Comput. Phys. 47 (1982) 412. 191 M.D. Feit, J.A. Fleck, J. Chem. Phys. 78 (1983) 301. [IO] N. Moiseyev, J. Chem. Phys. 101 (1994) 9716. [ I I] U. Peskin, N. Moiseyev, J. Chem. Phys. 99 (1993) 4590.

30

R. Chen, H. Guo/Computer

Physics

[ 121 J. Dai, J.Z.H. Zhang, .I. Chem. Phys. IO3 (1995) 1491. [ 131 R. Kosloff, J. Phys. Chem. 92 ( 1988) 2087. [ I4 I J. Manz, in: Femtochemistry and Femtobiology, V. Sundstrom, ed. (World Scientific, Singapore, 1997). 1151 C. Leforestier, R.H. Bessling, C. Cerjan, M.D. Feit, R. Friesner, A. Guldberg, A. Hammerich, Cl. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, 0. Roncero, R. Kosloff, J. Comput. Phys. 94 (1991) 59. [ 161 D. Kosloff, R. Kosloff, J. Comput. Phys. 52 (1983) 35. [ 171 T.J. Park, J.C. Light, J. Chem. Phys. 85 (1986) 5870. [ 181 H. Tal-Ezer, R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [ 191 R. Kosloff, Annu. Rev. Phys. Chem. 45 (1994) 145. [ 20 [ C. Lanczos, J. Res. Nat]. Bur. Stand. 45 (1950) 255. [ 21 I R. Kosloff, H. Tal-Ezer. Chem. Phys. Len. 127 (1986) 223. 122 I A.D. Hammerich, J.G. Muga, R. Kosloff, Israel J. Chem. 29 (1989) 461. 1231 B. Hartke, R. Kosloff, S. Ruhman, Chem. Phys. Lett. 158 (1989) 238. [ 24 [ Y. Huang, W. Zhu, D. Kouri, D.K. Hoffman, Chem. Phys. Lett. 206 ( 1993) 96. [ 25 J Y. Huang, D.J. Kouri, D.K. Hoffman, Chem. Phys. Lett. 225 (1994) 37. [ 26 I Y. Huang, D.J. Kouri, D.K. Hoffman, J. Chem. Phys. 101 (1994) 10493. 127 [ W. Zhu, Y. Huang. D.J. Kouri, C. Chandler, D.K. Hoffman, Chem. Phys. Lett. 217 (1994) 73. [28] D.J. Kouri, W. Zhu, G.A. Parker, D.K. Hoffman, Chem. Phys. Len. 238 ( 1995) 395. [ 291 Y. Huang, S.S. Iyengar, D.J. Kouri, D.K. Hoffman, J. Chem. Phys. 105 (1996) 927. [30[ V.A. Mandelshtam, H.S. Taylor, J. Chem. Phys. 103 (1995) 2903. [ 311 V.A. Mandelshtam, H.S. Taylor, J. Chem. Phys. 102 (1995) 7390. [32 [ V.A. Mandelshtam, in: Multiparticle Quantum Scattering with Applications to Nuclear, Atomic and Molecular Physics, D.G. Truhlar, B. Simon, eds. (Springer, New York, 1996). [ 331 L.-W. Wang, Phys. Rev. B 49 (1994) 10154. [ 341 T. Takatsuka, N. Hashimoto, J. Chem. Phys. 103 (1995) 6057. 135 [ R.N. Silver, H. Ridder, A.F. Voter. J.D. Kress, J. Comput. Phys. 124 (1996) 115. [36] H.-G. Yu, SC. Smith, J. Chem. Phys. 107 (1997) 9985. [37[ SK. Gray, G.G. Balint-Kurti, J. Chem. Phys. 108 (1998) 950. 1381 R. Chen, H. Guo. J. Chem. Phys. 105 (1996) 1311. 1391 R. Chen, H. Guo, J. Chem. Phys. 105 (1996) 3569. I40 [ R. Chen, H. Guo, Chem. Phys. Lett. 261 ( 1996) 605. I41 ] R. Chen, H. Guo, Chem. Phys. Lett. 279 (1997) 252. 1421 R. Chen, H. Guo, J. Chem. Phys. 108 (1998) 6068. [43] H. Guo, J. Chem. Phys. 108 (1998) 2466. 1441 H. Guo, Chem. Phys. Lett. 289 (1998) 396. [45] J.P. Boyd, Chebyshev and Fourier Spectral Methods (Springer, Berlin, 1989). 146 1 S.K. Gray, J. Chem. Phys. 96 ( 1992) 6543. [47] C. Lanczos, Applied Analysis (Prentice Hall, Englewood Cliffs, NJ, 1956).

Communications

119 (1999)

19-31

[48] R. Chen, H. Guo, Phys. Rev. E (1998) 7288. 1491 R. Chen, H. Guo, L. Liu, J. Muckerman, J. Chem. Phys. 109 (1998) 7128. [50] D.O. Harris, G.G. Engerholm. W.D. Gwinn, J. Chem. Phys. 43 (1965) 1515. [51] AS. Dickinson, P.R. Certain, J. Chem, Phys. 49 (1968) 4029. [52] D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods (SIAM, Philadelphia, PA, 1977). [53] Y. Huang, W. Zhu, D. Kouri, D.K. Hoffman, Chem. Phys. Lett. 214 (1993) 451. 1541 M. Berman, R. Kosloff, H. Tal-Ezer, J. Phys. A 25 ( 1992) 1283. [551 G. Ashkenazi, R. Kosloff, S. Ruhman, H. Tal-Ezer, J. Chem. Phys. 103 (1995) 10005. 1561 S.M. Auerback, C. Leforestier, Comput. Phys. Commun. 78 (1993) 55. [ 57 1 Y. Saad, Numerical Methods for Large Eigenvalue Problems (Manchester Univ. Press, Manchester, 1992) [ 581 B.N. Parlett, The Symmetric Eigenvalue Problem (PrenticeHall, Englewood Cliffs, NJ, 1980). [59] R.E. Wyatt, Phys. Rev. E 51 (1995) 3643. [60] H.O. Karlsson, J. Chem. Phys. 103 (1995) 4914. [61] W.H. Thompson, H.O. Karlsson, W.H. Miller, J. Chem. Phys. 105 (1996) 5387. [62] R. Chen, H. Guo, J. Comput. Phys. 136 (1997) 494. [63] H.-G. Yu, S.C. Smith, J. Chem. Sot. Faraday Trans. 93 (1997) 861. 1641 J.K. Cullum, R.A. Willoughby, Lanczos Algorithms for Large Symmetric Eigenvalue Computations (Birkhauser, Boston, 1985). [65] A. McNichols, T. Carrington, Chem. Phys. Lett. 202 ( 1993) 464. [66] M.J. Bramley, T. Carrington, J. Chem. Phys. 99 (1993) 8519. 1671 M.J. Bramley. J.W. Tromp, T. Carrington, G.C. Corey, J. Chem. Phys. 100 (1994) 6175. [68] P.-N. Roy, T. Canington, J. Chem. Phys. 103 (1995) 5600. [69] R. Chen, H. Guo, Chem. Phys. Lett. 277 (1997) 199. [70] X.T. Wu, E.F. Hayes, J. Chem. Phys. 107 (1997) 2705. [71] D. Neuhauser, J. Chem. Phys. 93 (1990) 2611. [72] D. Neuhauser, J. Chem. Phys. 95 ( 1991) 4927. 1731 D. Neuhauser, J. Chem. Phys. 100 (1994) 5076. [74] Y. Wang, T. Carrington, G.C. Corey, Chem. Phys. Len 228 (1994) 144. [ 751 V.A. Mandelshtam, T.P. Grozdanov, H.S. Taylor, J. Chem. Phys. 103 (1995) 10074. [76] V.A. Mandelshtam, H.S. Taylor, J. Chem. Sot. Faraday Trans. 93 (1997) 847. [77] H. Kono, Chem. Phys. Lett. 214 (1993) 137. [78] N. Moiseyev, P.R. Certain, F. Weinhild, Mol. Phys. 36 (1978) 1613. 1791 R.E. Wyatt, Adv. Chem. Phys. 53 (1989) 231. [ 801 M.R. Wall, D. Neuhauser, J. Chem. Phys. 102 ( 1995) 8011. [81] V.A. Mandelshtam, H.S. Taylor, J. Chem. Phys. 106 (1997) 5085. [ 821 V.A. Mandelshtam, H.S. Taylor, Phys. Rev. Lett. 78 ( 1997) 3274.

I?. Chen. H. Guo/Cumputer

Physics

[ 831 V.A. Mandelshtam, H.S. Taylor, .I. Chem. Phys. 107 (1997) 6756. [84] E Grossmann, V.A. Mandelshtam, H.S. Taylor, J.S. Briggs, Chem. Phys. Lett. 279 (1997) 355. ] 85 ] J. Main, V.A. Mandelshtam, H.S. Taylor, Phys. Rev. Lett. 79 (1997) 825. 1861 J.W. Pang, T. Dieckman, J. Feigon, D. Neuhauser, J. Chem. Phys. 108 (1998) 8360. [ 87 ] V.A. Mandelshtam, H.S. Taylor, J. Chem. Phys. 108 (1998) 9970. [88] J.W. Pang, D. Neuhauser, Chem. Phys. Lett. 252 ( 1996) 173. [ 891 K.C. Kulander, E.J. Heller, J. Chem. Phys. 69 ( 1978) 2439. [90] S.-Y. Lee, E.J. Heller, J. Chem. Phys. 71 (1979) 4777. [91] E.J. Heller, R.L. Sundberg, D. Tannor, J. Phys. Chem. 86 (1982) 1822. 19’21 D.J. Tannor, E.J. Heller, J. Chem. Phys. 77 (1982) 202. [ 931 H. Guo, T. Seideman, Phys. Chem. Chem. Phys.

Corwmnications

119 (1999)

19-31

31

I (1999) 1265. [94] H.A. Kramers, W. Heisenberg, Z. Phys. 31 (1925) 681. 1951 P.A.M. Dirac, Proc. Royal Sot. London I14 (1927) 710. 1961 S.C. Althorpe, D.J. Kouri, D.K. Hoffman, J.Z.H. Zhang, J. Chem. Sot., Faraday Trans. 93 (1997) 703. [97] S.C. Althorpe, D.J. Kouri, D.K. Hoffman, Chem. Phys. Lett. 275 ( 1997) 173. [98] D. Kouri, M. Arnold, D.F. Hoffman, Chem. Phys. Lett. 203 (1993) 166. [ 991 V.A. Mandelshtam, J. Chem. Phys. 108 ( 1998) 9999. [ 1001 D. Neuhauser, M. Baer. J. Chem. Phys. 90 (1989) 4351. [ 1011 R. Kosloff, D. Kosloff, J. Comput. Phys. 63 (1986) 363. [ 1021 T. Seideman, W.H. Miller, J. Chem. Phys. 96 (1992) 4412. [ 1031 G.-J. Kroes, D. Neuhauser, J. Chem. Phys. 105 (1996) 8690. [ 1041 G.-J. Kroes, M.R. Wall, J.W. Peng, D. Neuhauser, J. Chem. Phys. 106 ( 1997) 1800.