1507
J . Phys. Chem. 1984,88, 1507-1510 two reactions it is generally necessary to consider one radical species in addition to the major reactants CO, 02,and H 2 0 as effectively independent variables. (iii) Correlations of stable species with the intermediate radicals indicate the H concentration to be least correlated with the CO, 02,and H 2 0 concentrations. (iv) Large correlations between the various rate constants result when both k48 and k 1 5are considered among the dependent set of variables. (v) Correlations between species are generally largest at equilibrium and next largest during the induction period. In many instances, a minimum correlation was found during the kinetic period. (vi) To obtain a high correlation between [OH] or [HI and [CO], it is necessary to have another species containing hydrogen (such as H 2 0 ) in the independent set. (vii) According to second-order derived sensitivities, the correlation between reactions 48 and 15 is approximately linear, whereas the correlation between reactions 48 and 11 is highly nonlinear. Finally, a word about chemical intuition and sensitivity methodology should be included here. Chemical knowledge is certainly useful in the identification of some important aspects in a kinetic model. However, there are an enormous number of nonlinear effects that cannot be forseen prior to a sensitivity
calculation. Even for the simple mechanism considered here, the coupling and correlations that exist among the different variables are very complex and without the inspection of derived coefficients the magnitude and sense of these interdependences could not have been predicted. The application of chemical intuition is probably most useful in conjunction with a sensitivity study. In summary, derived sensitivities provide a valuable source of information for the elucidation of the interdependences of the variables of a kinetic mechanism. The kind of sensitivity analysis described above helps one to understand how the correlations among rate constants are filtered by the dynamics and ultimately reflected in the species concentrations predicted by the model under investigation. This estimation of system behavior is of particular importance when calculated species concentrations agree with experimental data, since in this case the sensitivity coefficients provide a tool for model evaluation.
Acknowledgment. We acknowledge the U S . Department of Energy and the Mobil Research and Development Corp. for support of this research. R.Y. acknowledges the support of a fellowship granted by the Exxon Educational Foundation. We also thank Dr. Britton Chang and Dr. Ken Brezinsky for helpful discussions. Registry No. Carbon monoxide, 630-08-0.
Correlation of Zero-Point Energy wlth Molecular Structure and Molecular Forces. 2. Approximation for n-Paraffins‘ Takao Oi and Takanobu Ishida* Department of Chemistry, State University of New York, Stony Brook, New York 1 1 794 (Received: June 17, 1983)
The approximation method for the zero-point energy (ZPE) developed previously on the basis of Lanczos’ T method has been applied to n-paraffins, C1-CI4, and to the sets of eigenvalues of n-paraffins in species A2, B2,A,, and B,. The best approximation to use for a given molecular system is predictable from the absolute error curve and distribution of eigenvalues. For the weighting function, w(X) 0: Xk, the best approximation for all n-paraffins tested is obtained by equating the range of expansion to the largest normalized eigenvalue and setting k = 0; the absolute value of the error is 2.5% for C, and falls below 1% for molecules larger than propane.
Introduction In our previous paper2 a second-order approximation method for the zero-point energy (ZPE) was developed on the basis of Lanczos’ ~ r n e t h o d . ~The approximation provides an explicit relation between the ZPE and the elements of F and G matrices through the traces of H (= FG) and H2in the following manner: h ZPE = -[bdX0,0 + b,(Xo,E)tr(H) + b&,F) tr(H2)l (1) 2 where bo = aoXo’/2E2/D
(2)
bl = alAo-1/2F/D
(3)
bz = a2X0-3/2/D (4) D = aoE2 + al< u2 (5) and, when the shifted Chebyshev polynomials of the first kind T n * ( x )are used for the orthogonal expansion, a, = -1, ul = -8,
+
and a2 = 8 / 3 . Xo and & are two variable parameters of the approximation which are determined as follows. Let the spectrum of vibrational eigenvalues, A, = 4?r2vZ2,of an N-atomic molecule be arranged in ascending order: 11,
0022-3654/84/2088-1507$01.50/0
Xi,
..*,A n
(=
Amax)
(6)
where n = 3N - 6 or 3 N - 5. Then, define a spectrum of normalized eigenvalues, x , % XJA,;
...)
..., xn
(= Xmax =
(7) where A, is a normalization reference point. E is then used as a range of x’s for the T-method expansion, which must satisfy the condition ~ 2 ,
XI,
Xmax/Q
Xmax 5 F (8) The best values of A, and $. were determined by minimizing the weighted sum, S, of the squares of the errors defined as
S = e w ( X , ) A I 2= e w ( X , ) [ b o+ bl A, I
(1) Research supported by the Office of Basic Energy Sciences, US. Department of Energy, Contract No. DE-AC02-80ER10612. (2) Oi, T.; Ishida, T. J . Phy. Chem. 1983,87, 1067. (3) Lanczos, C. “Applied Analysis”; Prentice-Hall: Englewod Cliffs, NJ, 1956.
...t
A29
i
+ b2X: - X,1/2]2
(9)
with respect to variations in Xo and t . The weighting function used was w(X,)= A i l k 0 1984 American Chemical Society
(10)
1508 The Journal of Physical Chemistry, Vol. 88, No. 8, 1984
Oi and Ishida
TABLE I: Best Approximations for n-Paraffins
approx A
B C
complete spectra parameters a1 ( k = 0 ) precision -3.1% (C,) to
out-of-planespectra a3 ( k = 0) -1.3% (C,) to -0.3% (Clo) -0.3% (C14)
parameters
b l ( k = O)a
b l ( k = O)a
precision
-2.510 (C,) to -0.03% (C,) c l ( k = 0.173)b
-3.1% (C,) to
parameters precision
t0.1% (C,) to t5.8% (C,,,)
-2.7% (C14) c l (k = 0.173)b t0.6% (C,) to 1.4%( C 1 4 )
a The generally most reliable approxiniation within the framework of the present weighting function. See text. k = 0.173 is the smallest value of k that yields at least one positive root for the optimization equations, aS/a A,, = ash g = 0. See ref 2.
and the optimum values of Xo and [ were obtained for various values of the weighting parameter k., Three cases were considered: case A, for which a constraint X,= ,A, was set; case B, another special case in which the range and case C in which 6 and of expansion 6 was equated to xmax; Q(= Xman/XO) were treated as independent variables. The optimization of S yielded the following sets of solutions. Case A: three positive roots for for every k value tested, ranging from k = 0 to lo6. Each of three roots, {(al), E(a2), and f(a3), leads to a distinct approximation. Case B: three k-dependent positive roots for f (and the corresponding values of Xo), which are designated as [(bl), E(b2), and f(b3), and two k-independent roots, [(b4) and t(b5). E(bl), 5(b2), and ((b3) yield an identical approximation, while [(b4) and l(b5) each give distinct approximations. Case C: two positive roots, [E(cl), ~ ( c l ) and ] [((c2), Q ( c ~ ) ] .They produce an identical approximation. In this paper we present an application of these approximations to the Z P E s of n-paraffins. Only six approximations, Le., those generated by the solutions a l , a2, a3, b l , b4, and c l , have been considered. The b5 approximation has been omitted, because it gave consistently and significantly poor results for the molecules we tested previously.2
Calculation Numerical approximations for the Z P E s of n-paraffins have been obtained by using six sets of approximation parameters and known eigenvalues of n-paraffins, and the results have been compared with the exact ZPE calculated from the same sets of eigenvalues. The vibrational eigenvalues used for C, (ethane) through C14 are those computed from the F matrices (VFF-IV) of Schachtschneider and Snydere4 In order to facilitate examination of the effect of changes in the frequency distributions, we also calculated the sum of all eigenvalues for the species A2, B,, A,,, and B, (VFF-II)4 except those for the torsional modes. For simplicity we will designate the sets of these (Az, B2, A,, BJ eigenvalues (but excluding the torsions) as the out-of-plane (OP) spectra and the complete sets of eigenvalues as the complete spectra (CS). We have included the results on methane which we previously reportedZfor comparison. The methane data were obtained from the F matrix of Jones and McDowell.’ Results and Discussion Typical results of the ZPE approximations expressed as relative errors have been deposited in the supplementary material, and the best approximations have been summarized in Table I. Among the case A approximations for the CS systems, a1 yields a precision which is least sensitive to variations in k, while all case A approximations for the O P are relatively insensitive to the changes in k. The insensitivity is a desirable property, since it (4) Schachtschneider,J. H.; Snyder, R. G. Spectrochim. Acta 1963, 19, 117. ( 5 ) Jones, L. H.; McDowell, R. S. J . Mol. Spectrosc. 1959, 3, 632.
5=
0.2
-0.1
t
k =
i
0.0
’
0.0
0.2
0.4
0.6
0.8
1.0
x ( m x/x,, Figure 1. Absolute error and distribution of eigenvalues A, as functions of x for case A (X,= Amx) at k = 0. The unit interval of the distribution histogram is Ax = 0.01, and the heights are normalized to unity for the
most populous interval. allows a wider margin of error in choosing a proper k value a priori. For the C S system a3 is most dependent on k . The k dependence of a2 is intermediate between those of a1 and a3. The best k values for the a2 approximation range from 18 (C,) to 0.5 (Cs-C,o) for the CS, while the range of the variation is relatively small for the O P spectra; -20 for C4 and 12 for C12-C14. For both the complete spectra and the O P spectra the best among the case B approximations is that of b l with k = 0 and the best of case C approximations is that of c l with k = 0. To facilitate an explanation for the differences in the approximations between C S and OP spectra, Figures 1-4 have been prepared. In the lower half of each figure the absolute error calculated as A(x) = bo + blx bzx2 - x1I2 (1 1)
+
has been plotted as a function of x (= A/X,). The range of x varies depending on whether the type of approximation is A, B, or C. In the upper half of each figure the distribution of eigenvalues of a few representative n-paraffins has been presented as histograms, with intervals of 0.01 in the unit of x. Each member of degenerate sets has been counted separately. The most populated interval is shown as a vertical line going across the full width of the spectrum box, and the heights of other lines have been drawn in proportion to the population relative to that of the most populous interval. It is seen that the eigenvalues are grouped into two distinct classes, Le., one of the large, C-H stretching, eigenvalues and the other of the smaller eigenvalues consisting of C-C stretching, various angle bending, wagging, rocking, and torsional modes. For paraffins there are no eigenvalues in the rather wide middle region. The similar situation holds with a large majority of compounds, and only those having the multiple bonds such as the carbonyl group have some population in the midrange. Consequently, a
The Journal of Physical Chemistry, Vol. 88, No. 8, 1984 1509
Approximation Method for Z P E s of n-Paraffins
C7(CS) [f =63]
c14 (OP) [f =42]
0.2 p:
2
E
0.1
-
I
!
( ( b l ) = 0.867171 k = 0.0
,~ I € ( c l ) = 1.438375
I I I I
-
s=
I
W
0.0
0.2
0.4
0.6
o.8T X,,
1.0 = 0.867
x ( E X/X,, Figure 2. Absolute error and distribution of eigenvalues Xi as functions of x for case B ([ = Xmsx/XO) at k = 0.
I 1
OL
2
E
((bl)
=
1.033384
7 ( c l ) = 1.443468 = 0.173
a
I l 1
-
0.1
L
-
4 -0.1
j 0.0
j 0.2
x
0.4
X/X,)
0.6
0.8
1.0
1.2
1.41
x m o x =1.443 I Figure 4. Absolute error and distribution of eigenvalues A, as functions of x for case C at k = 0.173. large error in an approximation in the midrange would not affect the precision. The lack of population in the low end of the O P spectra is due to the artificial exclusion of the torsional modes. From Figure 1 we see the reason why a1 ( k = 0) is the best approximation for the CS, while a3 ( k = 0) is best for the OP spectra. For the CS, the error curve of a1 ( k = 0) is such that the negative errors in approximating the ZPE contributions of the small eigenvalues (other than torsion) and the positive errors contributed by the torsional modes tend to cancel each other. For the O P spectra which lack the torsional modes, the cancellation of errors occurs between the lower and upper halves of the group of “small” eigenvalues, which is best achieved by a3 ( k = 0). In both a1 and a3 the C-H stretching eigenvalues do not significantly contribute to the total errors for the CS and OP spectra. In the b l ( k = 0) approximation (cf. Figure 2), the “large” eigenvalues and the torsional modes contribute positive errors, while all other eigenvalues yield negative errors. This leads to a more complete cancellation between the positive and negative errors in the CS than in the OP spectra, resulting in the generally better b l ( k = 0) approximation for the CS than for the OP spectra (cf. Table I). Figure 3 illustrates a consequence of giving an excessive weight to the “large” eigenvalues. It results in practically zero error for the C-H stretching modes and depletion of this source of positive errors which would otherwise have counteracted the negative errors due to the “small” eigenvalues. Figure 4 shows that the c l ( k = 0.173) approximation yields large positive error at the high end of the spectra, resulting in the overall positive errors for all molecules tested. Since the torsional modes contribute positive errors as usual, this method leads to the results which are better for the OP spectra than for the CS systems. The comparison of the complete and OP spectra has illustrated the importance of the contributions of small eigenvalues to the overall precision. This is a direct consequence of the present form (=
1510
J. Phys. Chem. 1984,88, 1510-1513
of the weighting function, eq 10. The error is at its highest at x = 0 due to the little weight given to this end when k I0.
give better results for a particular spectrum are predictable on the basis of the absolute error curves.
Conclusion Within the framework of the present weighting function the generally most reliable results are obtained for n-paraffins by the b l ( k =-0) approximation; for the complete spectra the erior is -2.5% for C1 and -1.5% for C3 and falls below 1% beyond C3 while, for the out-of-plane spectra, the error ranges from -3.1% for C3 to -2.7% for CI4. Other types of approximations which
Acknowledgment. We gratefully acknowledge use of computer time at the Computer Center of the State University of New York at Stony Brook, Supplementary Material Available: Two tables of typical results of the approximations for the ZPEs of n-paraffins (2 pages). Ordering information is given on any current masthead page.
Potential Surface Walking and Reaction Paths for C P vBe Be 4- 2H (’A’)
+ H2
- BeH,
Douglas O’Neal, Hugh Taylor, and Jack Simons* Department of Chemistry, University of Utah, Salt Lake City, Utah 84112 (Received: June 20, 1983)
- -
By combining the surface walking algorithm of Simons et al. with locally determined forces on an ab initio surface, the reaction paths for the model reactions Be + H2 BeH, Be + 2H (‘A,) were studied. This represents the first application of this algorithm to an ab initio surface which is generated locally as the walking proceeds.
Introduction The development of ab initio molecular gradient techniques’+ has made it much more feasible to examine theoretically the energetics of chemical reactions. In fact, it may even now be practical to use these methods within dynamical studies of molecular reactivity. Analytical expressions for molecular gradients (forces) and Hessians (force constant matrices) allow one to avoid tedious and inefficient finite difference calculations and provide efficient tools to use in locating both stable molecular geometries and reaction transition states. Stable molecular geometries have gradient components along all molecular distortions equal to zero and all Hessian eigenvalues (local-normal-mode frequencies) positive. Transition states have vanishing gradients and Hessians with only one negative eigenvalue. When combined with an efficient surface walking procedure,’,* ab initio analytical gradients can be used to proceed, on a single potential energy surface, from a reactant geometry, along various deformation directions, to transition-state geometries and onward to the product molecules. Simons et al.’ have recently developed and implemented, within the context of globally known two- and three-dimensional model potential energy surfaces, such an algorithm which we now utilize for the first time in an ab initio study
where the surface is not globally known but is actually generated as the walk proceeds. The model reactions
H2
+
Be
-
I
Be
+
2H
occur on the same (‘Al) potential energy surface in C, symmetry. The two product states (Be 2H and BeH2) arise through a bifurcation of the Be + H2 reaction path (defined below) which occurs after the transition state has been surmounted. In the present work the reaction Be H 2 BeH H is eliminated from consideration since the molecular deformations are (arbitrarily and artificially) restricted to retain C, symmetry. For this reason, we must consider our reaction study to be of a modelistic nature. It does, however, represent the first ab initio application of the surface walking algorithm of ref 7 to a locally generated energy surface which contains branching. The local description of the potential energy surface generated via our walk-driven ab initio quantum-chemical calculations also contains information about the forces and curvatures transverse to the direction of the reaction path. In a future work, we plan to make use of this information to carry out a dynamical study of the above reactions of Be H2 utilizing generalizations of the reaction path Hamiltonian method of Adams, Handy, and Millers9 We are particularly interested in developing an improvement of the reaction path dynamics method of ref 9 to permit the incorporation of channel branching.
+
+
-
+
+
(1) P. Pulay in “Modern Theoretical Chemistry”, H. F. Schaefer 111, Ed., Plenum Press, New York, 1977, Chapter 4. (2) J. D. Goddard, N. C. Handy, and H. F. Schaefer 111, J. Chem. Phys., 71, 1525 (1979); Y. Osamura, Y. Yamaguchi, and H. F. Schaefer 111, ibid., 75,2919 (1981); 77, 383 (1982); B. R. Brooks, W. D. Laidig, P. Saxe, J. D. Goddard, Y. Yamaguchi, and H. F. Schaefer 111, ibid., 72, 4652 (1980). (3) J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Brinkley, Int. J . Quantum Chem., 13, 225 (1979); R. Krishnan, H . B. Schlegel, and J. A. Pople, J . Chem. Phys., 72, 4654 (1980). (4) H. F. King and M. Dupuis, J. Comput. Phys., 21, 1 (1976); M. Dupuis, J. Rys, and H. F. King, J . Chem. Phys., 65, 111 (1976); M. Dupuis and H. F. King, ibid., 68, 3998 (1978). (5) A. Komornicki, K. Ishida, K. Morokuma, R. Ditchfield, and M. Conrad, Chem. Phys. Lett., 45,595 (1977); K. Ishida, K. Morokuma, and A. Komormichi, J . Chem. Phys., 66, 2153 (1977). (6) P. Jmrgensen and J. Simons, J . Chem. Phys., in press. (7) J. Simons, P. Jmrgensen, H. Taylor, and J. Ozment, J. Phys. Chem.,
The Surface Walking Algorithm To efficiently utilize local force and curvature information, the dependence of the electronic energy E on each of the molecule’s n internal distortion coordinates {Xi;i = 1, 2, ..., n) = X is Taylor series expanded’ about the geometryI0 X = 0
87, 2745, (1983). (8) C. J. Cerjan and W. M. Miller, J . Chem. Phys., 75, 2800 (1981).
(9) J. E. Adams, N. C. Handy, and W. H. Miller, J. Chem. Phys., 72.99 (1980).
0022-3654/84/2088-15 10$01.50/0
0 1984 American Chemical Society
