Ab initio Based DMBE Potential Energy Surface for the Ground

Feb 3, 2010 - Introduction. The ethynyl radical C2H is an important intermediate in hydrocarbon reaction systems in various chemical processes in...
0 downloads 0 Views 1MB Size
J. Phys. Chem. A 2010, 114, 2655–2664

2655

Ab initio Based DMBE Potential Energy Surface for the Ground Electronic State of the C2H Molecule S. Joseph and A. J. C. Varandas* Departamento de Quı´mica, UniVersidade de Coimbra, 3004-535 Coimbra, Portugal ReceiVed: October 27, 2009; ReVised Manuscript ReceiVed: January 11, 2010

We report a single-sheeted potential energy surface for the lowest state of C2H using double many-body expansion theory and accurate ab initio data points calculated at MRCI/AVQZ level of theory. The topographical features of the new global potential energy surface are analyzed in detail. 1. Introduction The ethynyl radical C2H is an important intermediate in hydrocarbon reaction systems in various chemical processes in the laboratory, such as combustion, discharges, and photolysis.1 It is supposed to be one of the most abundant polyatomic species in the interstellar medium,2 playing an important role in the formation and destruction of polyynes H(CC)n in comets and the interstellar media.1,3 Moreover, C2H serves as a precursor molecule for the formation of large carbon rich aliphatic and aromatic hydrocarbons in the pyrolysis reactions of hydrocarbons, particularly C2H2.4,5 In addition, C2H is a product of some reactions such as C2 + CH46 and H + C2H2,7 acting as a source of C2 and CH through photodissociation processes,8 mainly via the reactions:

C2H + hν f C2 + H

(1)

C2H + hν f C + CH

(2)

In fact, reaction 1 has also been used to explain properties of the emission of C2 in comets.9 Molecular structure makes ethynyl radical a specific one, because it is the simplest organic triatomic molecule having a triple C-C bond with a low-lying electronic state; see Figure 1 for a schematic diagram of some of its low-lying electronic states. Studies of C2H are therefore of great interest to both chemists and astrophysicists. In the laboratory, C2H was first observed in 1963 by Cochran and co-workers10 by electron spin resonance in a solid argon matrix at liquid helium temperatures, and later in 1973 in the gas phase.11 Subsequently, C2H as well as its isotopic analogues have been studied extensively with various experimental and spectroscopic methods, including ESR12,13 (electron spin resonance), LMR14,15 (laser magnetic resonance), microwave and millimeter-wave spectroscopy,16-18 matrix isolation IR (infrared) spectroscopy,19-22 color center laser spectroscopy,23-25 diode laser spectroscopy,26-29 FTIR (Fourier transform infrared) emission spectroscopy,30,31 and ultraviolet spectroscopy.32 Recently, C2H and its isotopic variant C2D radicals as well as the corresponding anions have been investigated by Zhou et al.33 using slow electron velocity-map imaging (SEVI). They attributed several higher-lying vibronic levels of the C2H molecule, with the spectra yielding new insights into the vibronic * Corresponding author. E-mail address: [email protected].

Figure 1. Dissociation scheme of C2H into C2 + H and C + CH. Solid lines correspond to 2Σ+ states and dotted lines represent 2Π states. (Adapted from ref 87.)

˜ and A ˜ states. In turn, Shi et al.34 studied coupling between the X various thermochemical properties, namely the reaction threshold energies and C-H bond dissociation energies of HC2n and HC2n- (n)1-3) linear molecular species, using guided ion beam tandem mass spectrometry (GIBMS). Moreover, these authors have reported the experimental C-H bond dissociation energy of the neutral ethynyl radical. At the same time, the experimental studies have been augmented by extensive theoretical ones,35-39 mutually promoting the development and understanding of the spectroscopy of C2H. Most such studies, both theoretical and experimental, have been devoted to the low lying A2Π r X2Σ+ transition. This transition lies in the IR region, and assignment of its spectrum has been complicated particularly for vibrationally excited levels due to conical intersections involving the A2Π state. For this reason, there is some uncertainty40 on the structure and energetics of C2H. For example “pure” gradient-corrected functionals such as BP86, BLYP, and PWP86 predict a bent structure41 while Hartree-Fock and post-Hartre-Fock methods predict (refs 42 and 43 and references therein) it to have a linear geometry instead. In recent years, a series of high-level ab initio ˜ and A ˜ states of C2H (12A′, studies has been performed on the X 2 2 2 A′, and 1 A′′) by Peyerimhoff and co-workers,39,44,45 and Mebel et al.46 However, the wealth of the papers by Peyerimhoff and co-workers39,44,45 concerns the vibronic energy levels and transition probabilities associated to them, while the one by Mebel et al.46 focuses on conical intersections. More recently, Tarroni et al.47,48 calculated the theoretical, dipole allowed,

