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+(X1Σ+) → C+(2P) + H2(X1Σ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.


Introduction
The C + + H 2 , 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 CH 2 + complex formed by the reaction is a crucial intermediate in interstellar matter and combustion process. 1 Hence the reaction C + + H 2 has been widely researched. The cation CH 2 + belongs to a class of molecules which are termed ''quasilinear''. From the spectroscopic side, the linear and bending problems of CH 2 + attracted the researchers' interest in the 1960s. 2 This issue was rst solved theoretically due to the lack of spectral data, Schaefer and Bender 3 predicted a bent equilibrium geometry for the electronic ground state (1 2 A 0 ) of CH 2 + in 1971. Later, the bent equilibrium geometry was conrmed by the Coulomb explosion imaging experiment. 4,5 Then, extensive research on CH 2 + 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 rst time, which tting a precise 3-D single-valued potential energy surface (PES). Utilizing an extensive multi-congurational 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( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S g + ) reaction. Recently, Schneider and Warmbier 9 constructed the CH 2 + ground state PES by tting internuclear distances polynomials with the multireference conguration 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 + + H 2 ( 1 S g + ) / CH + ( 1 S g + ) + H reaction for different initial rotational states of the H 2 (v ¼ 0) and H 2 (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 tting 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 + + H 2 /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-specic 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(1 2 A 0 ) of CH 2 + 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 singlepair and triple-pair (USTE). 19,20 For verifying this PES, the QCT has been performed on H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S 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 CH 2 + (1 2 A 0 ) 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 correlationconsistent basis sets have been applied during our work. This procedure involves 6 active orbitals (5A 0 + 1A 00 ), with a total of 443 (166A 0 + 177A 00 ) conguration state functions at AV5Z and AVQZ basis sets, respectively. In all 3255 ab initio grid points have been computed for C + -H 2 channels, the region was dened by 1.2 # R H 2 /a 0 # 4.4, 1.4 # r C + -H 2 /a 0 # 10 and 0 # g/deg # 90, while, for H-CH + , they cover geometries dened by 1.8 # R CH + /a 0 # 3.6, 1.4 # r H-CH + /a 0 # 10 and 0 # g/deg # 180, R, r and g 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 signies 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 limit where E CAS x is the energy when X / N and a ¼ 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 formula where A 5 is written as the auxiliary relation With a ¼ À3/8, c ¼ À1.17847713 and A 5 (0) ¼ 0.0037685459 are "universal-like" parameters. 19 Eqn (3) could be converted to (E N , A 3 ) two-parameter rule, which has access to the actual extrapolation process.
3. Analytical potential energy function of CH 2 The analytical function of CH 2 + (1 2 A 0 ) PES can be represented as where V A (1) is the isolated atomic energy, V AB (2) is a two-body energy term and V ABC (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 + ( 2 P) and H( 2 S) are all in their ground states. So, the onebody energy term V A (1) in eqn (5) are zero.

A. Two-body energy terms
The analytic energy function of the two-body terms V AB (2) for CH + (X 1 S + ) and H 2 (X 1 S 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 potentials which the potential energy of diatomic tends to innitely great when R AB / 0. The long-range potential is expressed as which the potential energy of diatom is equal to zero as R AB / N. The potential function in eqn (9) is truncated up to 8th power (n ¼ 8), we can obtain 11 parameters, including 2 nonlinear parameters b i (i ¼ 1, 2) and 9 linear parameters a i (i ¼ 0, 1,., 8) for CH + (X 1 S + ) and H 2 (X 1 S g + ) by the tting procedure.
All the tted parameters of the CH + (X 1 S + ) and H 2 (X 1 S g + ) are listed in Table 1 of the ESI. †

B. Three-body energy term
For the calculation of three-body term, the method of threebody distributed polynomial is adopted, 27,28 which was applied to calculate ground and excited states of FH 2 , 29 NH 2 (ref. [30][31][32] and NH 3 . 33,34 where 3) being the symmetric coordinates, g i (j) are determining parameters of the nonlinear range and R i,ref 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 tting procedure and tting result shows that the total root mean square derivation is rmsd ¼ 0.55 kcal mol À1 . All the tted parameters of the least square method are listed in Tables 2 and 3 of the ESI. †  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 + (X 1 S + ) and H 2 (X 1 S g + ). In order to evaluate the quality of the tting, we calculate the rmsd. The rmsd of CH + (X 1 S + ) and H 2 (X 1 S g + ) PECs are 0.05 kcal mol À1 and 0.005 kcal mol À1 , respectively. As a whole, it reveals a high quality tting 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.  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 Nh ), TS2(C 2v ) and TS3(C Nv ) barriers. Fig. 2 and 3 show the main topographical characteristics of the new CH 2 + PES computed in this work. Obviously, there is a correct and smooth behavior in the entire conguration 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 signicant characteristic of this map is that there is a TS1(D Nh ) linear transition state at R 2 ¼ R 3 ¼ 2.050a 0 with an energy of 878 cm À1 above the GM of CH 2 + but still 33 972 cm À1 below the energy of the C + + H 2 asymptote. This compares well with the MRCI(Q)/AV6Z PES, 11 where the transition state is computed to local at R 2 ¼ R 3 ¼ 2.064a 0 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(b) plots for the bond stretching in [H-C-H] + which the angle is xed at138.6 . It can be found from Fig. 2(b) that there is a deep well for CH 2 + PES, which is GM. Fig. 2(c) shows the contour plots to the insertion of C + + H 2 reaction. In this gure, the stationary points which are corresponding to TS1(D Nh ), TS2(C 2v ), the LM and the GM. As shown in Table 2, the LM is predicted to locate at R 1 ¼ 1.648a 0 and R 2 ¼ R 3 ¼ 2.757a 0 , 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 Panel (a) of Fig. 3 illustrates the contour plot of C + atom moving around H 2 (X 1 S g + ) which the bond length at equilibrium geometry R H 2 ¼ 1.401a 0 . Diatoms follow the X-axis centered on the origin. In addition, the corresponding contour plot of H atom moving around a xed CH + (X 1 S + ) is shown in (b) of the same gure, which the bond length is xed at equilibrium geometry R CH + ¼ 2.135a 0 , 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 CH 2 + PES can be also viewed in a relaxed triangular plot 49 using scaled hyperspherical coordinates (g* ¼ g/Q and b* ¼ b/Q), where the Q, g and b are written as 0 Clearly visible in Fig. 4 are all stationary points discussed above, which correspond to a GM, a LM and three transition states: TS1(D Nh ), TS2(C 2v ) and TS3(C Nv ) barriers. Fig. 5 shows the minimum energy paths (MEPs) for H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S g + ) reaction obtained from both the new PES and MRCI(Q)/AV6Z PES 11 for the collinear conguration :HC + H ¼ 180 . The MEPs indicate the potential energy of CH 2 + as a function for corresponding reaction coordinate of R CH +-R H 2 , R CH + and R H 2 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 R CH + ¼ 2.660a 0 and R H 2 ¼ 1.525a 0 and the little barrier locates at R CH + ¼ 2.210a 0 and R H 2 ¼ 2.972a 0 . 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 R CH + ¼ Fig. 1 Potential energy curves of CH + (X 1 S + ) and H 2 (X 1 S g + ).
The circles indicate the CBS(Q,5) energies. Table 2 Stationary points at the valence region for ground state of CH 2 + PES, with the unit of R 1 , R 2 and R 3 in a 0 , q HCH in deg, DV in kcal mol À1 , u 1 , u 2 and u 3 in cm À1 2.645a 0 and R H 2 ¼ 1.511a 0 and the little barrier is found to locate at R CH + ¼ 2.196a 0 and R H 2 ¼ 3.044a 0 . The well depth and the barrier height are computed to be 1.223 eV and 0.006 eV. Moreover, the reaction H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S g + )

Global minimum
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.

Dynamics of H + CH + reaction
On the new PES, QCT 50,51 calculation was performed for the H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S 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 A. The ICS is then written as where b max is the maximum impact parameter, N r is the number of trajectories that go into a certain reaction channel and N t 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 modied version of the ABC 52 quantum scattering code method have also been performed for the same reaction. We can nd 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 DCS is mainly used to study product and reagent relative velocity k-k 0 , which is the most common vector correlation. The global angular distributions of H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S 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 nd 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 Finally, rate constants for the reaction H( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S g + ) are computed over the temperature range 10-   where g e is the electronic degeneracy factor, we adopt g e ¼ 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.

Conclusions
In our research process, we have constructed a high quality global PES for the ground state of CH 2 + (1 2 A 0 ) 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( 2 S) + CH + (X 1 S + ) / C + ( 2 P) + H 2 (X 1 S 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 conicts to declare.  Rate constant for H( 2 S)+ CH + (X 1 S + ) / C + ( 2 P)+ H 2 (X 1 S 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.