Accurate global potential energy surface for the ground state of CH2 + by extrapolation to the complete basis set limit.

A full three-dimensional global potential energy surface is reported for the ground state of CH2 + by fitting accurate multireference configuration interaction energies calculated using aug-cc-pVQZ and aug-cc-pV5Z basis sets with extrapolation of the electron correlation energy to the complete basis set limit. The topographical characteristics have been compared in detail with a potential energy surface of the same type recently reported [J. Chem. Phys., 2015, 142, 124302] based on a least-squares fit to accurate high level ab initio MRCI(Q) energies, calculated using AV6Z basis set. The new three-dimensional global potential energy surface is then used in quasiclassical trajectory calculations for H(2S) + CH+(X 1Σ+) → C+(2P) + H2(X 1Σg +) reaction. The integral cross sections, differential cross sections and the rate coefficients have been computed. A comparison shows that our potential energy surface can be applied to any type of dynamic study. This journal is © The Royal Society of Chemistry.

Introduction

The C+ + H2, ion-molecule reaction has been the research core of extensive experimental and theoretical study owing to its important role in astrophysics and atmospheric chemistry. Particularly, the CH2+ complex formed by the reaction is a crucial intermediate in interstellar matter and combustion process.[1] Hence the reaction C+ + H2 has been widely researched.

The cation CH2+ belongs to a class of molecules which are termed ‘‘quasilinear’’. From the spectroscopic side, the linear and bending problems of CH2+ attracted the researchers' interest in the 1960s.[2] This issue was first solved theoretically due to the lack of spectral data, Schaefer and Bender[3] predicted a bent equilibrium geometry for the electronic ground state (12A′) of CH2+ in 1971. Later, the bent equilibrium geometry was confirmed by the Coulomb explosion imaging experiment.[4,5] Then, extensive research on CH2+ has been carried out theoretically.

Among a series of theoretical studies, Stoecklin and Halvick[6] reported the theoretical research result of the caption reaction for the first time, which fitting a precise 3-D single-valued potential energy surface (PES). Utilizing an extensive multiconfigurational wave function with the augmented Dunning correlation consistent basis set (aug-cc-pVQZ) to calculate large ab initio points. Then Halvick et al.[7] based on this PES to study quasiclassical trajectory (QCT) calculations, along with phase space theory and quantum rigid rotor calculations for H + CH+. And the PES was used by Zanchet et al.[8] for the state-to-state rate constant study with quantum wave packet analysis of the H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) reaction. Recently, Schneider and Warmbier[9] constructed the CH2+ ground state PES by fitting internuclear distances polynomials with the multireference configuration interaction (MRCI) and aug-cc-pVTZ ab initio points. And for verifying this PES, they have performed quantum scattering and QCT calculations. Additionally, Herráez–Aguilar et al.[10] have performed a dynamical study of the endothermic and barrierless C+ + H2(1Σg+) → CH+(1Σg+) + H reaction for different initial rotational states of the H2(v = 0) and H2(v = 1) manifolds, with the QCT and the Gaussian binning methodology on Schneider and Warmbier's[9] PES. Most recently, Li et al.[11] reported a new many-body expansion (MBE) PES by fitting MRCI/aug-cc-pV6Z ab initio energies. This PES is used by Guo et al.[12] to analyze the effect of isotopic substitution on three-dimensional dynamic properties of the reactions C+ + H2/HD/HT → CH+ + H/D/T. In addition, the time-dependent wave packet propagation approach was used to compute thermal rate constants and integral cross sections of the H + CH+ reaction in the coupled states approximation by Sundaram et al.[13] Faure et al. also presents a detailed theoretical study of state-to-state and initial-state-specific rate coefficients are computed in the kinetic temperature range 10–3000 K.[14]