10.1021/jp910269w  2010 American Chemical Society Published on Web 02/03/2010

2656

J. Phys. Chem. A, Vol. 114, No. 7, 2010

absorption spectra of C2H and its deuterated isotopomer C2D up to 10 000 cm-1 using ab initio diabatic energies and dipole moment surfaces. In turn, Szalay et al.49 obtained heats of formation of ethynyl radicals from high level ab initio calculations. These calculations include a nonrelativistic electronic energy from extrapolated coupled-cluster calculations (up to CCSDTQ), corrections for scalar relativistic and spin-orbit effects, and the diagonal Born-Oppenheimer correction. They also give the equilibrium geometry and harmonic and fundamental frequencies of C2H. Regarding global potential energy surfaces (PESs) involved in the reaction C + CH f C2H f C2 + H, Rayez and co-workers50 have performed a CASSCF (complete active space self-consistent field) study and later built diabatic PESs in collaboration with one of the authors51,52 for the three lowest states of C2H. Such a PES displays a linear geometry for equilibrium C2H. Very recently, C2H has gained additional attention since Mebel and co-workers53 proposed a new photoinduced mechanism for formation and growth of polycyclic aromatic hydrocarbons (PAH) in low-temperature environments via successive ethynyl radical additions, an alternative to high-temperature hydrogen-abstraction-C2H2-addition (HACA) sequences to form polycyclic aromatic hydrocarbon (PAH)-like species under low-temperature conditions in the interstellar medium or planetary atmospheres, such as on Titan. They have found from their ab initio calculations that the driving force of such a mechanism is the presence of ethynyl radicals. This prompted us to build a global PES for the lowest adiabatic sheet of 12A′ of the ethynyl radical, aiming at a further understanding of the structure, energetics, and spectroscopy of the title radical, and progressing toward a PES for larger hydrogen-carbon containing species. This will be done by utilizing DMBE54-58 theory. Specifically, we will use the single-sheeted DMBE formalism which takes proper account of long-range forces. The paper is organized as follows. Section 2 describes the ab initio calculations carried out in the present work. In turn, section 3 deals with the analytical representation of the PES, with subsection 3.1 focusing on the two-body energy terms, and section 3.2 on three-body ones. The main topographical

Joseph and Varandas features of the DMBE potential energy surface are discussed in section 4. The concluding remarks appear in section 5. 2. Ab Initio Electronic Structure Calculations The ab initio calculations have been carried out at the MRCI59 level using the full valence complete active space60 (FVCAS) wave function as reference. Such a FVCAS reference wave function, which involves 9 correlated electrons in 9 active orbitals (7a′ + 2a′′), amounts to a total of 2236 configuration state functions. The core orbitals have been kept frozen in all calculations; a comment on the core-correlation contribution will be presented later. For the basis set, we have selected the augcc-pVQZ (denoted for simplicity by AVQZ) one of Dunning61,62 with calculations being carried out using the Molpro63 package. A total of 1830 grid points were chosen to map the PES of H-C2 over the regions defined by 1.5 e rC2/a0 e 3.5, 2.0 e RH-C2/a0 e 10. 0 and 0 0.0 e γ/deg e 90. For the C-CH interactions, a grid defined by 1 0.5 e rCH/a0 e 3.5, 2.0 e RC-CH/a0 e 10.0, and 0.0 e γ/deg e 180 has been utilized. For both channels, r, R, and γ define atom-diatom Jacobi coordinates. All calculated ab initio energies have been subsequently corrected by the double many-body-expansion-scaled-externalcorrelation64 (DMBE-SEC) method to account for excitations beyond singles and doubles and, most importantly, for the incompleteness of the basis set. The total DMBE-SEC interaction energy is then written as

