344
J . Phys. Chem. 1986,90, 344-356
J. R. Barker, ”Large Molecular Energy Transfer: New Techniques and New Controversies” J. C. Stephenson, “Vibrational Relaxation of Chemical Bonds in Liquids and on Surfaces” J. I. Brauman, “Prediction of Rate Constants for Ionic Reactions in the Gas Phase and in Solution” J. A. Howard, “Measurement of Absolute Propagation and Termination Rate Constants for Alkylperoxyls in Solution by the Hydroperoxide Method” J. H. Expenson, “Kinetics and Mechanisms of Inorganic Reactions in Aqueous Solutions: Intermediates” T. H. Dunning, Jr., “Theoretical Characterization and Calculation of Potential Energy Surfaces for Chemical Reactions” J. Wolfrum, “Laser Stimulation and Observation of Bimolecular Reactions” Professor Frederick Kaufman was scheduled to give a talk entitled “Kinetics of Bimolecular Processes”. H e was unable to attend the meeting due to sudden illness. Tragically, Professor Kaufman died shortly after the conference ended. His loss removes not only one of the leading figures in gas kinetics in the world, but also a man deeply admired for his wit, dedication, and humanity. At some future date, a special issue of The Journal of Physical Chemistry will be dedicated to Professor Kaufman. The titles of the talks give a clear indication of the scope of the conference. The hour-long discussion periods concluding each session were especially valuable. Under the direction of the chairman, and with the invited lecturers providing guidance, the discussions provided an in-depth look at the current state of each area. A special feature of this meeting was the bringing together of kineticists from many disparate areas of chemistry. In all these
areas modern techniques allow us to focus on the fundamental interactions. At this level, the commonality of the effects are clear and we can all learn a great deal from each other. The new methodologies and broadened theoretical understanding have led to two distinct directions in experimental work: the study of truly fundamental processes in the form of state-testate chemistry, and the study of reactions involving molecules of increasing complexity and the measurement of rate constants with even greater accuracy. To be meaningful, the goal of all such studies is the verification and expansion of the theoretical base of understanding of reactive processes. In the case of state-to-state studies, the experiments offer a specific challenge to theory in that the observables, such as product state distributions, are extremely sensitive to the potential energy surface itself and the dynamics on such surfaces. Kinetic studies under thermal conditions, designed largely to ascertain such system characteristics as product identity and the absolute value of the rate constants as functions of temperature and pressure, are less sensitive to the details of the potential energy surface and thus are amenable to semiempirical approaches based on transition-state theory. Benson’s “Thermochemical Kinetics” has proven to be an extremely useful tool for explaining some observations and predicting others. Thus, theory may play an important role in this area, since current and developing theoretical methods may be more than adequate for finding the general structures and energetics of transition-state species required to describe thermal elementary steps. Ultimately, one expects that the concordance of key state-tostate experiments and the predictions from the best ab-initio theories will provide guidelines for the proper application of the semiempirical methods mentioned above. This will be of immense importance to our understanding of chemistry as a whole.
Theoretical Studies of the Energetics and Mechanisms of Chemical Reactions: Abstraction Reactions Thorn. H. Dunning, Jr.,* Lawrence B. Harding, Raymond A. Bair, Robert A. Eades, and Ron L. Shepard Theoretical Chemistry Group, Chemistry Division, Argonne National Laboratory, Argonne, Illinois 60439 (Received: July 26, 1985)
In this article we first review methods developed to locate and characterize critical points (minima and saddle points) on molecular potential energy surfaces. We then briefly discuss the computational advances which are currently driving the development of algorithms for calculating potential energy surfaces. Finally, we present a sampling of the studies of simple abstraction reactions carried out in the Theoretical Chemistry Group at ANL over the past five years. In this latter section we focus on a series of prototypical chemical reactions (0 + H2, H + HX, CH + H,, H + HCO, OH + H2, H + HCN, and H2 + C2H). Most of the above reactions involve direct processes; however, for the CH + H2 and H + HCO reactions indirect, addition-elimination pathways also play a role. For the homologous series of abstractions, H + HX (X = F-I), simple exponential relationships are found both between the barrier heights and the reaction exoergicities and between the two bond extensions at the saddle point. Finally, the last two reactions, H2 + CN and H, + C2H, are both highly exoergic and are predicted to have very early transition states. One consequence of this is the presence of unusually low, saddle point, vibrational frequencies leading to pronounced curvature in the Arrhenius plots for the rate constants of these two reactions.
Introduction Many chemical processes, including the oxidation of hydrocarbons, the chemistry of the atmosphere, and laser action in chemically reacting systems, are the result of complex chain reactions involving highly reactive species whose chemistry is difficult to characterize in the laboratory. Recent advances in the methodology of quantum chemistry, combined with continuing advances in computer technology, have opened up the possibility of determining the chemistry of these speciesfromfirst principles. Such studies provide a clear and compelling opportunity for a significant advance in our understanding of one of the most basic of chemical processes-the elementary chemical reaction.
The determination of the mechanism and calculation of the energetics and rates of chemical reactions requires, first of all, the potential energy surface for the composite molecular system of interest. The structure and energy of the transition state for a chemical reaction is of obvious importance. With this information, and corresponding information on the reactants, transition-state theory methods’ can be used to calculate the overall rate of the reaction. (1) P. Pechukas in “Dynamics of Molecular Collisions”, Part B, W. H. Miller, Ed., Plenum Press, New York, 1976, p 269; D. G. Truhlar, W . L. Hase, and J. T. Hynes, J . Phys. Chem., 87, 2664 (1983).
0022-3654/86/2090-0344$01 SO10 0 1986 American Chemical Society
Energetics and Mechanisms of Chemical Reactions
The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 345
To calculate the potential energy surface for a chemical reaction we must solve the electronic Schroedinger equation
He(R;r) $e(R;r) =
$e(R;r)
(1)
where (R} represents the coordinates of all of the nuclei and (r) refers to the electronic coordinates. The effective potential obtained from the solution of (1) and the mutual electrostatic interaction of the nuclei
defines the potential energy surface upon which the nuclei move. For bound species and transition states, solution of the nuclear Schroedinger equation
11. Characterization of Potential Energy Surfaces for Chemical Reactions: Location of Stationary Points Of the 3N nuclear coordinates in the nuclear Schroedinger equation, (3), three describe translation of the system as a whole while another three describe overall rotation of the system. Thus, the molecular potential energy surface, V(R), which describes the interaction energy of the atoms, is a function of only 3N - 6 nuclear coordinates. These coordinates, referred to as internal nuclear coordinates, can be chosen in a number of ways, e.g., for a triatomic system we might choose the distances between the three or the distances between two of the atoms atoms, (RI2,Rz3,Rl3), and,the included angle, (R12,R23,8123). The choice of which coordinate system to use is often just a matter of convenience. Specification of the molecular potential energy surface as a function of 3N - 6 internal coordinates
(3) determines the vibrational and rotational energy levels of the system. For transition states, one of the vibrational frequencies, that corresponding to the reaction coordinate, is imaginary. Equation 3 can be solved by using a variety of techniques, the simpliest of which is the method of G w i d (see also the FG method3). With this data, and information on the structures and electronic energies, reaction rates can be calculated from transition-state theory.’ While the methods used to solve the electronic Schroedinger equation have changed little in the past few years, great strides have been taken in improving the efficiency of the computational techniques and algorithms used. To date, only for the simpliest reaction, H H2, has a potential energy surface of chemical accuracy, Le., with errors less than 1 kcal/mol, been r e p ~ r t e d , ~ although the energetics of many, more complex chemical reactions have been calculated to an accuracy of a few kcal/mol. These latter calculations, while lacking absolute accuracy, have nonetheless provided a wealth of valuable detail on the energetics and mechanisms of the reactions studied. In the next few years, application of the new techniques and algorithms developed for solving the electronic Schroedinger equation over the past decade can be expected to greatly expand the number of chemical reactions for which quantitative or near-quantitative potential energy surfaces are known. The potential energy surface, (3), is a function of ( 3 N - 6) variables. Its characterization presents a challenge as formidable as that presented by the calculation of the energy itself. During the past few years, techniques for characterizing molecular potential energy surfaces have been rapidly evolving. Thus, reliable methods are now available for the location of minima (bound species) and saddle points (transition states).s These techniques provide a systematic means of characterizing the critical regions of molecular potential energy surfaces and their exploitation can be expected to greatly increase our understanding of the mechanisms, energetics, and dynamics of chemical reactions. In this paper we will first briefly discuss the methods developed to characterize minima and saddle points on molecular potential energy surfaces, comment on the computational advances which are currently driving the development of quantum chemical algorithms, summarize the technique used to calculate vibrational energies, and, finally, present a sampling of the studies of simple abstraction reactions carried out in the Theoretical Chemistry Group at Argonne National Laboratory over the past five years. In the latter section, we will focus on selected features of prototypical chemical reactions; in this way it is hoped that the conclusions drawn will extend beyond the specific systems considered.
+
(2) W. D. Gwinn, J . Chem. Phys., 55, 477 (1971). (3) E. B. Wilson, Jr., J. C. Decius, and P. C. Cross, “Molecular Vibrations,
The Theory of Infrared and Raman Vibrational Spectra’’, McGraw-Hill, New
York, 1955. (4) B. Liu, J . Chem. Phys., 58, 1925 (1973); P. Siegbahn and B. Liu, J . Chem. Phys., 68, 2457 (1978); M. R. A. Blomberg and B. Liu, J . Chem. Phys., 82, 1050 (1985). (5) For a review of these methods, see S. Bell and J. S. Crighton, J . Chem. Phys., 80, 2464 (1984).
presents a formidable practical problem. For example, even for a triatomic system the potential surface is a function of three variables and thus cannot be visually displayed. For a four-atom system it is a function of six (6) variables and a brute-force characterization of the surface is no longer possible, e.g., if only ten (10) points were required to represent the dependence of the potential on one of the coordinates then characterization of the entire surface would require lo6 calculations, an intractable task at present. Fortunately, not all regions of the potential energy surface are of equal importance. For chemical reactions the regions of greatest importance are the minima (representing the reactants, products, and any bound intermediates) and saddle points (representing transition states). Both minima and saddle points are stationary points on the potential energy surface and thus the gradient, Le., the first derivative of the potential energy with respect to the internal coordinates, vanishes at these points qi = av/aRi =
o
(i = 1, ..., 3~ - 6 )
(4)
For minima, which correspond to ualleys on the potential surface, the eigenvalues of the matrix of second derivatives, the curvature or Hessian matrix
k , = aZv/aRi aRj must be all positive, Le., the potential energy surface must be concave upward in all directions. For saddle points, which correspond to mountain passes on the potential energy surface, there can be one and only one negative eigenvalue of the Hessian matrix. The eigenvector corresponding to the negative eigenvalue specifies the direction of motion over the mountain pass, Le., it is the reaction coordinate at the saddle point. A number of methods have been developed in recent years to locate stationary points on potential energy surfaces. We will present a brief sketch of one such technique here, the NewtonRaphson methods6-’ This discussion should only be considered illustrative of the modern approaches that have been developed to locate stationary points on molecular potential energy surfaces; considerable work in this area continues. For a discussion of alternative techniques for locating stationary points, the reader is refered to the paper by Bell and C r i g h t ~ n . ~ Most modern techniques for locating stationary points on molecular potential energy surfaces are based on a quadratic representation of the surface in the vicinity of the stationary point of interest. Thus, if [R,] is an approximation to the desired stationary point, the energy in this region is written as E(R) = E,
+ cq,”AR, + cCkl,”AR, AR, I
(6)
I>J
with the gradients and Hessian elements evaluated at the point (6) P. Pualy, Mol. Phys., 17, 197 (1969). (7) J. W. McIver and A. Komornicki, Chem. Phys. Lett., 10, 303 (1971). (8) D. Poppinger, Chem. Phys. Lett., 35, 550 (1975).
346 The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 [R,,]. At the stationary point, [R,], the gradient must vanish and therefore from (6) R, = Ro - ko-lqo
(7)
In general, the gradient at the point predicted by (7) is not identically zero since the potential energy surface is not truly quadratic over the region of interest. If this is the case, the Hessian is recomputed at the new point and the process is repeated until the calculated gradient (or the predicted changes in the stationary point coordinates) is below some preset threshold. In our work we compute (6) from a grid of points; this requires the solution of (1) at least ( 3 N - 4)(3N - 5 ) / 2 times for each iteration in the search. Gradients can now be computed analytically for many types of electronic wave functions and, recently, analytical gradients have become available for the correlated wave functions needed for the study of chemical reactions.”* However, the Hessian matrix still cannot be efficiently computed for such correlated wave functions. To obtain information on the curvature of the potential energy surface, the Hessian matrix can be computed by finite differencing of the gradients. In fact, in most implementations of the Newton-Raphson method5 the Hessian matrix elements are simply estimated from the changes in the gradient in the iterative process outlined above. While the Newton-Raphson method converges in only one iteration if the surface is quadratic, instabilities are often encountered in molecular applications, Le., the method overshoots the stationary point, first in one direction and then another. This problem can be addressed by using line searches to adjust the length of the correction vector predicted by (7). In this procedure the energy or the square of the gradient norm, 1qI2,is calculated at selected points along the correction vector, -ko-]q0. The corrected stationary point is then taken to be that point at which these quantities are a minimum (energy minimization, of course, only works for locating minima). While the Newton-Raphson method and its many variants work well for the location of minima, difficulties can be encountered in the location of saddle points. For example, to locate a saddle point with this technique, the starting guess must be in the region of the potential energy surface where there is one and only one negative eigenvalue; otherwise the procedure diverges or tends toward the nearest stationary point. Further, line search techniques are of limited utility as the correction vector may lie along a line of vanishing curvature leading to one of the asymptotes. Recently, new techniques have been proposed which first isolate the direction of negative curvature and then minimize the energy in the remaining degrees of f r e e d ~ m . ~ These ~ ’ ~ . techniques ~~ which are based on conjugate gradient algorithms show considerable promise. Other methods are based on alternative expansions of the Further improvements in methods for locating saddle points are expected as work in this area continues. 111. Calculation of Potential Energy Surfaces for Chemical
Reactions
During the past few years there have been developments in both computer hardware and quantum chemistry software that have had a significant influence on the way theoretical chemists approach and solve problems. Although methods such as HF (Hartree-Fock), GVB (generalized valence bond), MCSCF (multiconfiguration self-consistent-field), and CI (configuration (9) M . Dupuis, J . Chem. Phys., 74, 5758 (1981). ( I O ) Y. Osamura. Y. Yamaguchi, and H. F. Schaefer 111, J . Chem. Phys., 77, 383 (1982). ( I I ) P. Jorgensen and J. Simons, J . Chem. Phys., 79, 334 (1983). (12) M. Page, P. Saxe, G. F. Adams, and B. H. Lengsfield 111, J . Chem. Phys., 81, 434 (1984). (13) S . Bell. .I. S. Crighton, and R . Fletcher, Chem. Phys. Left., 82, 122 (1981). (14) H. B. Schlegel, J . Compur. Chem., 3, 214 (1982). (15) C. J. Cerjan and W. H. Miller, J. Chem. Phys., 75, 2800 (1981). (16) A. Banerjee, N. Adams, J. Simons, and R. Shepard, J. Phys. Chem.. 89. 5 2 (1985).
Dunning et al. interaction) are employed now as they were years ago, there has evolved a qualitatively different approach to their use. We will discuss these trends here; for a discussion of the methods themselves the reader is referred to recent reviews on the GVB method,17 the MCSCF method,l* and the CI method.19 A . Computer Hardware Trends and Quantum Chemistry Computations. Until the later 1970’s, most theoretical chemists depended on access to a central computing facility for all of their computing needs. A typical central facility had to satisfy many different user communities, each of which had their special requirements and placed their own special demands on the facility. This inevitably lead to situations where the users of these central facilities had to use nonoptimal configurations to solve their computational problems. In the late 1970’s, however, there were two important disruptions of this status quo. First, the development of the supercomputer with processing speeds of 10 to 50 times the capacity of the older central facilities started a trend toward new kinds of central facilities, ones that were dedicated toward more specific user communities. Unfortunately, despite an abortive attempt at centralization in the late 1970’s (the short-lived National Resource for Computation in Chemistry), computational chemistry was not one of these select user communities, and it is only recently that chemists in this country have obtained widespread access to supercomputers. An important feature of supercomputers is that their maximum computation rates are achieved for only a few types of operations, typically vector operations such as dot products, vector-vector multiplies, and/or vector-vector adds. Another important feature of supercomputers is that the capacity of the CPU to perform floating point operations far exceeds the capacity of external memory devices such as disk drives to supply data to the CPU. This means that, for a balance to be achieved between CPU and 1 / 0 rates, each element transferred from the external device must be used in many floating point operations. Traditional algorithms therefore had to be modified, sometimes extensively, before codes showed significant performance improvements on these supercomputers.20 Second, super-minicomputers, such as the Digital Equipment Corporation VAX 11/780 introduced in 1977, became available that had adequate capacity to satisfy the computational needs of small groups of users. Further, memory prices for these superminicomputers were such that typical super-minicomputer configurations had ten times the amount of memory of the older central facilities. These minicomputers had slower CPU speeds than those in the central facilities but the fact that users had dedicated access to the machine more than compensated for this deficiency. Subsequent hardware improvements in CPU capacity have resulted in super-minicomputer and super-minicomputer/ attached processor systems that possess approximately the same capability as a large central facility of a decade ago. An important feature of the new attached processor systems, such as the Floating Point Systems, Inc. Model 164 scientific computer (FPS-164) used by our group, is that, like supercomputers, they also perform at their maximal rates for only a few types of operations; in addition, their CPU capacity also often outstrips their I/O capacity. It is ready access to these large-memory dedicated super-minicomputer and, more recently, super-minicomputer/attached processor systems that has, up to now, had the most significant impact on (17) F. W. Bobrowicz and W. A. Goddard 111 in “Methods of Electronic Structure Theory”, H. F. Schaefer 111, Ed., Plenum Press, New York, 1977, p 79. (18) A concise discussion of modern MCSCF techniques is not presently available. For a presentation of the relevant points, see H.-J. Werner and W . Meyer, J . Chem. Phys., 74, 5794 (1981); R . Shepard, I. Shavitt, and J . Somins, J . Chem. Phys., 76, 543 (1982); B. H. Lengsfield, 111, J . Chem. Phys., 77, 4073 (1982); J. Olsen, D. L. Yeager, and P. Jorgensen, Adcan. Chem. Phys., 54, 1 (1983); F. B. Brown, I . Shavitt, and R. Shepard, Chem. Phys. Lett., 104,363 (1984); P. E. M. Siegbahn, Chem. Phys. Lett., 109,417 (1984); P. J. Knowles and N. C. Handy, Chem. Phys. Lett., 111, 315 (1984). (19) 1. Shavitt in “Methods of Electronic Structure Theory”, H. F. Schaefer 111, Ed., Plenum Press, New York, 1977, p 189. (20) M. F. Guest and S. Wilson in “Supercomputers in Chemistry”, P. Lykos and I. Shavitt, Eds., American Chemical Society, Washington, D.C., 1981, ACS Symposium Series N o . 173, p I
Energetics and Mechanisms of Chemical Reactions both the solution of chemical problems and on the development of software for quantum chemistry computations. B. Computer Software Development for Quantum Chemistry Computations. Although the electronic structure methods used in our group (HF, GVB, MCSCF, and CI) have been used in various forms since the late 196O’s, the computational algorithms have changed significantly, particularly in the last few years. This has partially been in response to the access to dedicated computers by an increasing number of computational chemists and also from a natural evolution of the methodology and formalism. Central computer charging algorithms of a decade ago were prejudiced against algorithms requiring large memory and consequently most programs were designed to use a minimal amount of memory, even at the expense of increased 1/0 or in some cases increased floating point operations. In contrast, the cost of computing on dedicated computers is most appropriately measured by the elapsed time required for a given calculation. Minimization of elapsed times involves reduction of 1 / 0 times as well as reduction of CPU times, even at the expense of larger memory requirements. Many computational chemistry algorithms also involve the extensive use of iterative procedures. It has often proven useful to replace iterative procedures with noniterative procedures or with alternate iterative procedures with improved convergence characteristics in order to reduce the total 1/0 requirements (and thus the elapsed times). The following are a few examples of the changes in methodology that we have employed in our research group in order to maximize the performance of our codes on our FPS-164 and on the CRAY X-MP in the National Magnetic Fusion Energy Computer Center (Livermore, CA). The GVB method]’ is a contracted MCSCF method where the C S F (configuration state function) coefficients are nonlinear functions of a small number of variational parameters, typically one parameter for each correlated electron pair. The optimization of the orbital expansion coefficients and nonlinear variational parameters requires the construction of several Fock matrices, each of which may be regarded as a matrix-vector product where the matrix is an ordered list of two-electron integrals and the vector is an ordered set of effective one-particle density matrix elements. Older GVB programs typically required a separate reading of the integral list for each Fock matrix, thereby minimizing memory requirements. Newer programs not only construct all of the Fock matrices simultaneously with a single integral read, but also allow for the complete elimination of 1/0 in those cases where the integrals may be stored entirely in memory. A general MCSCF wave function is an arbitrary linear combination of CSFs where the orbitals are optimized along with the CSF expansion coefficients. The earliest molecular MCSCF wave functions were computed in a way similar to the GVB wave functions just discussed. The problems with this approach were that the effort involved increased very rapidly with the number of CSFs in the expansion and that solution of the equations was often only slowly convergent. Both of these problems have been solved with the development of the exponential operator MCSCF m e t h ~ d . ~ lWith - ~ ~ this approach, the orbital coefficients and CSF expansion coefficients are simultaneously optimized by using information about the first and second derivatives of the energy with respect to changes in these coefficients. If all the second derivative information is computed, an iterative procedure may be used which results in a drastic reduction in the number of iterations required. This is an example of using larger memories and changing basic algorithms to reduce both the 1/0 and the number of floating point operations required for a calculation. This full second-order scheme is practical for wave functions consisting of a few hundred C S F S ? ~ For larger CSF expansions, either a microiterative procedure may be used to solve for the second-order wave function corrections25-26or some of the sec(21) (22) (23) (24) (25)
B. Levy, Int. J. Quantum Chem., 4, 297 (1970). L. G . Yaffe and W. A. Goddard 111, Phys. Rev. A, 13, 1682 (1976). E. Dalgaard and P. Jorgensen, J . Chem. Phys., 69, 3833 (1978). R. Shepard, I. Shavitt, and J. Simons, J. Chem. Phys., 76,543 (1982). B. H. Lengsfield 111, J. Chem. Phys., 77, 4073 (1982).
The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 341 ond-derivative information may simply be discarded. The calculation of the required first- and second-derivative information in the exponential operator MCSCF method involves matrix and vector operations and results in efficient computational procedures on vector-oriented computers. The selection of configurations for MCSCF wave functions has also changed. Instead of the painstaking empirical selection of individual CSFs as was routine a decade ago, a more a priori approach is now standard. In the MCSCF implementation of the GVB method, e.g., the configurations are chosen to properly describe the dissociation of the molecular system into its constituent fragments. In the CASSCFZ7(complete active space SCF) and FORS2* (fully optimized reaction space) methods, two other currently popular approaches, the important (active) orbitals are determined from chemical intuition, aided perhaps by a few test calculations. The CSF set is then obtained by allowing all possible occupations of the electrons in the active orbitals. Since this CSF set grows very rapidly with the number of orbitals, it is sometimes necessary to restrict the C S F space with additional occupation or spin-coupling restrictions. One such approach is to divide the active orbitals and electrons into disjoint subsets, the division again being aided by chemical intuition and/or test calculations. Each subset may then be treated as is most appropriate, the possibilities ranging from allowing all possible distributions of the electrons among an orbital subset to allowing only double occupancy of the orbitals. The final wave function is formed from the product of the occupations of each of the disjoint sets. With these methods, MCSCF wave function expansions of several thousand CSFs are becoming routine, 100-1000 times larger than the expansions of a decade ago. A CI wave function is a linear combination of CSFs where only the expansion coefficients are optimized to give the lowest energy. Conventional C I methods19 consist of constructing the matrix representation of the Hamiltonian operator within the space of CSFs. The elements are written to an external storage device, such as a disk, as they are computed. Optimization of the expansion coefficients requires repetitive reading of these matrix elements from external storage to solve the matrix eigenvalue problem. Construction of the CI matrix using these techniques is very time consuming and for long CSF expansions requires large amounts of disk space. Conventional CI methods do not achieve a balance between the CPU and 1 / 0 times. These problems placed a premium on minimizing the C S F expansion length and prompted the development of elaborate, and unfortunately sometimes unreliable, CSF selection and extrapolation schemes. CSF expansions of several hundred to ten thousand were typically used with the conventional CI approach. An alternative approach to conventional CI is the direct CI method.29 In the direct CI method, the matrix-vector products are computed directly from a sorted integral list and a set of coupling constants. Although the direct CI method was initially very limited in the types of CSF expansions allowed, the development of the graphical unitary group approach30 (GUGA) generalized the method to allow arbitrary sets of CSFs; it is, however, most efficient for MR-SDCI (multiple-reference single and double replacement configuration interaction) calculations. The GUGA approach results in a very efficient computational scheme for the calculation of the required coupling coefficients, a time-consuming step in conventional CI, along with an effective organization and representation of the CSFs. The resulting method is appropriate for arbitrary spin states with any number (26) F. B. Brown, I. Shavitt, and R. Shepard, Chem. Phys. Lett., 105, 363 ( 1984).
(27) B. 0. Roos, P. R. Taylor, and P. E. M. Siegbahn, Chem. Phys., 48, 157 (1980). (28) K. Ruedenberg and K. R. Sundberg in “Quantum Science”, J.-L. Calais, 0. Goscinski, and Y. Ohm, Eds., Plenum Press, New York, 1976, p 505. (29) B. 0. Roos and P. E. M. Siegbahn in “Methods of Electronic Structure Theory”, H. F. Schaefer 111, Ed., Plenum Press, New York, 1977, p 277 (30) I. Shavitt in ‘The Unitary Group for the Evaluation of Electronic Energy Matrix Elements”, J. Hinze, Ed., Springer-Verlag, Berlin, 198 I , p 5 1.
Dunning et al.
348 The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 of open shells and is most efficient if the reference CSFs are constructed from a small set of internal orbitals while the MRSDCI CSF set is constructed from both the internal orbitals and a larger set of external orbital^.^' In the direct CI method efficient matrix-vector products are used to solve the eigenvalue equation. These products are constructed by reading blocks of integrals from external storage; each block of integrals is then used in a series of matrix-matrix product^,^*^^^ one matrix being a block of integrals and the other matrix being formed from a contiguous sequence of C S F coefficients. The high efficiency achieved by the direct CI procedure is at the expense of less flexible CSF selection schemes than those of the conventional CI method. Completely general CSF selection would result in the replacement of the efficient matrix-matrix multiply operations with a much less-efficient series of scalar operations. The C S F selection that is presently employed in the GUGA-based direct C I codes is only in the choice of the reference CSFs, usually those CSFs that provide a qualitative description of the relevant molecular deformations, e.g., the GVB configurations, and in the selection of the internal orbital components of the full MR-SDCI CSF expansion set.34 MR-SDCI expansions of several hundred thousand CSFs are not uncommon and for special cases have been as large as 1.5 million, roughly a factor of 1000 times larger than the expansion lengths of a decade ago. The fact that both the orbitals and C S F coefficients are optimized in MCSCF wave functions may be used to advantage in the calculation of the gradient of the electronic energy.35 The gradient may be calculated almost as efficiently for a general MCSCF wave function as for a S C F wave function, a special case of an MCSCF wave function with a single CSF. The only minor complication is that the two-particle density matrix is no longer a simple function of the one-particle density matrix. Although more difficult to evaluate than the gradient, the second-derivative of the energy with respect to a geometry deformation may also be calculated for MCSCF wave functions. Again, the fact that the orbitals and CSF coefficients are variationally determined is an advantage in the calculation of the Hessian matrix for MCSCF wave functions. The calculation of the gradient of a MR-SDCI wave function, however, poses additional difficulties which are just now being addressed.36
ij 1 u* 0
v
“i t ‘1.
x
3
W
W
a
1 I
1
--__ Reactants
Tran s i t ion Products State Figure 1. Definition of terms: energetics of chemical reactions (a) without and (b) with zero-point energy corrections.
where the sum runs over all of the 3N - 6 modes (3N - 7 modes for a saddle point). The enthalpy change at 0 K and the vibrationally adiabatic threshold for a chemical reaction are given by
+ (Ezpe@ - EZpereaCt) = AEb + (EzpesadPt Ezpereact)
AH, = AE,,,
AE,,,
-
(11) (12)
A”,is the energy change in the reaction at absolute zero; AE,,, IV. Calculation of Vibrational Energy Levels of Reactive Systems The energy defect and barrier height of a chemical reaction are simply
where Vsadptis the value of the potential at the saddle point and VrCactand Vprd are the values for the reactants and products, see Figure 1. The quantity AEb is sometimes referred to as the classical barrier to reaction. For an exoergic reaction Vprd will be more negative than V,,,,, and, thus, AEIxnwill be negative. Accurate prediction of the energetics of chemical reactions requires the addition of vibrational energies to the above electronic energies. In the harmonic approximation the zero point vibrational energy of the system at a stationary point is
(31) H. Lischka, R. Shepard, F. B. Brown, and I. Shavitt, Int. J . Quantum Chem., SlS, 91 (1981). (32) V. R. Saunders and J. H. van Lenthe, Mol. Phys., 48, 923 (1983). (33) R. Ahlrichs, H.-J. Bohm, C. Ehrhardt, P. Scharf, H. Schiffer, H. Lischka, and M. Schindler, J . Comput. Chem., 6, 200 (1985). (34) I. Shavitt in “Advanced Theories and Computational Approaches to the Electronic Structure of Molecule”, C. E. Dykstra, Ed., Reidel, Boston, 1984, p 185. (35) A. Banerjee, J. 0. Jensen, J. Simons, and R. Shepard, Chem. Phys., 87, 203 (1984). (36) D.J . Fox, Y. Osamura, M . R. Hoffmann, J. F. Gaw, G. Fitzgerald, Y . Yamaguchi, and H. F. Schaefer 111, Chem. Phys. Lett., 102, 17 (1983).
is the threshold for the reaction if the vibrational level occupancies are conserved along the reaction path. AH, and AE,,, are illustrated in Figure l . The vibrationally adiabatic threshold is expected to be less than the barrier height for the reactions of interest here. This is a result of the fact that the bonds are weaker at the saddle point than in the reactants and is so in spite of additional vibrational modes present in the composite system at the saddle point. The procedure for calculating the vibrational levels and corresponding normal modes used in our laboratory was put forward by Gwinn;2 see Harding and F ~ m l e r .This ~ ~ method employs the Hessian matrix, (5), in Cartesian coordinates, [XI. The Hessian matrix in Cartesian coordinates, which is of order 3N, is obtained from the Hessian matrix in internal coordinates, which is of order 3N - 6, by numerical differentiation. The Cartesian Hessian matrix is then transformed to mass-weighted coordinates, Le., if X A i is the ith Cartesian coordinate of atom A and X B jis the j t h Cartesian coordinate of atom B, then in mass-weighted coordinates the Hessian matrix elements are khe,/(mAmB)”2. Diagonalization of the mass-weighted Hessian matrix yields the fundamental frequencies and normal modes. For a nonlinear polyatomic molecule there are six zero (or nearly zero) eigenvalues of the mass-weighted Hessian matrix; these correspond to overall translation and rotation of the system. The 3N - 6 nonzero eigenvalues of the Hessian matrix are the squares of the vibrational frequencies (w,). For minima all of the nonzero eigenvalues are positive and the corresponding eigenvectors are the normal modes. Saddle points, on the other hand, have one (and only one) negative eigenvalue and, therefore, the asso(37) L. B. Harding and W. C. Ermler, J . Comput. Chem.. 6, 13 ( I 985).
Energetics and Mechanisms of Chemical Reactions REACTANTS
SADDLE POINT
The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 349
PRODUCTS
3 x 7 n
/
3"
Figure 2. Orbital correlation diagram for the 0 + H2reaction. On the oxygen atom the lobes represent the in-plane 2p orbitals while the circle represents the out-of-plane 2p orbital; the 2s orbital is not shown. On the hydrogen atoms the circles represent the 1s orbital.
ciated frequency will be imaginary (arxn). The eigenvector for the negative eigenvalue corresponds to the reaction coordinate at the saddle point.
V. Mechanisms and Energetics of Chemical Reactions: A Selection of Examples In this section we briefly present the major features of a selection of abstraction reactions which have been studied in the Theoretical Chemistry Group at Argonne National Laboratory over the past five years. As noted in the Introduction the examples have been chosen to be representative of general classes of reactions; in this way it is hoped that the conclusions drawn here will extend beyond the specific examples discussed. Because of space limitations and because this is a progress report on work carried out in the Theoretical Chemistry Group at the Argonne National Laboratory, we have limited the examples chosen to those actively studied in our Group. A . Orbital Description of Abstraction Reactions: q 3 P ) Hz.38 The changes in the electronic configuration of a molecular system resulting from chemical reaction can be qualitatively understood by determining how the orbitals of the reactants are transformed into those of the products. An orbital correlation diagram for the simple abstraction reaction 0 + H2 -OH + H (13)
+
+
is given in Figure 2. In (13) the reactants, O(3P) Hz('Zg'), give rise to two states of the collinear OH2 complex: a 311state and a 38-state. Of these, only the 311state directly correlates with the products, OH(*lI) H@). The 311state splits into two distinct states of 3A" and 3A' symmetry for noncollinear geometries. The GVB wave function for reaction 13 is (for brevity the oxygen 1s orbital is omitted)
+
where x is a general four-electron triplet spin function. In the reactants (la,3a,la) correlate with the 2s, 2pa, and 2pa orbitals of the oxygen atom and (2aa,2gb) correlate with the left and right Hz-bond orbitals; in the products (1 a,la,2a.&b) correlate with the l a , l and ~ left and right a-bond orbitals of the OH molecule, respectively, and 3a correlates with the 1s hydrogen orbital. The orbitals directly involved in the transfer of the bond pair from reactants to products in a chemical reaction comprise the active set. The remaining valence orbitals are referred to as the inactive (or semiactive) orbitals. The valence orbitals for (13) in the reactants, saddle point, and products regions are plotted in Figure 3. First let us consider the active orbitals, (2aa,2ub,3a). At the saddle point the 2aa orbital, which is the left (near) Hz-bond orbital in the reactants, has changed only slightly. However, the 2ab orbital, which correlates with the right (far) Hz-bond orbital in the reactants, has extensively delocalized onto the oxygen atom; at the saddle point it, in fact, has the qualitative character of a bonding combination of the right Hz-bond orbital and the oxygen ~
~~~
~
~~~~~~
(38) S. P. Walch, T. H. Dunning, Jr., and R. C. Raffenetti, J . Chem. Phys., 72, 406 (1980); T. H. Dunning, Jr., S. P. Walch, and A. F. Wagner in "Potential Energy Surfaces and Dynamics Calculations", D. G . Truhlar, Ed., Plenum Press, New York, 1981, p 329.
-
2pa orbital. In this way the ( 2 ~ & ? . 0bond b ) pair maintains a high overlap (S 0.85) as the bond is transferred from reactants to products and the barrier to reaction (- 12 kcal/mol) is only a small fraction of the energy required to break the H2 bond (109.5 kcal/mol) . The other active orbital is the 3a orbital which correlates with the oxygen 2po orbital in the reactants and the hydrogen 1s orbital in the products. In the saddle point region, the singly occupied 30 orbital builds in a nodal plane in the bonding region in order to remain orthogonal to the ( h a & , ) bond pair. The resulting change in the phase of the orbital between reactants and products does not affect the wave function, (14), with three active electrons; however, it does influence the reaction path when four electrons are actively involved;39see section VD. As can be seen in Figure 3, the inactive orbitals of the OH, complex are responsive to, but not directly involved in, the transfer of the electron pair from H2 to OH. Thus, the 1a orbital polarizes away from the near hydrogen atom, while the 1~ orbital polarizes slightly toward this atom, as a result of O H bond formation. The inactive orbitals do, however, play an important role in determining the relative energetics of abstraction and exchange reactions.4n The geometry of the saddle point is, of course, also determined by the interactions between the orbitals. Thus, proper alignment of the orbitals, Le., minimum disruption of the (20,,2Ub) bond pair, along with minimization of the repulsive interactions between the orbitals on the end atoms, especially the 3u and 2ab orbitals, dictate the collinear saddle point found for reaction 13. B. Trends in A Series of Abstraction Reactions: H HX (X = F-Z).40Since the mid-1930's there have been numerous attempts to relate the properties of the transition state to the overall energetics of the reaction. From this ~ o r k ~ l two - ~ 'propensity rules have emerged: (1) Very exoergic reactions have small activation energies, with the activation energy increasing with decreasing exoergicity (Evans-Polanyi-type relationship4'). (2) Very exoergic reactions have transition states which resemble the reactants, with the transition states becoming more like the products with decreasing exoergicity (Hammond postulate41). The reactions of hydrogen atoms with the hydrogen halides to yield molecular hydrogen and a halogen atom H HX H2 X (15)
+
+
+
+
is one of the simpliest series of reactions for which a consistent set of calculations have been reported for all of the members of the family X = F-L4" These calculations permit a detailed examination of the systematic trends observed in this prototypical series of reactions. The calculated energy defects for the reactions, (1 5), decrease monotonically from 24.9 kcal/mol in H H F to -29.6 kcal/mol in H + HI; the errors in the calculated defects decrease from 5.9 ( H + HF) to 1.3 kcal/mol ( H HI). As noted above, an inverse correlation is expected between the barrier height and the exoergicity of the reaction and such a trend is found here. Further, since an endoergic reaction is just an exoergic reaction in reverse, for large endoergicities we expect the barrier height to approach the endoergicity. A functional relationship between the barrier height and the reaction energy defect which reflects this symmetry is
+
+
where AEo(AErxn),the intrinsic barrier height, is a symmetric function of the reaction energy defect which approaches zero for large defects and H(AE,,,) is the Heaviside function, Le., H = (39) (40) (41) (42) (43) (44) (45) (46) (47)
W. A. Gooddard 111, J . Am. Chem. SOC.,94, 793 (1972).
T. H. Dunning, Jr., J . Phys. Chem., 88, 2469 (1984).
M. G. Evans and M. Polanyi, Trans. Faraday SOC.,34, 11 (1938). G. S. Hammond, J . Am. Chem. SOC.,77, 334 (1955). H. S.Johnston and C. Parr, J . Am. Chem. SOC.,85, 2544 (1963). R. A. Marcus, J . Phys. Chem., 72, 891 (1968). M. H. Mok and J. C. Polanyi, J . Chem. Phys., 51, 1451 (1969). N. Agmon and R. D. Levine, Chem. Phys. Leu., 52, 197 (1977). N. Agmon, J. Chem. SOC.,Faraday Trans. 2, 74, 388 (1978).
350
Dunning et al.
The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 ~' ~
ACTIVE ORBITALS
INACTIVE ORBITALS
la
3'0 L a b -
3a
177
r--,
-----I
1
7
7
I
~
I
~
u1 O H H
0
O H H
O H
O H H
H H
H
+
Figure 3. Contour plots of the GVB orbitals for the 0 H2 reaction for the reactants, saddle point, and product. The contour lines are evenly spaced with increments of 0.05 au; only the contours between -0.5 and +0.5 au have been plotted. Positive contours are traced in solid lines, negative contours in short dashed lines. and nodal lines in long dashed lines. Positions of the nuclei are denoted by an asterisk.
0.8
40
=
30
n
s
E" 1 '.
X
I
0.4
QL
a
20 Y
0.2-
W
n
w
H+HI
a 10
0.0 0.0
0.4
0.6
0.8
Figure 5. Relationship between the bond extensions, AR" and ARHX, for the H HX abstraction reactions. The circles refer to the calculated points, the solid line to eq 19.
+
0 0
0
-20
AE,,,
20
40
(kcal/mol)
Figure 4. Dependence of the barrier height on the reaction energy defect in the H HX abstraction reactions. The circles refer to the calculated points. The solid line includes H HF in the least-squares fit, eq 16; the dashed h e does not.
+
0.2
+
0 if AE,,, 10,and H = 1 if AE,,, > 0. Previous work suggested As shown an exponential relationship between AEo and AE,xn.40 in Figure 4,the calculated barrier heights are well described by
AEo(AErx,)= 8.97e4,03521AErxnI (17) In fact, if we delete the H + H F results, for which the error in the calculated energy defect is largest (5.9 kcal/mol), the remaining barrier heights are predicted to an accuracy of k O . 1 kcal/mol by (1 7). The large discrepancy evident in Figure 4 for the H + HF reaction may be due to the error in the calculated energy defect for this reaction and not to any deviation from the relationships 16 and 17; note, e.g., the excellent fit of the barrier heights reported in ref 40. Estimates of the barrier heights obtained from analysis of the available experimental data are -33 ( H + HF), -4 (H + HCl), -1.5 ( H + HBr), and -0 kcal/mol ( H + HI); see, e.g., ref 48 and 49. These estimates are to be compared to the calculated (48) J. T. Muckerman, Surface 5 , unpublished.
barrier heights of 3 1.4, 7.3, 3.8, and 2.3 kcal/mol, respectively. Similar calculations on the H H2 reaction, for which the barrier is known to an accuracy of f0.2 kcal/mol, predict a barrier which is too large by approximately 2 kcal/mol. As expected, the errors are somewhat larger in the H + H X abstraction reactions. If we define the bond extensions in the abstraction reactions by (18) AR, = R s p ( i )- R,(i)
+
where Rsp(i)is the length of bond i at the saddle point and R,(i) the bond length in the corresponding diatomic molecule, then for very exoergic reactions we would expect that AI?" would be large and ARHX would be small. For very endoergic reactions the opposite would be expected. This is indeed the case for the H HX reactions as can be seen in Figure 5. In fact, the bond extensions in reaction 15 are remarkably well represented by e - 3 . 5 7 5 A R ~ ~+ e-3.575AR~x = 1 (19)
+
The average deviation of the constant in (19) from unity is just 0.02. This equation is plotted as the solid line in Figure 5. A relationship of this form can be derived from BEBO (bond energy-bond order) consideration^,^^ although the constants in the (49) R. N. Porter, L. B. Sims, D. L. Thompson, and L. M. Raff, J . Chem. Phys., 58, 2855 (1973).
Energetics and Mechanisms of Chemical Reactions
l-----
The Journal of Physical Chemistry, Vol. 90, No. 3, 1986 351 a barrier of >75 kcal/mol for this pathway. Consider now the optimal approach in which the C H and H 2 bond axes are nearly parallel. The active orbitals for this pathway are plotted in Figure 7 . As is evident in the figure, the C H H, reaction may be regarded as an attack by one of the C H lone pair orbitals on the H, bond pair; compare, e.g., the behavior of the 2a