The purpose of present study is to build a high quality global PES for the ground state(12A′) of CH2+ from MRCI(Q)[15]ab initio energies based on the reference full valence complete active space (FVCAS)[16] wave function, the aug-cc-pV5Z(AV5Z) and aug-cc-pVQZ(AVQZ) basis sets of Dunning[17,18] have been applied. We all know that in order to get a highly precise PES, it usually requires a large basis sets. But, we did not use such a large basis sets in this work, instead we have extrapolated the total energy to the complete basis set (CBS) limit by using a uniform single-pair and triple-pair (USTE).[19,20] For verifying this PES, the QCT has been performed on H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+), the differential cross sections (DCSs), integral cross sections (ICSs) and the rate coefficients will be computed.

The paper is organized as follows. Section 2 introduces the theoretical approaches, such as ab initio calculations and the application of extrapolation. Section 3 introduces the analytic expression of CH2+(12A′) PES. Main topographical features are discussed in section 4. Section 5 describes the QCT calculations. Finally, the conclusion is presented in section 6.

Ab initio calculations and extrapolation scheme

The MRCI(Q)[15] approach is one of the best methods to obtain the precise PESs. All ab initio calculations are performed at the MRCI(Q) level with the FVCAS[16] as reference. MOLPRO 2012 (ref. 21) is a kind of program package about the quantum chemistry, in association with the Dunning et al.[17,18] correlation-consistent basis sets have been applied during our work. This procedure involves 6 active orbitals (5A′ + 1A′′), with a total of 443 (166A′ + 177A′′) configuration state functions at AV5Z and AVQZ basis sets, respectively. In all 3255 ab initio grid points have been computed for C+–H2 channels, the region was defined by 1.2 ≤ RH/a0 ≤ 4.4, 1.4 ≤ rC/a0 ≤ 10 and 0 ≤ γ/deg ≤ 90, while, for H–CH+, they cover geometries defined by 1.8 ≤ RCH/a0 ≤ 3.6, 1.4 ≤ rH–CH/a0 ≤ 10 and 0 ≤ γ/deg ≤ 180, R, r and γ are atom-diatom Jacobi coordinates for both channels. To get more precise energy points, the USTE[19,20] method is adopted. During the calculations, the core is frozen and ignoring the relativistic effect.

In order to carry out the extrapolation, electronic energy in the MRCI(Q) calculation is expressed by a sum of two terms[19]where the superscript CAS represents the complete-active space and the superscript dc represents the dynamical correlation energies, in addition the subscript X signifies that the electronic energy computed in the AVXZ basis set. The X = Q, 5 are used during the calculation.

Using a two-point extrapolation program suggested by Karton and Martin(KM),[22] the CAS energies are extrapolated to the CBS limitwhere ECASx is the energy when X → ∞ and α = 5.34 is an effective decay index.

The USTE protocol[19,20] has been triumphantly implemented to extrapolate the dynamical correlation energies in MRCI(Q) calculations, which is extrapolated by the formulawhere A5 is written as the auxiliary relation

With α = −3/8, c = −1.17847713 and A5(0) = 0.0037685459 are “universal-like” parameters.[19]Eqn (3) could be converted to (E∞, A3) two-parameter rule, which has access to the actual extrapolation process.

Analytical potential energy function of CH2+(12A′)

The analytical function of CH2+(12A′) PES can be represented as a MBE[23,24] formwhere VA(1) is the isolated atomic energy, VAB(2) is a two-body energy term and VABC(3) is the three-body energy term which is zero at all dissociation limits. In this work, the title system obeys the following dissociation scheme:where C+(2P) and H(2S) are all in their ground states. So, the one-body energy term VA(1) in eqn (5) are zero.

Two-body energy terms

The analytic energy function of the two-body terms VAB(2) for CH+(X1Σ+) and H2(X1Σg+) are imitated employing the Aguado and Paniagua[25,26] approach, which the function for title diatomic systems can be represented as summation of the short-range and long-range potentialswherewhich the potential energy of diatomic tends to infinitely great when RAB → 0. The long-range potential is expressed aswhich the potential energy of diatom is equal to zero as RAB → ∞. The potential function in eqn (9) is truncated up to 8th power (n = 8), we can obtain 11 parameters, including 2 nonlinear parameters β(i = 1, 2) and 9 linear parameters a(i = 0, 1,…, 8) for CH+(X1Σ+) and H2(X1Σg+) by the fitting procedure. All the fitted parameters of the CH+(X1Σ+) and H2(X1Σg+) are listed in Table 1 of the ESI.†

