Volume 137, number 6
CHEMICAL PHYSICS LETTERS
3 July 1987
QUANTUM REACTIVE SCA’ITERING IN THREE DIMENSIONS USING HYPERSPHERICAL (APH) COORDINATES. TESTS ON H+H, AND D+H, *
Gregory A. PARKER ‘, Russell T PACK, Billy J. ARCHER ’ and Robert B. WALKER Theoretical Division (T-12, MS J569), Los Alamos National Laboratory, Los Alamos, NM87545, USA Received 2 April 1987
Accurate CC calculations of 3D quantum reactive scattering using hyperspherical (APH) coordinates are reported for both zero and non-zero total angular momentum J and for both symmetric (H + HZ) and unsymmetric (D + HZ) mass combinations. Accurate 3D reactive scattering calculations are now possible for heavier systems and higher energies than ever before.
1. Introduction
We report the development and implementation of a general new quantum theory of reactive molecular collisions (rearrangement scattering) in three physical dimensions (3D) using adiabatically adjusting principal axis hyperspherical (APH) coordinates [ 11. These are the first accurate 3D calculations of reactive molecular scattering using any hyperspherical coordinate system to include non-zero total angular momentum J and to be applied to both symmetric ( H + H2) and unsymmetric ( D + HZ) mass combinations. Many calculations of reactive molecular scattering using hyperspherical coordinates have been reported for systems restricted to collinear geometries, but the only 3D hyperspherical calculation to our knowledge is the recent one of Kuppermann and Hipes [ 21 for the symmetric system H + H2 for J= 0. As we have shown elsewhere [ 31, the hyperspherical coordinate systems of Kuppermann [ 2,4] and of Delves [ 51 do not treat the arrangement channels equivalently, and they have a singular point right at one of the transition states for many systems; this makes accurate calculations with J#O difficult. In addition, virtually all [ 4,6-91 hyperspherical coordinate systems * Work performed under the auspices of the US Department of Energy.
’ Also at: Department of Physics and Astronomy, University of Oklahoma, Norman, OK 730 19, USA.
564
except the APH system require the use of half-integer angular momenta or (equivalently) discontinuous functions to construct off parity, odd J wavefunctions [ 3 1. The APH system [ 1] has the advantages relative to other hyperspherical coordinate systems that (1) it provides convenient mapping and visualization of potential energy surfaces and wavefunctions without requiring half-integer angular momenta or discontinuous functions; (2) all arrangement channels are treated completely equivalently, easing construction of accurate basis functions in all arrangements; (3) the body-frame (BF) axes follow the reaction smoothly and minimize coupling for systems with near linear transition states; and (4) it is directly applicable both to symmetric systems, such as H + HZ, and also to unsymmetric systems, such as Li+HF and F+ HZ. We believe it will allow benchmark coupled channel (CC) calculations for these heavier, chemically important systems. Relative to most theories that do not use hyperspherical coordinates, the present theory has the advantages that (1) all coordinate matching is by simple projection regardless of the atomic masses; (2) only regular solutions are required, and irregular solutions need never be generated, and (3) the most expensive parts of the calculation are energy independent, so that, once they are done, the scattering calculations can be rapidly carried out at the large number of energies needed to map out the reso-
0 009-2614/87/% 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
Volume
137, number
CHEMICAL
6
nances expected scattering.
to be important
PHYSICS
in reactive
2. Theory The APH coordinates, which were variationally derived elsewhere [ 11, turn out to be very similar to the coordinates of Smith [ 61 and Johnson [ 7) but differ from them in two significant ways. All three use BF axes tied to the instantaneous principal axes of inertia, but their [ 6,7] axes are appropriate near oblate symmetric top geometries, and the APH axes are appropriate near linear or prolate top geometries. All three systems use the same hyperradius p (O
=
_
‘=
wapp5
a
ificose -w2sin28
$
+T,+AJ:+BJ;+CJ:
J 2
(1)
Yax'
where the J, are the BF components of the total angular momentum, A= [pp2( 1 +sin@]-‘, B= (21(p2sin28)-‘, C= [~‘(l -sir@] -I, p is the threebody reduced mass, and T,, = - ( fi2/2,up2)
(
X-F---
4
a
sm28 ae
sin2Ba
ae
1
a=
+ -sin28 ax2 > *
{T,+t(A+B)fi=J(J+l)+[C-j(A+B)]fi2A2 + ( 15fi2/8pp2) + V(p, 8, x)}@$ x;P),
at each J and n (the BF z component of J) at many values of p (the centers of the propagation sectors) using finite element methods similar to those used by others [ 2, lo]. However, our modification [ 111 of the finite element engineering code ADINA [ 121 uses quadrilateral elements with 9 nodes per element and a very non-uniform mesh to make the mesh much denser in classically allowed than in classically forbidden regions. The full solution for a given total J, M (SF component of J), and parity p is expanded in the form, Y~JMP~ =4 s ~-5/=y;~(p) @fj(e, x;p)
Xm4(%
BY Y),
(4)
where the fi;rP, are normalized, good-parity Wigner rotation functions [ 31 of the Euler angles 0, /_I,and y. Eq. (4) is substituted into the Schrodinger equation [E-T-
V(p, e,x)] ‘PJMp”=O,
(5)
and a set (n = 1,2,...N) of radial solutions v/, regular at small p, are propagated from small p to large p using the VIVAS [ 131 method. Then, the solutions are projected onto functions of ordinary Delves [ 51 hyperspherical coordinates in each arrangement via a unitary transformation obtained by a simple 2D quadrature, and the boundary conditions, expressed in Delves coordinates, are applied directly to generate the scattering matrix S. The present calculations are full close-coupled or coupled channel (CC) calculations, and all matrix elements generated when eq. (4) is substituted into eq. (5) are kept. However, the APH equations strongly suggest a centrifugal sudden (CS) decoupling approximation [ 31 which we will explore in future publications.
(2)
The coordinates 8 and ,Y_ define the surface of the upper half of a sphere (the internal coordinate sphere or “hypersphere”). To get energy-independent basis functions @( 8, x;p) on the surface of this sphere, we solve a 2D equation,
= &J(P) @z(e,
3 July 1987
LETTERS
(3)
3. Test calculations and results To ensure that our programs are working correctly, test calculations were performed on the H+H,+H,+Hand D+H,+DH+H reactionsusing the Porter-Karplus (PK2) surface [ 141. Our future calculations will use better potential energy surfaces, but the present calculations were done to allow com565
Volume 137, number 6
CHEMICAL PHYSICS LETTERS
parison with existing accurate calculations [ 151. In the present computations, 60 surface functions 0 are determined at each of 60 to 120 values of p using 600 to 1200 finite elements. These functions include all open and several closed channels for all total energies below 1.5 eV. The CC calculations converge rapidly with respect to the number of surface functions used. The resulting S matrices are symmetric and, even at the highest energies, unitary to a part per thousand or better. In the present work, calculation of the surface functions @ is the expensive step; each 0 requires only 6 s of Cray-XMP time, but huge numbers of them are needed. Future optimization and vectorization will make that step much faster. Fig. 1 shows perspective plots of an H3 wavefunction on the surface of the “hypersphere”. The upper half surface of the sphere has been projected onto a plane, so that the center of these plots is the point 0=0 (the north pole); the edge of the plots is at 0 = n/2 (the equator). x = 0 is at the front of the plots, and x goes from -rt to IC (or 0 to 2~) as one.goes around the figure. Fig. la is the ground state wavefunction @?t( 8, x;p) at p= 5.0 a,. At this distance, the parts of the wavefunction in the different arrangement channels do not overlap. In each channel, motion along the tops of the ridges corresponds to rotation of HZ fragments. From the fact that the ridge heights vary little, one sees that rotation is only slightly hindered at this large distance. Motion across the ridges corresponds to vibration of H2 fragments and has the Gaussian appearance of ground state vibrational wavefunctions. At this distance, it is clear that the ground state surface function 018 consists mostly of a symmetric linear combination of H2( o=O, j=O) wavefunctions in each of the arrangement channels. There are 6 arrangement channels here instead of 3 to allow handling inversion easily as discussed earlier in this Letter. Parts of the wavefunction directly across the plot from one another are reflections of each other. Fig. lb is the same surface wavefunction 02 at the transition state distance, ~~3.166 ao. Now, the amplitude is not in the arrangement channels; instead, it is concentrated in the transition state regions between arrangement channels. It is this overlap between channels that allows flux from one arrangement into another in this APH approach, that 566
3 July 1987
a)
b)
Fig. 1. Perspective plots of the ground state surface eigenfunction 0% as a function of 13and ,y at two values of p. See the text for discussion. The values of p are (a) 5.0 and (b) 3.166 uO.
is, that allows reaction to occur. Note also here that rotation is very strongly hindered at this distance; the behavior of the wavefunction in both x (around the figure) and 8 (going away from the collinear edge) is basically vibrational (Gaussian). Sample results for total reaction probabilities into both final arrangements are shown in figs. 2 and 3. The three groups of curves in fig. 2a are for reactive transitions out of (o,, j,) = (0,O) into three final states,
3 July 1987
CHEMICAL PHYSICS LETTERS
Volume 137, number 6 0 24
:a)
O+H, -20
0
J=O
O”t Of (0.0)
A out of (0.2) -4 ii -
0
0 O”t Of (0.4)
-6.0 0
;
-6.0
-100
:b)
0 0.50
H+H2 J=l p=O * out Of (0.1)
z &
0.60
$ m 0
0.40
040
0 out
Of (0.5)
% 0.20
.bbSblW
1
0.60
0.70
ENERGY
0 50
ENERGY
A O”t of to,31
0.80
0.80
1.00
(eV)
Fig. 2. (a) Total reaction probabilities for specific transitions H+H,(o,,j,)-rH,(or,jr)+HforJ=Oand (o,,j,)=(O,O) plotted versus total energy in eV. The group of curves that are highest near the center of the figure (0.75 eV) are for (or, jr) = (0,l); the middle curves are for (0,2); and the lower curves are for (0,O). The present results for the three transitions are the open triangles, asterisks, and open diamonds, respectively; the results of Kuppermann and Hipes [ 21 are the open hexagons, the open squares, and the plus signs; and the results of Schatz and Kuppermann [ 151 are the dark squares. (b) The present total reaction probabilities from specific initial states into all possible product states versus total energy in eV for H + H,( o,,j,) +H2 + H with J= 1 and even inversion parity (p=O). The asterisks are for (o,, j, ) = (0,l) ; the triangles are for ( 0,3); and the diamonds are for (0,s). This case has odd reflection parity, and only A= 1 occurs, so that it is a measure of the contributions of non-linear geometries to the reaction.
(or, jr)=(O,O), (O,l), and (0,2), for H+H* reactions with J=O. In each of these three groups, there are three curves comparing the present results with earlier results [ 2,151. Overall, the agreement is good and within mutual uncertainties. The present results agree with those of Schatz and Kuppermann [ 151 slightly better than with those of Kuppermann and Hipes [ 21, differing from the results of ref. [ 151 by at most lOoh and usually by less than 5%. The pres-
0.60
0.7n
0.81 0.60
0.99
1.00 0.70
(eV)
Fig. 3. Log plot of the total reaction probabilities at low total energies from specific initial states into all possible product states for D+H*(n,, j,)+DH+H for J=O. The open hexagons are the present results for (u,, j,) = (0,O); the dark squares are the corresponding results of Schatz and Kuppermann [ 151. The open triangles are for (n,, j,) = (0,2), and the open squares for ( 0,4). The inset shows a linear plot of the present results for (n,, J, ) = (0,O) and includes higher energies than those given in ref. [ 151.
ent results are smoother with energy than both the earlier results, but, until we complete convergence tests with respect to all parameters, we cannot say that they are more accurate. Also, part of the small differences between the three calculations appears to be due to the use of different atomic masses for H in the different calculations [ 3,151. The dips in the curves near total energy E~0.97 eV are due to the onset of reactive scattering into u= 1 vibrationally excited products. Similar dips (not shown) due to v = 2 products are seen near 1.4 eV. Present results for total reaction probabilities into all final states from three initial states are shown in fig. 2b for non-zero total angular momentum (J= 1) and even inversion parity for the H+H2 reaction. This case includes only A = 1. Earlier results [ 15] (not shown) for this case were limited to E&0.6 eV and were nearly zero. We agree that they are very small for low E, but we note with interest that these reaction probabilities, which represent the contribution of non-linear transition state geometries to the reaction, rise to become fully as large as the J=O results at the higher energies shown. This behavior may provide a stringent test of reduced-dimensionality models of reactive scattering. Fig. 3 shows the total reaction probabilities into all 567
Volume 137, number 6
CHEMICAL PHYSICS LETTERS
final states out of three initial states for the unsymmetric-mass D + Hz+ DH + H reaction with .I= 0. The main figure shows a log plot of the probabilities at low energies, so that the present results can be compared for one of the transitions with the only other results available, those of Schatz and Kuppermann [ 151. The agreement is excellent. The inset shows a linear plot of the total reaction probability out of (Di,j,) = (0,O) to higher energies than were possible in the earlier calculations 1151.
4. Conclusions We conclude that 3D reactive scattering calculations using APH coordinates are viable and accurate and that they can be carried to higher energies and heavier systems than were formerly possible. Applications to reactions of genuine chemical interest are underway. We note that other groups [ 16,171 have published similar results in very recent months for isotopic hydrogen reactions using other potential surfaces and other new methods, and we believe that the present work, together with these other recent developments [ 2,16,17], signals the dawning of a new and brighter day in reactive scattering theory and in the understanding of chemical reactivity.
Acknowledgement
We thank George Schatz for sending us unpublished results, William A. Cook for much assistance with the finite element code, and Melvin L. Prueitt for assistance with computer graphics.
568
3 July 1987
References [ I] R.T Pack, Chem. Phys. Letters 108 (1984) 333. [ 21 A. Kuppermann and P.G. Hipes, J. Chem. Phys. 84 (1986) 5962. [ 31 R.T Pack and G.A. Parker, Quantum Reactive Scattering in Three Dimensions Using Hyperspherical (APH) Coordinates. Theory, to be published. [4] A. Kuppermann, Chem. Phys. Letters 32 (1975) 375. [5] L.M. Delves, Nucl. Phys. 9 (1959) 391; 20 (1960) 275. [6] ET. Smith, Phys. Rev. 120 (1960) 1058; J. Math. Phys. 3 (1962) 735; R.C. Whitten and F.T. Smith, J. Math. Phys. 9 (1968) 1103. [7] B.R. Johnson, J. Chem. Phys. 73 (1980) 5051; 79 (1983) 1906,1916. [ 81 C.A. Mead, Chem. Phys. 49 (1980) 23. [9] W. Zickendraht, Ann. Phys. (NY) 35 (1965) 18; Phys. Rev. 159 (1967) 448. [ lo] M. Mishra, J. Linderberg and Y. Ghm, Chem. Phys. Letters 111 (1984) 439; J. Linderberg and Y. ohm, Intern. J. Quantum Chem. 27 (1985) 273. [ 111 G.A. Parker, B.J. Archer, W.A. Cook and R. T Pack, unpublished. [ 121 K.J. Bathe, Massachusetts Institute of Technology, Acoustics and Vibration Laboratory report 82448-l (revision of December 1978). [ 131 G.A. Parker, T.G. Schmalz and J.C. Light, J. Chem. Phys. 73 (1980) 1757; G.A. Parker, J.C. Light and B.R. Johnson, Chem. Phys. Letters 73 (1980) 572. [ 141 R.N. Porter and M. Karplus, J. Chem. Phys. 40 (1964) 1105. [ 151 G.C. Schatz and A. Kuppermann, J. Chem. Phys. 65 (1976) 4668; private communications of unpublished work. [ 161 F. Webster and J.C. Light, J. Chem. Phys. 85 (1986) 4744. [ 171 K. Haug, D.W. Schwenke, Y. Shima, D.G. Truhlar, J. Zhang and D.J. Kouri, J. Phys. Chem. 90 (1986) 6757.