J. Phys. Chem. 1995,99, 7522-7529
7522
Influence of Atomic Fine Structure on Bimolecular Rate Constants: The C1(2P) -t- HCl Reaction George C . Schatzt Theoretical Chemistry Group, Chemistry Division, Argonne National Laboratory, Argonne, Illinois 60439 Received: August 17, 1994; In Final Form: December I , 1994@
In this paper we use quantum reactive scattering calculations to study the influence of atomic fine structure on the rate constants for atom-diatom reactions, using the C1(2P) HC1 reaction as an example. To do this, we calculate rate constants using two different methods, first a conventional single surface reactive scattering method based on propagating coupled channel equations using hyperspherical coordinates, and second a multiple surface method in which the electronic fine structure of the C1 is explicitly included in the calculation. In this second calculation, the Hamiltonian combines the nonrelativistic potential surfaces for the 2C and 211 states of ClHCl with a spin-orbit Hamiltonian in which the spin-orbit parameter is taken to be a constant and with Coriolis-induced electronic and rotational couplings. Rate constants are evaluated using the J-shifting approximation, using J = l/2 calculations as the reference. We find that the biggest effect of the fine structure is a static one which arises because the spin-orbit term in the Hamiltonian stabilizes the reagents and products differently than the 2Z1/2 saddle point. This causes the activation energy to increase by about 29% of the atomic 2P1/2 - zP3/2 energy difference. When this effect is factored out, the single and multiple surface rate constants are found to be in very good agreement (provided that the usual electronic statistical factors are applied to the single surface rate constant), with differences of less than 20% over the temperature range that we considered. We also study how the size of the spin-orbit constant influences this result, and we find that the overall rate is relatively insensitive to spin-orbit splitting (except for the shift in activation energy). Fine structure resolved rate constants for C1(2P1/2) are found to vary by many orders of magnitude as the spin-orbit splitting is varied between zero and the experimental value, corresponding to a change in the dynamics from statistical to adiabatic limits
+
+
I. Introduction The conventional approach' to the inclusion of electronic fine structure in the calculation of thermal rate constants involves multiplying a single surface rate constant that is calculated for the lowest electronically adiabatic potential energy surface by an electronic statistical factor that is simply the ratio of transition state to reagent electronic partition functions. This approach is based on transition state theory (TST), and usually it is used in TST calculations where the nuclear part of the rate constant is also expressed using ratios of partition functions. However, it has also been used in collision dynamics applications where the nuclear part of the rate constant is obtained from trajectory calculations or from quantum scattering calculations. Another way to treat the electronic fine structure contribution to thermal rate constants, which has occasionally been used: involves using the adiabatic correlations between reagent and transition state electronic states to determine the (thermally averaged) fraction of states which correlate with the reactive potential surface. This fraction is then multiplied by the single surface rate constant to determine the resulting rate constant. To go beyond the simple treatments of fine structure just discussed requires doing a dynamics calculation of some sort. There has been a significant amount of past work on fine structure effects in atom-molecule reactive collision^,^ but little of this has produced specific results that relate to thermal rate constants. Early work on the fine structure problem was based on semiclassical theory that was closely connected with studies of atom-atom collisions. Tully? for example, studied the +Permanent address and address for correspondence: Department of Chemistry, Northwestern University, Evanston, IL 60208-31 13. @Abstractpublished in Advance ACS Abstracts, May 1, 1995.
0022-365419512099-7522$09.00/0
F(2P1/2) H2 reaction using a combination of trajectories and time-dependent Schrodinger solutions and found that the cross section for this nonadiabatic reaction is about 10% of the adiabatic F(2P3/2 Hz)reaction. However, it is not clear from this study how the overall thermal rate constant is influenced, as Tully did not study the F(2P3/2) H2 reaction to determine if the m, = l/2 and 3/2 components have reactivity that is consistent with adiabatic correlations. Further semiclassical work on F(2P1/2) H2 was done by Lam and George: but again the overall thermal rate constant was not considered. Quantum scattering studies of fine structure effects in reactions like F H2 have been attempted on numerous occasions in the past, but the approximations used have prevented the determination of thermal rate constants. George, Zimmerman, and Baer6 studied the collinear dynamics of F Hz,including the 2 and ll potential surfaces that correlate to the reagent ground states and also including spin-orbit coupling, but Coriolis effects could not be described. Rebentrost and Lester7 studied F H2 in three dimensions using coupled channel methods that were exact except for restrictions on the potential surfaces that omitted the possibility of reaction. They found that the transition 2P1/2 2P3/2 occurs efficiently for certain (resonant) rotational transitions where spin-orbit splitting is balanced by rotational energy transfer. In this paper we present realistic quantum reactive scattering calculations that include reagent and product fine structure effects, with the aim of determining the influence of these fine structure effects on the thermal rate constants. The theory that we present is specifically appropriate for reactions of the type
+
+
+
+
+
+
-
A(2P)
+ BC - AB + C(2P)
where BC and AB are closed shell ('E) diatomics, and only the
0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 19, 1995 7523
Influence of Atomic Fine Structure on Rate Constants
22and 211electronic states of ABC that correlate to the ground terms of the reagent and product atoms are included. The computational approach is based on solving coupled channel equations for the reactive collision dynamics, using nuclear basis functions that are developed in hyperspherical coordinates. The Hamiltonian in this calculation includes all the relevant angular momentum couplings associated with electron spin, electron orbital, and nuclear orbital angular momentum. Except for the use of hyperspherical coordinates, the theory is very similar to that developed by Rebentrost and L e ~ t e r including ,~ the use of a diabatic electronic basis which is derived from the 2C and electronic states. We will specifically consider an application of this approach to the C1(2P) HCl reaction, but, in order to provide more general results, we include a study of how the thermal rate constant varies with the magnitude of the atomic spin-orbit splitting. The C1 HC1 reaction has been extensively studied8-10using quantum scattering calculations that are based on a single surface hyperspherical coordinate coupled channel method developed by Schatz.8 On the basis of these and other studies, it is known that a simple LEPS (London-EyringPolanyi-Sato) surface due to Bondi et ai." describes the thermal rate constants quite ac~urately,~ and, in addition, the single surface reactive scattering calculations require modest vibration-rotation basis sets for convergence (2-3 vibrations, 1220 rotations for the lowest vibration, and fewer for the higher vibrations). Further, it has recently been establishedlO that the centrifugal sudden8 and J-shift12 approximations are quite accurate, which means that, for single surface calculations, only the lowest partial wave (J = 0) needs to be considered to generate complete cross sections and rate constants. In the present application, we consider the full coupled channel (CC) description of the lowest partial wave (J = l/2) and then invoke the J-shift approximation to determine rate constants. Our multisurface calculations can obviously be used to determine much more information about the reaction dynamics than simply rate constants; this will be the subject of future publications. The basic theory used in these calculations is described in section I1 of this manuscript. Details of the application to C1 HC1, together with tests of convergence, are given in section 111, and the results of calculations are in section IV. Section V summarizes our results and conclusions.
+
+
+
11. Reactive Scattering Method A. Coordinates, Angular Momenta, and the Hamiltonian Operator. We use the notation of Rebentrost and Lester7 wherein the electronic orbital angular momentum vector of the atom A is denoted L, electron spin is denoted S, nuclear rotational angular momentum of the diatomic BC is denoted N,and the nuclear orbital angular momentum of A relative to BC is denoted 1. The total electronic angular momentum is denoted j, and the (electronic plus nuclear) total angular momentum is denoted J, so that j = L S and J = j N
+
+ +
1.
If we let R be the mass scaled13 atomic-diatom vector and r be the mass scaled diatom internuclear vector, then the Hamiltonian is given by
H =P/2p
+ 12/2pR2+ p2/2p + N2/2pr2 + He, + H,,
(1)
where p is the scaled reduced mass,13 P and p are the radial momenta associated with the distances R and r, respectively, Hel is the nonrelativistic electronic Hamiltonian, and H,, is the spin-orbit Hamiltonian. This Hamiltonian neglects the mass-
polarization terms in the electronic Hamiltonian? which are not likely to be of importance for the low energy processes that we consider. In the following treatment, we replace 1 by J - j - N in the centrifugal term in (1). Expanding the centrifugal term then leads to the expression 12/2uR2= (J2
+ j2+ N2)/2pR2 - (2J.j + 2J.N + 2N*j)/2pR2 (2)
The cross terms in (2) lead to three types of Coriolis coupling: orbital-electronic, orbital-rotational, and electronic-rotational. All of these terms will be evaluated accurately in the coupled channel expansion that we give below. The electronic Hamiltonian Hel in (1) will be assumed to be exactly given by a two state expansion in terms of the 2Z and 211 adiabatic electronic states that correlate to the ground states of the reagents and products. A diabatic representation is used to represent these states on the basis of a procedure that is similar to that of Rebentrost and Lester but somewhat simplified by using approaches that are familiar in atom-atom ~ c a t t e r i n g . ~ ~ , ~ ~ Details are given in sections 1I.B and 1II.B. The spin-orbit Hamiltonian H,,is assumed to have the form
H,, = AL*S= (A/2)('j2 - L2 - S2)
(3)
where the spin-orbit parameter A is assumed to be constant, independent of the internuclear distances, and equal to the value associated with the separated atom A. (If atoms A and C are different, then different d values would have to be used for the reagent and product regions.) For C1 HCl, d has the value -588 cm-', which is small but not negligible compared to the barrier height of roughly 3200 cm-'. The assumption that 1is constant is probably not very accurate for reactive collisions near the saddle point, where there is significant reorganization of the valence electronic structure. Accurate calculations of the geometry dependence of spin-orbit coupling16 typically show that [dl is significantly smaller in the region in which the Z-ll separation is larger than it is at long range. Fortunately, this problem is less important than might be imagined, as the spinorbit Hamiltonian has no effect to first order on the 2Z surface energy near the saddle point (see eq 25 below and subsequent discussion). Thus, this problem is probably not serious for determining thermal rate constants. The assumed spin-orbit Hamiltonian (eq 3) does produce one effect on the thermal rate constant that will be important to this paper. For a 2P atom, eq 3 predicts that the j = 3/2 state of the atom has an energy d/2, while j = l/2 has an energy -a. Since d is negative for a C1 atom, the j = 3/2 state is thus stabilized by 1/21dl relative to the nonrelativistic energy. In contrast to this, near the barrier top, the energy of the 221,2 potential is not altered (to first order) by the inclusion of H,,. Overall this means that the effective barrier to reaction is increased by1/21d (or equivalently by l/3 of the splitting between t h e j = 3/2 and j = '12 states of the atom) as a result of spinorbit interaction. For C1, this amounts to a 294-cm-' increase in the activation energy, which will have a sizeable effect on the thermal rate constant at room temperature. We will see later that treatment of spin-orbit coupling beyond first order leads to some reduction in the magnitude of this effect, but the basic idea is the same. This result is quite general for reactions involving open shell atoms, and it means that dynamics calculations which omit spin-orbit effects in the potential surface are subject to an error on the order of a substantial fraction (like 50%) of the spin-orbit constant.
+
7524 J. Phys. Chem., Vol. 99, No. 19, 1995
Schatz
B. Coupled Channel Expansion. We will use a body-fixed representation to develop the coupled channel equations. In this representation we let Q, be the projection of j along the body-fixed z axis, which we take to be the direction R. Similarly, A and Z are the projections of L and S, and QN is the projection of N. In this representation, the body-fixed electronic states @2j) are related to the spin and orbital parts of the electronic wave function by
(b) for the spin-orbit terms, we have
q;= +(j + 1) - L(L + 1) - S(S + l)]S, ti2
(11)
(c) for the Coriolis terms, we have
uc;
J(J =
+ 1) + F ( F + 1) - 2Q2
4,-
R2
(4)
where (llmll2m2ll3m3) is a Clebsch-Gordan coefficient and, in the pure precession limit, we take L = 1, S = l/2 everywhere. Now we couple the vectors j and N, calling the resultant vector F (F = j N). The resultant z projection of F is called Q. The resultant electron-nuclear wave functions associated with F and Q are written
+
INjFQ) =
C JNQN)@2j)(NQdQjlFQ)
(5)
QNQ,
where INQN) is the usual (nuclear) rotational wave function. Once we have the internal states INjFQ), it is not difficult to write down the coupled channel expansion for the wave function associated with each partial wave J and z projection M,as
where D is a rotation matrix that depends on the polar angles 8, @ associated with R, @ is an eigenfunction of the BC vibrational Hamiltonian, and g is an R-dependent coefficient that is determined numerically by solving the Schrodinger equation. For the present discussion, we assume that @ solves the Schrodinger equation associated with the isolated BC molecule as follows:
where v(r) is the diatomic internuclear potential. Substitution of eq 6 into the Schrodinger equation, using the Hamiltonian in eq 1, leads to the coupled channel equation d2g/dR2 = Ug
(8)
where g is a matrix of expansion coefficients and U is a matrix that may be expressed as follows:
Here, we use the collective index “t” to denote (vNjFR). In (9), the first term contains the total energy E, the second term includes the electrostatic potential couplings induced by the difference between Hel and v(r), the third describes spin-orbit effects, and the fourth describes Coriolis effects. Specific expressions for the second through fourth terms in eq 8 are as follows: (a) for the electrostatic terms, we write
G1 -- 2+DLnI(@vNI(NjFQlHelti
Q‘)6Q,Q’+116,Cf7‘SNAfSjj’Svv’
Q‘)
(I2)
+
where E*(J,Q) = [J(J 1) - Q(Q f 1)]1’2. Note that the Coriolis Hamiltonian in eq 2 has been simplified by use of F and the angular momentum eigenfunctions in eq 5 . Specific details of the electrostatic contribution to U for the case of C1 HCl are given in section II1.B. Parity decoupling may be introduced into the coupled channel equations using a generalization of what has previously been described for single surface ~cattering,’~ in which one uses parity adapted rotation matrices in eq 6. This leads to a factorization of the problem into two equal sized ones (for half-odd integer J values). C. Reactive Scattering and Hyperspherical Coordinates. The coupled channel equations given in eq 8 above are appropriate for the description of nonreactive atom-diatom scattering but must be modified for reactive collisions because the use of isolated BC vibration-rotation eigenstates for expanding the wave function is inappropriate. To treat reactive problems, we introduce the Delves hyperspherical coordinates,” following a prescription given previously for single surface reactions by Schatz8 and Koizumi and Schatz.ls In the following, we omit most details of the treatment that are the same as in the single surface problem. The Delves coordinates are defined using a hyperradius e equal to (R2 and an angle 7 equal to tan-’(r/R). The remaining hyperangles are the same as for Jacobi coordinates. Because of the strong tie with Jacobi coordinates, the Delves coordinates are strongly dependent on the arrangement channel. In the present application, we take advantage of this by defining diabatic electronic states in which the angular momentum quantization axis is the same as the R direction. We then include all of the arrangement channels a in the basis expansion. This leads to overlaps between electronic states defined with different quantization axes. However, the same electronic basis is used in each channel, so since the transformation between hyperspherical coordinates associated with different arrangement channels is simply a kinematic rotation,13 the overlap can be evaluated in terms of Wigner rotation matrices. The relevant overlap of basis functions is
+
+
where the primed functions refer to one arrangement channel and the unprimed to another. The overlap of rotation matrices is well-known to bel3
where d is a reduced rotation matrix and where A is the angle between V I ~ ~ ’ F ~ ’ ) I ~ ~ ~ ~ ) I the D R’s ~ ~for~ each ) arrangement channel (a function of the three internal coordinates of the ABC triat~mic’~).The overlap of electronic wave functions in (13) involves an (10)
Influence of Atomic Fine Structure on Rate Constants analogous formula:
J. Phys. Chem., Vol. 99, No. 19, 1995 7525 where got, is the rotational-electronic energy of the ClHCl complex for angular momentum J. In the present application, we approximate got, by BeI(1 l), where the correspondence between the total angular momentum J and the orbital angular momentum I is taken to be J = 1 l/z. (The conditions on the validity of this approximation are discussed below.) Be is the rotor constant associated with the 21: saddle point. Another complication that arises in single surface studies of C1 HC1 involves the inclusion of excited bends for which 152~1> 0, as these are not contained in the J = 0 calculation and would not be included via the J = 0 version of eq 20. In previous single surface studies,19 the contribution of these excited bends has been included through a "bend-shift'' approximation, which is similar to eq 20 but involves the bend energy rather than the rotational energy. We will use this approximation if calculating single surface rate constants. For multiple surface calculations, it tums out (see below) that the excited bends are already contained in the J = l/2 basis set, so bend-shifting is not done in this case. The accuracy of J-shifting and bend-shifting for single surface studies of C1 HCl has been studied in the past through comparisons with quantum scattering studies in which all partial wave probabilities were explicitly determined. l9 Comparison of rate constants indicates excellent agreement (few percent). In the present multiple surface studies we expect that the accuracy should be comparable, as the additional complications associated with electronic angular momentum should have little influence on the results given that the partial wave sum is dominated by large J terms (Le., J >> N and j ) , where J = I l/2 should be quite accurate. An important comparison that we consider in this paper is between the multiple surface rate constant that we have just defined and a single surface rate constant that is obtained by first doing scattering calculations for just the zI:surface and then correcting the resulting rate constant for multiple surface effects using a transition state theory expression that was mentioned in the Introduction. We denote Pm : as the single surface partial wave summed cumulative probability and kss(T) as the rate constant that is derived from this probability after correction for multiple surface effects. The detailed expression for this rate constant is
+
except that the normalization factor in (14) is omitted since the electronic wave functions are normalized. Substitution of eqs 14 and 15 into eq 13 leads to an overlap that now depends only on internal coordinates.
This integral is evaluated numerically. A similar procedure is applied to the evaluation of Hamiltonian matrix elements between basis functions that are defined relative to different arrangement channels. The coupled channel expansion in hyperspherical coordinates is given by eq 6 , except that the coefficient g is now a function of e and the eigenfunctions @ are functions of TI, with a reference Hamiltonian which needs to be chosen so as to keep the eigenfunctions localized to the desired arrangement channel region, as dicussed in ref 6. This leads to eigenfunctions that are nonorthogonal, resulting in overlaps between the nuclear basis functions that are the same as in the single surface problem. Propagation of the coupled channel equations, and matching to asymptotic boundary conditions, is the same as in the single surface calculations. The asymptotic solutions involve an obvious generalization of the single surface solutions in which the electronic states are factored in. This leads to a scattering matrix S which is labeled by initial and final values of the quantum numbers v, N,j , F, and s2 and the arrangement channel a. It is not difficult to use eq 5 to transform this into a representation which uses the quantum numbers v, N , QN, j , and Q j (useful for determining probabilities associated with specific rotational and electronic alignments and orientations), but, for the purpose of calculating the partial wave cumulative reaction probability that is needed in determining rate constants, one only needs to evaluate the sum
eUm
+
+
+
+
where the arrangement channel indices a and a' are chosen to be appropriate for the reaction of interest. D. Thermal Rate Constants. The multisurface thermal rate constant k(7') may be determined from the cumulative probability using the standard formula where Q:: is the transition state (ts) electronic partition function.
111. Application to C1 where Q&(n is the reagent nuclear partition function per unit volume (describing translational, vibrational, and rotational motions) and Q:Lag is the corresponding reagent electronic partition function. The cumulative reaction probability P,, in this case is the sum over partial waves of eU,:
In the present application, we have used the J = l/2 reaction probabilities to determine probabilities for the higher partial waves using the J-shifting a p p r o ~ i m a t i o n ~ ~ q ~ ~
+ HCl +
A. Potential Surface. In our application to C1 HC1, we have assumed that the zXsurface is given by the LEPS surface of Bondi et al." The parameters in this surface are such that single surface quantum scattering rate constants agree with experiment without including spin-orbit effects. This means that the present calculations which include spin-orbit effects will be in poorer agreement with experiment, but we will not attempt to correct for this here. The 211 surface is also represented as an extended LEPS function,20 with the Sat0 parameters taken to be zero for HC1 and Clz. These parameters give a barrier height of 0.872 eV for the 211surface, located at an HC1 distance of 2.860~0and having a linear ClHCl geometry, which is to be contrasted with a 22surface barrier of 0.371 eV at an HCl distance of 2 . 7 7 2 ~ . These two surfaces are in qualitative agreement with surfaces reported by Duggan and
7526 J. Phys. Chem., Vol. 99, No. 19, 1995
TABLE 1: Electronic Matrix Elements jQ, = +y2,+ 3 t 2 jQ, = + 3 t 2 , j’Q,’ j’Q,’ jQ,’ j’Q,’ j‘Q,’ j’Q,’
= V2, +V2
= +3t2,+ I t 2 = +’I2,
-It2
= +Y2, -Y2 = +V2,
+It2
= +V2,
-It2
Vw
- ‘15V20
Schatz
+It2
jQ, = +V2,
-112
0 0
0 0 0 VOO
(-fi/5)V20
0 0
0
(+&/~)VZO
0
vw + ‘I5v20
0 0 0 0 0
jQ, = + 3 t 2 ,
0 0
vw + ‘tsv20
Grice21 using diatomics-in-molecules methods. Previously, Yamashita and MorokumaZ2determined a limited portion of these two surfaces using ab initio calculations and found that the 2rI surface drops below 2Z at long range, corresponding the van der Waals well regions in the reagent and product valleys. The present surfaces match the Yamashita and Morokuma calculations in the repulsive regions but do not describe the van der Waals region. We have done calculations with a 211 surface in which a van der Waals well was added to match the Yamashita and Morokuma results, and we found that this curve crossing increases the rate constant at 300 K by about 25%. For the present purposes this is a small enough effect that we can ignore it. B. Electrostatic Coupling Matrix. With the assumption of just two potential surfaces, it is possible to evaluate the electrostatic coupling matrix in eq 10 quite simply as follows. We first introduce eq 5 into eq 10 so that matrix elements. involving the electronic states bQj) may be determined, and then nuclear matrix elements of these may be separately calculated. If we assume that the primed and unprimed variables in eq 10 refer to the same arrangement channel a,then the Wigner rotation matrices are orthogonal and eq 10 reduces to
He, = v,
+ V2$2(COS
6)
(23)
where 6 is the polar angle between the R axis and the electron coordinate z axis and where we take14
+
V, = l / g ( ~ x 213,)
(24a)
= 54(Ex - E,)
(24b)
v20
Equation 23 is designed so that diagonalization of it in a basis of p-orbitals will give the Z and ll energies Ex and En when the spin-orbit and Coriolis terms are neglected. Of course in the present case these terms are included, and this leads to coupling between the surfaces. If we evaluate the matrix elements of eq 23 using the electronic angular momentum eigenfunctions U s l j ) in eq 22, it is not difficult to determine the electrostatic coupling matrix given in Table 1. This shows that the states with lS2jl = 3/2 are not coupled to the others by electrostatic effects and correspond simply to the 2113/2 states throughout the reaction. By contrast, the two Q, = 1/2 states are coupled by the electrostatic
- ’/SV20
0 0
jQ, = + ‘ t 2 .
0 ( -&/5)V20 0 0 vw 0
+It2
jQj = + ‘ t 2 , - I t 2
0 0 (+d%)Vzo
0 0 VoO
Hamiltonian and upon diagonalization of this portion of the matrix in Table 1 yield the 2Z1~2and 2111/2 states. C. Effect of Spin-Orbit Hamiltonian on Born-Oppenheimer Potentials. If the spin-orbit Hamiltonian from eq 3 is included with the electrostatic terms in Table 1, the only additions are along the diagonals, with -A being added to the j = l/2 terms and A/2 being added to the j = 3/2 terms. If the sl, = l/2 portion of the resulting matrix is diagonalized, one obtains the roots
=v
E
--+-V2,&-V2, A 1 0 0 4 1 0
Zk
10
[
:4, ’4”
1+--+-7
;]1/2
(25)
The lowest energy root in eq 25 determines the spin-orbit modified 2Z potential surface. In the asymptotic limit where VZO= 0, this reduces to Ex A/2, while, in the limit in which VZOis large compared to A (approximately the case near the saddle point), it reduces to Ex (Le,, spin-orbit is quenched on the Z surface). This means that the barrier to reaction, namely, the difference between the transition state and asymptotic energies, is increased by 1A1/2 by the presence of the spinorbit Hamiltonian. A more accurate evaluation can be obtained by expanding eq 25 through order (A/V~O)~, yielding a barrier charge equal to 1A1/2 5A2/(6V20). Since V20 is negative at the barrier top, the second term in this expression serves to reduce the increase in the barrier to reaction arising from the spin-order Hamiltonian. For C1 HC1, using the LEPS surfaces described previously, the first order correction to the barrier is $294 cm-’, while the second order correction is -43 cm-’, so that the overall correction is +25 1 cm-‘. One can also evaluate eq 25 in infinite order, giving an overall correction of 255 cm-’, which is 43% of A. Since the spin-orbit correction to the Z potential surface varies slowly with nuclear coordinate, we expect that its primary effect on the reaction dynamics should be to shift the threshold for reaction up in energy by an amount equal to the correction at the saddle point. We will test this prediction in section IV below. D. Basis Sets and Numerical Parameters. The multisurface coupled channel calculations are extremely time consuming, so it is necessary to use as small a vibration-rotation basis set as possible. For the present calculations, we used single surface calculations to find a small basis that still gives accurate results for energies near the reactive threshold, and then we used the same basis for the multisurface calculations. Note that the complete set of 2P electronic states appropriate for J = It2 was used in all calculations, and a centrifugal sudden approximation was not used to reduce the number of projections explicitly coupled in the calculations. Most of our multisurface calculations were done using a basis of 172 states. This basis consists of rotational states N = 0-1 1 for vibrational state v = 0 and N = 0-3 for v = 1 for each of the two arrangement channels that are important for C1+ HC1. Table 2 gives the specific choices of rotational and electronic quantum numbers used in the basis. To test converge of the
+
+
If the primed and unprimed variables refer to different arrangement channels, then the overlap of rotation matrices yield d-matrices, as in eq 16, that should be inserted into the middle of eq 22. We next express the electronic Hamiltonian in terms of p-orbitals that are attached to the reagent (or product) C1 atom, using the expre~sion’~
-3t2
+
J. Phys. Chem., Vol. 99, No. 19, 1995 7527
Influence of Atomic Fine Structure on Rate Constants
TABLE 2: Rotational and Electronic Quantum NumbertP Associated with J = ‘/z N
i
0 0 0 0
312 312
1
312 312 312 312
512
V2 V2
112
112
I12
1 1 1 1 1 1
‘I2
1
I12
1 1 2
‘12
F 312 312 ‘I2 ‘I2 512
312
V2 ‘I2 312 312 ‘I2
8
‘12
1 -1 0 0
2
I12
’I2 lI2
‘I2 -‘I2
Y2
V2
2
312
’I2
2
Y2
Y2
2
312 312 ’I2 ‘I2 ’I2 ‘I2
‘12
-‘I2 ‘I2 -‘I2 ’I2 -‘I2
312
2 2 2
2 2
v 2
‘12
v2
- ‘I2
V2
‘I2
312
-112
’12
Y2 -V2
-112
- ‘I2 I12
312
-V2
38 86
0.0173 0.0107 0.0094
110
0.4
0.5
E(eV)
Figure 1. Cumulative reaction probability versus energy E (eV) for C1 + HC1. The solid curve shows the 80 state single surface results, while circles show 32 state single surface results. 5.0
‘12
.
I
-112
-112
‘I2 ‘12
-‘I2
TABLE 3: Cumulative Reaction Probabilities for C1 HCI, E = 0.40 eV basis 134 172
0.0 0.3
V2
The j , QN, and Qj quantum number combinations for N > 2 are analogous to those for N = 2. Note that QN and Q, are not in oneto-one correspondence to F and Q for a given Nand j .
em
a
‘I2 ’I2 -‘I2 -312
(I
basis
I
‘12
-1 1 0 0 1 -1 0 0
/n
- ‘I2 -’I2
-2 1 -1
112
‘I2
‘I2
-‘I2
0 0
‘12 -112
‘I2
‘12
-‘I2
-1 -1 1
-‘I2
V2 V2
2 2 2
0 0 0 0 1
‘12
-‘I2 ‘I2 -’I2
Q,
QNb
‘I2 -‘I2 ‘I2 - ‘I2 ’I2 - ‘I2
3.0
4.0
1
21 6 1 n 3.0
P
+
em 0.0099 0.0120
multisurface results with respect to basis, we also did calculations for basis sets of 38, 86, 110, and 134 states, each with just one vibration and rotations N = 0-3,O-7,0-9, and 0-1 1, respectively. Table 3 shows the cumulative reaction probabilities for these basis sets, as well as for the 172 state basis. All of these results refer to just one of the two equivalent parity subblocks associated with the *P states of C1. The spin-orbit parameter I for all these calculations is 588 cm-’, and the total energy E used in this table is 0.40 eV, which is close to the effective threshold for the reaction. Evidently, the cumulative probability is well converged with respect to the rotational basis at 134 states. The comparison of 134 and 172 states shows that the change in reaction probability associated with adding a vibrational state is relatively small at 0.40 eV. To study convergence properties for still larger numbers of basis functions, we have done single surface J = 0 calculations, one with 32 states, corresponding to the same nuclear basis as the 172 state multiple surface basis given above, and one with 80 states, corresponding to N = 0- 19 for v = 0, N = 0- 13 for v = 1, and N = 0-5 for v = 2. The comparison of cumulative reaction probabilities for these two basis sets for energies near the threshold for reaction is shown in Figure 1. This shows that the two results are in good agreement (few percent or better for energies at or below 0.50 eV, and they agree to within 25% for energies up to 0.60 eV. We will show later (section 1V.B) that this is adequate for making thermal rate constants for temperatures up to 400 K agree to within a few percent. This result refers to the single surface calculation, but, in view of the multisurface comparisons given in the previous paragraph,
V
1.0
L
p’.
I
‘ I
0.0 0.3
I
,//
/
,
-, /
0.4
0.5
0.6
EW)
Figure 2. Cumulative reaction probability versus energy E, showing 172 state multiple surface result for 1 = 0 (solid curve) and 1 = 588 cm-’ (dashed curve). Also shown are the 32 state single surface result (circles), and the Q N = Qhp = 0 component of the 1 = 0 probability (dotted). it seems likely that the 172 state basis will be adequate for determining multisurface rate constants within a few percent.
IV. Results and Discussion A. Cumulative Reaction Probabilities. Figure 2 presents multisurface cumulative probabilities (just one of the two identical parity components) for J = ‘/2 from our 172 state calculations. The solid curve presents the calculated probability for I = 0, while the dashed curve shows the corresponding probability for A = 588 cm-’. The circles show the single surface 32 state J = 0 results (the same as in Figure l), and we find that these are in good agreement with the corresponding il = 0 multisurface results right at threshold, but not at higher energy. This means that the multisurface reaction dynamics is dominated by single surface effects right at threshold. In addition we see the anticipated effect of making I nonzero, namely, that the reaction probability is uniformly shifted up in energy. The magnitude of the shift is about 0.03 eV, which is in quantitative agreement with the estimates made in section 1II.D based on adiabatic potentials. Note that the probabilities in Figure 2 increase from an effective threshold near E = 0.40 eV to a value close to unity and then level off for about 0.02 eV before increasing further
7528 J. Phys. Chem., Vol. 99, No. 19, 1995
Schatz
a
-15,
-I3.O
.16.0
,
,
.
,
I
I ’
-17.0 2.0
I 2.5
3.0 IOOOWT
3.5
.18
4.0
Figure 3. Arrhenius plot of rate constant (cm3/s)versus (1000 K)/T for C1 HC1. The solid line is the multiple surface result for 1 = 0, while the dashed line is the corresponding 1 = 588 cm-l result. The
+
dotted line shows single surface rate constants from the 32 state basis, and circles show corresponding 80 state single surface results. Also shown (squares) are the results of several experimental studies that are tabulated in ref 23, and quasiclassical trajectory results (solid line at bottom) from ref 24. at higher energy. This leveling off arises because only one state is energetically accessible at the transition state. This is easy to understand for the single surface J = 0 c a l c ~ l a t i o nwhere ,~~ the ground vibrationally adiabatic state is nondegenerate (corresponding to zero quanta in the bend and symmetric stretch modes), but it is also true to the multiple surface calculation provided that just one parity component is considered. (The two parity components would give an overall probability of two, corresponding to the doubly degenerate 21:state.) At energies above 0.50 eV, the single and A = 0 multiple surface results separate from one another, with the deviation being close to a factor of 2 at energies above 0.55 eV. While this effect could in part be due to a lack of convergence in the calculations, it seems likely as we now show that the dominant effect responsible for this is excited bends present in the multisurface calculations. These excited bends correspond to Q N = f l in Table 2, and they are missing in the J = 0 single surface calculation for which QN = 0 is required. The Q N = A1 states are expected to become accessible at approximately 0.08 eV above the reactive threshold, which is consistent with the energy where the single surface and multiple surface results show deviations. We have decomposed our multiple surface results into contributions from different QN’s, and Figure 2 shows that the QN = QN.= 0 result for A = 0 (dotted curve) is in excellent agreement with the single surface probability. In addition, note in Table 2 that for N 2 1, the number of QN = zk 1 states having j = 3/2, Qj = (the combination of electronic quantum numbers that adiabatically correlates to the reactive 22surface) is the same as those having ZZN = 0. Thus at high energies one expects the multiple surface probability to be about twice the single surface probability as observed. Thus, it seems clear that the excited bends are the source of the difference in reaction probabilities. Differences between single and multiple surface reaction probabilities could also arise if reaction on the 211 surface could occur, but this should not be important in the energy range being considered as the barrier on that surface is 0.87 eV and would not be shifted enough by spin-orbit effects to be relevant. B. Rate Constants. Figure 3 shows an Arrhenids plot of the thermal rate constants that have been obtained from the reaction probabilities in Figure 2, using J-shifting to perform the partial wave sum. The multiple surface results in this plot
I
I 200
100
0
300 400 Spin-orbit Parameter A (cm
500
600
Figure 4. Thermal rate constant (log scale) versus spin-orbit parameter 1 (cm-I), with the overall multiple surface result indicated by the solid curve upon which are circles that denote specific energies where calculations were done. Also shown is the contribution to the thermal rate constant arising from individual fine structure transitions, 2P3/2 2P3/2(dashed line), ’P3/2 2P1i2and 2P1/2 *P3/2(solid line), and 2P1,2 2P1/2 (dotted line). All rate constants have been multiplied by exp(0.4311k7) before plotting to remove the “static” contribution to the dependence on the spin-orbit parameter.
-
-
-
-
are based on A = 0 and 588 cm-’, and they show the expected decrease in k(T) with increasing A. The single surface results are based on both the 32 state basis (dotted) and the 80 state basis (circles), and the results are indistinguishable (agree to within a few percent) on the scale of the plot. For reference, we also include experimental results in Figure 3 as summarized in ref 24 and quasiclassical trajectory rate constants from ref 25 for the same *I:LEPS surface. Figure 3 shows that the single and multiple surface results are in excellent agreement with each other, with differences that are less than 20% over the temperature range considered, and much smaller than the uncertainty in the experiments. This means that the transition state approximation in eq 17 for the electronic degrees of freedom is extremely accurate. This result is, of course, completely consistent with our discussion of the cumulative reaction probabilities, where we noted good agreement of multiple surface and single surface results for energies where only one state is energetically accessible at the transition state. Note that the difference at high energy due to bending in Figure 2 is at least partially included via the bend-shift approximation in Figure 3. Otherwise the multiple surface rate constants would have been larger than their single surface counterparts. C. Variation of Rate Constants with Spin-Orbit Parameter. To further study the fine structure dependence of the thermal rate constants, in Figure 4 we present multiple surface rate constants that we have calculated at 300 K as a function of the spin-orbit parameter A. Included in the plot are five different kinds of rate constants: the overall thermal rate constant (solid curve with circles), and then the contribution to this rate constant from the possible combinations of reagent and product fine structure states, namely, 3/2 3 / ~ ,3/2 Il2, ‘12 3/2, and ’12 ‘12. In calculating these rate constants, we have shifted the energy scale down by 0.43A to remove the static A dependence from the overall thermal rate constant (as explained in section II1.C). Figure 4 shows that the overall thermal rate constant (with the energy shift) is nearly independent of the spin-orbit parameter, showing only a 13% change as A varies between zero and 588 cm-’. However the spin state resolved rate constants show significant variation with A, especially for small
-
- - -
Influence of Atomic Fine Structure on Rate Constants
J. Phys. Chem., Vol. 99, No. 19, 1995 7529
A.
For A = 0, all the fine structure states are asymptotically degenerate, so electrostatic coupling causes strong mixing of j = I12 and 3/2. As a result, we find that, to a very good approximation, each spin state contributes statistically to the overall cumulative probability, which means that the 3/2 3/2 rate constant is 4/9ths of the total, 3/2 l/2 and its reverse are each 2/9ths,and l/2 I12 contributes ‘19th. As A increases, the rate constants associated with the excited spin statej = I12 make less and less contribution to the thermal rate constant, and eventually one reaches the “adiabatic limit”, where the excited spin state rate constants decay exponentially with increasing A (Le., there is a linear dependence of the curves in Figure 4 on A). In this large A limit, the thermal rate constant is completely dominated by the 3/2 3/2 transition, and it is this case that is found to occur when A has the physically appropriate value for Cl(2P).
-
-
-
-
V. Conclusion In this study we have used a very advanced scattering calculation to settle a very old question; namely, how accurate is the transition state theory treatment of electronic fine structure in determining thermal rate constants for chemical reactions? Our results for C1 HCl indicate that this approximation is extremely accurate, with errors that are less than 20% over the temperature range where we are able to determine accurate thermal rate constants. In addition we find that C1 HCl is very much in the adiabatic limit with respect to the correlation between asymptotic and saddle point electronic states, with the overall rate constant being completely dominated by reaction between the ground electronic states of the reagents and products. We find that the barrier on the excited 21-1 surface is sufficiently higher than the zZ barrier, that is has essentially no influence on the reactivity for energies that are important to the thermal rate constant. Further, the reactivity on the *Z surface is completely dominated by a single vibrationally adiabatic barrier near threshold in both the single and multiple surface dynamics, so there is no mechanism for making the single and multiple surface rate constants different for the low temperatures we have considered. However, the J = ‘12 multiple surface reaction probability matches its J = 0 single surface counterpart only at low energies due to the participation of excited bend states at higher energy. Indeed the higher energy dynamics that we did not study extensively will likely show even more important differences between single and multiple surface results. One of the most important conclusions of this work concerns the effect of spin-orbit splitting. To a very good approximation, we find that the only effect of spin-orbit on the cumulative reaction probability is to shift the reactive threshold up in energy with increasing A. This results in an increase in the activation energy and a corresponding decrease in the rate constant, which is quite significant for C1 HCl. We developed a simple way to predict the magnitude of this shift based on an analysis of Born-Oppenheimer potentials that include the spin-orbit Hamiltonian, and, in the case of C1 HC1, we find that this shift is about 43% of A (or 29% of the ’P3/2 - *P1/2 splitting) in magnitude. One consequence of this shift is that potential surfaces that have been determined without including spinorbit carry an intrinsic error that should be corrected to be truly realistic. The present scattering calculations can be used in many ways other than the determination of thermal rate constants, and we plan to study other properties, including spin state resolved cross sections and electronic and nuclear alignment propensities, in
+
+
+
+
future studies. In addition, it should be apparent from the small nuclear basis sets that we were forced to use to generate multiple surface results that it will be important to develop angular momentum decoupling approximations that apply to the mixed electronic-nuclear problems in future studies.
Acknowledgment. It is a pleasure to celebrate Mostafa El-Sayed’s 60th birthday through this issue and to wish him many future birthdays as J . Phys. Chem. Editor. This research was supported by the Office of Basic Energy Sciences, Division of Chemical Sciences, US. Department of Energy, under Contract No. W-31-109-ENG-38. Preparation costs associated with this manuscript were supported by NSF Grant No. CHE9016490. Computations were done at the National Magnetic Fusion Energy Computer Center and at the National Center for Superconducting Applications. Portions of this work were initiated at a workshop on “Electronically Nonadiabatic Processes in Chemical Dynamics, Spectroscopy and Electronic Structure” sponsored by Argonne National Laboratory, July 2024, 1994. I thank D. G. Truhlar for helpful discussions. References and Notes (1) Truhlar, D. G. J . Chem. Phys. 1972, 56, 3189. (2) Johnson, B. R.; Winter, N. W. J . Chem. Phys. 1977, 66, 4116. (3) For reviews, see: Nakamura, H. Int. Rev. Phys. Chem. 1991, 10, 123. Baer, M. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. II, p 219. (4) Tully, J. C. J . Chem. Phys. 1974, 60, 3042. (5) Lam, K. S.; George, T. F. In Semiclassical Methods in Molecular Scartering and Specrroscopy; Child, M. S., Ed.; Reidel: Dordrecht, The Netherlands, 1979; p 179. (6) Zimmermann, I. H.; George, T. F. Chem. Phys. 1975, 7, 323; J . Chem. Phys. 1975,61,2468. Zimmermann, I. H.; Baer, M.; George, T. F. J. Chem. Phys. 1979, 71,4132. Baer, M. In Molecular Collision Dynamics; Bowman, J. M., Ed.; Springer: Berlin, 1983; p 117. (7) Rebentrost, F.; Lester, W. A., Jr. J . Chem. Phys. 1975, 63, 3737; 1976, 64, 3879; 1976, 64, 4223; 1977, 67, 3367. (8) Schatz, G. C. Chem. Phys. Lett. 1988, 150, 92. (9) For a recent summary of theoretical studies of C1 HC1, see: Schatz, G. C.; Sokolovski, D.; Connor, J. N. L. In Advances in Molecular Vibrations and Collision Dynamics; Bowman, J. M., Ed.; JAI Press: Greenwich, CT, 1994; Vol. KB, p 1. (10) Schatz, G. C.; Sokolovski, D.; Connor, J. N. L. Can. J . Chem. 1994, 72, 903. (11) Bondi, D. K.; Connor, J. N. L.; Manz, J.; Romelt, J. Mol. Phys. 1983, 50, 467. (12) Bowman, J. M. Adv. Chem. Phys. 1985, 61, 115. (13) Schatz, G. C.; Kuppermann, A. J . Chem. Phys. 1976, 65, 4642. Schatz, G. C.; Hubbard, L. M.; Dardi, P. S.; Miller, W. H. J . Chem. Phys. 1984, 81, 231. (14) Reid, R. H. G. J. Phys. B: At. Mol. Phys. 1973, 6, 2019. (15) Ma, Z.; Liu, K.; Harding, L. B.; Komotos, M.; Schatz, G.C. J . Chem. Phys. 1994, 100, 8026. (16) Langhoff, S. R. J . Chem. Phys. 1980, 73, 2379. (17) Delves, L. M. Nucl. Phys. 1959, 9, 391; 1960, 20, 275. (18) Koizumi, H.; Schatz, G. C. In Advances in Molecular Vibrations and Collision Dynamics; Bowman, J. M., Ed., JAI Press: Greenwich, CT, 1990; Vol. lA, p 139. (19) Sun, Q.; Bowman, J. M.; Schatz, G. C.; Sharp, J. R.; Connor, J. N. L. J . Chem. Phys. 1990, 92, 1677. Colton, M. C.; Schatz, G. C. Int. J . Chem. Kinet. 1986, 18, 961. (20) Kuntz, P. J.; Nemeth, E. M.; Polanyi, J. C.; Rosner, S. D.; Young, C. E. J . Chem. Phys. 1966, 44, 1168. (21) Duggan, J. J.; Grice, R.J . Chem. Soc., Faraday Trans. 2 1983,80, 139. (22) Yamashita, K.; Morokuma, K. J . Chem. Phys. 1990,93,3716; 1991, 94, 831. (23) Chatfield, D. C.; Friedman, R. S.; Lynch, G.; Truhlar, D. G. J . Phys. Chem. 1992, 96, 57. (24) Garrett, B. C.; Truhlar, D. G.; Wagner, A. F.; Dunning, T. H., Jr. J . Chem. Phys. 1983, 78, 4400. (25) Kornweitz, H.; Broida, M.; Persky, A. J . Phys. Chem. 1987, 91, 5496.
+
JP9422046