Three-body energy term

For the calculation of three-body term, the method of three-body distributed polynomial is adopted,[27,28] which was applied to calculate ground and excited states of FH2,[29] NH2 (ref. 30–32) and NH3.[33,34]where P((Q1,Q2,Q3) is the j-th (j = 1, 2, 3) polynomial, with Q (i = 1, 2, 3) being the symmetric coordinates, γ( are determining parameters of the nonlinear range and R( represent reference geometries. So, there are 150 linear coefficients, 9 reference bond distances and 9 nonlinear ones in all. In this work, a total of 3255 CBS points has been calculated in fitting procedure and fitting result shows that the total root mean square derivation is rmsd = 0.55 kcal mol−1. All the fitted parameters of the least square method are listed in Tables 2 and 3 of the ESI.†

Features of the PES

Table 1 presents the results of Re, De, ωe, ωeχe, αe and Be of CH+(X1Σ+) and H2(X1Σg+) together with the other theoretical[11,35-37,40-42] and experimental[38-42] data. The equilibrium internuclear Re of CH+(X1Σ+) to be 2.135a0, which is only 0.002a0 shorter than the experimental result[38] and 0.001a0 shorter than the latest result computed by Li et al.,[11] while for the case of H2, all the results are 1.401a0,[40,42−44] it shows a high precision. For the dissociation energies De, CH+(X1Σ+) differs from the theoretical and experimental results by less than 0.005 eV.[11,39] For the H2 (X1Σg+) has dissociation energy De = 4.749 eV, this compares well with the corresponding theoretical values are De = 4.748 eV (ref. 11) and corresponding experimental values De = 4.478 eV.[45] Overall, it can be concluded that the other spectroscopic constants are in good agreement with these literature results. Fig. 1 shows the potential energy curves (PECs) of CH+(X1Σ+) and H2(X1Σg+). In order to evaluate the quality of the fitting, we calculate the rmsd. The rmsd of CH+(X1Σ+) and H2(X1Σg+) PECs are 0.05 kcal mol−1 and 0.005 kcal mol−1, respectively. As a whole, it reveals a high quality fitting process. Shown in Fig. 1, the PECs at the CBS/USTE(Q,5) calculations nicely, showing accurate and smooth behavior both in short and long regions.

Spectroscopic constants of CH+ and H2 diatoms, with the unit of Re in a0, De in eV and ωe, ωeχe, αe and βe in cm−1

	R e	D e	ω e	ω e χ e	α e	β e
CH + (X 1 Σ + )
This work	2.135	4.257	2861.948	59.629	0.447	14.311
Theor.[11]	2.136	4.252	2853.027	58.515	0.489	14.201
Theor.[35]	2.136	4.244	2851	58.1	0.489	14.199
Theor.[36]	2.144	—	2849.03	66.448	0.490	14.094
Theor.[37]	2.127	4.660	3111.0	38.44	—	—
Exp.[38]	2.137	4.26 [39]	2857.56	59.32	0.495	14.178
H 2 (X 1 Σ g + )
This work	1.401	4.749	4404.613	126.636	2.233	60.861
Theor.[40]	1.401	4.748	4403.60	126.602	2.232	60.864
Theor.[41]	1.403	4.748	4395.22	126.119	2.221	60.735
Theor.[42]	1.401	4.711	4389.66	121.560	3.162	60.826
Exp.[43]	1.401	4.478 [45]	4401.21	121.33	3.062	60.853
Exp.[44]	1.401	4.476	4395.20	117.99	2.993	60.809

Fig. 1

Potential energy curves of CH+(X1Σ+) and H2(X1Σg+). The circles indicate the CBS(Q,5) energies.

Table 2 collects all known stationary points (geometry, energy and vibrational frequencies) of the new PES for ground state of CH2+. For better comparisons, the results of other theoretical and experimental work reports are also gathered in Table 2. A global minimum (GM) is found to local at θHC = 138.6°, R1 = 3.859a0 and R2 = R3 = 2.063a0 with the minimum energy −0.333Eh relative to the all dissociation asymptote (R1 represents the H2 interatomic separation, while R2 and R3 represent the two CH+ interatomic separations), its maximum deviation is only 0.004a0 for the CH+ bond length (R2 and R3) when compared to the data of the MRCI(Q)/AV6Z PES[11] with the θHC = 141.0°, R1 = 3.896a0 and R2 = R3 = 2.067a0. Comparing with the theoretical works of Stoecklin and Halvick,[6]R1 = 3.865a0 is 0.006a0 longer than our results. This compares well with the corresponding experimental values[2,46] of θHC = 139.8°. For the harmonic frequencies, the PES from this work computes values of 3048 cm−1, 3272 cm−1 and 921 cm−1, there is good consistency with the theoretical results calculated by Brinkmann et al.[47], which are 3011 cm−1, 3260 cm−1 and 965 cm−1, respectively. Table 2 also collects the attributes of a local minimum (LM) and three transition states: TS1(D∞h), TS2(C2v) and TS3(C∞v) barriers.

Stationary points at the valence region for ground state of CH2+ PES, with the unit of R1, R2 and R3 in a0, θHCH in deg, ΔV in kcal mol−1, ω1, ω2 and ω3 in cm−1

Feature	R 1	R 2	R 3	θ HCH	ΔVa	ω 1	ω 2	ω 3
Global minimum
This work	3.859	2.063	2.063	138.6	−99.64	3048	3272	921
Theor.[11]	3.896	2.067	2.067	141.0	−99.54	2983	3283	957
Theor.[47]	3.890	2.067	2.067	140.39	—	3011	3260	965
Theor.[6]	3.865	2.075	2.075	137.3	—	—	—	—
Exp.[2,46]	3.922	2.088	2.088	139.8	—	—	3133.4	—
Local minimum
LM(C2v)	1.648	2.757	2.757	34.8	−29.63	3254	1317	1220
Theor.[11]	1.666	2.645	2.645	36.7	−29.28	2937	1444	1307
Transition state
TS1(D∞h)	4.099	2.050	2.050	180.0	−97.13	3083	5489	1077i
Theor.[11]	4.127	2.064	2.064	180.0	−96.54	2939	5431	1319i
TS2(C2v)	2.106	2.362	2.362	52.9	−21.44	2265i	1766	1493
Theor.[11]	2.170	2.435	2.435	52.9	−23.41	2149i	2179	1745
TS3(C∞v)	1.525	2.660	4.185	0.0	−15.96	3557	664i	1527
Theor.[11]	1.511	2.645	4.156	0.0	−16.89	3676	172i	1769

Relative to the C+ + H2 asymptote.

Fig. 2 and 3 show the main topographical characteristics of the new CH2+ PES computed in this work. Obviously, there is a correct and smooth behavior in the entire configuration space. The salient features of these contour maps corresponding to several important stationary points for the main reaction. Fig. 2(a) illustrates a contour map for linear [H–C–H]+ stretch. The significant characteristic of this map is that there is a TS1(D∞h) linear transition state at R2 = R3 = 2.050a0 with an energy of 878 cm−1 above the GM of CH2+ but still 33 972 cm−1 below the energy of the C+ + H2 asymptote. This compares well with the MRCI(Q)/AV6Z PES,[11] where the transition state is computed to local at R2 = R3 = 2.064a0 with an energy of 1050 cm−1 above the GM and 33 765 cm−1 below the reactants asymptote. The corresponding infrared spectrum[48] result for this linear barrier is 1089 cm−1.

Fig. 2

(a) Contour plot for bond stretching for linear [H–C–H]+ configurations. (b) Contour plot for bond stretching in [H–C–H]+, keeping the included angle fixed at 138.6°. (c) Contour plot for the C2v insertion of the C+ into H2. (d) Contour plot for bond stretching in linear [H–H–C]+ configurations. Contours are equally spaced by 0.02Eh, starting at −0.329Eh for panels (a), −0.333Eh for panels (b), −0.332Eh for panels (c) and −0.199Eh for panels (d).

Fig. 3

(a) Contour plot for a C+ moving around a H2 molecule fixed at the equilibrium geometry RH = 1.401a0 and lying along the X-axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.002Eh, starting at −0.213Eh. Shown in dash are contours equally spaced by −0.0001Eh, starting at −0.174508Eh. (b) Contour plot for a H atom moving around a CH+ fixed at the equilibrium geometry RCH = 2.135a0 and lying along the X-axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.005Eh, starting at −0.33Eh. Shown in dash are contours equally spaced by −0.0002Eh, starting at −0.156468Eh.

Fig. 2(b) plots for the bond stretching in [H–C–H]+ which the angle is fixed at138.6°. It can be found from Fig. 2(b) that there is a deep well for CH2+ PES, which is GM. Fig. 2(c) shows the contour plots to the insertion of C+ + H2 reaction. In this figure, the stationary points which are corresponding to TS1(D∞h), TS2(C2v), the LM and the GM. As shown in Table 2, the LM is predicted to locate at R1 = 1.648a0 and R2 = R3 = 2.757a0, so agreeing with MRCI(Q)/AV6Z PES.[11] The main characteristics of the new PES for collinear [H–H–C]+ stretch are shown in the contour map of Fig. 2(d). The collinear TS3(C∞v) is found to locate at R1 = 1.525a0 and R2 = 2.660a0 with the energy of 15.96 kcal mol−1. This compares well with the R1 = 1.511a0 and R2 = 2.645a0 and 16.89 kcal mol−1 for the MRCI(Q)/AV6Z PES.[11]

Panel (a) of Fig. 3 illustrates the contour plot of C+ atom moving around H2(X1Σg+) which the bond length at equilibrium geometry RH = 1.401a0. Diatoms follow the X-axis centered on the origin. In addition, the corresponding contour plot of H atom moving around a fixed CH+(X1Σ+) is shown in (b) of the same figure, which the bond length is fixed at equilibrium geometry RCH = 2.135a0, which is in very good agreement with the MRCI(Q)/AV6Z PES.[11] The two plots clearly show that there is a smooth behavior both in the long and short range.

All the main topographical characteristics of CH2+ PES can be also viewed in a relaxed triangular plot[49] using scaled hyper-spherical coordinates (γ* = γ/Q and β* = β/Q), where the Q, γ and β are written as

Clearly visible in Fig. 4 are all stationary points discussed above, which correspond to a GM, a LM and three transition states: TS1(D∞h), TS2(C2v) and TS3(C∞v) barriers.

Fig. 4

Relaxed triangular plot of the new PES of the present work in hyperspherical coordinates. Contours are equally spaced by 0.003Eh, starting at −0.332Eh. Also indicated are the various stationary points.

Fig. 5 shows the minimum energy paths (MEPs) for H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) reaction obtained from both the new PES and MRCI(Q)/AV6Z PES[11] for the collinear configuration ∠HC+H = 180°. The MEPs indicate the potential energy of CH2+ as a function for corresponding reaction coordinate of RCH–RH, RCH and RH as the internuclear distance between C+–H and H–H, respectively. As shown Fig. 5, there is a little barrier connecting a deep well and a shallow well. For the new PES, the relatively deeper well is found to locate at RCH = 2.660a0 and RH = 1.525a0 and the little barrier locates at RCH = 2.210a0 and RH = 2.972a0. The well depth and the barrier height are computed to be 1.184 eV and 0.016 eV. Comparing with the MRCI(Q)/AV6Z PES,[11] the relatively deeper well locates at RCH = 2.645a0 and RH = 1.511a0 and the little barrier is found to locate at RCH = 2.196a0 and RH = 3.044a0. The well depth and the barrier height are computed to be 1.223 eV and 0.006 eV. Moreover, the reaction H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) is exoergic by 0.49 eV base on the both PESs. It can be seen from Fig. 5, the results of the two PESs are in good agreement.

Fig. 5

Minimum energy paths for H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) reaction obtained from both the new PES and MRCI(Q)/AV6Z PES[11] as a function of RCH–RH. Collinear configuration ∠HC+H = 180°.

