Finite element basis for the expansion of radial wavefunction in quantum scattering calculations

Volume 187, number I,2

CHEMICAL PHYSICS LETTERS

29 November 1991

Finite element basis for the expansion of radial wavefunction in quantum scattering calculations Woonglin Hwang, Yoon Sup Lee Department ofChemistry and CenterforMolecular

Science, Korea Advanced Institute ofScience and Technology,

P.O. Box 150. Cheongryang, Seoul. South Korea

and Seung C. Park Department of Chemistry and Institute ofBasic Science, Kangweon National University, Chunchon 200-701, South Korea

Received 1 May 1991; in final form 2.3August 1991

Radial wavefunctions in quantum scattering calculations are expanded in terms of two shape functions for each finite element. This approach is the R matrix version of Kohn’s variational method and also directly applicable to Smatrix in the log-derivative version. The linear algebra involved amounts to solving definite banded systems. In this basis set method, R matrix or log-derivative matrix is greatly simplified and the computational effort is linearly proportional to the number of radial basis functions, promising computational efficiencies for large scale calculations. Convergences for test vases are also reasonably rapid.

1. Introduction

There have been many studies of performing exact quanta1 scattering calculations in variational methods in recent years. Kouri, Truhlar and co-workers [ l-31 expanded the amplitude density with squareintegrable basis sets using their L*-amplitude density generalized Newton variational principles. McKay and co-workers [ 4,5] have used non-square integrable basis functions for electron-atom and electron-molecule scattering calculations to exactly solve the Schriidinger equation in the asymptotic region by the Schwinger variational principles. Using non-square integrable functions requires more computational efforts than pure L* calculation does. More recently, Miller and co-workers [6,7] applied the Kobn variational principle directly to S or T matrix rather than to the reactance matrix, or equivalently formulated the S matrix version of the Kohn variational principle, using a damping function and a complex basis function to enforce an outgoing boundary condition in each channel. Spurious solutions can be avoided in that method. 180

Manolopoulos et al. [ 8,9 ] have formulated the logderivative version of the Kohn variational principle in which the arbitrary damping function that arises in Miller and co-workers’ formulation [6,7] is not needed. They have used the composite basis set which is categorized into internal functions and boundary functions, both of which are contracted from a primitive radial basis set of Lobatto shape functions to an optimized radial basis. The log-derivative matrix is eliminated by directly constructing the S matrix with the stiffness matrix#’ in their formulation. This S matrix version has several desirable features, such as the automatic unitarity of S matrix and the useful capability to iteratively solve very large quanta1 scattering calculations for which transitions from only a few initial states are required. It is expected that the use of finite element type basis will lead to a simple and convenient way of constructing the Smatrix that Manolopoulos and co-workers have derived, be91 The K matrix in Manolopoulos’ paper is sometimes called the stiffness matrix in finite element method. See for example ref.

[‘Ol.

0009-2614/9 I/$ 03.50 0 199I Elsevier Science Publishers B.V. All rights reserved.

Volume 187,

number I,2

,29 November199 I

CHEMICALPHYSICSLETTERS

cause the local piecewise functions can be divided into internal and boundary functions. In L2-variational method, the matrix element calculation scales as M”, where M is the number of L*basis functions. As a result, the computational effort for large A4are dominated by the matrix operations which scale as M3. The advantage of a variational method is that the size of M required to achieve a given accuracy can be adjusted to be small, since it is essential in large scale calculations to reduce the size of basis set as much as possible. It is also important to devise a new computational method for which computational effort increases more slowly than M3 in L*-variational methods. We present here a new basis set method in quanta1 scattering calculations. Our approach is based on the R matrix version of Kahn variational principle, since our basis set is well suited to this variational method due to the inhomogeneous boundary conditions of the basis functions. Translational basis functions are expanded in terms of the piecewise polynomials that are used in ordinary finite element method. The stiffness matrix can be directly applied to S matrix, and the resulting expression of S matrix is greatly simplified and directly amenable to iterative solution. The most important advantage of our method is that both the matrix element calculation and the matrix operations scale as N, which is the number of L2-translational basis functions, and so does the computational labor. Section 2 describes the theoretical reformulation of the variational R matrix using the finite-elementtype piecewise functions in potential scattering. In section 3, the computational details and the numerical results applied to an elastic and a multichannel scattering are presented.

N- 1, with rN= r, and r, = r,,. All finite elements have equal lengths. In each element (r,+ ,, r;), !P(r) is expressed by a basis of real shape functions u2:s and uZi+ , ‘s localized in the intervals ) ri+ , , r,_ , (, as shown in fig. 1 and given by

y(r)= i [C2iu2r(r)+c2r+Iu2i+L (r)l ,

(1)

uzi(r)=(1-t)2(lt2t),

(2)

Uzi+I(r)=(l-t)‘sgn(r-r,)t,

(3)

i=O

where t= ~(r-r; I h ’

Sgn(r-ri)=-1

ifrr,,

sgn(r-r,)=t

and h is the length of each finite element. cZiand czl+ , are the values of Y’(r,) and its first derivative, respectively.

2.2. R matrix version of Kohn variational principle in elartic scattering We begin here with a simple case of s-wave potential scattering. Any method for solving the potential scattering which uses a global translational basis set is also directly applicable to both inelastic and reactive scattering problems [ 6-81. The Schriidinger equation for s-wave potential scattering is as follows, (H-0 =

p(r)

c-k$+V(r)-fk’

Y(r)=O.

>

(4)

The wavefunction Y(r) should be regular at the origin and has the boundary conditions

2. Theory 2.1. Expansion of radial wavefunction by a finite element basis The reaction coordinate wavefunction, y/(r), is expanded by the finite element functions, which are piecewise cubic polynomials used by Laloyaux et al. [ 1 I 1. The interaction region (I,, rob),where r,< r,, is subdivided into N finite elements (ri+,, r,), i=O to

Fig. 1. Twotypes of piecewise functions used in the expansion of

translaCona1wavefunction at every discretization mesh point.

181

Volume 187, numbqr I,2

Y(O)=O,

CHEMICAL PHYSICS LETTERS

Y(r-m)=I(r)-O(r)S,

(5)

where I(r) and O(r) are the flux normalized incoming and outgoing waves, respectively, and S is the scattering matrix, eris. The boundary condition at r=O leaves cl1to be zero. In this single channel case, the variational functional is m Y=-

Y(H-E) Fdr ,

s 0

Y(r)= f:

r=l

(6)

C,&(r).

(7)

The standard formula for the single channel R-matrix is well known [ 12] : (8)

where nl A,=

I

[ fujuJ+u,( I’-jk2)~j]

dr.

November1991

ment basis set is valued one, and all other basis functions are zero since they are localized in finite elements uO(ro)=l

and

uj(ro)=O, i#O.

(10)

When we substitute eq. ( 10) into eq. (8), the R matrix is simplified by R=;(A-I),,,

.

(11)

This simple form of R matrix can be exploited to a great computational efficiency since only one inverse matrix element need to be calculated. Additional and more significant advantage is in the bandedness of the A matrix. The computational effort required for the general matrix inversion scales as N3, but in the sparse matrix with a definite bandwidth, the required operations can be made proportional to N [ 161 as demonstrated by the numerical results of section 3. The present method has a great computational efficiency in comparison with other Lz-variational methods for which computational efforts increase as N 3.

(9)

II

One of the advantages of this numerical method against the analytical basis function expansion is the easiness of the integration of eq. (9), regardless of the shape of potential. The symmetric A matrix in eq. (8) corresponds to the stiffness matrix in the usual finite element method [ 131 #*. The Kohn’s derivation places no restriction on the basis functions ui( r) other than the linear independence and the completeness of the basis set. The inherent completeness of the present basis set in the interval 0 ~r
cinities of the eigenvalues as discussed in section III of ref. [141.

182

29

2.3. Multichannel scattering case In the multichannel case, the variational R matrix isgiven by [ 121 (12)

y(r)=

C

W(r0)

MY)

.

(13)

Similarly to the elastic scattering case of eq. ( 11) eq. ( 12) is greatly simplified by the relation R= L (A-‘)oo. 2P

(14)

As in the elastic scattering case, eq. ( 14) offers a significant computational saving other than that obtained from the bandedness of the matrix A&, since only a small portion of the inverse matrix elements of A are required in the calculation. The dimension of the matrix that should be actually inverted corresponds to the number of open channels. As was shown by Manolopoulos and co-workers’ formulation [ 141, the scattering matrix can be solved

Volume 187, number 1,2

iteratively for each column of the open channel in the log-derivative version using L_‘=R=(A_‘)oo=(Aoo-AosA~~~so)_’

)

Table 1 Fractional error in tans (accurate value is - 1.7449393)and CPU time at k=O. 15 for s-wave scattering by an exponential potential

(15) Node s,

where the index B is for all other basis functions except uo, and S=a+2iBT[A-Cl-‘B,

(16)

where u, B and C are identical to those in Manolopoulos and co-workers’ paper [ 141. From eq. ( 16), the symmetric property of A matrix guarantees the unitarity of the open channel scattering matrix So,.

CPU time (s)

Fractional CITOT

10 I5 20 25 30 60

gjr b1

Linpack c1

band d,

0.106 0.316 0.703 I .304 2.176 16.255

0.023 0.043 0.07 I 0.111 0.151 0.559

0.009 0.019 0.022 0.028 0.030 0.058

0.0010 0.0002 0.0001 0.0000 0.0000 o.ooocl

a) Node corresponds to the number of finite elements. b, pir is a Gauss-Jordan elimination algorithm. c) Linpack is DSIFA and DSISL in LINPACK. d, band is our modified algorithm.

3. Test calculations and discussion The A matrix is very sparse. For example, in the multichannel scattering case of 30 nodes of translational coordinate, the dimension of A matrix becomes 60n, where n is the number of internal basis functions. But most of the matrix elements are zero, and the upper bandwidth is 4n. In order to manipulate such a large banded matrix efficiently by eliminating unnecessary operations on zero elements, we have modified DSIFA and DSISL in LINPACK package [ 161. The diagonal pivoting of DSIFA is also used in our program so that the upper bandwidth is increased to 6n, and not to 8n. However our program performs operations only on nonzero elements to get a maximum efficiency such that the required operations are not increased by unnecessary operations on zero elements that may result from pivoting. The size of core memory to store the matrix elements is reduced substantially in band storage mode as s(n*N), where n and N are numbers of the internal and translational basis functions, respectively, and hereafter these notations will be used in this paper. Test calculations of both elastic and inelastic scatterings have been carried out on SUN3/ 180. 3. I. Elastic scattering As a test case, we consider the same potential as used in previous studies [ 6,17,18 ], V(r) = -e-’ .

29 November 199I

CHEMICAL PHYSICS LETTERS

(17)

Tables 1 and 2 show s-wave potential scattering results. The values of CPU time are averaged out from

Table 2 Fractional error in tan S (accurate value is 2.2003827) and CPU time at k=0.55 for s-wave scattering by an exponential potential Node ‘I

10 15 20 25 30

CPU time (s)

Fractional error

gir”’

Linpack ”

band d’

0.1 I1 0.317 0.7 1.3 2.187

0.025 0.047 0.075 0.114 0.154

0.012 0.016 0.02 0.023 0.026

O.OOL3 0.0002 0.0000 0.0000 0.0000

‘) Node corresponds to the number of finite elements. b’ gjr is a Gauss-Jordan elimination algorithm. ‘1 Linpack is DSIFA and DSISL in LINPACK. ” band is our modified algorithm.

20 repeated runs, since the time fluctuates slightly. We have used three matrix inversion subroutines in order to compare amounts of computational labor. gjr denotes a Gauss-Jordan elimination algorithm with partial pivoting to evaluate the whole inversematrix elements. We let Linpack denote DSIFA and DSISL that calculate only one inverse-matrix element according to eq. ( 11). band corresponds to our modified algorithm of calculating only one inversematrix element. The simple form of R matrix in eq. ( I1 ) affords a lot of time savings, as can be seen from the comparison of gjr’s and Linpack’s results. While gjr calculates the whole inverse-matrix elements, Linpack is made to calculate only one element. The fact that gjr and Linpack may differ slightly in computational 183

Volume 187, number 1,2

CHEMICAL PHYSlCS LETTERS

efficiency has no important bearing in our main stream of argument. More practical and significant savings are achieved using our band algorithm since the CPU time grows linearly as the number of translational basis functions that is actually twice the number of nodes. Practical usefulness of the present method in multichannel scattering case will become clear for large scale calculations, and one example is described in section 3.2. 3.2. Inelastic scattering The purpose of developing any basis set method is for the applications in inelastic and reactive scattering. Here we consider the standard benchmark calculation of Secrest-Johnson [ 6,191 model for collinear Het H2 vibrationally inelastic scattering. All quantities are given in the same units as those of Secrest and Johnson, as shown in table 3 and 4. We

29 November

have used our band algorithm for the matrix operations of eq. ( 14). The CPU time increases also very slowly as N, as will be shown in the results. The first application is for E=6 and the other quantities are given in table 3. There are three open channels, and a total of five vibrational states in this energy is sufficient to obtain a converged result similarly as in the previous study [ 61. The numerical results do not suffer from the numerical instability of exponential growth in closed channels, which is an important advantage of the finite element method [ 201. Table 3 shows three inelastic vibrational transition probabilities, Pa,, Pazand P,z. It can be found that the CPU time scales as N. Good convergence is obtained such that the transition probabilities are reasonably well converged with 40 translational basis functions (20 nodes), Computational accuracy to the limit of numerical results [6] can be easily achieved in the present method due to the N dependency of CPU

Table 3 Inelastic probabilities P,,, fez and P,,and CPU timeas a function ofthe number offiniteelements Vo= 12 and n= 5. Our band algorithm is used for the matrix operations

(node). E=6.0,

y=O.315, m=0.667,

CPU time

Node

PO,

PO2

IO I5 20 25 30 35 40 60

0.28185( -1) 0.30326 ( - 1) 0.30013(-l) 0.30033 ( - 1) 0.30047 ( - 1) 0.30053( -1) 0.30056( - 1) 0.30059( - 1)

0.96906( -5) 0.11377(-4) O.ll246(-4) 0.11259(-4) 0.11266( -4) 0.11269(-4) 0.11270(-4) 0.11272(-4)

0.12911(-2) 0.14719( -2) 0.14706( -2) 0.14715( -2) 0.14719( -2) 0.14720( -2) 0.14721(-2) 0.14721(-2)

008’

0.30058

0.11273(-4)

0.14723( -2)

( - 1)

1991

(s)

0.36 0.52 0.70 0.96 1.18 1.2 1.64 2.22

a) Numerical resultsof ref. [6].

Table 4 Inelastic probabilities PO,, P,4, Pz3 and Pzs and CPU time as a function of the number of finite elements ~~0.667, V,= 20.0 and n = 8. Our band algorithm is used for the matrix operations Node

PO,

Pi4

PZJ

P25

10 20 30 40 50 60 80 COB’

0.45559( -6) 0.34639 0.39660 0.39457 0.39433 0.39426 0.39423

0.49211(-6) 0.60587( -3) 0.38980( -3) 0.38988(-3) 0.38976( -3) 0.38973( -3) 0.38971(-3) 0.38971(-3)

0.29690 0.26740 0.23352 0.23344 0.23343 0.23343 0.23343 0.23345

0.26235( 0.22253( 0.19624( 0.19627(

,0.39423

‘) Numerical results of ref. [6].

184

(node). E= 12.0, y=O.3,

CPU time (s) -5) -5) -5) -5)

1.0 2.02 3.02 4.06

0.19628( -5) 0.19628( -5)

5.06 6.16 8.t2

0.19628( -5) 0.20137( -5)

Volume 187, number 1,2

CHEMICAL PHYSICS LETTERS

Miller and Jansen op de Haar [6] have obtained converged results with 18 translational basis functions. Although we are unable to directly compare computation time with other methods, we believe that the great reduction of linear algebraic operations in eq. ( 14) and the N dependency of computational effort of the present method will prove to be advantageous in many foregoing applications. The linear dependence of CPU time on N in our method is comparable to that on the number of sectors in close coupling methods, but close coupling method probably requires more involved calculations in each sector than our method. Table 4 summarizes the results for the higher energy E= 12.0. Results are fully converged with eight vibrational basis functions. Low order polynomials used in finite element method are poor basis functions for highly oscillating wavefunctions that may occur in the asymptotic region, but in the close-in and nonclassical region, they will show good convergence. The convergence rate in table 4 (E= 12.0) is slightlylower than that in table 3 (E~6.0) as expected, but it is satisfactory within the context of computational efficiency. time.

4. Conclusions The idea of using the finite element basis for the translational wavefunction in R matrix, or log-derivative version of the Kahn’s variational principle has several specific advantages that make it more convenient than other methods. ( 1) The use of finite element basis removes the complexity associated with choosing radial basis functions in molecular dynamics such as contracting primitive basis functions to an optimized radial basis set, and is well applicable to the various problems in molecular dynamics regardless of the forms of both the potential energy function and the kinetic operator. (2) The simple form of piecewise polynomial expansions for the translational function leaves the integration, or the matrix element calculation elementary, which leads to the easy programming. (3 ) The boundary condition of this basis given by eq. ( 10) permits eqs. (8) and ( 12) to be reduced in a simple and convenient manner to eqs. ( I 1) and ( 14), respectively. This allows us to calculate only quite a few elements

29 November

I99

I

in the inverse matrix of A, resulting in a great computational efficiency. (4) In multichannel scattering with a total of M basis functions, the computational effort of calculating matrix elements scales as M2, and that of the linear algebra scales as M3. Thus the matrix operations, and not the integrals, prevail the computational effort in conventional large scale L*variational methods. By employing the piecewise local basis functions, the time of the linear algebra scales as 0 (n ‘N) , where 12and N are the number of internal and radial basis functions, respectively. (5) The matrix to be inverted is real, energy independent and symmetric, which guarantees a unitarity of S matrix, and this method is also directly amenable to iterative solutions in the log-derivative version. This last point is not specific to the present method but the general consequence of the formalism utilized. This basis set method has good convergence in the interaction region, where the potential has some value. This fact will be clearly understood when we consider that finite element method aims to construct the global function which need not be as oscillatory [ 201. The basis set expansion may encounter difficulty in highly oscillating regions that occur in the asymptotic region. Hence, the most efficient procedure appears to be a hybrid algorithm that incorporates the basis set method and-the usual close coupling methods in different regions of space [ 2 11. An extension of the present method may be useful even in reactive scatterings, by employing the hyperspherical polar coordinate in exchange region, since exchange integral in Miller’s formulation [ 221 can be avoided due to the incorporation of channel orbitals under one hyperradius. The work in this direction is in progress.

Acknowledgement The research is in part supported by Korea Science and Engineering Foundation and by the Basic Science Research Institute program of the Ministry of Education of Korea. Authors are grateful to the referee for critical comments and one of us (SCP) would like to thank Dr. David Manolopoulos for valuable comments.

185

Volume 187, number 1,2

CHEMICAL PHYSICS LETTERS

References [ ] K. Haug, D.W. Schwenke, D.G. Truhlar, Y. Zhang, J.Z.H. Zhangand D.J. Kouri, J. Chem. Phys. 87 (1987) 1892. [2] D.W. Schwenke, K. Haug, M. Zhao, D.G. Truhlar, Y. Sun, J.Z.H. Zhang and D.J. Kouri, I. Phys. Chem. 92 (1988) 3202. [3] C. Yu, Y. Sun, D.J. Kouri, P. Halvick, D.G. Truhlar and D.W. Schwenke, J. Chem. Phys. 90 ( 1989) 7608. [4] R.R. Lucchese, K. Takatsuka and V. McKay, Phys. Rept. 131 (1986) 147. [5] R.R. Lucchese, D.K. Watson and V. McKay, Phys. Rev. A 22 (1980) 421. [6] W.H. Miller and B.M.D.D. Jansen op de Haar, J. Chem. Phys. 86 (1987) 6213. [7] J.Z.H. Zhang, S.-l Chu and W.H. Miller, J. Chem. Phys. 88 (1988) 6233. [S] D.E. Manofopoulos and R.E. Wyatt, Chem. Phys. Letters I52 (1988) 23. [9] D.E. Manolopoulos, M. D’Mello and R.E. Wyatt, J. Chem. Phys. 93 ( 1990) 403. [IO] 0. Axelsson and V.A. Barker, Finite element solution of boundary value problems (Academic Press, New York, 1984).

186

29 November 1991

[ I I] Th. Laloyaux, Ph. Lambin, J.-P. Vigneron and A.A. Lucas, J. Comput. Phys. 83 ( 1989) 398. [ 121R.K. Nesbet, Variational methods in electron-atom scattering theory (Plenum Press, New York, 1980). [ 131N. Sato and S. Iwata, J. Comput. Chem. 9 ( 1988) 222. [ 141D.E. Manolopoulos, M. D’Mello and R.E. Wyatt, J. Chem. Phys. 9 1 ( 1989) 6096. [ 151D.J. Zvijac and J.C. Light, Chem. Phys. 12 (1976) 237. [ 161J.J. Dongarra, C.B. Moler, J.R. Bunch and G.W. Stewart, LINPACK user’s guide (SIAM, Philadelphia, 1979) ch. 2. [ 171G. Staszewska and D.G. Truhlar, Chem. Phys. Letters I30 (1986) 341. [ 181C.W. McCurdy, T.N. Rescigno and B.L. Schneider, Phys. Rev. 36 (1987) 2061. [ 191D. Secrest and B.R. Johnson, J. Chem. Phys. I5 ( 1966) 4556. [20] A. Askar, A.S. Cakmak and H.A. Rabitz, Chem. Phys. 33 (1978) 267. [ 211 H. Rabitz, National resource for computation in chemistry, NRCC Proceedings, Vol. 1 (Argonne National Laboratory, Argonne, 1979) p. 144. [ 221 W.H. Miller, J. Chem. Phys. 50 ( 1969) 407.