V(R) ) VFVCAS(R) + VSEC(R)

(3)

where

VFVCAS(R) )

(2) (3) (RAB) + VABC,FVCAS (R) ∑ VAB,FVCAS AB

(4)

Figure 2. Ab initio and scaled scaled (DMBE-SEC) potential energy curve of CH(2Π). Also shown is the corresponding energy difference. Solid and dotted lines represent calculations without and with inclusion of core electrons; see text for details.

Potential Energy Surface for C2H

VSEC(R) )

J. Phys. Chem. A, Vol. 114, No. 7, 2010 2657

(2) (3) (RAB) + VABC,SEC (R) ∑ VAB,SEC

(5)

AB

and R ) {Ri} (i ) 1 - 3) defines the set of interatomic coordinates. In turn, the first two terms of the SEC series expansion assume the form

(2) VAB,SEC (RAB)

)

(2) (2) VAB,FVCAS/CISD (RAB) - VAB,FVCAS (RAB) (2) FAB

(6)

(3) VABC,SEC (R)

)

(3) (3) VABC,FVCAS/CISD (R) - VABC,FVCAS (R) (3) FABC

(7) (2) parameter in eq 6 has been Following previous work,64 the FAB chosen from the requirement that the corrected ab initio diatomic curves reproduce the experimental dissociation energies, while (3) has been estimated as the average of the two-body FABC F-factors. The ab initio potential energy curves for C2(X1Σg+) and CH(2Π) have then been suitably scaled using experimental dissociation energies, thus leading to small though important empirical corrections of the former. The (good) performance of the DMBE-SEC method has recently been examined in a study of the ground electronic state of SH2,65 with the reader being addressed to such a work for details. For C2(X1Σg+), the range of calculated values for the experimental dissociation energy (De) obtained is pretty large,66 varying from 50 504.6 ( 314.8 to 51 708.5 ( 174.9 cm-1 (see Table 3 of ref 66). Of them, we have chosen 51 708.5 cm-1 since this value agrees essentially within the error bars with the best available theoretical estimate reported in the above theoretical paper.66 For CH, we have in turn employed the value of67 29 358.74 cm-1. Using such values, the DMBE-SEC (2) (2) (3) ) 0.6483, FCH ) 0.9098, and FCCH ) procedure leads to FCC 0.8226. As noted above, the core has been kept frozen in all calculations. However, such effects have been accounted for

Figure 3. As in Figure 2 but for C2(1Σ+).

by the (small but important) empiricism introduced above. This can be seen from Figures 2 and 3 for CH and C2, respectively; the notation MRCI-Cx and SEC-Cx, with x ) 0, 2, or 4, indicating that the involved core electrons have been frozen during the MRCI calculations. The bottom panels of Figures 2 and 3 show the differences between the ab initio and DMBESEC energies. As shown from the bottom right-hand-side panels, the differences between the MRCI curves including corecorrelation effects68 and the DMBE-SEC ones are typically of up to a few mEh magnitude at the most relevant regions, although in highly repulsive regions they may be somewhat larger. Since the effect on relative energies in HC2 should be even smaller (especially if due mostly to three-body effects), core-correlation effects can likely be neglected without sacrificing significantly the accuracy of the PES. This will allow considerable savings of computational time, while being in alignment with our previous work65 on SH2. 3. Double Many-Body Expansion Potential Energy Surface 3.1. Two-Body Energy Terms. The potential energy curves of C2(X1Σg+) and CH(X2Π) have been extorted using the extended Hartree-Fock approximate correlation energy method for diatomic molecules, including the united atom limit69 (EHFACE2U), which shows the correct behavior at the asymptotic limits. They assume the general form69

V(R) ) VEHF(R) + Vdc(R)

(8)

where EHF denotes the extended Hartree-Fock type energy, and dc is the dynamical correlation energy. This is modeled semiempiricaly as described70

Vdc(R) ) -



