J. Phys. Chem. 1995,99, 8613-8616
8613
Reliability of the Single-Point Calculation Technique at Characteristic Points of the Potential Energy Surface J. Espinosa-Garcia” and J. C. Corchado Departamento de Quimica Fisica, Universidad de Extremadura, 06071 Badajoz, Spain Received: December 20, 1994; In Final Form: March 17, 199.5@
The errors associated with the application of the single-point calculation technique in the construction of the reaction path were analyzed in terms of both energy and entropy. Since the lower and higher levels yield curves that are shifted relative to each other, this technique has a destabilizing effect at the minimum (E’ error), while, at the saddle point, it has a stabilizing effect on the reaction coordinate (E” error) and a destabilizing effect on the remaining coordinates (Ekperror). Likewise, along the minimum energy path, this technique artificially locates the maximum away from the saddle point. To regard this shift away from the saddle point as a variational effect would be a conceptual error. A modification of the technique is proposed that avoids these computational and conceptual errors.
I. Introduction The accuracy of the dynamical description of a chemical reaction depends on the quality of the information about the potential energy surface (PES). Unfortunately, the complete construction of the PES is even today limited to a very few molecular systems, such as atom-atom or atom-diatom, and it becomes prohibitive as the number of degrees of freedom rises. For polyatomic systems, one has, therefore, to introduce approximations in order to calculate the rate constants. The conditions that will have to be satisfied by these approximations are (1) simplicity, so that only partial knowledge of the PES is necessary; (2) computational economy, so that it will be possible to use large basis sets and high correlation levels; (3) reproducibility of the results, to give an adequate approximation to the best theoretical and experimental values; and (4)conceptually free of error. An interesting and economical alternative to the more costly construction of an analytical PES is to compute the minimum energy path (MEP), which is the steepest descent from the saddle point toward reactants and products. This methodology requires the evaluation of the energy, gradient, and Hessian matrix at several points along the MEP, where a generalized normal-mode analysis is performed projecting out the path direction as well as over all rotations.’ A priori, this methodology satisfies the above four conditions. However, the energy and the energy derivative calculations are very costly if the level used (correlation energy and basis set) is to give sufficient accuracy to reproduce the most important feature on the reaction path. A commonly used technique2 to avoid this high cost is the single-point calculation which optimizes the geometry at one level of calculation (level A) and then makes a single-point calculation of energies at a higher level (level B) (Le., without optimization). We shall denote the resulting energy as level B//level A (B//A), where the left part (before the I/) refers to the electronic energy singlepoint calculation and the right part to the geometry optimization. The advantages of being useful and easy to use make this a very popular technique. However, it presents certain conceptual and computational problems which have not been rigorously analyzed. Thus, the standard single-point calculation does not
* Author to whom correspondence should be addressed. @
Abstract published in Advance ACS Abstracts, May 1, 1995.
satisfy condition 4. We shall here analyze these problems from the two points of view of energy and entropy. Firstly, in terms of energy, while this technique only poses a computational problem at the stationary points, in the MEP construction, it can artificially locate the energy maximum away from the saddle point. The single-point calculation technique thus presents a conceptual problem: what is simply a numerical defect can be mistaken for a variational effect. Secondly, we shall analyze the vibrational (entropy) effects. The frequencies obtained at the lower level A are added at each point along the MEP at the level B//A to construct the vibrationally adiabatic ground-state potential curve
where VMEPis the Bom-Oppenheimer potential along the MEP at the B//A level, tint is the local zero-point energy at level A, and s, the reaction coordinate, is defined as the signed distance from the saddle point along the MEP, with s > 0 referring to the product side. This technique also presents a conceptual problem: at the new energy maximum (B//A level), the frequencies calculated at the lower level A do not define a true saddle point. Obviously, these problems disappear if we construct the complete MEP with a sufficiently correlated level of calculation and a very large basis set. While such an ideal methodology would require no additional approximations, it would at present still be prohibitively costly for large polyatomic systems. The recent “direct ab initio” methodology of Truong3 is an advance in this sense, but even when calculating the reaction path at the quadratic configuration interaction with single and double excitations (QCISD) level of theory using the triple-5 plus polarization 6-31 lG(d,p) basis set, the potential energy had to be scaled along the MEP by a factor to match the barrier height at a higher level of calculation. In the present work, we analyze the possible sources of error in the MEP construction using the single-point calculation technique, and we propose a modification of the technique to avoid such errors. 11. Analysis of the Variations of Energy
As here we are only interested in the problems involved in the practical application of the MEP construction and not in
0022-3654195/2099-8613$09.00/0 0 1995 American Chemical Society
Espinosa-Garcia and Corchado
8614 J. Phys. Chem., Vol. 99, No. 21, 1995
TABLE 1: Geometric Parametere and Total Energiesbfor the NH3 Molecule d(N-H) L(HNH) E E' f
level A//A'
level BIBd
1.032 104.2 -55.455 42
1.012 106.1 -56.386 92
level B//Ae
-56.385 89 0.65
Distances in angstroms and angles in degrees. The experimental values are 1.012, 106.7 (ref 8). In hartrees. RHF/STO-3G//RHF/STO3G. MP2/6-31G(d,p)//MP2/6-3lG(d,p). e MF'2/6-3lG(d,p)//RHF/STO3G. f E ' = E(B//A) - E(B//B), in kcal mol-'. E E'.
Figure 1. Hypothetical two-dimensional surface at two levels of calculations: A//A and B//B. obtaining accurate results, we use as level A one of the simplest ab initio levels, (R or U)HF/STO-SG, and as level B we use electronic correlation and a larger basis set, (R or U)MP2/631G(d,p). The acronyms are standard. The tunneling effect is likewise not analyzed. Although the ideas and the concepts are general, for the sake of clarity we shall consider a numerical example. As test reaction we use NH3 H TS NH2 H2, which we have studied thoroughly at higher levels of calculation: geometric optimization at (R or U)MP2/6-3 1G(d,p) and correlation energy using the MP4, post-MP4, and QCI methods with a larger basis set.4 (a) Minima on the MEP. It is well-known that the geometric optimization is computationally more costly than the energy calculation. Thus, in several fields of chemistry, it is commonplace to make these optimizations at a given level (level A) and then to make a single-point calculation (without optimization) at a higher level (level B). We denote this technique as B//A. However, it presents the computational problem that we show schematically for a hypothetical potential curve (two-dimensional surface) in Figure 1. If we optimize at level A, we obtain a value of the equilibrium distance r (rtJ which minimizes the energy at the A//A level but not at the B//B level. Thus, the energy at the B//A level is slightly higher than the B//B energy. We shall designate this error as E', being the difference between the B level energies as calculated using the geometries of each level's minimum:
+
E' = E(B//A)
- -
+
- E(B//B)
This E' error has always a destabilizing effect. Obviously, the greater the difference between the A and B levels, the greater is E'. As one is usually interested in changes of energy, there will be compensation, or at least reduction, of this error, as will be the case with the heat of reaction, for example. In Table 1 we list the geometric parameters, the total energies, and the error E' for the NH3 molecule. The small value of E' shows that the single-point calculation technique (B//A) is adequate for the minimum. In view of the great differences between the two levels of calculation (B//B and A//A), this agreement is encouraging. (b) Saddle Point. In the case of a saddle point the situation is completely different. Let us recall that for a surface of k
A
c
B rap
r
Figure 2. Reaction coordinate at two levels of calculation: A//A and B//B.
dimensions this stationary point is a minimum in k - 1 directions of the space and a maximum in the remaining direction. In the k - 1 directions which are minima, the behavior of the single-point calculation technique is similar to the minimum on the MEP (see above the Figure 1). Thus, we obtain an Ekp error, which has also a destabilizing effect. Let us now look at the relationship between Ekp and the E' error in the reactants, E'R. In a polyatomic reaction, for the reactants we have (3N1 - 6) (3N2 - 6) intemal coordinates, while the saddle point is defined by ~ ( N I N2) - 7 (minima); that is, the saddle point presents five coordinates more. Therefore, Ekp > E'R usually, and the compensation of one error by another in the changes of energy (barrier height) is not as clear as before. Fortunately, these five coordinates are normally associated with low frequencies, Le., flat curves, and their contribution to the ESP error will be small. Thus, for the sake of simplicity, we shall assume in the rest of the paper an effective compensation, Ekp E'R, so that the barrier height is free of E' errors. However, in the remaining direction (reaction coordinate), whenever the A//A and B//B curves are shifted relative to each other, as is the usual case, the single-point calculation (BNA) gives mathematically a more stable energy than the B//B calculation and, therefore, a lower barrier height. We shall designate this error for the reaction coordinate as E". Figure 2 shows this behavior schematically. The E" error is independent of the position (above/below) of the A//A curve relative to the BIB, depending only on the relative shift of the two curves. In general, as the geometric parameters that are most closely related to the reaction are the most sensitive to the level of optimization, it seems logical to suppose that E" > ESP. Obviously, this separation of the reaction coordinate (E" error)
+
+
Single-Point Calculation Technique
Ek R
J. Phys. Chem., Vol. 99, No. 21, 1995 8615
I
Figure 3. Plot of the expected barrier at the B / B level (AI$) and the E’ and E” errors for a hypothetical reaction.
from the remaining coordinates (Eip error) is impossible in practical cases. The difference between the energies of the respective saddle points (s = 0), E(B//A) - E(B//B), will therefore include both errors. In kinetic studies, the barrier height plays a central role. In Figure 3, we plot this barrier and the E‘ and E” errors simultaneously for a hypothetical reaction. If we denote as AI$ the expected barrier (B//B level) and A,?$ the obtained barrier (B//A level), we get the following relationship:
@ = AI$
If we assume that E’, = ELp, the single-point calculation technique (B//A) yields a lower barrier height than the expected barrier (B//B) due to the E” error. However, as will be seen below, this E“ error can be removed, although care must be taken to avoid conceptual errors. (c) Minimum Energy Path. The foregoing would be the situation if we only analyzed the stationary points (reactants, products, and saddle point). However, in kinetic studies more information is required along the MEP, which is obtained following an intrinsic reaction coordinate (IRC).5 Normally, the IRC is obtained at a low level of calculation (A//A level), and the single-point calculation technique (B//A) is applied along these geometries. In Figure 4 we plot the IRC results, and in Table 2 we show the geometric parameters and changes of energy along the MEP for the NH3 H NH2 H2 reaction, our test reaction, at three levels of calculation: A//A, B//B, and B//A. Two effects are evident in the B//A technique. First, the maximum of the new curve (B//A) (s = 0.4, Table 2) is artificially shifted with respect to the maximum of the original curve (A//A) (s = 0). Note that now this new maximum is not a true saddle point. Second, this shifted new maximum (B//A) is in good quantitative agreement with the somewhat smaller B//B barrier height. If we take the energy of this new maximum on the reaction coordinate as the energy of the saddle point, we avoid the E” error (note that the ELp error is still present). This situation, however, poses a conceptual question: Should this shift in the single-point calculation be regarded as a variational effect? It seems clear that the answer should be no. Using this technique, one artificially shifts the maximum, and one would be regarding as a variational effect what is only a simple numerical error. We will analyze this effect again, when we introduce the frequencies to calculate the vibrationally adiabatic ground-state curve. In Table 2, we quantitatively analyze the situation. First, the maximum of the B//A curve appears at s = 0.40 amu”2bohr and is, therefore, shifted with respect to the original curve, A// A. The difference between the energies E(B//A) - E(B//B) at
-
rEX,,A
r
(S=0.4)A
Figure 4. Minimum energy path for the test reaction NH3 NH2
+ Hz, at three levels of calculation:
+
H A//A, BIB, and B//A.
-
TABLE 2: Geometric Parameters and Changes in Energy (kcal mol-’) along the MEP
+ Ekp - E“ - E’,
+
4 !*
( S d (S=O?
+
level N / A b
sa -0.20 -0.10 0.0 0.10 0.20 0.30 0.40
level B I B d
level B//Ac
rIe
rd
AE
1.082 1.124 1.161 1.200 1.249 1.301 1.355
1.141 1.067 1.009 0.953 0.885 0.821 0.766
12.30 13.98 14.58 13.80 10.68 6.12 1.86
AE
rl
1-2
18.18 20.39 22.20 22.99 23.34
1.202 1.255 1.305 1.351 1.404 1.453 1.500
0.978 0.919 0.869 0.832
AE
20.08 21.90 22.40 22.10 0.800 21.30 0.780 20.30 0.768 19.24
Reaction coordinate in amut’2bohr; s > 0 for products. (R or U)HF/STO-3G//(R or U)HF/STO-3G. (R or U)MP2/6-3 lG(d,p)//(R or U)HF/STO-3G. (R or U)MP2/6-31G(d,p)//(R or U)MP2/631G(d,p). e N-H distance, in angstroms. f H-H distance, in angstroms.
their respective saddle points (-56.855 16 and -56.849 45 hartrees, respectively) includes the E&, and E” errors. However, if we take the maximum energy in the curve B//A (s = 0.4) (-56.846 92), we eliminate the E” error, so that now the difference E(B//A)max- E(B//B),=o will only contain the ELp error, which is 1.59 kcal mol-’. Therefore, E“ will be E(B// A),=o - E(B//B),=o - ELp = -5.17 kcal mol-’, the same as E(B//A)+J- E(B//A)max. As one can see, this E” error leads to an underestimate of the barrier height and, therefore, an overestimate of the rate constant, obviously related to the difference between the A and B levels. As we foresaw, Eip in the saddle-point region is greater than E‘R in the reactant region (1.59/0.65), although they are both considerably compensated, given the great differences between A and B levels. In the practical application of the single-point calculation along the MEP, we propose therefore to take as the barrier height that calculated using the maximum energy of the new curve B//A, avoiding the E” error, and to move this maximum to its original position (s = 0) at the A//A level, thus avoiding taking as a variational effect what is only a numerical defect. This maximum is now only affected by the Ekp error, which is almost compensated by the E‘R error from the reactants. In conclusion, from the point of view of energy exclusively, the
Espinosa-Garcia and Corchado
8616 J. Phys. Chem., Vol. 99, No. 21, 1995 TABLE 3: Frequencies (cm-') and ZPE (kcal mol-') at Two Levels of Calculation level N/Aa
s=o -3014 764 786 1497 1851 1856 1993 3843 4027 ZPE
s = 0.4
level B I B b
769 788 1193 1327 1907 3440 3795 400 1
-2021 719 73 1 1146 1315 1610 2077 3489 3605 2 1.oo
24.62
23.76
UHFISTO-3GIIUHFISTO-3G. 31G(d,p).
UMP2/6-3 lG(d,p)//UMP2/6-
TABLE 4: Vibrationally Adiabatic Ground-State Potential, If(kcal mol-') level B//Ab AI? AZPE
t
level B/B"
approx 1'
approx 2d
approx 3'
22.40 -1.28 21.12
18.18 - 1.42 16.76
23.34 -0.56 22.78
23.34 - 1.42 21.92
a (R or U)MFW6-31G(d,p)//(R or U)MP2/6-31G(d,p). (R or U)MP2/ 6-31G(d,p)//(R or U)HF/STO-3G. Energy (s = 0) frequencies (s = 0) (see text). Energy (s = 0.4) frequencies (s = 0.4). e Energy (s = 0.4) frequencies (s = 0).
+
+
+
saddle point always represents the maximum of the MEP, and localizations far from the saddle point are due to physically meaningless mathematical artifacts. To regard this shift away from the saddle point as a variational effect would be an error.
111. Analysis of the Frequencies and the Vibrationally Adiabatic Potential In the variational transition state theory, the dividing surface is varied along the reaction path to minimize the rate constant. Thus, we obtain the generalized transition state (GTS) at the value s * . ~ Thermodynamically, the minimum rate constant criterion is equivalent to maximizing the generalized standardstate free energy of activation, AG"(T,s).~Therefore, from a thermodynamic point of view, one must consider the effects of th.e potential energy (analyzed in section 11), zero-point energy, and entropy on the location of this generalized transition state. The theoretical calculation of the frequencies requires the evaluation of the Hessian matrix at several points along the MEP, where a generalized normal-mode analysis is performed projecting out the path direction as well as over all rotations' and is therefore computationally costly. Usually, these frequencies are calculated at the same level at which the geometries are optimized, Le., the MIA level. Thus, the stationary points are true stationary points. In Table 3 we list the frequencies for our test reaction at several points along the MEP at two levels of calculation. Note that at this NIA level the point s = 0.4, which corresponds to the maximum of the curve BNA, is not a true saddle point. The frequencies at the B//B level correlate better with the frequencies at the A//A level at the point s = 0 (in general, as is ~ e l l - k n o w n.the , ~ frequencies at the AJIA level are overestimated with respct to the higher BID3 level). Finally, in Table 4 we show the vibrationally adiabatic ground-state potential, for the saddle point, with three approximations. In approximation 1, we take the B//A energies and the N I A frequencies, both calculated using the s = 0 geometry for the saddle point. This approximation underesti-
e,
mates the barrier height and due to the E" error, and only satisfies conditions 1 and 2. In approximation 2 , we use for the saddle point the BIIA energy at its maximum, located at s = 0.4, and the projected frequencies are also calculted using the geometry of this point. Taking the maximum BIIA energy along the MEP we avoid the E" error and we achieve quantitative agreement for the barrier height, but is overestimated due to the projected frequencies at the point s = 0.4 (see Table 3), which is not a true saddle point at the NIA level. which is AGO at T = 0 K, means that This overestimate of AGO is also overestimated and, therefore, erroneously minimizes the rate constant. This approximation only satisfies conditions 2 and 3. Finally, in approximation 3, we take the BIIA energy at its maximum (s = 0.4), avoiding the E" error, and we move it to the original maximum, s = 0, avoiding taking as a variational effect what is a calculational defect. The frequencies are now calculated at the true saddle point at the N I A level (s = 0). This approximation leads to the best agreement (for the barrier height and for with the higher BID3 level. Thus, approximation 3 satisfies the four conditions.
e,
Conclusions The partial knowledge of the PES through the MEP construction is a useful methodology in predicting rate constants directly from ab initio electronic calculations, without the intermediate step of a potential surface fit, but its practical application presents computational and conceptual problems if the singlepoint calculation technique is used, as is the habitual procedure. This technique consists in optimizing geometries at one level (A) and then calculating the energies (without optimization) at a higher level (B), BIIA. At the stationary points, this technique presents two computational problems related to the shifting of the original curves, MIA and BIB: the E' error at the minimum, usually almost cancelled out when changes of energy are being considered (heat of reaction and barrier height), and the F' error at the saddle point, which leads to underestimates of the barrier height. Along the MEP, this technique presents a conceptual problem: it artificially locates the energy maximum away from the saddle point. In the present communication, we have proposed a modification that avoids these problems. We take the maximum of the single-pointcalculation curve, B//A, thereby avoiding the E" error, and we move this maximum to its original position (s = 0 for the NIA level), thus avoiding taking as a variational effect what is only a numerical defect.
References and Notes (1) Miller, W. H.; Handy, N. C.; Adams, J. E. J. Chem. Phys. 1980, 72, 99. (2) Some recent examples: Koseki, S.; Gordon, M. S. J. Phys. Chem. 1989, 93, 118. Garrett, B. C.; Koszykowski, M. L.; Melius, C. F.; Page, M. J. Phys. Chem. 1990, 94, 7096. Durant, J. L.; Rohlfing, C. M. 1.Chem. Phys. 1993, 98, 8031. Dobbs, K. D.; Dixon, D. A,; Komomicki, A. J . Chem. Phys. 1993, 98, 8852. Dobbs, K. D.; Dixon, D. A,; J . Phys. Chem. 1994, 98, 5290. Fang, D.-C.; Fu, X.-Y.; lnt. J . Quantum Chem. 1994, 49, 3. Sicilia, E.; Di Maio, F. P.; Russo, N. Chem. Phys. Letr. 1994,225, 208. Truong, T. N.; Duncam, W. J . Chem. Phys. 1994, 101,7408. Finally, our own research group has used this technique in recent years. (3) Truong, T. N. J . Chem. Phys. 1994, 100,8017. (4) (a) Espinosa-Garcia, J.; Tolosa, S.; Corchado, J. C. J . Phys. Chem. 1994, 98, 2337. (b) Espinosa-Garcia, J.; Corchado, J. C. J . Chem. Phys. 1994, 101, 1333. ( 5 ) Fukui, K. J. Phys. Chem. 1970, 74,4161; Pure Appl. Chem. 1982, 54, 1825. (6) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. In The Theory of Chemical Reactions; Bear, M., Ed.; CRC: Boca Raton, FL, 1985; Vol. 4. ( 7 ) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (8) Chase, M . W., Jr.; Davies, C. A,; Downey, J. R., Jr.; Fmrip, D. J.; McDonald, R. A,; Syverud, A. N. JANAF Thermochemical Tables. J . Phys. Chem. Ref Data Suppl. 1985, 14. JP943359G
