- Email: [email protected]
Reduced-dimensionality quantum reactive scattering calculations of the C(3P) + C2H2 reaction on a new potential energy surface
Chemical Physics 312 (2005) 61–67 www.elsevier.com/locate/chemphys
Reduced-dimensionality quantum reactive scattering calculations of the C(3P) + C2H2 reaction on a new potential energy surface Toshiyuki Takayanagi
*
Department of Chemistry, Saitama University, 255 Shimo-Okubo, Sakura-ku, Saitama City, Saitama 338-8570, Japan Received 28 September 2004; accepted 16 November 2004 Available online 8 December 2004
Abstract The reaction of the ground-state carbon atom with acetylene has been theoretically studied by a reduced dimensionality method, where three degrees of freedom were taken into account in the quantum dynamics calculations. A three-dimensional potential energy surface has been developed on the basis of the B3LYP/6-31G(d,p) level electronic structure calculations. The total reactive cross sections have been calculated by the time-independent generalized R-matrix method with the use of the negative imaginary potential. The reactive cross sections obtained in the energy range between 0.1 and 20 kJ/mol show a qualitative trend that is consistent with the barrierless reaction but the energy dependence was found to be much weaker than the experimental results. It has also been found that excitation of the initial rotation does not affect the reactivity but excitation of the CC stretch significantly decreases the reaction cross section. 2004 Elsevier B.V. All rights reserved.
1. Introduction The reaction of the ground-state atomic carbon C(3P) with acetylene has recently attracted considerable interest due to its fundamental importance in the chemistry of dense interstellar clouds [1,2]. In particular, this reaction is thought as a key step for the synthesis of larger carbon-containing species since the overall process of the reaction of C(3P) with unsaturated hydrocarbons may be represented by C + CnHm ! Cn+1Hm1 + H. The importance of the reaction of C(3P) with acetylene in interstellar clouds has been pointed out by the rate constant measurements. The rate constant at room temperature for this reaction was reported to be very large, suggesting no barrier in the potential energy surface of the C(3P) + C2H2 reaction [3–5]. The kinetic measurement was then extended down to the temperature range of 300–15 K [6,7] and it has been found that *
Tel.: 48 858 9113; fax: 48 858 3700. E-mail address: [email protected].
0301-0104/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.11.026
the rate constant has a negative temperature dependence and that the absolute value of the rate constant is very large (>2 · 1010 cm3 molecule1 s1) even at 15 K. This is quite important from an astrophysical point of view since, in interstellar clouds, the density is extremely low and a typical temperature is in the range of 10–50 K. Although the previous kinetic studies mentioned above played a decisive role for understanding of the importance of the carbon–acetylene reaction in interstellar clouds, those studies cannot identify the products of the reaction. According to the ab initio electronic structure calculations[8–10], it is energetically possible that the C(3P) + C2H2 reaction leads to the production of two C3H isomers: cyclic-C3H and linear-C3H. Kaiser and coworkers [8–11] have performed an extensive series of crossed molecular beam experiments at relatively high collision energies with mass spectroscopic detection of the C3H molecule. Although the mass spectroscopic technique cannot distinguish these isomers in principle, they have concluded that the cyclic-C3H molecule is mainly produced with a minor contribution of the
62
T. Takayanagi / Chemical Physics 312 (2005) 61–67
linear-C3H product at higher collision energies from the obtained angular distributions. Very recently, the crossed molecular beam experiment [12] have also revealed that the C(3P) + C2H2 reaction leads to C3 ð1 Rþ g Þ þ H2 products through the spin-forbidden intersystem-crossing mechanism. Ab initio electronic structure calculations are usually quite useful since they can provide information about the energetic accessibility of reaction products for which there is no thermochemical data and about the potential energy barriers along the reaction pathways from reactants to possible products. In fact, extensive ab initio calculations for the C(3P) + C2H2 reaction at various levels of theory have been carried out so far [8–11,13– 16]. In addition, Buonomo and Clary [16] have recently reported a quantum dynamics study for the first time. They have treated the C(3P) + C2H2 reaction with a reduced dimensionality model, in which two degrees of freedom, the distance R between C(3P) and the center of mass of C2H2 and the angle R makes with the C2H2 bond axis, are taken into account in the dynamics. They have developed a two-dimensional potential energy surface by the ab initio CCSD(T) method with the other seven degrees of freedom being optimized. They have then performed two-dimensional time-dependent wave packet calculations to obtain reaction probabilities. Although the wave packet calculations were carried out only in the entrance region of the potential energy surface, they were able to obtain the information about the cyclic-C3H/linear-C3H product branching by changing the location of the flux absorbing potentials. They have found that linear-C3H is formed preferentially over the energy range considered and that the cyclic-C3H product contributes at higher collision energies. Clearly, this result contradicts an interpretation of the crossed molecular beam experiments [8–10]. The purpose of the present study is to further elucidate the dynamics of the C(3P) + C2H2 reaction from a theoretical point of view. We extend the reduceddimensional model of Buonomo and Clary [16] so as to include one more degree of freedom; three degrees of freedom are taken into account in our quantum dynamics model. We mainly focus on the collision energy dependence and on the initial state dependence of total reaction cross sections.
2. Computational procedure 2.1. Construction of the potential energy surface In this study, the potential energy surface for the C(3P) + C2H2 addition reaction has been calculated as a function of three coordinates, R, r and h, describing the distance between C(3P) and the midpoint of CC in C2H2, the CC internuclear distance, and the orientation
angle of the approaching carbon atom with respect to the CC axis, respectively. This means that the C(3P) + C2H2 reaction system is treated as a standard atom–diatom system. The other six degrees of freedom have been optimized with respect to the total electronic energy. Our method is a natural extension of the theoretical work of Buonomo and Clary [16], who have carried out the reduced dimensionality quantum dynamics calculations considering two degrees of freedom as active coordinates. They have excluded the C–C internuclear coordinate from their dynamics calculations. However, our three-dimensional model is rather realistic in the sense of the Born–Oppenheimer-like approximation where fast hydrogen atom motions were assumed to be adiabatic. The potential energy values were calculated using the hybrid density functional B3LYP method with the 631G(d,p) basis set. A total of about 2600 points were ˚, calculated in the ranges of 1.0 6 R 6 5.0 A ˚ 1.0 6 r 6 2.2 A, and 0 6 h 6 90 using Gaussian 98 set of programs [17]. In an early stage of our investigation, we have tried to calculate all three-dimensional grid points in order to apply the three-dimensional interpolation method; however, we have encountered a self-consistent field convergence problem at many geometrical configurations. In addition, in the region where r is large, we often had difficulty in obtaining smooth potential energy values. From these reasons, we have fitted the calculated potential energies to the following many-body type analytical function: V ðR1 ,R2 ,R3 Þ ¼ V 2 ðR1 Þ þ
M X i,j,k
ci,j,k qi1 qj2 qk3 ,
with the constrains i + j + k 6 M. Here, Ri is the C–C diatomic internuclear distance. R1 is the C–C distance in the reactant acetylene molecule and the first term V2 is the two-body term, for which we have employed the standard Morse potential. The second term is the three-body term proposed by Aguado and Paniagua [18]. Notice that the variable qi ¼ Ri eai Ri vanishes for Ri = 0 and Ri ! 1. In our fit, we chose M = 7. Also, we have taken the permutation symmetry with respect the Jacobi angle, h, into account. All the parameters in the above analytical function were determined by the non-linear least-square method. The resulting three-dimensional potential energy surface is plotted in Fig. 1. For both h = 0 and h = 30 approaches of C(3P) toward the acetylene molecule, it is seen that there exist saddle points in the entrance region. On the other hand, no barrier is seen for both h = 60 and h = 90 approaches. In order to understand the reaction mechanism more clearly, the two-dimensional potential energy surface is plotted in Fig. 2 as a function of R and h, where the coordinate r was optimized with respect to the energy. An overall feature of our potential
T. Takayanagi / Chemical Physics 312 (2005) 61–67
63
energy surface is quite similar to the two-dimensional surface of Buonomo and Clary [16]. For a T-shaped approach (h = 90), we can see a shallow well around ˚ corresponding to the cyclic-C3H2 intermediR = 1.1 A ate. The path then reaches the linear-C3H2 intermediate ˚ . Two after passing the saddle point around R 0.95 A ˚ potential wells are also seen around R 2 A for both h = 0 and h = 180, which also leads to the linearC3H2 intermediate formation. In addition, barriers are ˚ , for small seen in the region between R = 2.5 4.2 A angles h. At this point, it is informative to compare the energy levels of intermediates obtained at the B3LYP/631G(d,p) level of theory to the result obtained at a more accurate level of theory. Cyclic-C3H2 (cyclopropenylidene) has an energy of 262 kJ/mol with respect to the reactants while linear-C3H2 (propargylene) has an energy of 432 kJ/mol at the B3LYP/6-31G(d,p) level. These values are somewhat smaller than the corresponding CCSD(T)/QZ2P level values [10], 225 and 386 kJ/mol, respectively. 2.2. Quantum dynamics calculations
Fig. 1. Contour plots of the potential energy surface for the C(3P) + C2H2 reaction as a function of Jacobi coordinates, R, r and h. Zero of the energy is taken to be the reactant C(3P) + C2H2 potential minimum. Contours are spaced by 20 kJ/mol; solid curves are used for energies that are positive relative to C(3P) + C2H2, dashed curves are used for energies which are negative, and bold lines denote the zero contours.
We have carried out quantum reactive scattering calculations by the generalized R-matrix propagation method [19,20], which has been developed by Aguilar and coworkers [21–23] to solve the time-independent close-coupling equations with the use of a complex absorbing potential. If the absorbing potential is properly located in the product region and efficiently absorbs the reactive flux, we can easily obtain stateto-all reaction probabilities without the information of the product side. The great advantage of this method is that it can treat a reactive scattering problem in a formulation very similar to an inelastic scattering problem. As mentioned previously, the reaction system has been treated as a standard atom–diatom system. The Hamiltonian is then given by H ¼
Fig. 2. Contour plots of the potential energy surface for the C(3P) + C2H2 reaction as a function of R and h. Contours are spaced by 20 kJ/mol.
h2 o2 h2 o2 ðJ jÞ2 h2 j2 h2 þ þ þ V ðR,r,hÞ, 2lR oR2 2lr or2 2lr r2 2lR R2
where lR is the reduced mass between C and C2H2, J is the total angular momentum operator, j is the rotational angular momentum operator of C2H2, and lr is the reduced mass of (CH)–(CH). The calculations for the total angular momentum quantum number J > 0 were carried out using the standard centrifugal sudden (CS) approximation [24]; The operator term, (J j)2, was replaced into [J(J + 1) 2X2 + j2], where X is the projection on the body-fixed Z axis. Thus, Coriolis terms that couple states with different X have been neglected. This means that the projection quantum number X is assumed to be conserved.
64
T. Takayanagi / Chemical Physics 312 (2005) 61–67
The total scattering wave function can be expanded using the product of vibrational–rotational basis sets and usual close-coupling equations can be obtained although the matrix elements are complex. A twodimensional eigenvalue problem at a given value of the Jacobi scattering coordinate R was solved using the discrete-variable-representation (DVR) method [25]. We employed the particle-in-box basis for the C–C stretch coordinate r and associated Legendre basis for the angular coordinate h. For the DVR calculations, 200 basis sets in r and 100 basis sets in h were employed. The calculations can be decoupled into subprograms by exploiting the permutation symmetry of the two identical C atoms. The scattering coordinate R was first divided into ˚ . We have em600 sectors in a range of R = 1.0–5.0 A ployed the linear Neuhauser–Baer [26,27] negative imag˚ . The inary potential, located in the range of 1 < R < 2 A absorbing strength parameter was determined so as that the imaginary potential effectively absorbs the reactive flux. Note that the present reactive probability includes both the cyclic-C3H and linear-C3H production channels unlike the two-dimensional wave packet calculations of Buonomo and Clary [16]. Calculations have been carried out in the range of total angular momentum between 0 and 150, with a step DJ being 1–2, depending on the total energy. The maximum value of the projection quantum number X was taken to be 4. Thus, we were able to calculate the cross sections for the C(3P) + C2H2(j) reaction up to j = 4.
3. Results and discussion The C(3P) + C2H2(v = 0, j = 0 20) reaction probabilities for the total angular momentum J = 0 are shown in Fig. 3 as a function of the total energy. Here, v is the vibrational quantum number of the CC stretch. The total energy is measured from the asymptotic energy level of the C(3P) + C2H2(v = j = 0) channel. Although many sharp resonances are seen, the reaction probabilities show a sudden increase just after the reaction threshold and can approximately be represented by a step function. This result is quite typical for barrierless reactions. Also, it is seen that values of average probabilities are quite large and range from 0.8 to 1. The reaction probability for C + C2H2(j = 0) is almost unity above the reaction threshold. This is highly in contrast with the wave packet result of Buonomo and Clary [16]; their calculation shows that the reaction probabilities for j = 0 range from 0 to 0.8 and are somewhat smaller than our result. Plots of cumulative reaction probabilities, which correspond to total reaction probabilities summed over all initial states, are given for different numbers of J(X = 0) in Fig. 4. It is seen that the cumulative reaction probability for J = 0 shows less resonance structures
Fig. 3. Reaction probabilities for the C(3P) + C2H2(v = 0, j) reaction for the zero total angular momentum (J = 0) as a function of the total energy. The zero of the total energy is measured from the asymptotic C(3P) + C2H2(v = j = 0) level.
Fig. 4. Plots of cumulative reaction probabilities for the C(3P) + C2H2 reaction for different number of J as a function of the total energy. The zero of the total energy is measured from the asymptotic C(3P) + C2H2(v = j = 0) level.
compared to the initial state-selected reaction probabilities presented in Fig. 3. This indicates that the increases and decreases in the initial state-selected reaction prob-
T. Takayanagi / Chemical Physics 312 (2005) 61–67
abilities at resonance energies are compensated each other in the cumulate reaction probability. Although resonances become less important as the value of J increases, a typical J-shifting behavior is clearly seen; i.e., the reaction threshold energy gradually increases as J increases. The number of open channels is also plotted in Fig. 4 as a dotted line. It is interesting to note that the sudden step feature can be seen only for the cumulative reaction probability for J = 0 at low energies. The total reaction cross sections obtained in this study are presented in Fig. 5. It is found that the cross sections for v = 0 and j = 0, 2, and 4 give almost identical results in the energy range considered. This behavior is somewhat different from the result of adiabatic capture centrifugal sudden approximation (ACCSA) calculations by Buonomo and Clary [16]. The ACCSA calculations show that the initial state selected reaction rate constant decreases as the initial rotational state increases, although their ACCSA calculations have been done up to j = 12. Since our quantum scattering calculations have been done for Xmax = 4, we could not obtain the converged cross section for such a large value of j. Also, unfortunately, Buonomo and Clary did not perform the wave packet calculations for j > 0. It should be important to compare the present cross sections with the 2D wave packet result by Buonomo and Clary [16]. Their cross section for j = 0 is also plotted in Fig. 5 as a bold line. It can be seen that our cross section is larger than the 2D result by a factor of about 2–2.5. It is suggested that the source of this disagreement mainly comes from the difference between the potential energy surfaces employed. Buonomo and Clary [16] have developed their surface on the basis of the ab initio calculations at the RCCSD(T)/cc-pVDZ level, while our
Fig. 5. Reactive cross sections for the C(3P) + C2H2(v = 0,1 and j = 0,2,4) reaction obtained from the present 3D quantum calculations as a function of the translational energy. Open circles show the cross sections for C(3P) + C2H2(j = 0) obtained from the 2D quantum calculations (see text). Also plotted as a bold line is the cross sections of the 2D wave packet calculations by Buonomo and Clary taken from [16].
65
potential surface is based on the B3LYP/6-31G(d,p) level calculations. In order to confirm this, we have carried out two-dimensional quantum scattering calculations using the same Hamiltonian employed by Buonomo and Clary but using our B3LYP/6-31G(d,p) potential surface. The calculated cross sections for the C(3P) + C2H2(j = 0) reaction are plotted also in Fig. 5. It is seen that the calculated 2D cross sections agree with the 3D results. On the other hand, we notice that Buonomo and Clary also carried out the ACCSA calculations on the same potential energy surface used in the 2D wave packet calculations. Interestingly, the ACCSA cross sections are much larger than their 2D wave packet results by a factor of about 2 and are quite close to the present 3D and 2D results (see Fig. 9 of Ref. [16]). As mentioned previously and also discussed in [16], this disagreement mainly comes from the difference in the magnitude of the reaction probability. Therefore, it can be concluded that the 2D potential energy surface of Buonomo and Clary yields somewhat large nonreactive probabilities. In other words, nonreactive reflection is very important for the 2D surface of Buonomo and Clary, while the so-called capture theory works well for our potential energy surface. This is interesting since the overall feature of the two potential energy surfaces are quite similar. Since the present reduced-dimensionality model includes the CC stretching, it should be interesting to see the effect of the CC vibrational excitation on the reactive cross section. The cross section for v = 1 and j = 0 is plotted in Fig. 5. It is noted that the cross section has a broad maximum around the translational energy of 5 kJ/mol. This is mainly due to the fact that the vibrationally adiabatic potential curve for v = 1 is slightly different from that for v = 0; the former has a very small barrier but the latter is completely attractive. A more important point which should be mentioned is that the cross section for v = 1 is smaller than that for v = 0 by a factor of 1.5–2. This result indicates that vibrational excitation of the CC stretching does not enhance the reactivity but rather suppresses the reaction. After detailed analyses of the computational results, we have found that this behavior is due to the existence of the vibrational quenching processes for the v = 1 collisions. Notice that the reaction for v = 0 occurs with the probability being nearly unity (see Fig. 3). On the other hand, there exists the nonreactive, vibrational quenching process, v = 1 ! v = 0, for the C(3P) + C2H2(v = 1) collision. As a result, the reaction probability for C(3P) + C2H2(v = 1) is always smaller than that for C(3P) + C2H2(v = 0). Hence, the cross section for v = 1 is smaller than that for v = 0. Finally, Fig. 6 compares the energy dependence of the calculated reactive cross sections with the experimental results for the C(3P) + C2H2 reaction. It can be seen from Fig. 6 that the agreement between theory and
66
T. Takayanagi / Chemical Physics 312 (2005) 61–67
Fig. 6. Plot of logarithm of the integral cross sections for the C(3P) + C2H2 reaction obtained by the present 3D quantum scattering calculations plotted against the logarithm of the translational energy. The bold straight line indicates the E0.8 dependence obtain from the crossed molecular beam experiment of [12].
experiment is less good. The crossed molecular beam experiment shows the E0.8 dependence in the translational energy range of 0.3–20 kJ/mol [12], although only relative cross sections have been measured. Our cross sections approximately show the E0.5 dependence above the total energy of 2 kJ/mol. In a smaller energy region, theory gives a much smaller slope. Interestingly, Clary et al. [1] reports the very good agreement between experiment and the 2D wave packet result in the energy range of 1–5 kJ/mol for the C(3P) + C2D2 reaction. It has been reported that the energy dependence of the cross section for C(3P) + C2D2 is slightly different from that for C(3P) + C2H2 (E0.71 vs E0.80). Although whether this slight difference comes from the experimental uncertainty or isotope effects is not known, it is highly desirable that absolute values of the reactive cross sections should be measured from the experimental side. Also, further theoretical studies should definitely be performed in the near future in order to understand all dynamical aspects of the C(3P) + C2H2 reaction. Those would include the improvement of the potential energy surface, the understanding of the electronically nonadiabatic effect due to spin–orbit interactions, and theoretical determination of reaction products using an ab initio direct dynamics technique. Investigation along this line is currently undertaken in our laboratory.
4. Summary The reaction between carbon atoms C(3P) and acetylene has been studied theoretically using a reduced dimensionality model, in which hydrogen atom motions were treated in an adiabatic approximation and the
three degrees of freedom describing the C–C2 relative motions were explicitly treated in the dynamics. A three-dimensional potential energy surface has been obtained based on the electronic structure calculations at the B3LYP/6-31G(d,p) level and an analytical function were then developed for a general use. We have performed time-independent quantum reactive scattering calculations using the generalized R-matrix propagation method combined with an appropriate negative imaginary potential, which absorbs the reactive flux efficiently. The calculated reaction probabilities showed a typical result for a barrierless reaction but we have found that the probabilities calculated by the present three-dimensional model are larger than those obtained by the two-dimensional model of Buonomo and Clary [16]. This discrepancy is due to the difference in the potential energy surfaces employed; nonreactive reflection is more important for the surface of Buonomo and Clary than the present B3LYP surface. Our total reactive cross sections were thus larger than those of Buonomo and Clary [16], but were close to the ACCSA cross sections on the same potential surface of Buonomo and Clary. It has also been found that excitation of the initial rotation does not enhance the reactivity, at least up to j = 4. However, it is found that excitation of the C–C stretch considerably decreases the reaction cross section. This is simply because of the opening of the additional vibrational inelastic channel at high energies. Comparison of the energy dependence of the computed total reactive cross sections with the experimental results shows significant disagreement; our cross sections revealed a much weaker energy dependence. This may be due to the inaccuracy of the potential energy surface mainly around the long range region.
References [1] D.C. Clary, E. Buonomo, I.R. Sims, I.W.M. Smith, W.D. Geppert, C. Naulin, M. Costes, L. Cartechini, P. Casavecchia, J. Phys. Chem. A 106 (2002) 5541. [2] I.W.M. Smith, Chem. Soc. Rev. 31 (2002) 137. [3] N. Haider, D. Husain, Z. Phys. Chem. (Munich) 176 (1992) 133. [4] N. Haider, D. Husain, J. Chem. Soc., Faraday Trans. 89 (1993) 7. [5] D.C. Clary, N. Haider, D. Husain, M. Kabir, Astrophys. J. 422 (1994) 416. [6] D. Chastaing, P.L. James, S.D. Le Picard, I.R. Sims, I.W.M. Smith, Phys. Chem. Chem. Phys. 1 (1999) 2247. [7] D. Chastaing, S.D. Le Picard, I.R. Sims, I.W.M. Smith, Astron. Astrophys. 365 (2001) 241. [8] R.I. Kaiser, D. Stranges, Y.T. Lee, A.G. Suits, Astrophys.J. 477 (1997) 982. [9] R.I. Kaiser, C. Ochsenfeld, M. Head-Gordon, Y.T. Lee, A.G. Suits, J. Chem. Phys. 106 (1997) 1729. [10] C. Ochsenfeld, R.I. Kaiser, Y.T. Lee, A.G. Suits, M. HeadGordon, J. Chem. Phys. 106 (1997) 4141. [11] R.I. Kaiser, A.M. Mebel, Y.T. Lee, J. Chem. Phys. 114 (2001) 231.
T. Takayanagi / Chemical Physics 312 (2005) 61–67 [12] L. Cartechini, A. Bergeat, G. Capozza, P. Casavecchia, G.G. Volpi, W.D. Geppert, C. Naulin, M. Costes, J. Chem. Phys. 116 (2002) 5603. [13] J. Takahashi, K. Yamashita, J. Chem. Phys. 104 (1996) 6614. [14] R. Guadagnini, G.C. Schatz, S.P. Walch, J. Phys. Chem. A 102 (1998) 5857. [15] A.M. Mebel, W.M. Jackson, A.H.H. Chang, S.H. Lin, J. Am. Chem. Soc. 120 (1998) 5751. [16] E. Buonomo, D.C. Clary, J. Phys. Chem. A 105 (2001) 2694. [17] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, A.G. Baboul, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B.
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
67
Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian 98, Revision A.9, Gaussian, Inc., Pittsburgh, PA, 1998. A. Aguado, M. Paniagua, J. Chem. Phys. 96 (1992) 1265. T. Takayanagi, Y. Kurosaki, J. Chem. Phys. 113 (2000) 7158. K. Yagi, T. Takayanagi, T. Taketsugu, K. Hirao, J. Chem. Phys. 120 (2004) 10395. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, M. Baer, Chem. Phys. Lett. 291 (1998) 346. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, J. Chem. Phys. 109 (1998) 5761. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, J. Chem. Phys. 111 (1999) 1979. P. McGuire, D.J. Kouri, J. Chem. Phys. 60 (1974) 2488. J.C. Light, I.P. Hamilton, J.V. Lill, J. Chem. Phys. 82 (1985) 1400. D. Neuhauser, M. Baer, J. Chem. Phys. 90 (1989) 4351. D. Neuhauser, M. Baer, J. Chem. Phys. 92 (1990) 3419.
Reduced-dimensionality quantum reactive scattering calculations of the C(3P) + C2H2 reaction on a new potential energy surface Toshiyuki Takayanagi
*
Department of Chemistry, Saitama University, 255 Shimo-Okubo, Sakura-ku, Saitama City, Saitama 338-8570, Japan Received 28 September 2004; accepted 16 November 2004 Available online 8 December 2004
Abstract The reaction of the ground-state carbon atom with acetylene has been theoretically studied by a reduced dimensionality method, where three degrees of freedom were taken into account in the quantum dynamics calculations. A three-dimensional potential energy surface has been developed on the basis of the B3LYP/6-31G(d,p) level electronic structure calculations. The total reactive cross sections have been calculated by the time-independent generalized R-matrix method with the use of the negative imaginary potential. The reactive cross sections obtained in the energy range between 0.1 and 20 kJ/mol show a qualitative trend that is consistent with the barrierless reaction but the energy dependence was found to be much weaker than the experimental results. It has also been found that excitation of the initial rotation does not affect the reactivity but excitation of the CC stretch significantly decreases the reaction cross section. 2004 Elsevier B.V. All rights reserved.
1. Introduction The reaction of the ground-state atomic carbon C(3P) with acetylene has recently attracted considerable interest due to its fundamental importance in the chemistry of dense interstellar clouds [1,2]. In particular, this reaction is thought as a key step for the synthesis of larger carbon-containing species since the overall process of the reaction of C(3P) with unsaturated hydrocarbons may be represented by C + CnHm ! Cn+1Hm1 + H. The importance of the reaction of C(3P) with acetylene in interstellar clouds has been pointed out by the rate constant measurements. The rate constant at room temperature for this reaction was reported to be very large, suggesting no barrier in the potential energy surface of the C(3P) + C2H2 reaction [3–5]. The kinetic measurement was then extended down to the temperature range of 300–15 K [6,7] and it has been found that *
Tel.: 48 858 9113; fax: 48 858 3700. E-mail address: [email protected].
0301-0104/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.11.026
the rate constant has a negative temperature dependence and that the absolute value of the rate constant is very large (>2 · 1010 cm3 molecule1 s1) even at 15 K. This is quite important from an astrophysical point of view since, in interstellar clouds, the density is extremely low and a typical temperature is in the range of 10–50 K. Although the previous kinetic studies mentioned above played a decisive role for understanding of the importance of the carbon–acetylene reaction in interstellar clouds, those studies cannot identify the products of the reaction. According to the ab initio electronic structure calculations[8–10], it is energetically possible that the C(3P) + C2H2 reaction leads to the production of two C3H isomers: cyclic-C3H and linear-C3H. Kaiser and coworkers [8–11] have performed an extensive series of crossed molecular beam experiments at relatively high collision energies with mass spectroscopic detection of the C3H molecule. Although the mass spectroscopic technique cannot distinguish these isomers in principle, they have concluded that the cyclic-C3H molecule is mainly produced with a minor contribution of the
62
T. Takayanagi / Chemical Physics 312 (2005) 61–67
linear-C3H product at higher collision energies from the obtained angular distributions. Very recently, the crossed molecular beam experiment [12] have also revealed that the C(3P) + C2H2 reaction leads to C3 ð1 Rþ g Þ þ H2 products through the spin-forbidden intersystem-crossing mechanism. Ab initio electronic structure calculations are usually quite useful since they can provide information about the energetic accessibility of reaction products for which there is no thermochemical data and about the potential energy barriers along the reaction pathways from reactants to possible products. In fact, extensive ab initio calculations for the C(3P) + C2H2 reaction at various levels of theory have been carried out so far [8–11,13– 16]. In addition, Buonomo and Clary [16] have recently reported a quantum dynamics study for the first time. They have treated the C(3P) + C2H2 reaction with a reduced dimensionality model, in which two degrees of freedom, the distance R between C(3P) and the center of mass of C2H2 and the angle R makes with the C2H2 bond axis, are taken into account in the dynamics. They have developed a two-dimensional potential energy surface by the ab initio CCSD(T) method with the other seven degrees of freedom being optimized. They have then performed two-dimensional time-dependent wave packet calculations to obtain reaction probabilities. Although the wave packet calculations were carried out only in the entrance region of the potential energy surface, they were able to obtain the information about the cyclic-C3H/linear-C3H product branching by changing the location of the flux absorbing potentials. They have found that linear-C3H is formed preferentially over the energy range considered and that the cyclic-C3H product contributes at higher collision energies. Clearly, this result contradicts an interpretation of the crossed molecular beam experiments [8–10]. The purpose of the present study is to further elucidate the dynamics of the C(3P) + C2H2 reaction from a theoretical point of view. We extend the reduceddimensional model of Buonomo and Clary [16] so as to include one more degree of freedom; three degrees of freedom are taken into account in our quantum dynamics model. We mainly focus on the collision energy dependence and on the initial state dependence of total reaction cross sections.
2. Computational procedure 2.1. Construction of the potential energy surface In this study, the potential energy surface for the C(3P) + C2H2 addition reaction has been calculated as a function of three coordinates, R, r and h, describing the distance between C(3P) and the midpoint of CC in C2H2, the CC internuclear distance, and the orientation
angle of the approaching carbon atom with respect to the CC axis, respectively. This means that the C(3P) + C2H2 reaction system is treated as a standard atom–diatom system. The other six degrees of freedom have been optimized with respect to the total electronic energy. Our method is a natural extension of the theoretical work of Buonomo and Clary [16], who have carried out the reduced dimensionality quantum dynamics calculations considering two degrees of freedom as active coordinates. They have excluded the C–C internuclear coordinate from their dynamics calculations. However, our three-dimensional model is rather realistic in the sense of the Born–Oppenheimer-like approximation where fast hydrogen atom motions were assumed to be adiabatic. The potential energy values were calculated using the hybrid density functional B3LYP method with the 631G(d,p) basis set. A total of about 2600 points were ˚, calculated in the ranges of 1.0 6 R 6 5.0 A ˚ 1.0 6 r 6 2.2 A, and 0 6 h 6 90 using Gaussian 98 set of programs [17]. In an early stage of our investigation, we have tried to calculate all three-dimensional grid points in order to apply the three-dimensional interpolation method; however, we have encountered a self-consistent field convergence problem at many geometrical configurations. In addition, in the region where r is large, we often had difficulty in obtaining smooth potential energy values. From these reasons, we have fitted the calculated potential energies to the following many-body type analytical function: V ðR1 ,R2 ,R3 Þ ¼ V 2 ðR1 Þ þ
M X i,j,k
ci,j,k qi1 qj2 qk3 ,
with the constrains i + j + k 6 M. Here, Ri is the C–C diatomic internuclear distance. R1 is the C–C distance in the reactant acetylene molecule and the first term V2 is the two-body term, for which we have employed the standard Morse potential. The second term is the three-body term proposed by Aguado and Paniagua [18]. Notice that the variable qi ¼ Ri eai Ri vanishes for Ri = 0 and Ri ! 1. In our fit, we chose M = 7. Also, we have taken the permutation symmetry with respect the Jacobi angle, h, into account. All the parameters in the above analytical function were determined by the non-linear least-square method. The resulting three-dimensional potential energy surface is plotted in Fig. 1. For both h = 0 and h = 30 approaches of C(3P) toward the acetylene molecule, it is seen that there exist saddle points in the entrance region. On the other hand, no barrier is seen for both h = 60 and h = 90 approaches. In order to understand the reaction mechanism more clearly, the two-dimensional potential energy surface is plotted in Fig. 2 as a function of R and h, where the coordinate r was optimized with respect to the energy. An overall feature of our potential
T. Takayanagi / Chemical Physics 312 (2005) 61–67
63
energy surface is quite similar to the two-dimensional surface of Buonomo and Clary [16]. For a T-shaped approach (h = 90), we can see a shallow well around ˚ corresponding to the cyclic-C3H2 intermediR = 1.1 A ate. The path then reaches the linear-C3H2 intermediate ˚ . Two after passing the saddle point around R 0.95 A ˚ potential wells are also seen around R 2 A for both h = 0 and h = 180, which also leads to the linearC3H2 intermediate formation. In addition, barriers are ˚ , for small seen in the region between R = 2.5 4.2 A angles h. At this point, it is informative to compare the energy levels of intermediates obtained at the B3LYP/631G(d,p) level of theory to the result obtained at a more accurate level of theory. Cyclic-C3H2 (cyclopropenylidene) has an energy of 262 kJ/mol with respect to the reactants while linear-C3H2 (propargylene) has an energy of 432 kJ/mol at the B3LYP/6-31G(d,p) level. These values are somewhat smaller than the corresponding CCSD(T)/QZ2P level values [10], 225 and 386 kJ/mol, respectively. 2.2. Quantum dynamics calculations
Fig. 1. Contour plots of the potential energy surface for the C(3P) + C2H2 reaction as a function of Jacobi coordinates, R, r and h. Zero of the energy is taken to be the reactant C(3P) + C2H2 potential minimum. Contours are spaced by 20 kJ/mol; solid curves are used for energies that are positive relative to C(3P) + C2H2, dashed curves are used for energies which are negative, and bold lines denote the zero contours.
We have carried out quantum reactive scattering calculations by the generalized R-matrix propagation method [19,20], which has been developed by Aguilar and coworkers [21–23] to solve the time-independent close-coupling equations with the use of a complex absorbing potential. If the absorbing potential is properly located in the product region and efficiently absorbs the reactive flux, we can easily obtain stateto-all reaction probabilities without the information of the product side. The great advantage of this method is that it can treat a reactive scattering problem in a formulation very similar to an inelastic scattering problem. As mentioned previously, the reaction system has been treated as a standard atom–diatom system. The Hamiltonian is then given by H ¼
Fig. 2. Contour plots of the potential energy surface for the C(3P) + C2H2 reaction as a function of R and h. Contours are spaced by 20 kJ/mol.
h2 o2 h2 o2 ðJ jÞ2 h2 j2 h2 þ þ þ V ðR,r,hÞ, 2lR oR2 2lr or2 2lr r2 2lR R2
where lR is the reduced mass between C and C2H2, J is the total angular momentum operator, j is the rotational angular momentum operator of C2H2, and lr is the reduced mass of (CH)–(CH). The calculations for the total angular momentum quantum number J > 0 were carried out using the standard centrifugal sudden (CS) approximation [24]; The operator term, (J j)2, was replaced into [J(J + 1) 2X2 + j2], where X is the projection on the body-fixed Z axis. Thus, Coriolis terms that couple states with different X have been neglected. This means that the projection quantum number X is assumed to be conserved.
64
T. Takayanagi / Chemical Physics 312 (2005) 61–67
The total scattering wave function can be expanded using the product of vibrational–rotational basis sets and usual close-coupling equations can be obtained although the matrix elements are complex. A twodimensional eigenvalue problem at a given value of the Jacobi scattering coordinate R was solved using the discrete-variable-representation (DVR) method [25]. We employed the particle-in-box basis for the C–C stretch coordinate r and associated Legendre basis for the angular coordinate h. For the DVR calculations, 200 basis sets in r and 100 basis sets in h were employed. The calculations can be decoupled into subprograms by exploiting the permutation symmetry of the two identical C atoms. The scattering coordinate R was first divided into ˚ . We have em600 sectors in a range of R = 1.0–5.0 A ployed the linear Neuhauser–Baer [26,27] negative imag˚ . The inary potential, located in the range of 1 < R < 2 A absorbing strength parameter was determined so as that the imaginary potential effectively absorbs the reactive flux. Note that the present reactive probability includes both the cyclic-C3H and linear-C3H production channels unlike the two-dimensional wave packet calculations of Buonomo and Clary [16]. Calculations have been carried out in the range of total angular momentum between 0 and 150, with a step DJ being 1–2, depending on the total energy. The maximum value of the projection quantum number X was taken to be 4. Thus, we were able to calculate the cross sections for the C(3P) + C2H2(j) reaction up to j = 4.
3. Results and discussion The C(3P) + C2H2(v = 0, j = 0 20) reaction probabilities for the total angular momentum J = 0 are shown in Fig. 3 as a function of the total energy. Here, v is the vibrational quantum number of the CC stretch. The total energy is measured from the asymptotic energy level of the C(3P) + C2H2(v = j = 0) channel. Although many sharp resonances are seen, the reaction probabilities show a sudden increase just after the reaction threshold and can approximately be represented by a step function. This result is quite typical for barrierless reactions. Also, it is seen that values of average probabilities are quite large and range from 0.8 to 1. The reaction probability for C + C2H2(j = 0) is almost unity above the reaction threshold. This is highly in contrast with the wave packet result of Buonomo and Clary [16]; their calculation shows that the reaction probabilities for j = 0 range from 0 to 0.8 and are somewhat smaller than our result. Plots of cumulative reaction probabilities, which correspond to total reaction probabilities summed over all initial states, are given for different numbers of J(X = 0) in Fig. 4. It is seen that the cumulative reaction probability for J = 0 shows less resonance structures
Fig. 3. Reaction probabilities for the C(3P) + C2H2(v = 0, j) reaction for the zero total angular momentum (J = 0) as a function of the total energy. The zero of the total energy is measured from the asymptotic C(3P) + C2H2(v = j = 0) level.
Fig. 4. Plots of cumulative reaction probabilities for the C(3P) + C2H2 reaction for different number of J as a function of the total energy. The zero of the total energy is measured from the asymptotic C(3P) + C2H2(v = j = 0) level.
compared to the initial state-selected reaction probabilities presented in Fig. 3. This indicates that the increases and decreases in the initial state-selected reaction prob-
T. Takayanagi / Chemical Physics 312 (2005) 61–67
abilities at resonance energies are compensated each other in the cumulate reaction probability. Although resonances become less important as the value of J increases, a typical J-shifting behavior is clearly seen; i.e., the reaction threshold energy gradually increases as J increases. The number of open channels is also plotted in Fig. 4 as a dotted line. It is interesting to note that the sudden step feature can be seen only for the cumulative reaction probability for J = 0 at low energies. The total reaction cross sections obtained in this study are presented in Fig. 5. It is found that the cross sections for v = 0 and j = 0, 2, and 4 give almost identical results in the energy range considered. This behavior is somewhat different from the result of adiabatic capture centrifugal sudden approximation (ACCSA) calculations by Buonomo and Clary [16]. The ACCSA calculations show that the initial state selected reaction rate constant decreases as the initial rotational state increases, although their ACCSA calculations have been done up to j = 12. Since our quantum scattering calculations have been done for Xmax = 4, we could not obtain the converged cross section for such a large value of j. Also, unfortunately, Buonomo and Clary did not perform the wave packet calculations for j > 0. It should be important to compare the present cross sections with the 2D wave packet result by Buonomo and Clary [16]. Their cross section for j = 0 is also plotted in Fig. 5 as a bold line. It can be seen that our cross section is larger than the 2D result by a factor of about 2–2.5. It is suggested that the source of this disagreement mainly comes from the difference between the potential energy surfaces employed. Buonomo and Clary [16] have developed their surface on the basis of the ab initio calculations at the RCCSD(T)/cc-pVDZ level, while our
Fig. 5. Reactive cross sections for the C(3P) + C2H2(v = 0,1 and j = 0,2,4) reaction obtained from the present 3D quantum calculations as a function of the translational energy. Open circles show the cross sections for C(3P) + C2H2(j = 0) obtained from the 2D quantum calculations (see text). Also plotted as a bold line is the cross sections of the 2D wave packet calculations by Buonomo and Clary taken from [16].
65
potential surface is based on the B3LYP/6-31G(d,p) level calculations. In order to confirm this, we have carried out two-dimensional quantum scattering calculations using the same Hamiltonian employed by Buonomo and Clary but using our B3LYP/6-31G(d,p) potential surface. The calculated cross sections for the C(3P) + C2H2(j = 0) reaction are plotted also in Fig. 5. It is seen that the calculated 2D cross sections agree with the 3D results. On the other hand, we notice that Buonomo and Clary also carried out the ACCSA calculations on the same potential energy surface used in the 2D wave packet calculations. Interestingly, the ACCSA cross sections are much larger than their 2D wave packet results by a factor of about 2 and are quite close to the present 3D and 2D results (see Fig. 9 of Ref. [16]). As mentioned previously and also discussed in [16], this disagreement mainly comes from the difference in the magnitude of the reaction probability. Therefore, it can be concluded that the 2D potential energy surface of Buonomo and Clary yields somewhat large nonreactive probabilities. In other words, nonreactive reflection is very important for the 2D surface of Buonomo and Clary, while the so-called capture theory works well for our potential energy surface. This is interesting since the overall feature of the two potential energy surfaces are quite similar. Since the present reduced-dimensionality model includes the CC stretching, it should be interesting to see the effect of the CC vibrational excitation on the reactive cross section. The cross section for v = 1 and j = 0 is plotted in Fig. 5. It is noted that the cross section has a broad maximum around the translational energy of 5 kJ/mol. This is mainly due to the fact that the vibrationally adiabatic potential curve for v = 1 is slightly different from that for v = 0; the former has a very small barrier but the latter is completely attractive. A more important point which should be mentioned is that the cross section for v = 1 is smaller than that for v = 0 by a factor of 1.5–2. This result indicates that vibrational excitation of the CC stretching does not enhance the reactivity but rather suppresses the reaction. After detailed analyses of the computational results, we have found that this behavior is due to the existence of the vibrational quenching processes for the v = 1 collisions. Notice that the reaction for v = 0 occurs with the probability being nearly unity (see Fig. 3). On the other hand, there exists the nonreactive, vibrational quenching process, v = 1 ! v = 0, for the C(3P) + C2H2(v = 1) collision. As a result, the reaction probability for C(3P) + C2H2(v = 1) is always smaller than that for C(3P) + C2H2(v = 0). Hence, the cross section for v = 1 is smaller than that for v = 0. Finally, Fig. 6 compares the energy dependence of the calculated reactive cross sections with the experimental results for the C(3P) + C2H2 reaction. It can be seen from Fig. 6 that the agreement between theory and
66
T. Takayanagi / Chemical Physics 312 (2005) 61–67
Fig. 6. Plot of logarithm of the integral cross sections for the C(3P) + C2H2 reaction obtained by the present 3D quantum scattering calculations plotted against the logarithm of the translational energy. The bold straight line indicates the E0.8 dependence obtain from the crossed molecular beam experiment of [12].
experiment is less good. The crossed molecular beam experiment shows the E0.8 dependence in the translational energy range of 0.3–20 kJ/mol [12], although only relative cross sections have been measured. Our cross sections approximately show the E0.5 dependence above the total energy of 2 kJ/mol. In a smaller energy region, theory gives a much smaller slope. Interestingly, Clary et al. [1] reports the very good agreement between experiment and the 2D wave packet result in the energy range of 1–5 kJ/mol for the C(3P) + C2D2 reaction. It has been reported that the energy dependence of the cross section for C(3P) + C2D2 is slightly different from that for C(3P) + C2H2 (E0.71 vs E0.80). Although whether this slight difference comes from the experimental uncertainty or isotope effects is not known, it is highly desirable that absolute values of the reactive cross sections should be measured from the experimental side. Also, further theoretical studies should definitely be performed in the near future in order to understand all dynamical aspects of the C(3P) + C2H2 reaction. Those would include the improvement of the potential energy surface, the understanding of the electronically nonadiabatic effect due to spin–orbit interactions, and theoretical determination of reaction products using an ab initio direct dynamics technique. Investigation along this line is currently undertaken in our laboratory.
4. Summary The reaction between carbon atoms C(3P) and acetylene has been studied theoretically using a reduced dimensionality model, in which hydrogen atom motions were treated in an adiabatic approximation and the
three degrees of freedom describing the C–C2 relative motions were explicitly treated in the dynamics. A three-dimensional potential energy surface has been obtained based on the electronic structure calculations at the B3LYP/6-31G(d,p) level and an analytical function were then developed for a general use. We have performed time-independent quantum reactive scattering calculations using the generalized R-matrix propagation method combined with an appropriate negative imaginary potential, which absorbs the reactive flux efficiently. The calculated reaction probabilities showed a typical result for a barrierless reaction but we have found that the probabilities calculated by the present three-dimensional model are larger than those obtained by the two-dimensional model of Buonomo and Clary [16]. This discrepancy is due to the difference in the potential energy surfaces employed; nonreactive reflection is more important for the surface of Buonomo and Clary than the present B3LYP surface. Our total reactive cross sections were thus larger than those of Buonomo and Clary [16], but were close to the ACCSA cross sections on the same potential surface of Buonomo and Clary. It has also been found that excitation of the initial rotation does not enhance the reactivity, at least up to j = 4. However, it is found that excitation of the C–C stretch considerably decreases the reaction cross section. This is simply because of the opening of the additional vibrational inelastic channel at high energies. Comparison of the energy dependence of the computed total reactive cross sections with the experimental results shows significant disagreement; our cross sections revealed a much weaker energy dependence. This may be due to the inaccuracy of the potential energy surface mainly around the long range region.
References [1] D.C. Clary, E. Buonomo, I.R. Sims, I.W.M. Smith, W.D. Geppert, C. Naulin, M. Costes, L. Cartechini, P. Casavecchia, J. Phys. Chem. A 106 (2002) 5541. [2] I.W.M. Smith, Chem. Soc. Rev. 31 (2002) 137. [3] N. Haider, D. Husain, Z. Phys. Chem. (Munich) 176 (1992) 133. [4] N. Haider, D. Husain, J. Chem. Soc., Faraday Trans. 89 (1993) 7. [5] D.C. Clary, N. Haider, D. Husain, M. Kabir, Astrophys. J. 422 (1994) 416. [6] D. Chastaing, P.L. James, S.D. Le Picard, I.R. Sims, I.W.M. Smith, Phys. Chem. Chem. Phys. 1 (1999) 2247. [7] D. Chastaing, S.D. Le Picard, I.R. Sims, I.W.M. Smith, Astron. Astrophys. 365 (2001) 241. [8] R.I. Kaiser, D. Stranges, Y.T. Lee, A.G. Suits, Astrophys.J. 477 (1997) 982. [9] R.I. Kaiser, C. Ochsenfeld, M. Head-Gordon, Y.T. Lee, A.G. Suits, J. Chem. Phys. 106 (1997) 1729. [10] C. Ochsenfeld, R.I. Kaiser, Y.T. Lee, A.G. Suits, M. HeadGordon, J. Chem. Phys. 106 (1997) 4141. [11] R.I. Kaiser, A.M. Mebel, Y.T. Lee, J. Chem. Phys. 114 (2001) 231.
T. Takayanagi / Chemical Physics 312 (2005) 61–67 [12] L. Cartechini, A. Bergeat, G. Capozza, P. Casavecchia, G.G. Volpi, W.D. Geppert, C. Naulin, M. Costes, J. Chem. Phys. 116 (2002) 5603. [13] J. Takahashi, K. Yamashita, J. Chem. Phys. 104 (1996) 6614. [14] R. Guadagnini, G.C. Schatz, S.P. Walch, J. Phys. Chem. A 102 (1998) 5857. [15] A.M. Mebel, W.M. Jackson, A.H.H. Chang, S.H. Lin, J. Am. Chem. Soc. 120 (1998) 5751. [16] E. Buonomo, D.C. Clary, J. Phys. Chem. A 105 (2001) 2694. [17] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, A.G. Baboul, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B.
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
67
Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian 98, Revision A.9, Gaussian, Inc., Pittsburgh, PA, 1998. A. Aguado, M. Paniagua, J. Chem. Phys. 96 (1992) 1265. T. Takayanagi, Y. Kurosaki, J. Chem. Phys. 113 (2000) 7158. K. Yagi, T. Takayanagi, T. Taketsugu, K. Hirao, J. Chem. Phys. 120 (2004) 10395. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, M. Baer, Chem. Phys. Lett. 291 (1998) 346. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, J. Chem. Phys. 109 (1998) 5761. F. Huarte-Larran˜aga, X. Gime´nez, A. Aguilar, J. Chem. Phys. 111 (1999) 1979. P. McGuire, D.J. Kouri, J. Chem. Phys. 60 (1974) 2488. J.C. Light, I.P. Hamilton, J.V. Lill, J. Chem. Phys. 82 (1985) 1400. D. Neuhauser, M. Baer, J. Chem. Phys. 90 (1989) 4351. D. Neuhauser, M. Baer, J. Chem. Phys. 92 (1990) 3419.