Dynamics of H + CH+ reaction

On the new PES, QCT[50,51] calculation was performed for the H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) reaction. In this work, we computed the ICSs, DCSs and rate constants. A total of 10 000 trajectories have been run for each of the collision energy. The time integration step is chosen to be 0.1 fs of classical motion equations, with H atom and the center of mass of the CH+ initially separated by 15.0 Å. The ICS is then written aswhere bmax is the maximum impact parameter, Nr is the number of trajectories that go into a certain reaction channel and Nt is the total number of trajectories.

As shown Fig. 6, the ICSs are expressed as a collision energy function. For comparison our results with the quantum wave packet calculations[11] and a modified version of the ABC[52] quantum scattering code method have also been performed for the same reaction. We can find that the ICS based on our surface smaller than the results based on MRCI(Q)/AV6Z PES[11] when the collision energy is less than about 40 meV. But when the collision energy is larger than about 40 meV, our results are consistency with quantum wave packet results. Overall, our results are reasonably good consistency with previous results.[11,52]

Fig. 6

Integral cross sections of reaction of H(2S) + CH+(X1Σ+) → C+(2P)+ H2(X1Σg+) as a function of collision energy. The QM and ABC from Li[11] and Schneider[9] are shown as well.

DCS is mainly used to study product and reagent relative velocity , which is the most common vector correlation. The global angular distributions of H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) reaction at collision energies of 10, 20, 30 and 40 kcal mol−1 based on the new PES are shown in Fig. 7. We can find that with the collision energies increase, the backward scattering phenomenon becomes more and more obvious. The enhanced phenomenon of backward scattering may be caused by an insertion reaction mechanism proposed by Pino et al.[53]