n)6,8,10

χn(R)

Cn Rn

(9)

with the damping functions for the dispersion coefficients assuming the form

2658

J. Phys. Chem. A, Vol. 114, No. 7, 2010

[

(

χn(R) ) 1 - exp -An

R R2 - Bn 2 F F

Joseph and Varandas

)]

n

(10)

In eq 10, An and Bn are auxiliary functions54,55 defined by

An ) R0n-R1

(11)

Bn ) β0 exp(-β1n)

(12)

where R0, R1, β0, and β1 are universal dimensionless parameters for all isotropic interactions: R0 ) 16.366 06, R1 ) 0.701 72, β0 ) 17.193 38, and β1 ) 0.095 74. In turn, F is a scaling distance defined (in atomic units) by F ) 5.5 + 1.25R0, where R0 ) 2(〈rA2〉1/2 + 〈rB2〉1/2) is the LeRoy71 parameter for onset the undamped R-n expansion, and 〈rX2〉 is the expectation value of the square radius for the outermost electrons of atom X (X ) A, B). Note that only the value of the C6 dispersion coefficient is known for C2(1Σ+), and hence C8 and C10 have been estimated using the rule55

Cn ) C6κnRa(n-6)/2 0

Figure 4. EHFACE2U potential energy curves for CH(2Π) and C2(1Σ+ g ). The circles correspond to the ab initio points calculated at the present level of theory. Dashed lines represent lowest adiabatic curve of C2(1Σ+) obtained72 by diagonalizing a 2 × 2 potential matrix.

(13)

where κ8 ) 1, κ10 ) 1.13, and a ) 1.54 (this expression has been misprinted,66,72 although the calculations there reported are correct). Since there is some uncertainty in the C6 value that is distinct but rather similar for both states, we have utilized here the commonly accepted value73 of 40Eha06. In turn, we have chosen C5 ) 0 for the electrostatic coefficient as this correlates with the excited state involved in the C2 Σ+ avoided intersection.66,72 The exponentially decaying part of the EHF type energy is represented by the general form55

been adopted for consistency and simplicity. The potential energy curve for CH(2Π) is obtained by fitting the ab initio points calculated elsewhere.74 Figure 4 shows the both diatomic curves and the values of all parameters for both diatomic potentials are given in Table 1 of the Supporting Information. 3.2. Three-Body Energy Terms. 3.2.1. Three-Body Dynamical Correlation Energy. The three-body dynamical correlation energy assumes the usual form75 3

D VEHF(R) ) - (1 + R

V(3) dc )

n

∑ air ) exp(-γr) i

(14)

∑ ∑ fi(R) C(i)n (Ri,θi) χn(ri)ri-n i)1

(16)

n

i)1

where

γ ) γ0[1 + γ1 tanh(γ2r)]

(15)

and r ) R - Re is the displacement from the equilibrium diatomic geometry. In turn, D, ai, and γi (i ) 1, ..., n) are adjustable parameters.55,69 For the ground state C2(X1Σ+g ) we have used the MRCI/AVQZ points calculated in this work to model the potential energy curves. Figure 4 shows the calculated and scaled energies, in addition to the fitted analytical EHFACE2U69 (extended Hartree-Fock approximate correlation energy including the united-atom limit of the diatomic) curves so obtained. Also shown for comparison is the lowest adiabatic curve of C2(1Σ+) that has been obtained72 by diagonalizing a 2 × 2 potential matrix. Not surprisingly, the agreement is seen to be good except close to the avoided crossing region, which is rather difficult to model via a single curve fit. There is also some deviation at long-range distances probably due to having utilized singlestate energies in the present work. The deviation at the very repulsive wing are due to having imposed here the proper Coulombic behavior at the collapsed diatomic limit for vanishingly small interatomic separations, Rf0. However, the differences are small, and hence the single-sheeted procedure has

where i labels the I-JK channel (I, J, and K label the atoms) associated with each atom-diatom interaction, Ri is the diatomic internuclear distance, ri is the separation of atom I from a certain point located in the JK diatomic bond distance, and θi is the included angle between the two associated vectors. In turn, Cn(Ri,θi) are atom-diatom dispersion coefficients to be defined later, and χn(ri) are the usual dispersion damping functions.70 Moreover, the three-dimensional switching function assumes the form

