- Email: [email protected]
Some observations on the computed forms of the bond rupture surfaces in H2, CH4 and C2H6
Volume 13 1,number 6
CHEMICAL PHYSICS LETTERS
2 1 November 1986
SOME OBSERVATIONS ON THE COMPUTED FORMS OF THE BOND RUPTURE SURFACES IN H,, CH, AND C,H,
E.M. EVLETH ’ and E. KASSAB Dynamique des Interactions Mokulaires, 4, Place Jussieu, 75230 Paris, France
ER 271, Tour 22, Uniuersitt Pierre et Marie Curie,
Received 15 August 1986; in final form 12 September 1986
The potential surfaces for bond ruptures in Ha, CH4, and C2He are examined in their reducedenergy forms. For CH SCF level are nearly identical to bond rupture in methane, the 6-31G** reduced energies found at the two-confguration previously published higher-level treatments. It is shown that reduced forms of the potentialenergy surface are linear in logarithmic space along considerable portions of the bond rupture coordinates in the three structures examined. This linearity enables a comparison of different computational methodologies to be made and allows assignment of regions of possible artifactual behavior. In addition, it enables computational procedures to be formulated in which the surfaces for radical-radical recombinations can be reliably explored at reduced cost.
1. Introduction
Dynamic modelling [l-4] of radical-radical recombination reactions is sensitive to the surface forms in bond distance regions (2-5 A) where the fragment-fragment interaction energies are small. These, as well as analysis [S] based on transitionstate theory, assume zero activation energies for the J = 0, non-zero energy corrected surfaces for radicalradical recombinations. At intermediate distances (i.e. 2-5 A), the generally assumed forms of radicalradical recombination surfaces are exponential in character [5,6]. These simple forms are derived for bonding interactions even though rM6 and higherorder terms will be necessary to describe the longrange interactions (e.g., see ref. [7]). The actual analytical form necessary to describe adequately the intermediate region for bond ruptures is largely untested by direct quantum-mechanical calculation. There is also some question as to what quantummechanical methodology to use in this region. For bond ruptures, the intermediate-to-long distance surface regions are conceptually open-shell in character, and closed-shell single-configuration SCF methods used in the vast majority of electronic struc-
ture calculations are invalid. Of the broadly available methods, the MP4/UHF(6-3 lG**) and MP4/ RHF(6-31G**) surfaces [8,9] for the CH bond rupture in CH4 have drawn considerable comment [ l4, 8-131, This particular surface is artifactual iu the region of dynamical interest iu comparison with other ab initio surfaces [ 10,121. Using this surface to model a tight Morse equation, the computed recombination rate constants are too low [2,3,11]. The surface appears to be improved by removing some of the energetic effects due to spin contamination [ 133. However, the single-reference MP method has convergence problems at intermediate distances for bond ruptures [2,10,14]. and the use of a multireference space built from a MC SCF MO basis is recommended. Nevertheless, we will show that calculations at high levels of theory may not always be necessary and, in favorable cases, surface forms can be obtained at a lower level.
2. Technical aspects The calculations presented here used the programs MONSTERGAUSS [15] and GAMESS *. The * We thank M. Dupuis, IBM Corporation, for a 1984 IBM
’ To whom correspondence should be addressed.
0 009-26 14/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
version of this program.
475
CHEMICAL PHYSICS LETTERS
Volume 13 I, number 6
MONSTERGAUSS MC SCF involves only doublereplacement configurations. In some cases this is a serious limitation and artifacts can be encountered. We have confnmed that this does not seem to be the case for the systems treated here by (i) comparing computations obtained at the double t single replacement level with GAMESS and (ii) verifying that the 2C SCF energies at very large distances equalled the sum of the RI-IF open-shell energies of the fragments. The two-configuration and larger MC SCF doublereplacement (2C SCF, 1OCSCF, and 17C SCF) computations were generally initiated from single-configuration SCF MOs obtained at about lSR,. Careful selection was then made of the orbital space undergoing rupture. This space was used in the 2C SCF computations. In obtaining our potential energy curves (figs. l-4), we used the MC SCF MOs obtained at the previous geometry as the initial guess for the next point. This procedure speeded subsequent MC SCF convergence. The most time consuming computations reported here were the gradient optimizations at the 2C SCF/3-21G level for ethane. These optimizations took between 2 and 5 min each on a NAS 9080 (8 MIPS). Based on these figures, computing the Hessian matrix components at the DZP small MC SCF level along the reaction pathway seems to be computationally feasible.
3. Results and comments 3.1. Some expectations on the forms of the surfaces The Morse equation is one of the most generally exploited surface forms for bond ruptures [ 11, E(r) = De {l - exp[-B(r
- r,)] }2 ,
(1)
where E(r) is the bonding energy at distance r, D, the bond dissociation energy and B = (k/m,)1 12, where k is the force constant. Eq. (1) can be revised to define a reduced energy, ER, which is unity (or -1) [lo] at re and zero at infinity: E, = 1 - E(r)/D, = 2 exp [-B(r - r,)] - exp [-2L?(r - re)] .
476,
(2)
21 November 1986
At large values of r the slope, d ln ER/dr, approaches the value -B. A region of linear behavior of In E, occurs at large r. However, the Morse equation yields energies which are too low compared to ab initio results (e.g. see ref. [ lo]) at intermediate r values. In fact, a number of energy and force-constant-bonddistance empirical relationships are linear in log space [6]. However, ab initio bond rupture calculations have not been generally tested for such behavior. In the case of the CH rupture in methane [8-10,12,13], we know of no In ER tests to see if linear regions occur. However, Cobos and Troe [ 1l] called attention to the problem of the distance dependence of an effective Morse B parameter. If In ER linear regions do
Table 1 6-31G** full CI a) 1OC SCF surface for H2 HH distance
Energy
(A)
(au)
Reduced energy b)
In ER
ER 0.6 0.7 0.74186 c) 0.75 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.6 1.7 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.25 3.5 4.0 10.0, 20.0
-1.146808 -1.164044 -1.165149 -1.165068 -1.162863 -1.152214 -1.136981 -1.119922 -1.102641 -1.086109 -1.070912 -1.045601 -1.035590 -1.027192 -1.014869 -1.007071 -1.002428 -0.999772 -0.998290 -0.997473 -0.996946 -0.996694 -0.996513 -0.996466
0.89 0.99 1.00 1.00 0.99 0.92 0.83 0.73 0.63 0.53 0.44 0.29 0.23 0.18 0.11 0.063 0.053 0.019 0.011 5.9 x 2.8 X 1.4 x 2.8x -
1o-3 1O-3 1o-3 1O-4
-0.12 -0.01 0.00 -0.00 -0.02 -0.08 -0.18 -0.31 -0.46 -0.63 -0.82 -1.23 -1.46 -1.70 -2.22 -2.77 -3.34 -3.95 -4.53 -5.12 -5.86 -6.61 -8.19 -
___a) 1OC SCF is equivalent to a full CI, see text. bl ER= (ED - ZEH)/(Ee - ~EH), where Ee is the energy at 0.74186 A, 2EH = 0.996466 au and ED is the energy at the distance given;Ee - EH = 0.168683 au (105.8 kcal/ mole). c) Ref. [ 171, optimized at the MP4 level.
Volume 131,number 6
CHEMICAL PHYSICS LETTERS
both (A) energy versus r(HH) and (B) In ER versus r(HH). The linearity of the In ER plot (fig. IB) extends well into the large H-H separation region, where any energy change is impossible to visually perceive (fig. 1A). The effective Morse B value in this linear region (~2.9 A-l) is about 50% higher than the Morse value [6]. This is also the case for methane and ethane. A m(energy) test on the 72-term variational ground state surface [ 181 of H2 also showed linear behavior in the 2-4 A region. We have not examined the very long distance region, dominated by r-6, r- 8, r--lo terms [7], to see where the departure from linear exponential behavior occurs. We also ob-
generally exist, the surface form of ab initio results can be anticipated and modelling schemes can be simplified [ 11. 3.2. Behavior of the H2 surface Table 1 shows the energy-distance behavior of the 1OCSCF/6-31G** H-H surface. This MC SCF level is equivalent to a full CI [ 161 (see comparative calculations at 0.74186 A shown in fig. 1). The MP4 calculation [ 171 does not represent a full CI even for so small a system. Also shown in table 1 are the reduced energies, ER, as a function of distance. Fig. 1 shows
-I
2 1 November 1986
oo-
-106. -108. -IlO2 ^ -112. z LK -114. iii W -I
l6*
-I
18.
SCF = -I 131272 MW = -1.164555 SDCI = -I 165149 IOCSCF = -I 165149
(A)
-7-6-
05
07
09
II
I5 DISTANCE,
20 H-H
25
30
35
A
Fig. 1. (A) Computed energy at the 6-31G** level for two- and ten-configuration MC SCF. The latter calculation is equivalent to a full CI as shown by reference calculations at 0.74186 A. See ref. [ 171 for MP4 calculations. (B) bi,#?R versus H-H distance. See table 1 for the energy details and method of computing ER. 477
21 November 1986
CHEMICAL PHYSICS LETTERS
Volume I3 1, number 6
served that HF and OH 2C SCF/6-3 lG* ruptures in HF and Hz0 showed linear In E, behavior in the 23.75 A region.
[ 131 projected UMP4/6-3 lG** values for the last three distances (0,27,0.088, and 0.023, respectively) are marginally different from the MRCD CI or 2C SCF values but significantly better than the unprojected values [9] (0.16,0.037, and 0.008). Obviously, what does change in all these calculations is the computed D, value for the CH bond. This value is difficult to precisely compute [ 191 and one can avoid doing so by using the experimentally known value to calculate the surface features from the theoretical reduced energies. For methane, the 2C SCF/6-3 lG** reduced energies serve equally well for scaling purposes as the MRD CI values for the CH region up to 3.5 A. Whether the progressive departure of the 2C SCF energies from those obtained at the MRD CI level above 3.0 A is of importance in dynamical modelling will require comparative computations. The 2C SCF procedure describes the first-order bonding interactions while ignoring other interfragment correlation effects. The 2C SCF and MRD CI
3.3. Behavior of the CH4 = CF3 + H surface Recently, Hirst [ 121 has published both the tetrahedral and optimized C-H bond rupture surfaces for CH, at the MRD CI/2C SCF(6-3 lG**) level. Table 2 and fig. 2 show a comparison of Hirst’s tetrahedral MRD CI computations with our 2C SCF computations. The E, values for both sets of computations are within the circles shown in the energy plot (fig. 2A) and are only perceptively different on the log scale (fig. 2B) above 3.5 A. The E, values obtained for the four points generated (O-73,0.3 1,O. 10 and 0.029 at R = 1.5,2.0,2.5, and 3.0 A, respectively) in the large basis set MRD CI calculations of Brown and Truhlar [lo] were not significantly different from either Hirst’s or our own values. Schlegel’s recent
Table 2 6-31G* * “tetrahedral surface” calculation for CH4 = CHs + H CH3-H distance
MRD CI energy a)
(A)
(au)
ERb)
In ER
-40.391645
1.00 d)
0.00
1.12 1.136
-40.390221
0.99
0.00
-40.379019 -40.338125 -40.296228 -40.262673 -40.239320 -40.224704 -40.216210 -40.211683
0.93 0.71 0.48 0.30 0.18 0.096 0.051 0.027
-0.07 -0.34 -0.73 -1.20 -1.74 -2.33 -2.98 -3.63
-40.207303 -40.206766
2.9 x 1o-3 -
-5.84 -
-40.206773
0.00 d)
20.0 100.0
ER
In ER
(au)
1.086
1.25 1.5 1.75 2.0 2.25 2.5 2.75 3.0 3.25 3.50 3.75 4.0 5.0
2C SCF energy c)
NC -40.217069
1.00 e)
-40.207665 -40.171087 -40.132646 -40.101602 -40.079687 -40.065886 -40.058043 -40.053915 -40.051837 -40.050809 -40.050162 -40.049909
0.94 0.73 0.50 0.31 0.18 0.096 0.049 0.024 0.012 5.9 x 1o-3 2.1 x 1o-3 5.4 x 10-4
40.049818
0.00 e)
0.00 -0.06 -0.32 -0.70 -1.17 -1.72 -2.34 -3.01 -3.71 -4.42 -5.13 -6.19 -7.52
a) Ref. [ 121, CH = 1.086. b, ER = (ED - E,)/(Ee - 2E,), where Ee is the energy at 1.086 A, E, the largest distance computed and ED the energy at the distance given. c, 2C SCF, HOM02/LUM02 pair. d) E 1.o86 - Eloo = 0.184872 au (116.0 kcal/mole). e, Et.12 - Ezo = 0.167525 au (105.1 kcal/moIe).
478
Volume I3 1, number 6
21 November 1986
CHEMICAL PHYSICS LETTERS
0
v
”
/d+-
a
v
-3
G E
(A)
a-8 C
6 0.5.
/
E 3
if! a
x
1.0.:l
/
/
/
,e I
/
, “, ’ ,m ’
-5 -4 3 -3 3
0 /m O’ZCSCF @ = MRCVPCSCF
_;,-G
x
x I -4-6 ,w /*
(B)
--2
@ =MRCI/2CSCF(OPTIMIZED) 0 = MP4AJHF
1.0
1.5 C-H
2.0
DISTANCE,
2.5
3.0
--I
3.5
,O 4.0
A
Fig. 2. (A) Reducedenergy curve for the tetrahedral CH rupture in methane at the 6-31G** 2C SCF and MRD CI/ZC SCF levels [ 8,9]. All values fall within the circles shown. (B) Comparison of lnER versus C-H distance. The MRCD optimized results are reported in refs. [ 8,9]. The MP4/UHF values are from the UHF/MP4 SDQ results reported in ref. [ 81.
calculations show sufficient superimposable linearity of the ln ER plots to permit a linear extrapolation into the longer-distance regions where deviation from linearity occurs at the 2C SCF level. Moreover, geometry optimization [ 121 along the reaction path has the effect of improving ln ER linearity at long distances while maintaining superimposability at shorter distances (fig. 2). It is also worth noting that the MP4 In ER curve shown in fig. 2 enables us to detect the point (2 A) where artifactual behavior begins to occur. However, the projected UMP4 In ER plot was found to be linear and nearly superimposable on the MC SCF curves. The fact that such behavior was obtained indicates to us that a major portion of the artifactual behavior was removed by the projection procedure. The presence of a linear region in the In ER plot at large r means that an effective constant Morse B value can be assumed for the region of dynamical interest. In the cases of methane or ethane (see below), these effective B values (=2.6-2.8 A-1) are larger than those used for their Morse curves (al .81.9A-l) [l]. The major question posed by the favorable com-
parison of the computations shown in table 2 and fig. 2 is whether one can generally obtain satisfactory curve forms at lower levels of theory. Our response is both guarded and optimistic. First of all, one assumes that the ln ER plots will eventually deviate from linearity at long distances. Secondly, an extended region of linearity may not even be encountered for interacting polar or polarizable fragments. Thirdly, one anticipates both a basis set and MC SCF level dependence of the effective B values. Without presenting the details of the computations, this latter effect is shown ln fig. 3 for the tetrahedral C-H bond rupture in CH4 at the 2C SCF/STO-3G and 17C SCF/ STO-3G levels and at the 17C SCF/3-21G level. All these In ER surfaces show regular behavior, the first with curvature. The effective B values in the 2.53.0 A region vary both with basis set and MC SCF level. At the STO-3G level, increasing the MC SCF level from two to seventeen configurations improves the linearity of the In ER plot and reduces the effective B value from 3.6 to 2.4 A-l. The In ER behavior of the 17C SCF/STO-3G surface is fortuitously close to the 6-31G** (B = 2.6-2.7) surfaces (fig. 2) up to 3.5 A. However, the 17C SCF/3-21G slope is 479
Volume 13 1, number 6
CHEMICAL PHYSICS LETTERS
21 November 1986
7
-6.
1.5
2.0 C-H DISTANCE,
2.5
3.0
A
3.5
Fig. 3. Comparison of lnE~ versus C-H distance in a tetrahedral scan for different basis sets and MC SCF levels. The values of De: 2C SCF/STO-3G, 131 kcal/mole; 17C SCF/STO-3G, 147 kcal/mole, and 17C SCF/3-21G, 119.8 k&/mole. R(CH) = 1.09 A.
smaller (1.9 A-l) than the 6-3 lG** values. Due to size consistency correlation artifacts, the computed bond energies (caption, fig. 3) vary with the MC SCF level. Therefore, one should not generally expect the superimposability of 2C SCF and MRD CI In ER curves that was found in fig. 2. The use of these small basis sets magnifies the problems but also illustrates some possibilities for exploitation. For small systems, one will normally use larger basis sets capable [lo]
of describing the intermediate-distance bonding regions. In this case, the necessary surface information for dynamic modelling may be extractable at the 2C SCF:level but ln ER versus r slopes should be verifled by performing calculations at the MRD CI level. If one can then find a smaller basis set, small MC SCF calculation which mimics the large basis set results at the ER level, one can risk treating a larger analogous system (e.g. R-CH2-H rupture in a larger system)
Table 3 2C SCF/3-21G optimized surface for the reaction CaHe = CHs + CHs
energy
CH, -CH3 distance (A)
2C SCF (au)
1.572 1.75 2.00 2.25 2.50 3.00 3.25 3.50 4.00 20.00
-78.808295 -78.800025 -78.711443 -78.740802 -78.716136 -78.689655 -78.683494 -78.681649 -78.679395 -78.678781
b,
ER a)
1 .oo 0.94 0.72 0.4% 0.29 0.084 0.036 0.021 4.7 x 1o-3 0.000
ERR
0.00 -0.07 -0.33 -0.74 -1.24 -2.48 -3.31 -3.87 -5.35 a)
Optimized parameters R(C-H)
L CCH
1.084 1.081 1.077 1.075 1.073 1.071 -1.071 1.071 1.070 1.073
110.4 108.3 105.7 103.3 101.0 96.9 95.0 93.1 89.4 c) 90.3 c)
a) E 1.572- E2o.o = 0.129514 au (81.3 kcal/mole). b) 2C SCF CC harmonic force constant = 3.80 mdyn A-‘. c) One expects 90” for these large distances and the values shown are probably artifacts of the optimization procedure in a region where the interaction energies are small.
480
CHEMICAL PHYSICS LETTERS
Volume 131, number 6
where neither large basis sets nor MRD CI comparisons are technically or financially feasible. 3.4. Behavior of the C,H, = 2CH3 surface Table 3 shows the optimized reaction surface for the CC rupture in ethane at the 2C SCF/3-21G level. Fig. 4 shows both the energy and In ER behavior as a function of CC distance. The In ER surface shows linear behavior in the 2.5-4.0 A region, similar to the CH rupture in methane. The 2C SCF/3-2lG level D, (table 3,81.3 kcal/mole) is reasonably close to the experimental reaction enthalpy (88 kcal/mole, less ZPE) [ 51. These computations were sufficiently inexpensive to suggest that, for dynamical applications, a 2C SCF/6-3 lG** level treatment can be envisaged. A few additional intermediate-distance MRD CI calculations will determine whether the effective B value can be obtained at the lower 2C SCF level.
2 I November 1986
We finally note that MNDO CI calculations [20] indicate that barriers exist for both CH, + H and CH, + CH, as well as for some other recombinations. Based on the results reported here and other work in progress, we view these barriers as artifacts. We are informed that the AM1 parameterization yields barrier-free surfaces for those radical-radical recombinations that have been reinvestigated so far [21].
4. Conclusions The surface forms for the bond ruptures in H,, CH4, and C2H6 show linear behavior in log(reduced energy) space for the region in the early stages of bond formation. This behavior can be used for diagnostic purposes. Comparison of the methane 2C SCF/ 6-31G** calculations with MRD CI results indicates that the reduced-energy surfaces are nearly identical. We conclude that the possibility exists of obtaining the bond-rupture surface of high-level calculations at lower levels.
Acknowledgement
(A)
CH-CH 3
3
-
CH3
+
CH3
We acknowledge a number of discussions with Professor J.J. Dannenberg. We also thank the CNRS computing center at CIRCE, Orsay, France for use of their facilities. This research project was financed in part by research contract ATP-CNRS-920042.
References
-6.
[l] C.J. Cobos and J. Troe, J. Chem. Phys. 83 (1985) 1010, and references therein.
-5’
[ 21 R.J. Duchovic and W.L. Hase, J. Chem. Phys. 82 (1985)
-4
P
1.5
20 c-c
3.0 &CE,
3.5
A
Fig. 4. Optimized 3-21G 2C SCF level CC bond rupture in ethane. Data taken from table 3. (A) Energy versus CC distanCe. (B) hI ER plot Similar t0 figs. 1-3.
3599. [3] W.L. Hase and R.J. Duchovic, J. Chem. Phys. 83 (1985) 3448, and references therein. [4] J.F. LeBlanc and P.D. Pacey, J. Chem. Phys. 83 (1985) 4511, and references therein. (5) S.W. Benson, Thermochemical kinetics (Wiley, New York, 1976) pp. 86-100,160-166. (61 H.S. Johnston, Gas phase reaction rate theory (Ronald Press, New York, 1966) chs. 4,11. [7] W. Koiix, Intern. J. Quantum Chem. 1 (1967) 169. (81 R.J. Duchovic, W.L. Hase, H.B. Schlegel, M.J. Frisch and K. Raghavachari, Chem. Phys. Letters 89 (1982) 120.
481
Volume I3 1, number 6
CHEMICAL PHYSICS LETTERS
19) R.J. Duchovic, W.L. Hase and H.B. Schlegel, J. Phys. Chem. 88 (1984) 1347. [ 101 F.B. Brown and D.G. TruhIsr, Chem. Phys. Letters 113 (1985) 441. [ 111 C.J. Cobos and J. Troe, Chem. Phys. Letters 113 (1985) 419. [12] D.M. Hirst, Chem. Phys. Letters 122 (1985) 225. [13] H.B. Schlegel, J. Chem. Phys. 84 (1986) 4530. [14] N.C. Handy, Faraday Symp. Chem. Sot. 19 (1984) [15] :A. Peterson and R.A. Poirier, MONSTERGAUSS, Department of Chemistry, University of Toronto.
482
2 1November 1986
1161A.C. Wahl and G. Das, in: Methods of electronic strut-
ture theory, Vol. 3, ed. H.F. Schaefer III (Plenum Press,
New York, 1977) pp. 51-78, and references therein. (171 M.J. Frisch, J.A. Pople and J.S. Binkley, J. Chem. Phys. 80 (1984) 3265. [181 W. Kolos and L. WoIniewicz, J. Mol. Spectry. 54 (1975) 303. 1191P.E.M. Siegbahn,Chem. Phys. Letters 119 (1985) 515. [201J.J. Dannenberg and K. Tanaka, J. Am. Chem. Sot. 107 (1985) 671; J.J. Darmenberg, J.C. Rayez, M.T. Rayez-Meaume and P. Halvich, J. Mol. Struct. THEOCHEM 123 (1985) 343. [211 J.J. Dannenberg, private communication.
CHEMICAL PHYSICS LETTERS
2 1 November 1986
SOME OBSERVATIONS ON THE COMPUTED FORMS OF THE BOND RUPTURE SURFACES IN H,, CH, AND C,H,
E.M. EVLETH ’ and E. KASSAB Dynamique des Interactions Mokulaires, 4, Place Jussieu, 75230 Paris, France
ER 271, Tour 22, Uniuersitt Pierre et Marie Curie,
Received 15 August 1986; in final form 12 September 1986
The potential surfaces for bond ruptures in Ha, CH4, and C2He are examined in their reducedenergy forms. For CH SCF level are nearly identical to bond rupture in methane, the 6-31G** reduced energies found at the two-confguration previously published higher-level treatments. It is shown that reduced forms of the potentialenergy surface are linear in logarithmic space along considerable portions of the bond rupture coordinates in the three structures examined. This linearity enables a comparison of different computational methodologies to be made and allows assignment of regions of possible artifactual behavior. In addition, it enables computational procedures to be formulated in which the surfaces for radical-radical recombinations can be reliably explored at reduced cost.
1. Introduction
Dynamic modelling [l-4] of radical-radical recombination reactions is sensitive to the surface forms in bond distance regions (2-5 A) where the fragment-fragment interaction energies are small. These, as well as analysis [S] based on transitionstate theory, assume zero activation energies for the J = 0, non-zero energy corrected surfaces for radicalradical recombinations. At intermediate distances (i.e. 2-5 A), the generally assumed forms of radicalradical recombination surfaces are exponential in character [5,6]. These simple forms are derived for bonding interactions even though rM6 and higherorder terms will be necessary to describe the longrange interactions (e.g., see ref. [7]). The actual analytical form necessary to describe adequately the intermediate region for bond ruptures is largely untested by direct quantum-mechanical calculation. There is also some question as to what quantummechanical methodology to use in this region. For bond ruptures, the intermediate-to-long distance surface regions are conceptually open-shell in character, and closed-shell single-configuration SCF methods used in the vast majority of electronic struc-
ture calculations are invalid. Of the broadly available methods, the MP4/UHF(6-3 lG**) and MP4/ RHF(6-31G**) surfaces [8,9] for the CH bond rupture in CH4 have drawn considerable comment [ l4, 8-131, This particular surface is artifactual iu the region of dynamical interest iu comparison with other ab initio surfaces [ 10,121. Using this surface to model a tight Morse equation, the computed recombination rate constants are too low [2,3,11]. The surface appears to be improved by removing some of the energetic effects due to spin contamination [ 133. However, the single-reference MP method has convergence problems at intermediate distances for bond ruptures [2,10,14]. and the use of a multireference space built from a MC SCF MO basis is recommended. Nevertheless, we will show that calculations at high levels of theory may not always be necessary and, in favorable cases, surface forms can be obtained at a lower level.
2. Technical aspects The calculations presented here used the programs MONSTERGAUSS [15] and GAMESS *. The * We thank M. Dupuis, IBM Corporation, for a 1984 IBM
’ To whom correspondence should be addressed.
0 009-26 14/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
version of this program.
475
CHEMICAL PHYSICS LETTERS
Volume 13 I, number 6
MONSTERGAUSS MC SCF involves only doublereplacement configurations. In some cases this is a serious limitation and artifacts can be encountered. We have confnmed that this does not seem to be the case for the systems treated here by (i) comparing computations obtained at the double t single replacement level with GAMESS and (ii) verifying that the 2C SCF energies at very large distances equalled the sum of the RI-IF open-shell energies of the fragments. The two-configuration and larger MC SCF doublereplacement (2C SCF, 1OCSCF, and 17C SCF) computations were generally initiated from single-configuration SCF MOs obtained at about lSR,. Careful selection was then made of the orbital space undergoing rupture. This space was used in the 2C SCF computations. In obtaining our potential energy curves (figs. l-4), we used the MC SCF MOs obtained at the previous geometry as the initial guess for the next point. This procedure speeded subsequent MC SCF convergence. The most time consuming computations reported here were the gradient optimizations at the 2C SCF/3-21G level for ethane. These optimizations took between 2 and 5 min each on a NAS 9080 (8 MIPS). Based on these figures, computing the Hessian matrix components at the DZP small MC SCF level along the reaction pathway seems to be computationally feasible.
3. Results and comments 3.1. Some expectations on the forms of the surfaces The Morse equation is one of the most generally exploited surface forms for bond ruptures [ 11, E(r) = De {l - exp[-B(r
- r,)] }2 ,
(1)
where E(r) is the bonding energy at distance r, D, the bond dissociation energy and B = (k/m,)1 12, where k is the force constant. Eq. (1) can be revised to define a reduced energy, ER, which is unity (or -1) [lo] at re and zero at infinity: E, = 1 - E(r)/D, = 2 exp [-B(r - r,)] - exp [-2L?(r - re)] .
476,
(2)
21 November 1986
At large values of r the slope, d ln ER/dr, approaches the value -B. A region of linear behavior of In E, occurs at large r. However, the Morse equation yields energies which are too low compared to ab initio results (e.g. see ref. [ lo]) at intermediate r values. In fact, a number of energy and force-constant-bonddistance empirical relationships are linear in log space [6]. However, ab initio bond rupture calculations have not been generally tested for such behavior. In the case of the CH rupture in methane [8-10,12,13], we know of no In ER tests to see if linear regions occur. However, Cobos and Troe [ 1l] called attention to the problem of the distance dependence of an effective Morse B parameter. If In ER linear regions do
Table 1 6-31G** full CI a) 1OC SCF surface for H2 HH distance
Energy
(A)
(au)
Reduced energy b)
In ER
ER 0.6 0.7 0.74186 c) 0.75 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.6 1.7 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.25 3.5 4.0 10.0, 20.0
-1.146808 -1.164044 -1.165149 -1.165068 -1.162863 -1.152214 -1.136981 -1.119922 -1.102641 -1.086109 -1.070912 -1.045601 -1.035590 -1.027192 -1.014869 -1.007071 -1.002428 -0.999772 -0.998290 -0.997473 -0.996946 -0.996694 -0.996513 -0.996466
0.89 0.99 1.00 1.00 0.99 0.92 0.83 0.73 0.63 0.53 0.44 0.29 0.23 0.18 0.11 0.063 0.053 0.019 0.011 5.9 x 2.8 X 1.4 x 2.8x -
1o-3 1O-3 1o-3 1O-4
-0.12 -0.01 0.00 -0.00 -0.02 -0.08 -0.18 -0.31 -0.46 -0.63 -0.82 -1.23 -1.46 -1.70 -2.22 -2.77 -3.34 -3.95 -4.53 -5.12 -5.86 -6.61 -8.19 -
___a) 1OC SCF is equivalent to a full CI, see text. bl ER= (ED - ZEH)/(Ee - ~EH), where Ee is the energy at 0.74186 A, 2EH = 0.996466 au and ED is the energy at the distance given;Ee - EH = 0.168683 au (105.8 kcal/ mole). c) Ref. [ 171, optimized at the MP4 level.
Volume 131,number 6
CHEMICAL PHYSICS LETTERS
both (A) energy versus r(HH) and (B) In ER versus r(HH). The linearity of the In ER plot (fig. IB) extends well into the large H-H separation region, where any energy change is impossible to visually perceive (fig. 1A). The effective Morse B value in this linear region (~2.9 A-l) is about 50% higher than the Morse value [6]. This is also the case for methane and ethane. A m(energy) test on the 72-term variational ground state surface [ 181 of H2 also showed linear behavior in the 2-4 A region. We have not examined the very long distance region, dominated by r-6, r- 8, r--lo terms [7], to see where the departure from linear exponential behavior occurs. We also ob-
generally exist, the surface form of ab initio results can be anticipated and modelling schemes can be simplified [ 11. 3.2. Behavior of the H2 surface Table 1 shows the energy-distance behavior of the 1OCSCF/6-31G** H-H surface. This MC SCF level is equivalent to a full CI [ 161 (see comparative calculations at 0.74186 A shown in fig. 1). The MP4 calculation [ 171 does not represent a full CI even for so small a system. Also shown in table 1 are the reduced energies, ER, as a function of distance. Fig. 1 shows
-I
2 1 November 1986
oo-
-106. -108. -IlO2 ^ -112. z LK -114. iii W -I
l6*
-I
18.
SCF = -I 131272 MW = -1.164555 SDCI = -I 165149 IOCSCF = -I 165149
(A)
-7-6-
05
07
09
II
I5 DISTANCE,
20 H-H
25
30
35
A
Fig. 1. (A) Computed energy at the 6-31G** level for two- and ten-configuration MC SCF. The latter calculation is equivalent to a full CI as shown by reference calculations at 0.74186 A. See ref. [ 171 for MP4 calculations. (B) bi,#?R versus H-H distance. See table 1 for the energy details and method of computing ER. 477
21 November 1986
CHEMICAL PHYSICS LETTERS
Volume I3 1, number 6
served that HF and OH 2C SCF/6-3 lG* ruptures in HF and Hz0 showed linear In E, behavior in the 23.75 A region.
[ 131 projected UMP4/6-3 lG** values for the last three distances (0,27,0.088, and 0.023, respectively) are marginally different from the MRCD CI or 2C SCF values but significantly better than the unprojected values [9] (0.16,0.037, and 0.008). Obviously, what does change in all these calculations is the computed D, value for the CH bond. This value is difficult to precisely compute [ 191 and one can avoid doing so by using the experimentally known value to calculate the surface features from the theoretical reduced energies. For methane, the 2C SCF/6-3 lG** reduced energies serve equally well for scaling purposes as the MRD CI values for the CH region up to 3.5 A. Whether the progressive departure of the 2C SCF energies from those obtained at the MRD CI level above 3.0 A is of importance in dynamical modelling will require comparative computations. The 2C SCF procedure describes the first-order bonding interactions while ignoring other interfragment correlation effects. The 2C SCF and MRD CI
3.3. Behavior of the CH4 = CF3 + H surface Recently, Hirst [ 121 has published both the tetrahedral and optimized C-H bond rupture surfaces for CH, at the MRD CI/2C SCF(6-3 lG**) level. Table 2 and fig. 2 show a comparison of Hirst’s tetrahedral MRD CI computations with our 2C SCF computations. The E, values for both sets of computations are within the circles shown in the energy plot (fig. 2A) and are only perceptively different on the log scale (fig. 2B) above 3.5 A. The E, values obtained for the four points generated (O-73,0.3 1,O. 10 and 0.029 at R = 1.5,2.0,2.5, and 3.0 A, respectively) in the large basis set MRD CI calculations of Brown and Truhlar [lo] were not significantly different from either Hirst’s or our own values. Schlegel’s recent
Table 2 6-31G* * “tetrahedral surface” calculation for CH4 = CHs + H CH3-H distance
MRD CI energy a)
(A)
(au)
ERb)
In ER
-40.391645
1.00 d)
0.00
1.12 1.136
-40.390221
0.99
0.00
-40.379019 -40.338125 -40.296228 -40.262673 -40.239320 -40.224704 -40.216210 -40.211683
0.93 0.71 0.48 0.30 0.18 0.096 0.051 0.027
-0.07 -0.34 -0.73 -1.20 -1.74 -2.33 -2.98 -3.63
-40.207303 -40.206766
2.9 x 1o-3 -
-5.84 -
-40.206773
0.00 d)
20.0 100.0
ER
In ER
(au)
1.086
1.25 1.5 1.75 2.0 2.25 2.5 2.75 3.0 3.25 3.50 3.75 4.0 5.0
2C SCF energy c)
NC -40.217069
1.00 e)
-40.207665 -40.171087 -40.132646 -40.101602 -40.079687 -40.065886 -40.058043 -40.053915 -40.051837 -40.050809 -40.050162 -40.049909
0.94 0.73 0.50 0.31 0.18 0.096 0.049 0.024 0.012 5.9 x 1o-3 2.1 x 1o-3 5.4 x 10-4
40.049818
0.00 e)
0.00 -0.06 -0.32 -0.70 -1.17 -1.72 -2.34 -3.01 -3.71 -4.42 -5.13 -6.19 -7.52
a) Ref. [ 121, CH = 1.086. b, ER = (ED - E,)/(Ee - 2E,), where Ee is the energy at 1.086 A, E, the largest distance computed and ED the energy at the distance given. c, 2C SCF, HOM02/LUM02 pair. d) E 1.o86 - Eloo = 0.184872 au (116.0 kcal/mole). e, Et.12 - Ezo = 0.167525 au (105.1 kcal/moIe).
478
Volume I3 1, number 6
21 November 1986
CHEMICAL PHYSICS LETTERS
0
v
”
/d+-
a
v
-3
G E
(A)
a-8 C
6 0.5.
/
E 3
if! a
x
1.0.:l
/
/
/
,e I
/
, “, ’ ,m ’
-5 -4 3 -3 3
0 /m O’ZCSCF @ = MRCVPCSCF
_;,-G
x
x I -4-6 ,w /*
(B)
--2
@ =MRCI/2CSCF(OPTIMIZED) 0 = MP4AJHF
1.0
1.5 C-H
2.0
DISTANCE,
2.5
3.0
--I
3.5
,O 4.0
A
Fig. 2. (A) Reducedenergy curve for the tetrahedral CH rupture in methane at the 6-31G** 2C SCF and MRD CI/ZC SCF levels [ 8,9]. All values fall within the circles shown. (B) Comparison of lnER versus C-H distance. The MRCD optimized results are reported in refs. [ 8,9]. The MP4/UHF values are from the UHF/MP4 SDQ results reported in ref. [ 81.
calculations show sufficient superimposable linearity of the ln ER plots to permit a linear extrapolation into the longer-distance regions where deviation from linearity occurs at the 2C SCF level. Moreover, geometry optimization [ 121 along the reaction path has the effect of improving ln ER linearity at long distances while maintaining superimposability at shorter distances (fig. 2). It is also worth noting that the MP4 In ER curve shown in fig. 2 enables us to detect the point (2 A) where artifactual behavior begins to occur. However, the projected UMP4 In ER plot was found to be linear and nearly superimposable on the MC SCF curves. The fact that such behavior was obtained indicates to us that a major portion of the artifactual behavior was removed by the projection procedure. The presence of a linear region in the In ER plot at large r means that an effective constant Morse B value can be assumed for the region of dynamical interest. In the cases of methane or ethane (see below), these effective B values (=2.6-2.8 A-1) are larger than those used for their Morse curves (al .81.9A-l) [l]. The major question posed by the favorable com-
parison of the computations shown in table 2 and fig. 2 is whether one can generally obtain satisfactory curve forms at lower levels of theory. Our response is both guarded and optimistic. First of all, one assumes that the ln ER plots will eventually deviate from linearity at long distances. Secondly, an extended region of linearity may not even be encountered for interacting polar or polarizable fragments. Thirdly, one anticipates both a basis set and MC SCF level dependence of the effective B values. Without presenting the details of the computations, this latter effect is shown ln fig. 3 for the tetrahedral C-H bond rupture in CH4 at the 2C SCF/STO-3G and 17C SCF/ STO-3G levels and at the 17C SCF/3-21G level. All these In ER surfaces show regular behavior, the first with curvature. The effective B values in the 2.53.0 A region vary both with basis set and MC SCF level. At the STO-3G level, increasing the MC SCF level from two to seventeen configurations improves the linearity of the In ER plot and reduces the effective B value from 3.6 to 2.4 A-l. The In ER behavior of the 17C SCF/STO-3G surface is fortuitously close to the 6-31G** (B = 2.6-2.7) surfaces (fig. 2) up to 3.5 A. However, the 17C SCF/3-21G slope is 479
Volume 13 1, number 6
CHEMICAL PHYSICS LETTERS
21 November 1986
7
-6.
1.5
2.0 C-H DISTANCE,
2.5
3.0
A
3.5
Fig. 3. Comparison of lnE~ versus C-H distance in a tetrahedral scan for different basis sets and MC SCF levels. The values of De: 2C SCF/STO-3G, 131 kcal/mole; 17C SCF/STO-3G, 147 kcal/mole, and 17C SCF/3-21G, 119.8 k&/mole. R(CH) = 1.09 A.
smaller (1.9 A-l) than the 6-3 lG** values. Due to size consistency correlation artifacts, the computed bond energies (caption, fig. 3) vary with the MC SCF level. Therefore, one should not generally expect the superimposability of 2C SCF and MRD CI In ER curves that was found in fig. 2. The use of these small basis sets magnifies the problems but also illustrates some possibilities for exploitation. For small systems, one will normally use larger basis sets capable [lo]
of describing the intermediate-distance bonding regions. In this case, the necessary surface information for dynamic modelling may be extractable at the 2C SCF:level but ln ER versus r slopes should be verifled by performing calculations at the MRD CI level. If one can then find a smaller basis set, small MC SCF calculation which mimics the large basis set results at the ER level, one can risk treating a larger analogous system (e.g. R-CH2-H rupture in a larger system)
Table 3 2C SCF/3-21G optimized surface for the reaction CaHe = CHs + CHs
energy
CH, -CH3 distance (A)
2C SCF (au)
1.572 1.75 2.00 2.25 2.50 3.00 3.25 3.50 4.00 20.00
-78.808295 -78.800025 -78.711443 -78.740802 -78.716136 -78.689655 -78.683494 -78.681649 -78.679395 -78.678781
b,
ER a)
1 .oo 0.94 0.72 0.4% 0.29 0.084 0.036 0.021 4.7 x 1o-3 0.000
ERR
0.00 -0.07 -0.33 -0.74 -1.24 -2.48 -3.31 -3.87 -5.35 a)
Optimized parameters R(C-H)
L CCH
1.084 1.081 1.077 1.075 1.073 1.071 -1.071 1.071 1.070 1.073
110.4 108.3 105.7 103.3 101.0 96.9 95.0 93.1 89.4 c) 90.3 c)
a) E 1.572- E2o.o = 0.129514 au (81.3 kcal/mole). b) 2C SCF CC harmonic force constant = 3.80 mdyn A-‘. c) One expects 90” for these large distances and the values shown are probably artifacts of the optimization procedure in a region where the interaction energies are small.
480
CHEMICAL PHYSICS LETTERS
Volume 131, number 6
where neither large basis sets nor MRD CI comparisons are technically or financially feasible. 3.4. Behavior of the C,H, = 2CH3 surface Table 3 shows the optimized reaction surface for the CC rupture in ethane at the 2C SCF/3-21G level. Fig. 4 shows both the energy and In ER behavior as a function of CC distance. The In ER surface shows linear behavior in the 2.5-4.0 A region, similar to the CH rupture in methane. The 2C SCF/3-2lG level D, (table 3,81.3 kcal/mole) is reasonably close to the experimental reaction enthalpy (88 kcal/mole, less ZPE) [ 51. These computations were sufficiently inexpensive to suggest that, for dynamical applications, a 2C SCF/6-3 lG** level treatment can be envisaged. A few additional intermediate-distance MRD CI calculations will determine whether the effective B value can be obtained at the lower 2C SCF level.
2 I November 1986
We finally note that MNDO CI calculations [20] indicate that barriers exist for both CH, + H and CH, + CH, as well as for some other recombinations. Based on the results reported here and other work in progress, we view these barriers as artifacts. We are informed that the AM1 parameterization yields barrier-free surfaces for those radical-radical recombinations that have been reinvestigated so far [21].
4. Conclusions The surface forms for the bond ruptures in H,, CH4, and C2H6 show linear behavior in log(reduced energy) space for the region in the early stages of bond formation. This behavior can be used for diagnostic purposes. Comparison of the methane 2C SCF/ 6-31G** calculations with MRD CI results indicates that the reduced-energy surfaces are nearly identical. We conclude that the possibility exists of obtaining the bond-rupture surface of high-level calculations at lower levels.
Acknowledgement
(A)
CH-CH 3
3
-
CH3
+
CH3
We acknowledge a number of discussions with Professor J.J. Dannenberg. We also thank the CNRS computing center at CIRCE, Orsay, France for use of their facilities. This research project was financed in part by research contract ATP-CNRS-920042.
References
-6.
[l] C.J. Cobos and J. Troe, J. Chem. Phys. 83 (1985) 1010, and references therein.
-5’
[ 21 R.J. Duchovic and W.L. Hase, J. Chem. Phys. 82 (1985)
-4
P
1.5
20 c-c
3.0 &CE,
3.5
A
Fig. 4. Optimized 3-21G 2C SCF level CC bond rupture in ethane. Data taken from table 3. (A) Energy versus CC distanCe. (B) hI ER plot Similar t0 figs. 1-3.
3599. [3] W.L. Hase and R.J. Duchovic, J. Chem. Phys. 83 (1985) 3448, and references therein. [4] J.F. LeBlanc and P.D. Pacey, J. Chem. Phys. 83 (1985) 4511, and references therein. (5) S.W. Benson, Thermochemical kinetics (Wiley, New York, 1976) pp. 86-100,160-166. (61 H.S. Johnston, Gas phase reaction rate theory (Ronald Press, New York, 1966) chs. 4,11. [7] W. Koiix, Intern. J. Quantum Chem. 1 (1967) 169. (81 R.J. Duchovic, W.L. Hase, H.B. Schlegel, M.J. Frisch and K. Raghavachari, Chem. Phys. Letters 89 (1982) 120.
481
Volume I3 1, number 6
CHEMICAL PHYSICS LETTERS
19) R.J. Duchovic, W.L. Hase and H.B. Schlegel, J. Phys. Chem. 88 (1984) 1347. [ 101 F.B. Brown and D.G. TruhIsr, Chem. Phys. Letters 113 (1985) 441. [ 111 C.J. Cobos and J. Troe, Chem. Phys. Letters 113 (1985) 419. [12] D.M. Hirst, Chem. Phys. Letters 122 (1985) 225. [13] H.B. Schlegel, J. Chem. Phys. 84 (1986) 4530. [14] N.C. Handy, Faraday Symp. Chem. Sot. 19 (1984) [15] :A. Peterson and R.A. Poirier, MONSTERGAUSS, Department of Chemistry, University of Toronto.
482
2 1November 1986
1161A.C. Wahl and G. Das, in: Methods of electronic strut-
ture theory, Vol. 3, ed. H.F. Schaefer III (Plenum Press,
New York, 1977) pp. 51-78, and references therein. (171 M.J. Frisch, J.A. Pople and J.S. Binkley, J. Chem. Phys. 80 (1984) 3265. [181 W. Kolos and L. WoIniewicz, J. Mol. Spectry. 54 (1975) 303. 1191P.E.M. Siegbahn,Chem. Phys. Letters 119 (1985) 515. [201J.J. Dannenberg and K. Tanaka, J. Am. Chem. Sot. 107 (1985) 671; J.J. Darmenberg, J.C. Rayez, M.T. Rayez-Meaume and P. Halvich, J. Mol. Struct. THEOCHEM 123 (1985) 343. [211 J.J. Dannenberg, private communication.