Fig. 7

Differential cross section as a function of the scattering angle θt for the title reaction at the collision energy of 10, 20, 30 and 40 kcal mol−1.

Finally, rate constants for the reaction H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+) are computed over the temperature range 10–1000 K by running QCT on the new PES of this work. By supposing a Maxwell–Boltzmann distribution on the collision energies, the rate constant is written as[54]where ge is the electronic degeneracy factor, we adopt ge = 1 in the present work. By comparing with the results of various theoretical[7,9,11] and experimental[55,56] studies, the results are shown in the Fig. 8. It can be seen from Fig. 8 that our result (QCT) is less than others work for the whole temperature range. The QCT calculation is defect, the error will be larger when the collision energy is low, so the error of rate constant is also very large when the temperature is very low. However, as the temperature increases, the images get closer and closer to the other results, especially at high temperatures it agrees well with Li et al.[11] and Federer et al.[56] So, it turns out that our new PES can be applied to any type of dynamic study.

Fig. 8

Rate constant for H(2S)+ CH+(X1Σ+) → C+(2P)+ H2(X1Σg+) reaction calculated in this work. Other theoretical works of Li,[11] Warmbier and Schneider[9] and Halvick et al.,[7] as well as experimental data from Luca et al.[55] and Federer et al.[56] are also included.

Conclusions

In our research process, we have constructed a high quality global PES for the ground state of CH2+(12A′) from MRCI ab initio energies based on the reference FVCAS wave function, both AVQZ and AV5Z basis sets subsequently extrapolated to the CBS limit. All known stationary points including geometries, energies and vibrational frequencies can be obtained, and all the results are consistency with the corresponding theoretical and experimental values. The consistency and accuracy of the CBS method have also been affirmed by comparing the MRCI(Q)/AV6Z PES.[11] Finally, QCT calculation has been performed on H(2S) + CH+(X1Σ+) → C+(2P) + H2(X1Σg+), the ICS, DCS and the rate coefficients are computed in detail and compared to the MRCI(Q)/AV6Z and other PESs, as well as experimental values in the literature. In summary, the new PES built here can be used for any type of dynamic study.

Conflicts of interest

There are no conflicts to declare.