1 fi ) {1 - tanh[ξ(ηRi - Rj - Rk)]} 2

(17)

Following other work,76 the numerical parameters have been fixed at η ) 6 and ξ ) 1.0a0-1. In addition, the switching function fi has been chosen from the requirement that its value must be +1 at Ri ) Rie and ri f ∞, and 0 when Ri f ∞. As discussed elsewhere,75 we have multiplied the two-body dynamical correlation energy for the ith pair by Πj*i(1 - fj) with corresponding factors for channels j and k, thus ensuring that the only two-body contribution at the ith channel is that of JK diatom pair. The atom-diatom dispersion coefficients in eq 16 are given by

Potential Energy Surface for C2H

C(i) n (Ri,θi) )

J. Phys. Chem. A, Vol. 114, No. 7, 2010 2659

∑ CLn(R) PL(cos θi)

(18)

L

where PL(cos θi) denotes the Lth Legendre polynomial. As usual, the expansion in eq 18 has been truncated by considering only the coefficients C06, C26, C08, C28, C48, and C010; all other coefficients have been supposed as making insignificant contributions and hence neglected. To estimate the dispersion coefficients, we have employed the generalized Slater-Kirkwood formula,77 with the dipolar polarizabilities calculated in the present work at the MRCI/AVQZ level of theory. Following the common procedure, the atom-diatom dispersion coefficients calculated for a set of internuclear distances were then fitted to the form 3

CL,A-BC (R) ) CL,AB + CL,AC + DM(1 + n n n

∑ airi)

×

i)1

3

exp(-a1r -

∑ biri)

(19)

i)2

where CnL,AB, for L ) 0, are the corresponding atom-diatom isotropic dispersion coefficients (see Table 2 of the Supporting Information; for L * 0, we have assumed CnL,AB ) 0). All parameters and coefficients in eq 19 have been explained elsewhere.76 The parameters that resulted from such fits are given in Table 2 of the Supporting Information, and the internuclear dependences of the dispersion coefficients are shown in Figure 5. 3.2.2. Three-Body Extended Hartree-Fock Energy. By removing, for a given triatomic geometry, the sum of the onebody and two-body energy terms from the corresponding DMBE-SEC interaction energies in eq 3, which was defined with respect to the infinitely separated ground-state atoms, one obtains the total three-body energy. By further subtracting the three-body dynamical correlation contribution, eq 16, from the total three-body energy that is calculated in that way, one obtains the three-body extended Hartree-Fock energy. This will be represented by the following three-body distributed-polynomial78 form

(3) VEHF )

6

3

j)1

i)1

∑ Pj(Q1,Q2,Q3) ∏ {1 - tanh[γij(Ri - Rij,ref)]} (20)

where Pj(Q1,Q2,Q3) is the jth polynomial, to the fifth order, defined as

Pj(Q1,Q2,Q3) ) c1 + c2Q1 + c3Q3 + c4Q12 + c5S2a2 + c6Q1Q3 + c7S2b2 + c8Q13 + c9Q1S2a2 + c10S33 + c11Q12Q3 + c12Q1S2b2 + c13Q3S2a2 + c14Q14 + c15Q12S2a2 + c16S2a4 + c17Q1S33 + c18Q13Q3 + c19Q12S2b2 + c20Q1Q3S2a2 + c21Q3S33 + c22S2a2S2b2 + c23Q15 + c24Q13S2a2 + c25Q1S2a4 + c26Q12S33 + c27S2a2S33 + c28Q14Q3 + c29Q13S2b2 + c30Q12Q3S2a2 + c31Q1Q3S33 + c32Q1S2a2S2b2 + c33Q3S2a4 + c34S2b2S33 (21)

Figure 5. Dispersion coefficients for the atom-diatom asymptotic channels as a function of the corresponding internuclear distance of the diatom.

Figure 6. Reference geometries used in the present work for the threebody EHF part of the potential energy surface. See text for details.

