J . Phys. Chem. 1988, 92, 3202-3216
3202
Quantum Mechanical Algebraic Variational Methods for Inelastic and Reactive Molecular Colllsions David W. Schwenke, Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455, and NASA Ames Research Center, Moffett Field, California 94035
Kenneth Haug, Meishan Zhao, Donald G. Truhlar,* Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455
Yan Sun, John Z. H. Zhang, and Donald J. Kouri Departments of Chemistry and Physics, University of Houston, Houston, Texas 77004 (Received: August 17, 1987)
The quantum mechanical problem of reactive or nonreactive scattering of atoms and molecules is formulated in terms of square-integrable(L2) basis sets with variational expressions for the reactance matrix. We present and test several formulations, involving expansions of the wave function (the Schwinger variational principle) or amplitude density (a generalization of the Newton variational principle), single-channel or multichannel distortion potentials, and primitive or contracted basis functions. The test results, for inelastic and reactive atom-diatom collisions, are very encouraging, and we anticipate that the methods will be very useful for a variety of collision calculations and that they will allow the accurate quantal treatment of systems for which other available methods would be prohibitively expensive.
1. Introduction Variational approaches are the most powerful general techniques for computational applications of quantum mechanics. For bound-state problems the expansion of the wave function in an L2 (Le., square integrable) basis leads to a familiar matrix eigenvalue problem which may be solved in two distinct steps: evaluation of the matrix elements and diagonalization of the matrix. Specialized techniques for both steps have been developed and improved by many workers over a long period of time. Expansion of scattering wave functions in a basis set also leads to two steps but also to new difficulties. One class of variational principles for scattering problems is based directly on the differential Schroedinger equation with nonhomogeneous boundary conditions. The most straightforward of these variational principles, due originally to Hulthtn and Kohn,’ requires non-square-integrable basis functions, leading to more difficult integrals in the first step; also, because of the nonhomogeneous boundary conditions, the second step is transformed from an eigenvalue problem to a set of linear equations that may have spurious solutions when certain of the matrices are singular. These problems and other aspects of the Kohn and related approaches are reviewed elsewhere.* Very recently the Kohn variational method has been formulated using a complex basis function to enforce an outgoing wave boundary condition in each channel, and it has been shown that spurious solutions can be avoided in this way.3 Another class of variational methods for scattering problems is based on integral equations. Consider, e.g., the Schwinger variational m e t h ~ d . ~In~this . ~ method ~ the scattering boundary (1) Kohn, W. Phys. Rev. 1948, 74, 1763. Hulthen, L. Ark. Mat. Astron. Fys. 1948, 35AI25.. (2) (a) Truhlar, D. G.; Abdallah, J., Jr.; Smith, R. L. Adu. Chem. Phys. 1974, 25, 21 1. (b) Callaway, J. Phvs. Reo,. 1978, 45, 89. (c) Nesbet. R. K.
Variational Methods in Electron-Atom Scattering Theory; Plenum: New York, 1980. (3) Miller, W. H.; Jansen op de Haar, B. M. D. D. J . Chem. Phys. 1987, 86, 6213. (4) Schwinger, J. Phys. Rev.1947, 7 2 , 742. Lippmann, B. A,; Schwinger, J. Phys. Rev. 1950, 79, 469. Kato, T. Prog. Theor. Phys. 1951, 6, 295. Lippmann, B.-A. Phys. Rec. 1956, 102, 264.
0022-3654 I88 12092-3202$01.50/0
conditions are enforced by Green’s functions and correct results can be obtained by expanding the wave function in a real, L2basis. Furthermore the approximate reactance matrix obtained from a Schwinger variational calculation may be interpreted as the exact reactance matrix for an approximate potential which does not yield spurious ~ o l u t i o n s . ~ The ~ ~Schwinger method has been applied in nuclear physics,5 for electron-atom, positron-atom, and electron-molecule scattering: and for collinear atom-diatom scattering,’ but not for three-dimensional atom-molecule scattering. Since the computational advantages of this kind of approach deserve exploration, in the present paper we report a formulation of this method for three-dimensional atom-molecule collisions, and we test its efficiency by numerical applications. We also consider a related variational method, which is a generalization of a variational principle first given by Newton for single-channel problem^,^ and which has t h e same advantagesI0 as mentioned for the Schwinger method but which has been shown in model problems to be more rapidly convergent.l’ The generalized ( 5 ) Blatt, J. M.; Jackson, J. D. Phys. Rev. 1949, 76, 1-8. Gerjuoy, E.; Saxon, D. S. Phys. Rev. 1954, 94, 478. Joachain, C. J. Nucl. Phys. 1965, 64, 529, 548. Hiifner, J.; Shakin, C. M. Phys. Reu. 1968, 175, 1350. Hiifner, J. Lemmer, R. A. Phys. Rev. 1968, 175, 1394. Zubarev, A. Sou. J . Part. Nucl.
(Engl. Transl.) 1976, 7 , 215. (6) Altshuler, S. Phys. Rev. 1953, 89, 1278. Zubarev, A. L. Sou: J . Nucl. Phys. (Engl. Transl.) 1978, 28, 287. Belyaev, V. B.; Podkopaev, A. P.; Wrzecionko, J.; Zubarev, A. L. J. Phys. B. 1979, 12, 1225. Lucchese, R. R.; McKoy, V . Phys. Rev. A 1980, 21, 112. Watson, D. K.; Lucchese, R. R.; McKoy, V.; Rescigno, T. N.; Phys. Rev. A 1980, 2/, 738. Maleki, N.; Macek, J. Phys. Rev. A 1980,21, 1403. Takatsuka, K.; McKoy, V. Phys. Rev. Lett. 1980.45, 1396. Phys. Rev. A 1981,24,2473. 1984,30, 1734. Lima, M. A. P.; Gibson, T. L.; Takatsuka, K.; McKoy, V. Phys. Rev. A 1984, 30, 1741. Gibson, T. L.; Lima, M. A. P.; Takatsuka, K.; McKoy, V. Phys. Rev. A 1984, 30, 3005. (7) (a) Miller, W. H. J . Chem. Phys. 1969, 50, 407. (b) Adams, J. E.; Miller, W. H. J . Chem. Phys. 1979, 83, 1505. (c) Hermann, M. R.; Miller, W. H. Chem. Phys. 1986, 109, 163. (8) Adhikari, S . K.; Sloan, I. H. Phys. Reu. C 1975, / I , 1133. Zubarev, A. L. Sou. J . Nucl. Phys. (Engl. Trans/.)1976, 23, 40. (9) Newton, R. G. Scattering Theory of Particles and Waues;McGrawHill: New York, 1966; Section 11.3. (10) Sloan, I. H.; Brady, T. J. Phys. Rev. C 1972, 6, 701. Brady, T. J.; Sloan, I . H.; Phys. Reu. C 1974, 9, 4. ( 1 1) Staszewska, G.; Truhlar, D. G.Chem. Phys. Lett. 1986,130, 341. J . Chem. Phys. 1987, 86, 2793.
0 1988 American Chemical Societv
The Journal of Physical Chemistry,.Vol. 92, No. I I, I988 3203
Variational Methods for Molecular Collisions Newton method is formulated here as a direct extension of our recent applicationi2 to inelastic and reactive molecular collisions of the method of moments for the translational radial amplitude density. The generalized Newton method and the Schwinger method are applied to the same problem, and the Newton method is also used to compare the convergence characteristics of two possible choices for the distortion potential. Section 2 presents the various methods for inelastic atomdiatom scattering; and section 3 shows how they may be extended to treat reactive scattering. Section 4 presents computational details, section 5 presents test results, and section 6 contains discussion, including a comparison to methods based on differential equations. Section 7 is a summary.
2. Theory: Inelastic Scattering In this section we consider nonreactive collisions A the Hamiltonian H =T
+ BC with
+ VIb+ Pnt
+F
(2)
and the distorted-wave Hamiltonian is H D = T + V i b +r"
(3)
Then the full Hamiltonian may be written
H=H~+P
Substituting the expansion of eq 12 into eq 11 yields
(5)
where p is the reduced mass for relative translational motion. (Note: p will be redefined for reactive scattering in section 3.) Then the full Schroedinger equation with its boundary conditions may be replaced by the Lippmann-Schwinger equation \kQ = $Q
+ GD'U\kno,
no = 1, ..., N
(6)
where no denotes the initial channel, $Qis the regular standingwave solution of the distorted-wave problem
( H D- E ) P = 0
where agnis a real coefficient and 9, is a real square-integrable ( L 2 )basis function. Note that since eq 11 always includes \kn or \k"o multiplied by 'U, the basis functions 9, need be nonzero only when U\knor U\k" is nonzero. This can be used to advantage since, by definition, 'U differs significantly from zero only when ptdiffers significantly from r".We determine the coefficients apnby requiring that when the expansion of eq 12 is substituted into eq 11 the results are stationary with respect to variations in the coefficients; Le., we require
(4)
The distorted-wave principal-value Green's function is
GD = (-h2/2p)P(E - P)-'
This yields the exact reactance matrix if \k" and \ k n o are exact, and it is stationary with respect to small errors in these wave functions if they are not exact.9 To proceed we expand the wave function as
(1)
where T is the kinetic energy, VVib is the vibrational potential of BC, and ptis the interaction energy of A with BC. The interaction potential is partitioned into a distortion potential and a coupling potential p by
V"t = r"
matrix form rather than reactance matrix form. This has the disadvantage that real algebra is replaced by complex algebra but the advantage that one can calculate physical observables from individual elements or columns of the transition matrix, whereas this requires the whole reactance matrix. 2.1. Schwinger Method. The Schwinger variational principle2c,4-9*11 for the reactance matrix can be written in bilinear form as
(7)
and
$ = bTa + aTb - aTca
(14)
where a superscript T denotes a matrix transpose,
and c,p' =
(@,plp,y)- (9,1'UGDUpp')
(16)
Invoking the stationary requirement of eq 13 yields a = c-lb
(17)
Substituting eq 17 into eq 14 and cancelling terms gives
'U = (-2p/h2)P
(8)
The reactance matrix K for the full problem may be calculated from Knn, = 'Knna + Xnno
(9)
where n denotes the final channel, OK is the reactance matrix due to the distortion potential, and Xnn, = =
(Vl~l*"O)
(*nl'UIP)
(loa)
(1Ob)
The simplicity of the relations 10a and 10b is a consequence of the normalizations used for various quantities in eq 6. These will be specified in detail in section 4. By use of different boundary conditions on the Green's function, the equations in this article may easily be written in transition (12) (a) Haug, K.; Schwenke, D. W.; Truhlar, D. G.; Zhang, J.; Kouri, D. J. J . Phys. Chem. 1986, 90, 6757. (b) Zhang, J.; Kouri, D. J.; Haug, K.; Schwenke, D. W.; Truhlar, D. G. J. Ckem. Pkys. 1988,88,2492. (c) Haug, K.; Schwenke, D. W.;Truhlar, D. G.; Zhang, Y.; Zhang, J. Z. H.; Kouri, D. J. J. Chem. Phys. 1987.87, 1892. (d) Zhang, J. Z. H.; Zhang, Y.; Kouri, D. J.; Garrett, B. C.; Haug, K.; Schwenke, D. W.; Truhlar, D. G. Faroday Discuss. Ckem. SOC.,in press.
$ = bTc-'b
(18)
The row and column orders of the matrices in eq 18 are N X N for 8,M X N for b, and M X M for c, where N is the number channds included in the calculation and M is the number of L2 basis functions. Once the matrices b and c have been formed, the task demanding the most computational resources will be the formation of c-lb since this step will scale as M3 for large enough M . Furthermore the work in forming b and c scales as M N and A@,respectively, so the formation of c-'b will always dominate in the limit of large M . Note that the requirement to solve linear equations of order M is not a consequence of using a variational principle-such operations are also required for nonvariational calculations involving expansion in a basis of size M . The advantage of a variational method is that it usually reduces the size of M required for a given accuracy in 8 or K for fixed N . The disadvantage is that the integrals involved in the formation of the required matrices, b and c in this case, are more complicated in the variational formu1ations.l' However, as already pointed out, the matrix element calculation scales only as A@,and so for large enough M the time required for this step should become negligible compared to the other operations involved, and variational methods should always be more efficient for large enough M . Note also that the matrix c in (16) is explicitly symmetric; this fact can also be utilized to save computer time and storage resources in the
3204
Schwenke et ai.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
calculations. In addition this symmetry ensures that K is symmetric and hence that the scattering matrix is unitary, so that the probability flux is conserved and microscopic reversibility is ensured at all stages of convergence. (In practice, of course, the symmetry may be lost if numerical errors are introduced in the Green's function or the integrations.) We note that symmetry of K was not guaranteed in the nonvariational methods we employed in previous papers.12 2.2. Generalized Newton Variational Principle. The variational principle introduced by Newton for the transition (T) matrix has been discussed in various contexts by many author^.^-",'^-^^ In the present article we use it in reactance matrix form with a distorted-wave Green's function as a multichannel variational principle for the amplitude density. For nonreactive scattering the amplitude density is defined by'*
1"" = U9"O
3. Theory: Reactive Scattering In this section we consider the extension of the theory of section 2 to rearrangement collisions. We will present the extension explicitly for the case of the formulation of section 2.2 and three arrangement channels, but the method of section 2.1 may be extended in precisely the same way, if desired, and the case of two arrangements is an obvious subcase. We begin by defining an arrangement channel quantum number a which takes on the values a = 1 for A + BC, a = 2 for B + CA, and a = 3 for C + AB. For each arrangement we let R, denote the mass-scaled radial separation between the atom and the center of mass of the diatom. Each arrangement also has its own V " l b , Vnt, p, fl,p ,and GD which we denote by adding a subscript a. Then the Hamiltonian is conveniently expressed as
(19)
We shall directly expand the amplitude density
-
H=H:+$
(28)
where is a coupling potential that vanishes when R, a. Each channel will be specified by a unique collective quantum number n; Le., n runs 1 to ii, for a = 1, R , + 1 to f l , + tiz for a = 2, and ii, ii2 + 1 to iil + ii2 + ii3 for a = 3. Since each channel is assigned a unique n, the specification of a is redundant. However, we will include (Y except for quantities involved in matrix computations because its specification clarifies the nature of the arrangement coupling. We expand the full scattering wave function for initial arrangement a. and initial channel no as
+
where A,, like is real, and as a consequence, we need not deal directly with the wave function. To obtain an integral equation for the amplitude density we multiply the Lippmann-Schwinger equation (6) by 21,which yields @ff,
= U p + UGDp
(21)
A variationally correct expression for the reactance matrix in terms
of the amplitude density may then be written as9-I1 Xnn, = Xfin, +
(VIUGDIP) + (PIGDuIP)- (PIGDIP) + (PIGDUGDIP)(22)
where \ k p is an arrangement component. We use mass-scaled coordinates in which all kinetic energy terms are associated with the same reduced mass p, which is independent of arrangement.l* Using this reduced mass, we write the Schodinger equation as ( E - H:)\k/o"o
where
= (-h2/2~)XU,,\kzV'0 a'
%fino
= (VIWIVo)
(23)
where ?fa,
Substituting (20) into (22) yields
g = p + BTA + ATB - ATCA
(24)
= (-2p/h2)[H - E
(25)
+ 6,,,(E - HE)]
= CU,,'*/V'O
(32)
n'
We also define pnas the solution of (H: - E)#', = 0
where
Bo, = (@plGD%lV) Copt = ( @BIGDI@o,)- ( @pIGDUGDI@,y)
(26) (27)
As with the Schwinger variational method the linear-equations coefficient matrix, in this case C, is symmetric so that g is formally symmetric even for finite basis sets. In practice, as discussed with respect to the Schwinger method above, to the extent that the Green's function or the integrations are approximate, the symmetry of eq 25 may be lost.
(33)
with initial arrangement a and initial channel n and G : as the principal-value Green's function of this Schmedinger equation normalized as in (5) except with the new definition of p. The solution of (30) can then be written \ k y o = 6naolpno
+ c;po
(34)
and the full reactance matrix may be written Knn, = 6 a a ~ K ~ n+o X n n ,
(35)
where Xnn, =
(13) Sloan, I. H.; Adhikari, S. K. Nucl. Phys. A 1974, ,4235, 352. 1975, A241, 429. 1975, A254, 545(E). (14) Rabitz, H.; Conn, R. Phys. Rev. A 1973, I , 577. Conn,R.; Rabitz, H. J. Chem. Phys. 1974,61,600. Rabitz, S.; Rabitz, H. J . Chem. Phys. 1977, 67 - , 2964 -- - .. (15) Kouri, D. J. J. Chem. Phys. 1973, 58, 1914. Kouri, D. J.; Levin, F. S. Phys. Rev. C 1975, 11, 352. (16) Takatsuka, Ki., McKoy, V. Phys. Rev. A 1981, 23, 2352. (17) Kuruoglu, 2.C.; Micha, D. A. J . Chem. Phys. 1980, 72, 3327. (18) Johnson, B. R.; Secrest, D. J. Math. Phys. (N.Y.)1966, 7 , 2187. Secrest, D. In Methods In Computational Physics; Alder, B.; Fernbach, S.; Rotenberg, M., Eds.; Academic: New York, 1971; Vol. 10, p 243. See also: Kouri, D. J. J . Chem. Phys. 1969.51, 5204. Baer, M.: Kouri, D. J. Phvs. Rec. A 1971, 4, 1924. J . Chem. Phys. 1972, 56, 1758
(31)
and we define a reactive amplitude density in arrangement channel a for initial channel no asI2
Analogously to eq (1 8), the variationally correct expression for $ is now
x- = x* - + BVB
(30)
(PnICYo?
(36)
As in section 2, the fact that the coefficient of X,,, in (35) is unity follows from the boundary conditions chosen for the various functions. A variational principle analogous to Newton's for the nonreactive case is given by Xnno
= X:no -
+
C($anIUna,G,41~po) + C( a'
a'
C(C?l@(C*) + a'
C~IG~,Q~~~~~PO~
( ~ ~ l G E , , ? f a , , a , C t ,(l 3~7~) o ? a'"'
where XnBno=
(vnl?faaolv~n~) (38)
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3205
Variational Methods for Molecular Collisions
As in the nonreactive case, this expression is symmetric. We note that 2r h2
'Un&@"') = - -(K - H:t)IP'")
(39a)
and
= 'U*'G,4 - 1
+ ,6
(39b)
Therefore, instead of ?faat, defined by (31), we may use an arrangement-dependent analogue of the 'U of section 2. This analogue is defined by
'U" = -(2r/hZ)(H -
where A is a constant basis set parameter. For a given A the Gaussian width parameter is determined by the prescription of Hamilton and Light;20Le., the widths are adjusted such that the overlap of two neighboring Gaussians is equal to exp(-$/2), where c is a numerical parameter chosen for good convergence. Since the grid is evenly spaced this yields uniform Gaussian widths. We will also perform calculations with functions defined as linear combinations of primitive Gaussian basis functions. We denote these as A,,(R), and we call them contracted basis functions? We expand these as
e)
(40)
X,,(R) = CZ;,,,,X,,(R)
There are many possible contraction schemes. In the present calculations we determine the contraction coefficients E",, by the following procedure. The conventional close coupling trial function is
and the resulting variational principle is
Q" =
x,, = x;, + C[(J/""l'U"G,4Icp"o) + (c?IG,4'Uqpflo)] a'
C [ (e.?lC:,tl@"O)- (~.?lG~,'U"'G:,(c~o)] (41)
R'CFn,(R)
Fnno(R) = CGnnpaa,Xmpp(R) B
x;, = ( ~ a " I ' U q p * )
(42)
(49)
+,(XI
n
and comparing this to eq 12 and 46 with t , replaced by
a ' " '
where (38) now becomes
(48)
m'
x,
yields (50)
Thus, in order for the $, to be an efficient basis they should be able to represent the F,,, as well as possible. But the radial functions F,, solve the standard close-coupling equations
The equivalence of expressions 37-38 and 41-42 follows immediately from (31), (33), and (39)-(40). We now make the expansion
CO=
(43)
B
where AB, is a real coefficient, and 9, denotes a real L2basis function. Note that, just as each channel has a unique n, each basis function has a unique p, so the sum in eq 43 runs from 1 to M ,for a = 1, from M I 1 to MI M2 for a = 2, and from M I M2 1 to M for a = 3. Thus each /3 is associated with a unique a, which we call aB. The expansions (43) yield an equation of the form (24) with the complication that some integrals now involve functions expressed in two different coordinate systems:
+
+
+
+
B , , = ( 9,IG:p"IJI"'"')
(44)
where 1, is the orbital angular momentum in channel n, k,Z is the wave vector squared given by
k,Z = 2 r ( E - t , J n ) / h 2
(52)
tunin is the internal eigenenergy of channel n, which has vibrational quantum number v, and rotational quantum number j,, and unnf =
-(2P/h2)Vnn,(R)
(53)
where V,, = $dx +,*(x)I/int(R,x)+,,,(x)
(54)
We then generate functions Xmn and associated eigenvalues K by diagonalizing the operator
~ :
The stationary reactance matrix is again given by eq 25.
4. Basis Functions and Computational Details In this section we make our methods more concrete by specifying in detail our computational methods. 4.1. Nonreactive Scattering. 4.1 .I. Basis Functions. First consider the basis functions 9, for the Schwinger variational calculations. We have purposely written the formalism of sections 2 and 3 in a very general way with regard to basis sets because there are many possible choices for these basis functions, and it may be critical for success in demanding applications to optimize this choice. Here we will consider only basis functions of the general form * ~ ( R J )= R-'tm8np(R)+ n 8 ( X )
(46)
where R denotes the mass-scaled radial separation between A and the center of mass of BC, x denotes the collection of other coordinates, & ( x ) is a standard vibrational-rotational-orbital function,lg and t,,(R) is a translational radial basis function. For the nonreactive scattering example in this paper, we will consider two choices of t,,(R). First, as before,12 we take t,,(R) to be a Gaussian centered at point RZ. This is called a primitive basis function and is denoted A,(R); the n subscript is omitted because we use a channel-independent set of evenly distributed R , values and uniform widths. The R: grid is a given by RZ = RB
+ ( m - 1)A
(47)
in the X, basis, To the extent permitted by the primitive basis, the functions A, decay approximately exponentially in the classically forbidden region, and for open channels they are oscillating waves of wavelength 2a/_~,, in the near-asymptotic region. We choose as our new basis the (Afin) eigenvectors with K~~ closest to k:. That is, we use the functions which are closest in energy to the function we are expanding. When the contracted basis functions are defined this way, we call them energy-adapted basis functions.22 If we neglect the coupling terms in (51) and look for solutions with homogeneous boundary conditions, such solutions satisfy an eigenvalue equation with the same operator as the one diagonalized to obtain the energy-adapted functions. Therefore the energy-adapted functions are expected to be very efficient for expanding F,, and hence for expanding the wave function, as required in the Schwinger method. In the Newton method, however, we expand the amplitude density, and, even when coupling and nonhomogeneous boundary conditions are neglected, this does not satisfy a corresponding eigenvalue equation. Hence (19) See,e.g.: Arthurs, A. M.; Dalgarno, A. Proc. R. SOC.London A 1960, 256, 540. Lester, W. A., Adu. Quantum Chem. 1975, 9, 199. (20) Hamilton, I. P.; Light, J. C. J . Chem. Phys. 1986, 84, 306. (21) Abdallah, J., Jr.; Truhlar, D.G. J . Chem. Phys. 1974, 61, 30. (22) Staszewska, G.; Truhlar, D. G. J . Chem. Phys. 1986, 86, 1646.
3206 The Journal of Physical Chemistry, Vol. 92, No. 11. 1988 we would use a different procedure to obtain efficient basis functions. This is not attempted here but is discussed in section 6. The internal basis functions &(x) are the same as used in ref 12. 4.1.2. Distortion Potential. We next consider the distortion potential used in eq 2. The motivation for using a distortion potential to generate the distorted solutions ryl used in the integral equations is that the equations involve the product 'UW so that H , is sensitive only to regions of space where both and JIIn are nonzero. By introducing a distortion potential we can minimize the size of this region and hence the number of basis functions required. In addition, if 'U were small enough, the first term would dominate the second in eq 9 for many elements of K. Again this should mean that only a few translational radial basis functions are required. In this paper we consider two types of distortion potentials, single channel and multichannel. The single-channel distortion potential is defined by
ps= C I+n n
) Vnn( +nI
(55)
In the multichannel distortion potential, we group together the most strongly coupled channels, namely the channels which differ only in their rotational and orbital quantum numbers. Thus we write
pM = C&njIAt)Vnd(+n,I n,n'
(56)
where And is zero or one depending on whether or not channels n and n' have different vibrational quantum numbers. For use in later equ_ationsit is convenient to define the symbol And which is equal to Ad for the multichannel distortion potential and equal to S,, for the single-channel distortion potential. 4.1.3. Distorted Waves and Green's Functions. We expand the distorted waves and Green's functions as
The functions ('?fnd and (')fnn. solve
(60) where Ad,, is defined below eq 56. The boundary conditions on (60) are ('lfnn'
i6ndk,,-'/2 sin (k,R
-
R--0
0
- lna/2)
Schwenke et al. We determine the radial functions (rlfnd and ("fnn, [as well as see eq 69 and 701 by the 9-point the auxiliary functions and finite difference boundary value method (FDBVM). The FDBVM has been described previously,'2bs23but in some of the calculations presented here we employed a new scheme for the finite difference grid. To explain this we must first recall a basic strategy used in both the old and new schemes; namely, when the finite difference approximation employed in the central part of the grid involves more than three points, it is systematically reduced in order near the ends of the grid so that the approximation to the second derivative never involves a point more than one step off an end of the grid. (The points one step off each end are the boundary values that determine which particular solution is approximated.) P r e v i ~ u s l y lwe ~ * used ~ ~ a grid which had a fixed stepsize, except that when the finite difference approximation to the second derivative involved more than three points in the central part of the grid, the stepsize was decreased by several orders of magnitude at the ends of the grid to minimize the error produced by the lower order methods required at the ends of the grid. We have now introduced two improvements. First we make the central part of the grid only approximately, but not exactly, evenly spaced. In particular, the grid points have been displaced so that a subset of the grid points are located on Gauss-Legendre quadrature nodes. Then repeated Gauss-Legendre quadrature can be and is used to evaluate the integrals in the R coordinate [which are given below, e.g., in eq 71-73, 118, 125, and 1261. In practice we divide the finite difference grid into about eight sections, and each section is spanned by the Gauss-Legendre nodes for 12-point Gaussian integration. Interdispersed are 24 other points to make the stepsize approximately constant. To check the convergence of the finite difference solution, we added extra points uniformly between these 36 points. The choice of 12-point quadrature was the result of rough optimization to determine the fewest points required to give sufficient accuracy in the integrations. The second change we make in the FDBVM concerns the strategy for reducing the stepsize at the ends of the grid. We found that making this change suddenly may cause numerical instabilities (a well-known problem in the field of numerical solutions of differential equations). Thus in the present calculations we introduced the change more gradually by decreasing the stepsize near the end of the grid at a slow fixed rate. A method that worked well is to decrease the last 20 stepsizes by a factor of 0.7 each. This leads to a stepsize at the end which is about times the central stepsize. This particular choice works well but has not been optimized. 4.1.4. Matrix Elements, Once the basis functions and distortion potentials are specified the next step in the calculation is the determination of the quantities required to evaluate the integrals in eq 15, 16, 23, 26, and 27. These integrals involve the distorted the Green's functions GD,the coupling potentials 'U, waves PO, and the basis functions. The x integrations are carried out by the methods discussed previously,'2b and the rest of the present section is devoted to the radial integrals. The radial quadratures required for the Schwinger method are #@,
(61a)
+ where
The radial integrals for the generalized Newton method are where Rf is a constant chosen to avoid numerical difficulties. It is only for the particular choice of boundary conditions given in eq 62 that eq 9 is valid without extra factors.
(23) Truhlar, D.G.; Kuppermann, A. J. Am. Chem. SOC.1971, 93, 1840.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3207
Variational Methods for Molecular Collisions
The remaining integrals in eq 63, 66, and and 71-73 are evaluated by using the extended trapezoidal rule as describedIzb previously or (repeated) Gauss-Legendre quadrature. 4.1.5. Linear Equations. For computational efficiency the final step of the calculation is written as a linear equation solution followed by a matrix multiply; e.g., eq 254s replaced by
3 - = 3- 8
+ BTA
(76)
where
CA = B It is convenient to define half-integrated Green's functions as
(77)
Then the scattering matrix S is calculated by S = ( 1 - iK)-'(l
+ iK)
(78)
Note that C is M X M whereas ( 1 - iK) is N X N . Since M >> N in the present calculations, the computational effort in (78) is negligible. 4.2. Reactive Scattering. 4.2.1. Basis Functions. For the reactive case we take
= Rag-ltmgng(Rorg)~,gng(Xag)
Once the half-integrated Green's functions are known, the calculation of c involves integrals of the same form as required to calculate b, and the calculation of B and C involve integrals of the same form as required to calculate We now consider methods to calculate these functions. For simplicity we only One way to calculate this is to determine (r)Jnd and consider ()'J,, numerically, for example, by the finite difference boundary value method discussed previously, to form g n , by eq 59, to substitute this into eq 70, and to evaluate the integral by quadrature. To facilitate this, as discussed previouslyIzband in section 4.1.3, we directly generate the functions and ('?fnd on a grid so that they need not be interpolated for the integrations over the R coordinate. Nevertheless, we found that such quadratures are difficult to converge, in part because of the discontinuous derivative of the integrand when R = R'. A second approach is to calculate directly. This can be accomplished by noting that it solves the (coupled) differential equation(s)
(79)
where 0 labels a basis function, cyB is the channel in whose coordinates this function is expressed, R, denotes the mass-scaled radial separation between A and the center of mass of BC, R, and R3 denote the corresponding atom-diatom distances for cy = 2 or 3, tmBng is a translational radial basis function as described in section 4.1 . l , x, denotes the collections of other coordinates (Le., all except R,), and )& x,J is again a standard vibrational-rotational-orbital function. For the reactive scattering examples presented in this paper, the translational radial basis functions are taken either as primitive Gaussians as described in section 4.1.1 or as sine functionsz4defined by
zB.
R,
(0,
< Ro or R, > R
,
EB,.
otherwise (80) where AP = RN- Ro and Ro and R, are parameters. In practice Ro is close to the smallest classical radial turning point, and RN is chosen such that H essentially vanishes for R > RN In the symmetric system considered here, Ro and R N are the same for all three arrangements. Physically the sine-type translational basis function of (80) is the eigenstate of a square well. Sine-type basis sets are easy to use since they have only two parameters (& and RN), but they make the quadratures more difficult than Gaussians because they are d e l o c a l i ~ e d . ~ ~ 4.2.2. Distortion Potential and Green's Function. As in the nonreactive case we consider distortion potentials of the form
e
Here A&, is unity if channel n and channel n'belong to the same distortion potential block and zero otherwise, and
subject to the boundary conditions
2P
u&t= - -Jdxa hZ k,"*
COS
(k,,R - 1,1r/2)
(2~kn~)-11* exp[-lknl(R
- Rf)]
>0 k,2 < 0 k,2
(75b)
It is an easy matter to solve eq 74 subject to the boundary conditions of eq 75 using the finite difference boundary value method. We call this the method of half-integrated Green's functions (HIGF). We found we could obtain especially high accuracy with the HIGF method. In addition, we note that the HIGF method avoids the calculation of the irregular functions ("fnn. This second point is significant since these are more difficult to work with numerically.
$an*
(x,) Vntsn(Ra,x,)4 a n 4 X a )
(82)
In the case of a single-channel distortion potential, which we denote we replace Ain, by a Kronecker delta; and for what we call a rotationally coupled distortion potential, we set A;,, equal to unity only for channel pairs which differ at most by rotational-orbital quantum numbers (for other pairs it is zero). In the reactive case one can also consider setting A&, equal to unity for all n and n' so that all the channels of a given arrangement are coupled in distortion potential. This is called full intraarrangement distortion.
e',
(24) Zvijac, D. J.; Light, J. C. Chem. Phys. 1976, 12, 237. Shima, Y . ; Baer, M. Chem. Phys. Lett. 1982, 91, 43. Shima, Y . ;Baer, M. J . Phys. B 1983, 16, 2169. (25) Schwenke, D. W.; Haug, K.; Truhlar, D. G.; Schweitzer, R.; Zhang, J. Z. H.; Sun,Y . ;Kouri, D. J. Theor. Chim. Acta 1987, 72, 237.
3208 The Journal of Physical Chemistry, Vol. 92, No. 1I , 1988 It was considered in a preliminary communication26on this research, and we will not repeat most of those results here. The radial functions determined for a given choice of distortion potential will be defined and denoted as in section 4.1.3 except that all quantities now have an additional label a, e+, the distorted wave is
= (R,,x,IV‘) = R;’?C(R,)
v”(R,,x,)
dan(X,) (83)
for a single-channel distortion potential, and eq 57 is 1CiaoYR,,x,) = R , ; l ~ A ~ ~ o & o , , ( ~ ,(r!C:o(R,o) o)
Schwenke et al. plicated by nonorthogonality of the states. For the generalized Newton variational method, we have three matrices to calculate namely (42), (44), and (45). Various approaches and combinations of approaches are possible for evaluating these integrals. We will consider two approaches in the following subsections. 4.2.3.1. Single-Channel Distortion Potential and Numerical Integration of Explicit Green’s Functions with an Insertion Basis. First consider Bod. According to (44), (85), (86), and (91), we may write this integral as
= ( anmlG:Ua’(a’n’k,,)
B,,
(84)
n
for a multichannel distortion potential. We will require additional details of the case of a single-channel distortion potential. In this case it is convenient to rewrite (83) in ket form
Iv”) = lank,)
= I4an)lakn)
(85)
where I&,,) is a vibrational-rotational-orbital eigenstate, and lak,) is an associated translational state. In the same way, the basis function @p can be expressed as
1%)
= I“,ngm,)
=
(86)
14,8np)l~pq9)
where lam) is a translational basis function in arrangement a. When H: acts on Id,,,), its vibrational-rotational-orbital operators are simply replaced by their eigenvalue, leaving the radial translational derivative unchanged, i.e., H:I&fl) = h:flld,n)
= ( amlE(d,,121”la’n’kn,)
Inserting unit operators yields Bp”, =
The solution of (89) is expressed in the coordinate representation as (R,lakfl) = (1 /Ra)(r)mR,)
(90)
Similarly,
(a’dIa’k,,O (98)
where law’) is a translational insertion basis set which is not necessarily identical with lam). In the present article though we assume these bases are the same, Le., (R,lV) = ( 1 / R a ) w L )
(99)
The summation over p‘ can be truncated as long as the results are converged to the desired accuracy. Similarly, using (45), (86), and (91) we may write C,, as
C,, = (anmlG:(l - Uu’G:,)Ja’n’m’) = (amlE(d,,la’n’m’) - ( a m l b( 4,,,lUa’E:la’n’m’)
R,2(amlElR,) (RaI(4,&’n’m’)
CS d R ,
where the internal eigenenergy now has an arrangement subscript added. The translational states satisfy the equation (89)
lap’) (ap’l
e’
= JdR,
(E - h:,,)lak,,) = 0
and
c S d R , R,Z(amlEIR,) (R,I (4,,Iw”la’n’d)
(87)
where h:, involves only one derivative operator and is given by
1; dR,R21R,) (R,1
(97)
( I00)
-
Ra2(amlElR,) (R,I ( d,,,I~“Ia’n ’1’) ( a’~’lE:Ia’m’)
r’
=
1% Ra2(amlEIR,)
[ (R,I (danla’n’m’) -
C (R,I (d,,,121Q’la’n’w’)( a ’ ~ ’ l ~ : l a ’ m ’(101 )] e’
The third matrix is
RE,,, = (ank,JUm’Ia’n’kd)= e‘
S d R , R,Z(ak,,IR,) (R,I(d,,I”U*’la’n’d) (a’dIa’k,,O (102)
Using
El@,,,)(91) Clearly
(aknlR,) = (l/R,)(r%(R,)
(103)
gives ( - 2 p / h 2 ) ( E - h:,,)%
= 1
(92)
In the coordinate representation, (92) reads
(93) (94)
we get (-2p/h2)(E - hE,,)R,-’&(R,,R,’)
= X S d R , R , (r%(R,)(R,l(d,nl~~’la ’n ’p ’) (a’yla’kd) e‘
(104)
Defining
&(KA‘)= ~3,’(~,lZlK9
%En,
Expressions 98, 101, and 104 have a common matrix element ( R , ~ ( & , J W d ~ a f nwhich ~ ) , is independent of energy. It is defined in ref 12 as ep’(R,,) and is discussed in detail in that reference. The above integrals also involve the matrix elements (a’p’la’k,,); (amlElR,), (ap’lg$la’m’), and (R,I(&,la’n‘m’). These can be put in computational form by using eq 94, 99, and 103 as follows: (a’p’la’k”,) =
= R,-’6(R, - RY)
(95)
Solving (95), we obtainI2 G(RaJL’) = (r)fi(R3(i%XR3 (96) where R,< = min (R,,RY), RZ = max (R,,RJ), and (i)fi(R,) is an irregular solution of (89). 4.2.3. Integrals for Reactive Scattering. For reactive scattering, the matrix elements involving different arrangements are com(26) Schwenke, D. W.; Haug, K.; Truhlar, D. G.; Sun, Y.; Zhang, J. 2. H.; Kouri, D. J. J . Phys. Chem. 1987, 91, 6080.
(amlEIR,) =JdR,‘
SdR,, Xe,(R,,) (r)flt’(R,,)
R a ” ( W K “ (RdIEIR,) = (1 /Rm)JdRa‘
Am(Ra‘)&(KA’) (106)
(a’p’1E:la’m’) = S d R , . RarXe,(R,)(a’mlE:IR,.)
(R,l#&’n’m’)
(105)
= Jdx,
(107)
1
d,n*(x,) dadXal) --m,(Ra,) Rd (108)
It is understood that x,, and Ref are functions of x, and R,.
The Journal of Physical Chemistry, Vol. 92, No. 1 I , 1988 3209
Variational Methods for Molecular Collisions An alternative method for calculating the integrals is to insert the identity C,Iap)(apl into the energy-independent matrix ep'(R,), i.e.
ep'(R,) = (R,I (&,lU"'Ia'nk') = 1 x--X,(R,)(anp(UdJa'n'p') Ir Ra
(109)
r, is the mass-scaled'2b bond length of the diatom in arrangement a, (...I...) is a Clebsch-Gordan coefficient, qmis a spherical harmonic, 7,is the angle defined by y, = arccos (i,.A,) (114) a rotation matrix element with Euler angles +,,B,,x,, dim, is a reduced rotation matrix element, and A,," = arccos (R,.R,,) (115)
@Mm
1
The tradeoff between the two methods of calculating the energy-independent matrix is discussed in ref 12b and 25. Note: Since the "insertion basis" l p ) is used to resolve the identity, we often need to use a large number of functions tb achieve convergence. The total number p,,, of functions I p ) will not affect the dimensions of the B, C, and XBmatrices, which are determined by the number of variation translational basis functions Im). The Appendix contains a proof that when the size of the insertion basis is the same as the size of the variation basis, Newton's variational principle reduces to the method of moments for the amplitude density. The matrix element (amlg(R,) can also be calculated directly by solving the inhomogeneous differential equations introduced in section 4.1.3 (HIGF method), although this was not done for the reactive scattering calculations in this article. 4.2.3.2. Multichannel Distortion Potential and Half-Integrated Green's Function without Insertion. The quantities we require OK&,, and Ea. We for the final approach are &,, (r)&, reiterate from section 4.2.2 that these quantities are identical with those required in the nonreactive formalism except that now they have an extra index. As with the nonreactive case, the introduction of the half-integrated Green's functions Epreduces the integrals for the matrices %:%, B,,, and CsFto a common form. This simplication remains possible for the reactive case because in eq 44-45 the Green's function always appears with a basis function of the same arrangement channel. Furthermore, the same half-integrated Green's functions are involved in the reactive and nonreactive case. First consider the calculation of %fl,,. The integration methods we use are very similar to those employed recently in a paper discussing the Born approximation for reactive ~ c a t t e r i n g . ~ ~ Introducing the expansion of (84) into (42) produces
c,,,
%n"n,
=
n hd
A:,,.A:gdSd7
R;'
c#I,,,*(x,)
X
''X,n(Roi)UoloR,,-'~olono'(X,,)!rX3n,(Ruo) ( 110)
where d7 is a 6-dimensional volume element. We need consider only the case of a # ao. [Otherwise (66) should be used.] Of the six dimensions, three can be integrated analytically by exin terms pressing the vibrational-rotational-orbital functions of body-frame coordinates. Thus we writeIzb 1 + a n ( X a ) = -xaudn(rJ X ra
-
See ref 12b for a detailed discussion of the sign conventions we employ. Since ?la,is independent of the Euler angles +,,B,,X,, the integrals over these quantities can be performed analytically. The remaining three degrees of freedom are most conveniently integrated by using the coordinates R,, R,,, and A,,,. Thus we have
~:;P(R,,R,,)
+
[(21, 1)(21, + 1)]1/2 =2 p ( ~ 9 3 ~ a ~ , , X 25 + 1
1
-~m,p;j;(ru,)Y , ; m ~ ( y , , , O ) d ~ ~ ~ ( A (1 , , O18) )
rue
where 5 3 3 is the same as 6'3 except -2p~"'."/A2 is replaced by unity, and J P o is a mass factor defined below. We emphasize that eq 118 applies only for a # ao. In terms of the independent variables R,, R,,, and A,,,, the other coordinates are given by
+ (kaaoRJ2 + 2.&uuoR,R,0
ra = Jnaao[R,:
cos A,ao]112 (119)
cos y, = (-l)p-~dPO(kOaoR, + R,, cos A,,,)/r, rao = Jnuao[R2
+ (J%~'QR,,)~+ 2k""oR,R,,
cos A,ao]112 (121)
cos y,, = ( - l ) p ~ ~ o ~ a a ~ ( k m+ a ~R,R ,cos o A,,o)/r,, J p o
= [ ( M - m,)(M &uao
- m,,)/(M
- m, - m & V l ' ' 2
= [m,m,,/(M - m,)(M -
(120)
(122) (123) (124)
where Pa* is the parity of the permutation from (12) to (aao), M is the total mass of the three atoms, and m, is the mass of the lone atom in arrangement a. With these definitions, it is easy to give the remaining required expressions: Le., (44) becomes
and (44) becomes where
xauj is a vibrational function solving
(27) Schatz, G. C.; Hubbard, L. M.; Dardi, P. S.; Miller, W. H. J . Chem. Phys. 1984, 81, 231.
for a # cyo, (72) and (73) for a = ao. The cos A,", integral in eq 118 is performed by using GaussLegendre quadrature, and the R , and R,, integrals in eq 72-73,
3210 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Schwenke et al.
TABLE I: Transition Probabilities for Nonreactive Scattering of I by HIat E = 0.045Eh MMAD~ RMP" 6' 0.26 5.51e 6.5Ie 0.8
'nmax
A
R? RLx
n
n'
I/
1' 2 3 4 2 3 4 3 4 4
1
1 1 2 2 2 3 3 4
6.25 (-I)$ 3.75 (-1) 3.58 (-6) 4.47 (-7) 6.25 (-1) 1.06 (-5) 1.87 (-6) 9.49 (-1) 5.14 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 2.18 (-4) 5.01 (-5) 6.25 (-1) 6.50 (-4) 1.66 (-4) 9.49 (-1) 5.11 (-2) 9.49 (-1)
16 0.2 5.31 8.31 1.8 6.26 (-1) 3.74 (-1) 7.05 (-6) 2.03 (-7) 6.26 (-1) 1.73 (-5) 3.77 (-7) 9.49 (-1) 5.11 (-2) 9.49 (-1)
31 0.12 5.31 8.91 1.6 6.25 3.75 4.04 1.06 6.25 2.42 4.68 9.49 5.11 9.49
(-1) (-1) (-6) (-6) (-1)
(-5) (-6) (-1) (-2) (-1)
31 0.2 5.31 11.31 1.8
46 0.2 5.31 14.31 1 .8
61 0.2 5.31 17.31 1.8
101 0.12 5.31 17.35 1.6
6.25 (-1) 3.75 (-1) 3.72 (-6) 6.71 (-7) 6.25 (-1) 1.13 (-5) 2.82 (-6) 9.49 (-1) 5.09 (-2) 9.45 (-1)
6.26 (-1) 3.74 (-1) 3.96 (-6) 5.82 (-7) 6.26 (-1) 1.12 (-5) 2.43 (-6) 9.49 (-1) 5.11 (-2) 9.49 (-1)
6.25 (-1) 3.74 (-1) 3.79 (-6) 4.99 (-7) 6.25 (-1) 1.17 (-5) 2.18 (-6) 9.49 (-1) 5.11 ( - 2 ) 9.49 (-1)
6.25 (-1) 3.75 (-1) 3.61 (-6) 4.60 (-7) 6.25 (-1) 1.08 (-5) 1.94 (-6) 9.49 (-1) 5.11 (-2) 9.49 (-1)
R matrix propagation. Method of moments for amplitude density with multichannel distortion potential. CNumberof distributed Gaussian basis functions. dSpacing between Gaussian centers (in mass-scaled bohr). e Position of first and last Gaussian (in mass-scaled bohr). /The channel quantum numbers are 1 (D = 0, j = 0); 2 (0, 2); 3 (1, 0); 4 ( I , 2). $6.25 (-1) 6.25 X IO-'.
TABLE 11: Transition Probabilities for Nonreactive Scattering of I by H2at E = 0.045Eb SVPb RMP" 20' 0.26 4.4Ie 8.21' 1.4
"ax
A R? RLx n 1/ 1 1 1 2 2 2 3 3 4
20 0.2 5.51 9.31 0.8
20 0.1 4.01 5.91 0.8
30 0.3 4.71 13.41 1.4
6.25 (-1) 3.75 (-1) 2.70 (-6) 1.67 (-6) 6.25 (-1) 1.23 (-5) 8.72 (-6) 9.49 (-1) 5.11 (-2) 9.49 (-1)
6.28 (-1) 3.70 (-1) 1.19 (-3) 6.12 (-4) 6.24 (-1) 1.99 (-3) 3.63 (-3) 9.50 (-1) 4.67 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 8.48 (-6) 2.81 (-7) 6.25 (-1) 2.47 (-5) 4.56 (-7) 9.49 (-1) 5.11 (-2) 9.49 (-1)
35 0.3 4.71 16.91 1.4
40 0.1 4.01 7.91 0.8
60 0.1 4.01 9.91 0.8
n' 11 2 3 4 2 3 4 3 4 4
6.25 (-l)g 3.75 (-1) 3.58 (-6) 4.47 (-7) 6.25 (-1) 1.06 (-5) 1.87 (-6) 9.49 (-I) 5.14 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 6.23 (-6) 8.59 (-7) 6.25 (-1) 1.51 (-5) 4.36 (-8) 9.49 (-1) 5.11 (-2) 9.49 (-1)
6.25 3.75 9.17 4.94 6.25 2.25 5.61 9.49 5.11 9.49
6.25 (-1) 3.75 (-1) 2.59 (-6) 7.31 (-7) 6.25 (-1) 1.90 (-5) 2.48 (-7) 9.49 (-1) 5.11 (-2) 9.49 (-1)
(-1) (-1) (-6) (-7) (-1) (-5) (-7) (-1) (-2) (-1)
6.25 (-1) 3.75 (-1) 3.96 (-6) 7.01 (-8) 6.25 (-1) 1.67 (-5) 3.44 (-7) 9.49 (-1) 5.11 (-2) 9.49 (-1)
R matrix propagation. Schwinger variational method with multichannel distortion potential. 'Number of distributed Gaussian basis functions. dSpacingbetween Gaussian centers (in mass-scaled bohr). eCenter of first and last Gaussian (in mass-scaled bohr). 'Channel quantum numbers are 1 (zi = 0, j = 0), 2 (0, 2), 3 ( I , 0), 4 ( I , 2). 86.25 (-1) = 6.25 X lo-'.
116, and 125-126 are performed by using repeated Gauss-Legendre quadrature. 5. Calculations
We will illustrate the above methods with calculations on four test problems. For the first two, in which we treat nonreactive I H2collisions on the modified Parr potential energy surface12b,28 and reactiveH H2 collisions on the Porter-Karplus no. 2 (PK2) surface,29we restrict ourselves to a fixed-length close-coupling expansion, and we study only the rate of convergence with respect to increasing the size of the translational basis. For the third test problem we demonstrate complete convergence of H + H2 transition probabilities for the accurate Liu-Siegbahn-Truhlar-Horowitz (LSTH) surface30 which has also been studied in previous work.i2s3i-34 For the fourth test problem, we return to reactive
+
+
(28) Parr, C. A. Ph.D. Thesis, California Institute of Technology, Pasadena, CA, 1969. (29) Porter, R. N.; Karplus, M. J . Chem. Phys. 1964, 40, 1105. (30) Liu, E.J . Chem. Phys. 1973,58,1925. Siegbahn, P.; Liu, B. J . Chem. Phys. 1978.68.2457. Truhlar, D. G.; Horowitz, C. J. J . Chem. Phys. 1978, 68,2466. 1979, 71,1514(E). (31) Walker, R. B.: Stechel, E. B.; Light, J. C. J . Chem. Phys. 1978, 69, 2922. (32) (a) Colton, M. C.; Schatz, G . C. Chem. Phys. Lett. 1986, 124, 256. (b) Schatz, G. C., private communication, July 1987. . (33) Webster, F.; Light, J. C. J . Chem. Phys. 1986, 85, 4744. (34) Linderberg, J.; Vessal, B. Znr. J . Quantum Chem. 1987, 31, 65.
H + H2 on the PK2 surface, and we test the rate of convergence with respect to the internal basis as well as the translational basis. 5.1. Nonreactive Examples. For a nonreactive test case we consider the problem
I
+ H~(uJ')
----*
I
+ H~(u'J'')
on the potential energy surface of Parr,2sslightly modified.i2b All calculations reported here for this system are for total angular momentum zero and for four channels at a total energy 0.045Eh (28.238 kcal/mol), with respect to the I + H2potential asymptote; this energy is between the thresholds for u' = 1 and u' = 2 . The channels included are u = 0, j = 0, 2 and u = 1, j = 0, 2. These channels define a model space for testing convergence of the radial translational basis; they are not intended to represent a converged state space. The basis functions are products of asymptotic vibrational eigenfunctions, Arthurs-Dalgarno rotational-orbital functions, and distributed Gaussians in the translational coordinate. The results are given in Tables I-V. Each table gives the results of several L2calculations a_nd compares them to the results obtained by a numerically accurate propagation calculation employing R matrix propagation (RMP) for the same four-channel problem. The RMP results are used for comparison to test whether the L2 calculations are converged with respect to the translational radial basis. Therefore in Tables I-V we report only those numerical parameters that are associated with the radial translational basis. We also report only a selection of the tests
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3211
Variational Methods for Molecular Collisions
TABLE III: Transition Probabilities for Nonreactive Scattering of I by Hzat 0.045Eh SVP'
n
n'
RMP"
4c
1 1 1 1 2 2 2 3 3 4
1 2 3 4 2 3 4 3 4 4
6.25 (-1) 3.75 (-1) 3.58 (-6) 4.47 (-7) 6.25 (-1) 1.06 (-5) 1.87 (-6) 9.49 (-1) 5.14 (-2) 9.49 (-1)
5.93 (-1) 4.07 (-1) 1.08 ( - 5 ) 1.07 (-5) 5.93 (-1) 5.32 ( - 5 ) 6.72 (-6) 9.54 (-1) 4.64 (-2) 9.54 (-1)
" R matrix propagation.
10 6.00 4.00 9.46 5.77 6.00 2.58 1.77 9.51 4.86 9.51
(-1) (-1) (-6) (-5) (-1) (-5) (-4) (-1) (-2) (-1)
20
60
2+2
2+4
4+4
6.24 (-1) 3.76 (-1) 4.04 (-6) 4.33 (-7) 6.23 (-1) 1.63 ( - 5 ) 7.17 (-7) 9.48 (-1) 5.11 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 3.96 (-6) 7.01 (-8) 6.25 (-1) 1.67 (-5) 3.44 (-7) 9.49 (-1) 5.11 (-2) 9.49 (-1)
6.01 (-1) 3.99 (-1) 1.43 ( - 5 ) 4.17 ( - 5 ) 6.00 (-1) 5.07 (-5) 2.17 (-4) 9.50 (-1) 4.96 (-2) 9.50 (-1)
6.32 (-1) 3.67 (-1) 7.64 (-5) 1.14 (-4) 6.32 (-1) 3.88 ( - 5 ) 2.86 (-4) 9.55 (-1) 4.49 (-2) 9.55 (-1)
5.74 (-1) 4.26 (-1) 2.58 (-6) 6.17 (-6) 5.74 (-1) 1.79 (-5) 2.81 ( - 5 ) 9.49 (-1) 5.06 (-2) 9.49 (-1)
bSchwinger variational method with multichannel distortion potential with Gaussian spacing A = 0.1 and starting point
RY = 4.01 in mass-scaled bohr. The Gaussian overlap parameter is c = 0.8. 'A single number means number of energy adapted basis functions. n
+ m means the first n energy adapted basis functions plus the m primitive Gaussians which lie closest to the translational turning point for that channel. TABLE IV: Transition Probabilities for Nonreactive Scattering of I by H2at E = 0.045Eh: Comparison of Variational Calculations with Multichannel Distortion Potential to Calculations Employing Conventional Close Coupling GNVP'
RMP" 3c 0.3d 5.61' 6.41'
4 0.3 5.61 6.51
6.21 (-1) 3.79 (-1) 3.43 (-6) 3.73 (-7) 6.21 (-1) 9.88 (-6) 1.42 (-6) 9.49 (-1) 5.09 (-2) 9.49 (-1)
6.26 (-1) 3.74 (-1) 3.51 (-6) 4.73 (-7) 6.25 (-1) 1.06 (-5) 1.97 (-6) 9.48 (-1) 5.20 (-2) 9.48 (-1)
"ar
A RP
REm,, n' 1f 1f 1 2 1 3 1 4 2 2 2 3 2 4 3 3 3 4 4 4
5 0.3 5.61 6.81
7 0.3 5.61 7.41
8 0.3 4.71 6.81
15 0.3 4.7 1 8.91
20 0.2 4.71 8.51
n
6.25 3.75 3.58 4.47 6.25 1.06 1.87 9.49 5.14 9.49
(-1)f
(-1) (-6) (-7) (-1) (-5) (-6) (-1) (-2) (-1)
6.25 (-1) 3.75 (-1) 3.61 (-6) 4.57 (-7) 6.25 (-1) 1.07 (-5) 1.93 (-6) 9.49 (-1) 5.11 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 3.60 (-6) 4.50 (-7) 6.25 (-1) 1.06 (-5) 1.89 (-6) 9.49 (-1) 5.12 (-2) 9.49 (-1)
6.25 3.75 3.60 4.56 6.25 1.06 1.92 9.49 5.11 9.49
(-1) (-1) (-6) (-7) (-1) (-5) (-6) (-1) (-2) (-1)
6.25 (-1) 3.75 (-1) 3.60 (-6) 4.56 (-7) 6.25 (-1) 1.06 (-5) 1.89 (-6) 9.49 (-1) 5.12 (-2) 9.49 (-1)
6.25 (-1) 3.75 (-1) 3.60 (-6) 4.50 (-7) 6.25 (-1) 1.06 (-5) 1.89 (-6) 9.49 (-1) 5.12 (-2) 9.49 (-1)
'
" R matrix propagation. Generalized Newton variational method with multichannel distortion potential. Number of distributed Gaussian basis functions (c = 1.4). dSpacing between Gaussian centers (in mass-scaled bohr). 'Position of first and last Gaussian (in mass-scaled bohr). /The channel quantum numbers are 1 ( u = 0, j = 0); 2 (0, 2); 3 (1, 0); 4 (1, 2). 86.25 (-1) = 6.25 X lo-'. TABLE V Transition Probabilities for Nonreactive Scattering of I by Hz at E = 0.045Eh: Comparison of Variational Calculations with Single-Channel Distortion Potential to Conventional Close Coupling RMP" 15c 0.4d 4.71e 1O.3lc
Mmax
A RP
R:,, n
n'
18 1 1 1 2 2 2 3 3 4
18 2 3 4 2 3 4 3 4 4
GNVPb
6.25 (-l)h 3.75 (-1) 3.58 (-6) 4.47 (-7) 6.25 (-1) 1.06 (-5) 1.87 (-6) 9.49 (-1) 5.14 (-2) 9.49 (-1)
6.10 (-1) 3.90 (-1) 2.14 (-6) 5.88 (-7) 6.10 (-1) 1.55 (-5) 1.40 (-6) 9.96 (-1) 3.97 (-3) 9.96 (-1)
15, 1 6 0.3 4.71, 4.41 8.91 6.14 (-1) 3.86 (-1) 3.59 (-6) 4.34 (-4) 6.14 (-1) 1.03 (-5) 1.85 (-6) 9.49 (-1) 5.14 (-2) 9.49 (-1)
5.83 4.17 2.93 1.69 5.83 7.90 2.97 9.96 4.05 9.96
15 0.2 4.71 7.51 (-1) (-1) (-7) (-9) (-1) (-7) (-9) (-1) (-3) (-1)
20 0.2 5.01 8.81 6.26 (-1) 3.74 (-1) 3.44 (-6) 4.26 (-7) 6.26 (-1) 1.08 (-5) 1.91 (-6) 9.49 (-1) 5.14 (-2) 9.49 (-1)
25 0.3 5.01 9.81
27 0.3 4.71 12.51
6.25 (-1) 3.75 (-1) 3.61 (-6) 4.63 (-7) 6.25 (-1) 1.05 (-5) 1.95 (-6) 9.49 (-1) 5.13 (-2) 9.49 (-1)
6.86 (-1) 3.14 (-1) 4.17 (-6) 5.40 (-7) 6.86 (-1) 7.47 (-6) 2.37 (-6) 9.49 (-1) 5.12 (-2) 9.49 (-1)
27 0.2 4.71 9.91 6.25 (-1) 3.75 (-1) 3.59 (-6) 4.80 (-7) 6.25 (-1) 1.09 ( - 5 ) 2.04 (-6) 9.49 (-1) 5.12 (-2) 9.49 (-1)
27 0.2 5.01 10.21 6.25 3.75 3.59 4.44 6.25 1.07 1.87 9.49 5.12 9.49
(-1) (-1) (-6) (-7) (-1) (-5) (-6) (-1) (-2) (-1)
30 0.2 4.71 10.51 6.25 3.75 3.59 4.49 6.25 1.07 1.90 9.49 5.12 9.49
(-1) (-1) (-6) (-7) (-1) (-5) (-6) (-1) (-2) (-1)
'R matrix propagation. 'Generalized Newton variational method with single-channel distortion potential. 'Number of distributed Gaussian basis functions (c = 1.4). dSpacing between Gaussian centers (in mass-scaled bohr). ePosition of first and last Gaussian (in mass-scaled bohr). [These two bases give the same results to three significant figures. gThe channel quantum numbers are 1 (v = 0, j = 0); 2 (0, 2); 3 (1, 0); 4 (1, 2). h6.25 (-1) = 6.25 X lo-'. of the various numerical parameters. For example, in tests not reported here we find that the best convergence is achieved when the Gaussian overlap parameter c (introduced in section 4.1.1) is assigned a value in the range 0.8-1.8, and the results for the present test case are not very sensitive to varying it in this range. We performed the radial integrals by Gauss-Hermite-like quadratures with displaced nodes as discussed in ref 12b, and the results in the tables are based on 9- to 19-point quadratures. Table I gives results obtained by the nonvariational method of moments
for the amplitude density with multichannel distortion potentials. Tables I1 and 111give results obtained by the Schwinger variational principle again with multichannel distortion potentials. Table I1 is for primitive basis functions, and Table I11 is for energy-adapted contracted ones. Tables IV and V give results for the Newton variational principle, Table IV for multichannel distortion potentials and Table V for single-channel ones. For this simple four-channel problem all distortion potentials are rotationally coupled; i.e., they are based on two coupled ro-
3212 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 TABLE VI: Parameters for Reactive Scattering Calculations with Sine Bases for Translation
parameter
footnote a
b c
d e e
f g g g
h i
.i
i i
set 1 9 45-1 80 5-20 60 1.84 5.74 0.06 0.80 10.82 0.02 9 11 0.6 3.2 7
set 2
run 1
run 2
run 3
9 135 15 70 1.66 6.06 0.04 0.70 10.92 0.01 11 13 0.5 3.4 9
66 1980 30 64 1.60 5.72 0.02 0.70 9.71 0.01 9 30 0.6 3.2 I1
87 2610 30 64 1.60 5.72 0.02 0.70 9.71 0.01 9 30 0.6 3.2
66 1320 20 64 1.60 5.72 0.02 0.70 9.71 0.01 9 30 0.6 3.2
11
11
"umber of channels. bNumber of basis functions to expand all amplitude density components. Number of radial translational basis functions for each channel. dNumber of sine functions in insertion basis for each channel. CBoundariesof interval over which sine bases are defined by eq 80b. /Quadrature step size for radial integrals. 8Boundaries and step size of finite difference grids for calculation of the Green's functions. *Number of quadrature points for angular integrals. 'Number of quadrature points for vibrational integrals with a # a' and a = a'. 'Interval for trapezoidal vibrational quadratures. tational channels in the same vibrational manifold. The integrals over the Green's functions were evaluated by the HIGF method of section 4.1.4. For these calculations eq 74 or their analogues for the Schwinger method were solved by the finite difference boundary value method of ref 12b with a variably spaced grid and grid spacing as given in eq 132b of ref 12b with A(F) equal to 0.05~~. Other numerical details are similar to those of ref 12b. The inelastic scattering calculations employing the L2methods are compared to calculations employing the R matrix propagation to demonstrate convergence. The methods used for the R matrix propagation calculations are explained 5.2. Applications to Reactive Scattering. We performed three sets of calculations for reactive scattering. For the first set of reactive test calculations we used the Newton method with a single-channel distortion potential for the H + H2 system on potential energy surface no. 2 of Porter and K a r p l u ~ . ' ~ (Although this potential energy surface is not very accurate, it defines a widely employed test system.) We consider zero total angular momentum and a total energy of 0.60 eV with respect to the H H2 potential asymptote; at this energy only the ground vibrational state is open. We included three channels ( u = 0; j = 0-2) in each arrangement, for a total of nine channels. This is not a converged internal basis; rather it defines a model space for examining convergence with respect to the radial translational basis. First we performed a series of calculations in which we monitored the rate at which the results for this channel selection converge to the radial-translational basis function limit as the number of sine functions is increased. In this series of calculations we calculated the Green's functions by the finite difference ~ boundary value method with boundaries at 0.80 and 1 0 . 8 2 ~and a step size of 0.02ao,and we performed all radial integrations by the trapezoidal rule with a step size of 0 . 0 6 ~ ~The . energy-independent exchange matrix was calculated by the insertion method of eq 109 with 60 sine functions. For both the variable-size basis used to expand the reactive amplitude density and the 60-function basis used for the insertion basis, the range parameters were Ro = 1 . 8 4 and ~ ~ R , = 5 . 7 4 ~ ~The . other parameters are given in Table VI, where they are called set 1. Then we performed another
+
(35) Light, J. C.; Walker, R. B. J . Chem. Phys. 1976, 65, 4272. ( 3 6 ) Mullaney, N. A,; Truhlar, D. G. Chem. Phys. 1979, 39, 91. (37) Truhlar, D. G.; Harvey, N. M.; Onda, K.; Brandt, M. A. In Algorithms and Computer Codes for Atomic and Molecular Quantum Scattering Theory; Thomas, L., Ed.; Lawrence Berkeley Laboratory: Berkeley, CA, 1979; Vol. I , Chapter 14.
Schwenke et al. TABLE VII: Parameters for Reactive Scattering Calculations with Gaussian Bases for Translation Jmax
parameter v," runs 1, 3 run 2 run 4 run 5 run 6, 7 run 8 run 9 run 10 run 11 runs 12, 13 run 14 run 1 5 run 16
u = Ob 0 0 0 0 1 1 1 1 2 2 2 2 3
7 7 7 7 7 7 7 9 7 9 9 9 9
D
L 1'
7 7 7 9 7 7 7 7 7
Nd N(F)' R?(a,Y A(a# 24 344 2.20 0.367 24 668 2.20 0.55 24 668 2.20 0.367 24 344 1.83 0.367 48 344 2.65 0.38 48 344 2.55 0.37 48 344 2.20 0.35 60 344 2.65 0.38 72 344 2.65 0.38 78 344 2.65 0.38 78 344 2.46 0.30 78 344 2.46 0.35 102 344 2.65 0.38
"Maximum vibrational quantum number in basis. *Maximum rotational quantum number in basis for ground vibrational level. 'Same for excited vibrational levels. dNumber of channels. eNumber of points in finite difference grid. fCenter of first Gaussian. g Spacing between Gaussians. calculation with a new set of parameters, keeping the translational basis fixed at 15 functions. The new parameters are given in Table VI as set 2. Next we consider H H 2 on the LSTH potential, and again we employ the Newton method with a single-channel distortion potential and a sine basis for the radial translational degree of freedom. We again consider zero total angular momentum, and results are obtained at four energies. We present two runs with different numbers of channels. In run 1, we included two runs with different numbers of channels. In run 1, we included three vibrational manifolds in each arrangement with 10, 8, and 4 rotational states, respectively, for a total of 66 channels. In run 2, we included 4 vibrational manifolds in each arrangement, with 12, 10, 6, and 1 rotational state, for a total of 87 channels. In run 3, we have used the same basis set and quadrature parameters as in run 1 except that the number of translational basis functions is reduced from 30 to 20. The numerical parameters for these runs are given in the last two columns of Table VI. For the third set of reactive scattering calculations we employed the Newton variational method with multichannel distortion potentials and primitive (Le., uncontracted) distributed Gaussians for the radial translational functions. We considered zero total angular momentum, and 0.60 eV total energy. We also considered both restricted and converged internal-state expansions. We performed over 40 calculations to test the effect of changing various parameters. We will present results from 16 of these calculations, and several parameters of these runs given in Table VII. For all these runs, the finite difference grid for the Green's functions extended from 0.361 to 13.872a0, the angular quadratures were performed with 30-point integrations, the vibrational quadratures were performed with 10-point integration, the vibrational wave functions were obtained in a basis of 30 harmonic oscillator functions, and the neighbor-Gaussian overlap parameterz0 c was set equal to 1.4. (The latter parameter determines the width of the Gaussians.) In all cases the distortion potential was rotationally coupled. All reactive scattering probabilities in this paper include the factor of 2 for identical product rearrangements.
+
6. Discussion Tables I and I1 show the results of Lzcalculations with primitive Gaussian basis sets and multichannel distortion potentials. These tables show that the method of moments for the amplitude density and the Schwinger variational method yield similar results as the R matrix propagation method for the I + H2 test case. Table I shows that, for the method of moments with 6 Gaussians per channel, only the probabilities greater than 0.05 are well converged, while the others are sometimes 2 orders of magnitude in error. When Gaussians are added further into the asymptotic region, the small probabilities slowly converge to a value in close
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3213
Variational Methods for Molecular Collisions
-
curate results with only 4 primitive Gaussians per channel. In fact small transition probabilities obtained with 4 primitive Gaussians per channel in a variational calculation in Table IV are much better converged than those obtained by the method of moments with 61 functions per channel in Table I. Thus, for the same degree of convergence, the use of the generalized Newton variational method for the amplitude density decreases the mmax required for convergence by more than a factor of 15 as compared to the method of moments with the same distortion potential. Table V shows that repeating the calculations of Table IV with a single-channel distortion potential also yields much less accuracy and less stable results for a given basis, and about 27 Gaussians are required to stably converge all the probabilities. Thus the multichannel distortion potential decreases the mmaxrequired for convergence by about a factor of 7 . In summary we find that the Newton variational formulation converges much better than the Schwinger formulation or the method of moments, when all are applied to the same test case with a multichannel Green's function. Furthermore, when the Newton variational method with multichannel Green's function is compared to the Newton method with single-channel Green's function, the multichannel Green's function shows much better accuracy for a given basis. The reactive scattering calculations are summarized in Tables VIII-XII. Table VI11 demonstrates, by comparing variational calculations with two sets of numerical parameters, that the calculations are very stable with respect to simultaneous variations in all numerical parameters. Table IX compares the method of moments for the amplitude density to the Newton variational principle for the amplitude density, when both are applied to the same reactive scattering test case with a single-channel Green's function and a sine basis, and set 1 is used for the other parameters in the variational calculations. The variational results converge to three significant figures with 9 sine functions, whereas the
TABLE VIII: Comparison of Nonreactive and Reactive Transition H2(v = O j ? H on the PK2 Probabilities for H H2(Y = O j =O) Surface with Restricted Internal State Expansion by Variational Calculations i i' set 1 set 2
+
nonreactive
0 2 0 1 2 1 1 2 2 2
0 0 0 0
reactive
o
nonreactive reactive
1 1 1 2 2
nonreactive reactive
+
2.30 X 7.43 X 5.70 X 4.41 X 3.63 x 9.76 X 4.69 x 2.77 X 2.40 X 2.26 x
lo-' 10-I 10-3 lo-' 10-3 low3 lo-' 10-3
2.31 7.42 5.63 4.42 3.58 9.76 4.71 2.75 2.41 2.20
X X X X
lo-' lo-'
IO-'
x 10-3 X
lo-'
x 10-3 lo-' lo-' x 10-3 X X
a E = 0.60 eV; 15 radial translational basis functions. The two sets of other parameters are explained in Table VI.
agreement with the RMP result. Table I1 shows that with 20-60 Gaussians per channel in the Schwinger method only the probabilities greater than 0.05 are well converged. Table 111 shows we obtain somewhat faster convergence with the Schwinger method when we employ energy-adapted basis functions. Tables IV and V are for the variational principle for the amplitude density, with multichannel and single-channel distortion potentials, respectively. Based on extensive convergence checks we found that the centers of the Gaussians must cover the range 5 . 0 - 9 . 9 ~when ~ a single-channel distortion potential is used but only 5 . 6 - 6 . 5 ~when ~ a multichannel distortion potential is used. The reduced coordinate space that must be spanned by the L2 basis in the latter case is a very dramatic finding, and as a result Table IV shows that the generalized Newton variational method combined with a multichannel distortion potential can yield ac-
-
TABLE IX: Comparison of Nonreactive and Reactive Transition Probabilitiesfor H + H2(v=Oj) H + H 2 ( v = O j ? on the PK2 Surface with Restricted Internal State Exoansion bv Variational (GNVP) and Nonvariational (MMAD) Methods" nonreactive j=O, j k o
N(G)~
MMAD
GNVP
MMAD
5 7 9
1.53c 2.39 2.32 2.30 2.30 2.30 2.30
2.36 2.31 2.30 2.30 2.30 2.30 d
8.00 7.25 7.42 7.40 7.41 7.43 1.43
10
12 15 20
MMAD 8.14' 5.00 5.21 6.03 5.75 5.71 5.69
j=1, j k l
GNVP
MMAD
9.36 7.37 7.42 9.78 9.78 7.43 7.43 9.75 9.76 1.43 7.43 9.76 9.76 d reactive ( ~ 1 0 ~ )
GNVP 5.43 5.56 5.69 5.69 5.69 5.70 d
MMAD 8.73 4.09 4.1 1 4.61 4.43 4.42 4.40
GNVP 4.35 4.36 4.41 4.41 4.41 4.41 d
j=2, j k 2 GNVP 9.77 9.17 9.76 9.76 9.76 9.16 d
MMAD 1.78 2.48 2.41 2.41 2.40 2.40 2.40
GNVP 2.46 2.41 2.40 2.40 2.40 2.40 d
j=1, j k 2
j=O, j'=1
j=O, j k o 5 7 9 10 12 15 20
(XlO)
j=O, j k 2
MMAD 4.86 3.23 3.58 3.67 3.61 3.64 3.63
GNVP 3.38 3.55 3.63 3.63 3.63 3.63 d
MMAD 13.1 4.21 4.37 4.90 4.72 4.70 4.69
GNVP 4.57 4.63 4.69 4.69 4.69 4.69 d
MMAD 9.06 2.47 2.73 2.77 2.75 2.77 2.71
j=2, j k 2
GNVP 2.69 2.74 2.11 2.17 2.77 2.77 d
MMAD 1.94 2.07 2.26 2.25 2.26 2.26 2.26
GNVP 2.07 2.22 2.26 2.26 2.26 2.26 d
' E = 0.60 eV. b N ( G ) is the number of translational basis functions used. CThismeans 1.53 X lo-'. dThis calculation was not done. eThis means x 10-3.
8.14
+
TABLE X: H HZGround-State Reaction Probabilities on LSTH Surface, Summed over Final Rotational States, for Several Energies Compared to Previous Values energy, eV 0.500 0.600 0.625 0.650
WSL (ref 31) 3.6 X 1.2 X lo-' 2.0 X lo-'
cs (ref 32a) 1.0 x 10-4 2.5 X lo-'
S (ref 32b) 2.4 x 10-4 3.8 X lo-*
WL (ref 33)
LV (ref 34)
5.0 X 1.3 X lo-' 2.5 x 10-I
2.3 X lo-' 4.6 x 10-1
MMAD" (ref 12b) 1.8 X lo4 3.6 X 1.0 X 10-I 2.0 x 10-1
Method of moments for the reactive amplitude density with single-channel distortion potential. density with single-channel distortion potential. See Table VI.
run l C 1.8 X lo4 3.6 X 9.8 X 2.0 x lo-'
GNVPb run ZC 1.8 X 3.7 X 1.0 X lo-' 2.1 x 10-I
run 3' 1.8 X 3.6 X 9.8 X 2.0 x 10-I
Variational method for the reactive amplitude
3214 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Schwenke et al.
TABLE XI: Nonreactive Transition Probabilities P,, ( v = v' = 0, E = 0.6 eV) for H vmax
0
1
2
run 1 2 3-5 6 7 8 9 10
3
16
5
b
Po0 0.142 0.152 0.151
Po2 0.764 0.772 0.772
p04
p22
PI3
PI5
0.0165 0.0161 0.0161
0.134 0.141 0.141
0.340 0.342 0.342
0.220 0.222 0.222
1.54 x 10-4 1.57 X 1.58 x 10-4
0.0671 0.0675 0.067 1 0.0669 0.0674
0.747 0.747 0.746 0.746 0.747
0.0190 0.0191 0.0191 0.0191 0.0192
0.0737 0.0739 0.0736 0.0735 0.0736
0.0307 0.0310 0.0310 0.0310 0.031 1
0.230 0.230 0.230 0.230 0.230
2.00 2.02 2.06 2.06 2.03
10
0.0542 0.0544 0.0545 0.0544 0.0544
0.732 0.732 0.732 0.732 0.732
0.0199 0.0201 0.0201 0.0201 0.0201
0.0643 0.0643 0.0642 0.0641 0.0642
0.0302 0.0305 0.0305 0.0305 0.0305
0.229 0.229 0.229 0.229 0.229
2.19 x 2.23 x 2.23 x 2.23 X
5
0.0540
0.73 1
0.0202
0.0639
0.0305
0.229
2.20 x 10-4
0.0538
0.741" 0.739"
0.0196" 0.0195"
0.0690
0.0304' 0.0303"
0.226
2.22 x lo4' 2.15 x 10-40
"ax
0 2 3, 4 4 5 9 12 5
11 12 13 14 15
+ H, on the PK2 Surface
4 5 6 10
'When the results for PlJ8and PJ9differ to the number of significant figures tabulated, both TABLE XII: Reaction Milliprobabilities 103P,,,( v = umax
0
run 1 2 3 4 5
"ax
0 2 3 4 4
Y'
= 0, E = 0.6 eV) for H
p24
x 10-4 x 10-4 X
lo4
x 10-4 x 10-4 2.16 x 10-4 10-4 10-4 10-4 lo4
Reference 38.
are given.
+ H2on the PK2 Surface
mPm 16.2 7.89 7.96 7.96 7.96
mpo, 7.39 13.7 13.8 13.8 13.8
mP02
mP03
14.5 7.03 7.10 7.10 7.11
0.742 1.34 1.36 1.36 1.36
mP12 6.51 12.0 12.1 12.1 12.1
mp,, 2.86 2.23 2.26 2.26 2.26
0.063 1 0.115 0.117 0.117 0.117
mP14
1
6 7 8 9 10
4 5 9 12 5
21.7 21.7 21.7 21.8 21.7
37.8 37.8 37.9 38.1 37.8
19.8 19.8 19.8 19.9 19.8
3.88 3.92 3.93 3.94 3.93
33.9 34.0 34.2 34.3 34.0
6.48 6.58 6.59 6.61 6.59
0.358 0.366 0.368 0.368 0.372
2
11 12 13 14 15
4 5 6 10 10
25.0 25.0 25.0 25.1 25.1
43.6 43.7 43.7 43.9 43.8
22.9 22.9 23.0 23.1 23.1
4.52 4.59 4.60 4.63 4.62
39.3 39.5 39.6 39.8 39.7
7.58 7.71 7.74 7.76 7.76
0.422 0.439 0.439 0.446 0.446
3
16
5
25.1
43.9
23.1
4.62
39.7
7.73
0.443
5
b
24.9
41.5" 42.2"
22.0" 21.9"
4.21" 4.25n
36.8" 36.1'
6.99" 6.94'
0.41 1' O.41On
'When the results for PJ1,differ from those for Pi:, to the number of significant figures tabulated, nonvariational method requires more than 15. With 5 basis functions, the variational results are converged to 9% or better, and with 7 basis functions to better than 2l/,%. In contrast the nonvariational calculations with these bases show errors as large as 227% and 12%, respectively. Thus the variational formulation converges more rapidly than the nonvariational for reactive as well as nonreactive scattering. Table X provides a near-state-of-the-art application of the variational method. As shown in columns 2-6 of this table, our previous -C2 results12bfor H Hz in the threshold region are in poor agreement with all previously reported results except those of ref 31, and even comparing to ref 31 they show a 14-19% disagreement at one energy. The last three columns of this table demonstrate excellent convergence of the present variational method with a sine basis, and comparison to the fourth last column shows excellent agreement, to 1 in the second significant figure of previous nonvariational resultslzb with a primitive Gaussian basis. We are very encouraged that two such different methods not only converge stably but to the same answer, and we believe this demonstrates that the results in the last three columns are converged to f l in the second significant figure. As mentioned in section 4, one can also consider more general choices of basis functions. For example, since the basis functions need not satisfy scattering boundary conditions, some or all of them could be expressed in other coordinate systems. Alternatively, we could imitate some of the effects achieved in other contexts with adiabatic basis functions by using in eq 46, instead of &&x), a vibrational-rotational-orbital function which is different for each translational basis function. A third alternative
+
*
both are shown. bReference 38
is to base contraction coefficients on the results of calculations involving fewer channels. This would work for both variational methods discussed here. Another possibility for future work is to use energy-adapted basis functions for the amplitude density. Note that multiplying \ k n o by R&*% and integrating over all x produces zn,no(R)= Cen,n(R)Fnn,(R)
(127)
Po= (l/R)C@n(X)znno(R)
(128)
n
where fl
It is the quantity znnowhich we expand in Newton's variational method. Thus if we have good approximations for F,,,,, we can use them to approximate the zm. Now, if we consider only open channels for simplicity, Fnno has the asymptotic form
-
Fnno
R--
k,,-1/2[6n,o sin (k,R
- l,,r/2)
+
KnnoCOS (k,R - l , r / 2 ) ]
(130)
thus we need to know Kw However, we can approximate it with R:nowhich must be calculated anyway. -Our scheme is then to take linear combinations of pairs of the Xfifl with K ~ : closest to k,2 to produce functions which behave approximately as sine and cosine for large R . We then produce an approximate Fw by using eq 130 for all R with Rtno.This is then substituted into eq 127 to get an approximate zn4 which is a linear combination of the
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3215
Variational Methods for Molecular Collisions TABLE XIII: Reaction Milliprobabilities 103qf ( Y = ~
a
method ref 38 MMAD GNVP GNVP
P
a
uncoupled rot, coupled fully coupled
N 117 78 84 117 117
Y'
= 0, E = 0.6 eV) for H + H2on the PK2 Surface R?(Un) A(4 Pno POI Po2 pal
mmax A(uo-~)
31 10 4 4 6
100.0 21.78 7.42 7.42 14.55
0.967 2.46 2.2 2.2 2.0
0.14 0.30 0.367 0.367 0.367
41.9" 43.2 43.9 44.1 44.1 44.1
22.0 23.5 23.1 23.3 23.3 23.3
4.23" 4.79 4.62 4.65 4.66 4.66
PII 71.3 75.8 76.3 76.6 76.7 76.7
P12
PI^
36.4" 40.8 39.7 40.0 40.0 40.0
6.96" 8.11 7.73 7.80 7.81 7.81
p14 0.410 0.399 0.443 0.450 0.451 0.451
Average of P,y and Pyi.
A,, but with R-dependent coefficients (this is because enddepends on R). We then reexpand this approximate z,, in the A, basis to get constant coefficients which then are the coefficients of the energy-adapted basis functions for the amplitude density. It would be interesting to test this scheme. The results for the final set of reactive scattering calculations are presented in Tables X I and XII. These tables refer to H + H2 reactive scattering on the PK2 surface, and we consider the convergence of the radial translational basis in Newton-type variational calculations with rotationally coupled distortion potentials. First consider the bottom two lines of Table XI and the bottom three lines of Table XI1 where we compare our converged results to close coupling calculations3*of Schatz and Kuppermann. The results agree within a few percent (the calculations of ref 38 are only claimed to be converged within about 5%). Next consider the convergence with respect to mmax,which is the number of Gaussians per channel, for various values of u,,. Runs 1-5 show convergence of 2% or better in the state-to-state =1 probabilities for u,,, = 0 and mmax= 2 , runs 6-9 for u, show convergence of 1% in the state-testate probabilities greater than 3 X lo-" with mmax= 4 and convergence of 2% or better in all tabulated probabilities with M,,, = 5, and runs 11-15 for u,, = 2 show 2% or better convergence for tabulated transition probabilities greater than with mmax= 4 and for all tabulated transition probabilities with mmax= 5. These convergence rates of 2-5 radial translational basis functions per channel are extremely encouraging compared to the 20-50 radial translational basis functions per channel mentioned in the Introduction. Since the slow step in big calculations scales as @, these improvements could lead to a savings of about 3 orders of magnitude in largescale applications. This confirms the promise of earlier tests" of the Newton variational method for potential scattering and justifies the present effort to generalize it to reactive scattering. Further examination of runs 9-16 shows that all tabulated transition probabilities greater than lo-' are converged to better than 3% not only with respect to the radial translational basis but also with respect to the internal basis. Finally Table XIIT compares a well-converged result obtained by using different methods. The final results given for the generalized Newton variational principle (GNVP) are believed converged to about f l in the third significant figure with respect to all numerical parameters (finite difference grid, quadratures, internal basis, and Gaussian basis). With 31 Gaussians per channel, the method of moments for the amplitude density (MMAD) agrees within 4%for reaction probabilities greater than The GNVP results with two kinds of distortion potential agree to 1% or better for these probabilities. Comparison of these well-converged results to those in Table XITI confirms that we can achieve 2% convergence for the rotationally coupled distortion with m, = 5, whereas with the fully potential probabilities > coupled distortion potential we can achieve 1% convergence for these probabilities with mmax= 4. We note that not only do the present results show rapid convergence with respect to the radial translation basis, they also show good convergence with respect to internal states, especially as compared to calculations involving the channel-permuting array.39 This is probably a consequence of including basis functions in all arrangements and fully coupling them-a procedure with roots ~~~
24.9 24.4 25.1 25.2 25.2 25.2
~
(38) Schatz, G. C.; Kuppermann, A. J . Chem. Phys. 1976, 65, 4668. (39) Kouri, D. J.; Baer, M. In The Theory of Chemical Reaction Dynamics; Clary, D. C., Ed.; Reidel: Dordrecht, Netherlands 1986; p 359.
in the Hartree-Fock method2S4 and first suggested for chemical reactions by Micha4' and Miller.' (This is called the Fock coupling s ~ h e m e . ' ~ .Miller ~ ~ ) and Jansen op de Haar3 recently proposed a Fock-coupled Kohn method for reactive scattering, and it is very encouraging that this method has now been used43to reproduce our results (last four columns of Table X) for H H2 reaction probabilities on the LSTH surface.
+
7. Summary New practical schemes for quantum mechanical calculations of transition probabilities in atom-molecule or molecule-molecule collisions are presented. We also discuss new procedures for generating the required integrals over the Green's functions and some computational considerations in choosing among basis sets for variational approaches. The new methods presented here are applicable to both inelastic and reactive collisions. Each method involves a two-potential formulation involving a conveniently chosen distortion potential; the rest of the interaction is called the coupling potential or reactive potential, and the difference between the full wave function or amplitude density and the wave function or amplitude density for scattering by the distortion potential is called the residual wave function or reactive amplitude density. We consider four combinations of approaches, in each of which the scattering by the distortion potential is included by integrals over a Green's function, and the scattering by the residual potential is treated variationally by expansion of the residual wave function or a related quantity in a square-integrable basis set. The integrals over the Green's function may be generated directly by a matrix finite difference boundary value method without explicit construction of irregular solutions for the distortion potential, and the variational problem for the residual interaction is reduced to a single solution of a large set of linear equations. The general method involves two prominent. choices: between single-channel or multichannel distortion potentials and between the Schwinger variational principle for the expansion of the residual wave function or the Newton variational principle for the expansion of reactive amplitude densities. We present numerical calculations testing all four of the possible combinations. In particular all four methods are applied to a nonreactive scattering problem, and they are all found to enhance the convergence with respect to nonvariational formulations presented previously. Further enhancement of the rate of convergence with respect to increasing the basis set is obtained by using energy-adapted contracted basis sets. On the basis of these calculations and comparisons we conclude that the generalized Newton variational principle with a multichannel distortion potential is the most efficient method. The combination of a single-channel distortion potential and the Newton variational principle is also applied to reactive scattering using the Fock scheme to couple the arrangement channels, and again the variational principle leads to faster con(40) See, e.g.: Seaton, M. J. Philos. Trans. R . SOC.London, A 1953,245, 469. Smith, K.; Henry, R. J. W.; Burke, P. G. Phys. Rev. 1966, 147, 21. Saraph, H.; Seaton, M. J. Shemming, J. Philos. Trans. R . SOC.London, A 1969, 246, 77. Harris, F. E.; Michels, H. H. Methods Comput. Phys. 1971, 10, 143. Nesbet, R.K. Adv. Quantum Chem. 1975,9,215. Truhlar, D. G. In Semiempirical Methods of Electronic Structure Calculation, Part B Segal, G . A,, Ed.; Plenum: New York, 1977, p 247. (41) Micha, D. A. Ark. Fys. 1965, 30, 411. (42) Schwenke, D. W.; Truhlar, D. G.; Kouri, D. J. J. Chem. Phys. 1987, 86, 2712. (43) Zhang, J. Z. H.; Miller, W. H., unpublished results, July 1987.