where S2a2 ) Q22 + Q32, S2b2 ) Q22 - Q32 and S33 ) Q33 3Q22Q3, with Qi (i ) 1-3) being symmetry cordinates that are defined,55,79 for the jth polynomial, by

()

Q1 Q2 ) Q3

(

13  13  13 0  12 - 12  23 - 16 - 16

)

( ) R1 - Rj,ref 1

R2 - Rj,ref 2 R3 -

(22)

Rj,ref 3

where the superscript “ref” means reference geometry. We should emphasize that although we have utilized a physically motivated formalism to represent the dynamical correlation, the resulting DMBE potential energy surface is purely ab initio except for the fact that the diatomic curves have been scaled slightly to reproduce the corresponding spectroscopic well depths (see section 2). Figure 6 displays the reference geometries that define the displacement coordinates used to write the EHF-type polynomiare obtained als in eq 20. As usual, the reference geometries Rj,ref i by first assuming their values to coincide with bond distances

2660

J. Phys. Chem. A, Vol. 114, No. 7, 2010

Joseph and Varandas

TABLE 1: Accumulated (acc) and Stratum (strat) Root-Mean-Square Deviations (kcal mol-1) of the DMBE Potential Energy Surface Na

energy

max. devb

N>rmsdc

rmsd

acc

strat

acc

strat

acc

strat

acc

strat

acc

strat

10 20 40 60 80 100 150 200 250 500 1000

0-10 10-20 20-40 40-60 60-80 80-100 100-150 150-200 200-250 250-500 500-1000

115 180 425 493 540 616 897 1529 1645 1705 1833

115 65 245 68 47 76 281 632 116 60 128

0.218 0.670 2.267 3.268 3.268 3.268 3.288 4.124 4.124 4.124 4.124

0.218 0.670 2.267 3.268 3.034 2.206 3.288 4.124 3.293 2.665 3.639

0.135 0.151 0.321 0.490 0.539 0.594 0.830 0.903 0.901 0.910 0.945

0.135 0.177 0.402 0.947 0.905 0.896 0.994 0.998 0.871 1.136 1.322

51 53 42 79 97 126 221 382 416 446 459

51 8 36 18 11 25 89 160 33 27 32

a Number of calculated DMBE-SEC points up to the indicated energy range. b Maximum deviation up to indicated energy range. c Number of calculated DMBE-SEC points with an energy deviation larger than the root-mean-square deviation (rmsd).

TABLE 2: Stationary Points Properties of the C2H(12A′) Molecule property

ab initioa

DMBEb

DMBEc

others

exp

global minimum

RCC/a0 RCH/a0 θCCH/deg ω1 (C-H(str))/cm-1 ω2 (bend)/cm-1 ω3 ((C-C)(str)/cm-1

2.291 2.011 180.0 3444 411 2010

2.307 2.007 180.0

2.285 (2.291d) 2.009 (2.011d) 180.0 (180.0d) 3436 398 1932

2.312,e 2.309,f 2.298,g 2.269,h 2.283i 2.032,e, 2.029,f 1. 999,g 1.995,h 2.005-2.009i 180.0,f 180.0,g 180.0h 3346,e 3353,f 3600,g 3521,h 3328,k 3550,l 3458m 570,e 489,f 549,g 512,h 275,k 370,l 452m 1975,e 1975,f 1848,g 2109,h 2017,k 1875,l 2069m

2.298j 1.969j 180.0j 3612,n,s 3299o,s 375,p,s 372q,s 1848,n,s 1840r,s

T-shaped minimum

R1/a0 R2/a0 R3/a0 ω1/cm-1 ω2/cm-1 ω3/cm-1

2. 424t 2. 406t 2. 406t 2417t 461t 1768t

2.423 (2.422d,t) 2.408 (2.404d,t) 2.408 (2.404d, t) 2085 (2410d,t) 525 (516d,t) 1748 (1754d,t)

R1/a0 R2/a0 R3/a0 ω1/cm-1 ω2/cm-1 ω3/cm-1

2. 486t 2. 747t 2. 256t 2564t 685it 1684t

2.422 (2.421d,t) 2.705 (2.704d,t) 2.256 (2.255d,t) 2326 (2590d,t) 563i (329id,t) 1748 (2227d,t)

feature

saddle point linking T-shaped and global minima

a This work, MRCI/AVQZ. b Reference 51. c This work. d DMBE-SEC. e Reference 86, CASSCF. f Reference 50. CASSCF/VTZ. g Reference 84. HF/6-31G**. h Reference 87. EOMIP/VTZ i Reference 43. Best estimate from CCSD(T), CCSDT, MR-AQCC, and full CI calculations. j Reference 79. k Reference 88. l Reference 45. m Reference 49. CCSD(T)/VQZ n Reference 19-21. o Reference 25. p Reference 85. q Reference 29. r Reference 26. s Fundamental vibrational frequencies ν1, ν2, ν3. t Values obtained from a local polynomial fit; see text for details.

of the associated stationary points. Subsequently, we relax this condition via a trial-and-error least-squares fitting procedure. Similarly, the nonlinear range-determining parameters γji have been optimized in this way. The complete set of parameters amounts to a total of 183 linear coefficients ci, 18 nonlinear coefficients γji, and 6 reference geometries Rj,ref i . Their optimal numerical values are collected in Tables 3 and 4 of the Supporting Information. A total of 1830 points has been used for the calibration procedure, with the energies covering a range up to 1000 kcalmol-1 above the C2H global minimum. Table 1 shows the stratified root-mean-squared deviation (rmsd) values of the final potential energy surface with respect to all the fitted ab initio energies. As shown, the fit shows a maximum rmsd of 0.94 kcal mol-1 up to the highest repulsive energy region. 4. Features of the DMBE Potential Energy Surface As expected, the DMBE potential energy surface obtained in this work displays a linear equilibrium geometry for 12A′ state of C2H. Moreover, it shows no barrier for the reaction C(3P) + CH(2Π) f C2H, and an isomerization reaction path with a T-shaped minimum at the top. Before and after this minimum there are saddle points (transition states) connecting such a T-shaped minimum with the two equivalent global minima.

The bond lengths of the global minimum are R(C-C) ) 2.2850a0, R(C-H) ) 2.0088a0 and the corresponding experimental80 values are 2.2997a0 and 1.9691a0. These calculated DMBE values show good agreement with the values from previously reported PES51 of 2.3073a0 and 2.0069a0. Recently, in a series of computations to determine the equilibrium geometry of the ethynyl radical, Szalay et al.43 give as their best theoretical estimates the values of 2.2827a0 for R(C-C) and 2.0049-2.0087a0 for R(C-H) as obtained from CCSD(T), CCSDT, MR-AQCC, and full CI calculations with basis sets up to core-polarized pentuple-ζ quality. From their mixed theoretical-experimental approach, they obtained43 empirical equilibrium geometrical parameters of 2.2808a0 and 2.020a0 for the C-C and C-H bond lengths, respectively. In turn, Sherrill and co-workers42 gave values of 2.2676a0 and 2.0031a0 from CCSD(T)/VXZ (X ) D, T, Q) levels of theory. Moreover, Santiago et al.81 performed B3LYP/6-31G(d) calculations, having obtained equilibrium bond distances of 2.3173a0 and 2.0220a0 in the above order. Our ab initio calculations predict a linear equilibrium geometry with bond lengths of 2.2912a0 and 2.0111a0, thus differing by 0.0085a0 and 0.0062a0 from the theoretical values of Szalay et al., possibly the most accurate result reported thus far. Table 2 shows that the harmonic frequencies of the DMBE PES are in good agreement with the

Potential Energy Surface for C2H

Figure 7. Isomerization minimum-energy path of C2H molecule as a function of angle θ calculated at AVQZ level of theory. The solid circles correspond to the 12A′′ state, and circles denote the 12A’ state; the global minimum of the 12A′ energy is taken as the zero of energy.

Figure 8. Isomerization minimum-energy path of C2H molecule as a function of angle θ, which is calculated at different levels of theory (AVTZ, AVQZ, AV5Z). The zero of energy is taken as the bottom of the AV5Z PES curve.

calculated raw MRCI/AVQZ values of 3444, 411, and 2010 cm-1 (CH stretch, CC stretch, CCH bending, respectively), with the corresponding DMBE-SEC values being 3436, 398, and 1932 cm-1. Szalay et al.43 have also calculated harmonic vibrational frequencies of the C2H molecule at the CCSD(T)/ VXZ (X ) T, Q) level, using both UHF and ROHF reference functions. The UHF-CCSD(T)/VQZ frequencies are 3458, 452, and 2069 cm-1, while the ROHF-CCSD(T)/VQZ ones are 3452, 381, and 2027 cm-1. Our ab initio results are consistent with their predictions, with the differences being (14, 41, 59) and (8, 30, 17) cm-1 in the above order. In addition, the DMBE PES predicts values of 5.045 and 7.821 eV for the dissociation energies De(CC-H) and De(C-CH), in very good agreement with the corresponding experimental values of 5.182 and 7.88 eV.83 The agreement with other reported values51 of 4.93 and 7.51 eV is also quite satisfactory. Ab initio calculations of the isomerization minimum-energy paths were carried out for the 12A′ and 12A′′ states of C2H (see Figure 7). For this, we have utilized Jacobi coordinates, {R, r, θ}, where r is the C-C distance, R the distance from the H atom to the middle of the C-C axis, and θ is the angle between the directions defined by r and R. It is clear that these two electronic states possess similar topographies, except when θ approaches 90° (see the insets of Figures 7). For each surface,

J. Phys. Chem. A, Vol. 114, No. 7, 2010 2661

Figure 9. DMBE PES curve for isomerization minimum-energy path of C2H molecule as a function of angle θ. The circles correspond to the DMBE-SEC points calculated at the present level of theory.

Figure 10. Contour plot for a C2V insertion of an H atom into diatomic C2. Contours are equally spaced by 0.01Eh, starting at -0.39Eh.

Figure 11. Contour plot for bond stretching C-C-H in linear configuration. Contours are equally spaced by 0.015Eh, starting at -0.41Eh.

there are two large potential wells corresponding to the linear structures (for both cases, the global minima) H-C-C and C-C-H, and a C2V saddle point (T-shaped, as also noted by Boye´-Pe´ronne et al.84 and Tarroni and co-workers47) separating

2662

J. Phys. Chem. A, Vol. 114, No. 7, 2010

Joseph and Varandas

Figure 12. Contour plot for an H atom moving around a fixed C2 diatomic in equilibrium geometry R1 ) 2.358a0, which lies along the X-axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.0045Eh, starting at -0.415Eh. The dashed lies are contours equally spaced by 0.0045Eh, starting at -0.2356Eh.

Figure 13. Contour plot for an H atom moving around a fixed CH diatomic RCH ) 2.12a0, which lies along the X-axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.05Eh, starting at -0.38Eh. The dashed lines are contours equally spaced by 0.005Eh, starting at -0.1336Eh.

the wells. The energy gap between the surfaces is almost always below 1 eV for all nuclear configurations. The DMBE surface predicts a C2V T-shaped minimum in its isomerization reaction path and a transition state that links the T-shaped and global minima (see Figures 8 and 9 for MRCI/AVTZ, AVQZ, and AV5Z calculations as well as the curve from the DMBE PES). Note that we have calculated a denser grid of energies in the vicinity of the stationary point. In fact, about 95 ab initio grid points have been calculated for the local fit, over a region defined by 2.3 e R2/a0 e 2.5, 2.3 e R3/a0 e 2. 5, with the angle ∠CHC varying 55 0.0 e θ/deg e 65.0. These points were then fitted to a local polynomial with a root-mean squared deviation less

than 2 cm-1. Both the local polynomial and our DMBE surface show a T-shaped minimum. This C2V stationary point has been properly characterized. We found that this C2V minimum lies 0.949 69 eV above the global minimum of C2H, with equilibrium bond lengths of R1 ) 2.423a0 and R2 ) R3 ) 2.408a0, and frequencies of ω1(C-H) ) 2085 cm-1 (2410), ω2(C-H) ) 525 cm-1 (516), and ω3(C-C) ) 1748 cm-1 (1754) (the values in parentheses are obtained from the local polynomial fit where the rmsd is