An ab initio analytical potential energy surface for the O( P) S( P) reaction useful for kinetic and dynamical studies

A total of 816 ab initio points at the PUMP4/6-311G(2 d ) level were used to derive an analytical expression for the lowest 3 A (cid:56) adiabatic potential energy surface (cid:126) PES (cid:33) of the reaction O (cid:126) 3 P (cid:33) (cid:49) CS (cid:126) X 1 (cid:83) (cid:49) (cid:33) ! CO (cid:126) X 1 (cid:83) (cid:49) (cid:33) (cid:49) S (cid:126) 3 P (cid:33) . Thermal rate constants calculated using the variational transition state theory and semiclassical tunneling correction were used as a tool to determine the optimum analytical surface. This was done by comparing the calculated rate constant at room and lower temperatures with the experimental values. The best analytical surface (cid:126) PES 3 (cid:33) reproduces the rate constant at low temperatures well. However, it has not been possible to obtain an analytical PES capable of reproducing both the rate constant at 300 K and the activation energy (cid:126) 150–300 K range (cid:33) . At higher temperatures, the contribution of the lowest 3 A (cid:57) adiabatic potential energy surface to the rate constant seems to be important to reproduce the experimental data. At present, the PES 3 is the most suitable analytical surface to be used for kinetic and dynamical single surface


I. INTRODUCTION
The reaction between atomic oxygen and carbon monosulfide, O͑ 3 P ͒ϩCS͑ X 1 ⌺ ϩ ͒→CO͑ X 1 ⌺ ϩ ͒ϩS͑ 3 P ͒ ⌬H 0 0 ϭϪ86.7Ϯ6.1 kcal/mol, has been the object of considerable interest because it is the primary source of vibrationally excited CO in the CO chemical laser. This laser was first detected by photolyzing a CS 2 /O 2 mixture in a laser cavity. 1 Measurements of the CO vibrational distribution from this reaction at thermal reactant energies have been reported by several authors using a variety of techniques. [2][3][4][5][6][7][8][9][10] These experiments indicate the existence of a peak at the twelfth or thirteenth vibrational level, with levels populated up to the thermodynamical limit.
The rate constant at room temperature ͑300 K͒ has been measured by several authors with good agreement among the reported values ͓͑12.9Ϯ3.1͒ϫ10 12 cm 3 mol Ϫ1 s Ϫ1 ͔. 11 The temperature dependence of the rate constant is not as well established, as it only relies on a single experimental work that considers the 150-300 K temperature interval, 12 a range where a curvature in the Arrhenius plot was followed. On the basis of these measurements, recommended values of the preexponential factor ͓Aϭ͑16.3Ϯ6.8͒ϫ10 13 cm 3 mol Ϫ1 s Ϫ1 ͔ and the activation energy ͑E a ϭ1.5Ϯ0.5 kcal/mol͒ are reported in a review of kinetic data. 11 The relative CO rotational populations at different vibrational levels of the molecule have recently been measured, [13][14][15] using laser induced fluorescence detection, in experiments which used thermal and hot ͑hyperthermal͒ oxygen atoms. Analyses of Doppler profiles have also permitted the study of the correlations between the relative velocity vector of reactants ͑k͒, the relative velocity vector of products ͑kЈ͒, and the rotational angular momentum vector of products ͑JЈ͒, showing only a negative kЈϪJЈ angular correlation. This means that the CO molecule is produced with JЈ preferentially perpendicular to kЈ.
There are only a few theoretical studies on this reaction. In an earlier paper 16 we presented a quasiclassical trajectory ͑QCT͒ study on a semiempirical analytical potential energy surface ͑PES͒ which fitted MNDO/CI stationary points of the 3 AЈ PES. In that work a nonlinear OCS minimum and a collinear transition state ͑TS͒ were adjusted. The QCT CO vibrational distribution was in agreement with the experimental one, but the rate constant at room temperature was very low ͓͑0.4Ϯ0.1͒ϫ10 12 cm 3 mol Ϫ1 s Ϫ1 ͔ in comparison with the experimental value. In a recent paper 17 we have carried out an ab initio study on the lowest adiabatic 3 AЈ and 3 AЉ PES of this system. This work has allowed us to characterize the most relevant properties of the stationary points of both surfaces, and to select a suitable method and basis set to produce dense grids of ab initio points, accurate enough to be used in the construction of analytical representations of the lowest 3 AЈ and 3 AЉ PES. The first one of these analytical PES is presented in this work.

A. Ab initio calculations
Three adiabatic PES, one of 3 AЈ and two of 3 AЉ symmetry, are involved in this reaction. The lowest energy 3 AЈ and 3 AЉ adiabatic surfaces on the reactant's side come from avoided crossings between the diabatic surfaces arising from the O͑ 3 P͒ϩCS͑X 1 ⌺ ϩ ͒ and O͑ 3 P͒ϩCS͑a 3 ⌸͒ states, the ground and excited states of reactants, respectively. Simi-a͒ Authors to whom correspondence should be addressed. larly, on the product's side they come from the ground S͑ 3 P͒ϩCO͑X 1 ⌺ ϩ ͒ and excited S͑ 3 P͒ϩCO͑a 3 ⌸͒ states of products. Different ab initio methods ͓MP4, MCSCF, CIPSI-3, CCDϩST4͑CCD͒, etc.͔ and standard basis sets have been used to characterize the stationary points on both surfaces. 17 From these calculations, it appears that: ͑a͒ the 3 AЈ PES always lies energetically below the 3 AЉ PES although not far away from it; ͑b͒ nonlinear OCS transition states have been found on both surfaces; ͑c͒ a very shallow nonlinear OCS minimum ͑not bound if the zero point vibrational energy is included͒ has been found on the 3 AЈ PES; ͑d͒ the wave function is mostly single determinantal in nature at all the stationary points of both PES, and the unrestricted fourth order Mo "ller-Plesset ͑UMP4͒ method is accurate enough to be used for the calculation of the large number of ab initio points required in the fitting procedure to obtain an analytical expression for the PES.
Over 990 points have been calculated to obtain a dense grid at all relevant portions of the lowest 3 AЈ PES at the projected unrestricted fourth order Mo "ller-Plesset ͑PUMP4͒ level using the standard 6-311G(2d) basis set. The PUMP4 method has been chosen to eliminate the small spin contaminations which have appeared in some regions of this surface. The GAUSSIAN 92 suite of programs 18 has been used in all the computations. The calculated points span the range 0.9рR CO р2.5 Å, 1.25рR CS р2.65 Å, ЄOCSϭ165°, 145°, 125°, 105°, 85°. Moreover, 31 points for the CO energy curve and 38 points for the CS energy curve have also been computed to obtain a fully ab initio surface for this system. The calculation of each triatomic point takes around 42 min on a CRAY Y-MP 232 computer and 2 h on an IBM RS/ 6000 3BT workstation. Table I gives the properties of the stationary points at the PUMP4 level. The PUMP4 geometry of the transition state has been obtained by interpolation within a grid of ab initio points calculated around the region where the TS is located. Its geometry differs only slightly from the UMP4 one, 17 although the PUMP4 energy barrier ͑2.26 kcal/mol͒ is much lower than the UMP4 value ͑5.57 kcal/mol͒.
As shown in our previous paper, 17 the barrier height for this reaction is quite dependent on the method and basis set used. In addition, the barrier shows a clear tendency to diminish as the quality of the ab initio method and the basis set is increased ͑Table II͒. A much better search of the TS at the largest basis set, by optimizing the geometry, would probably lead to a lower energy barrier. The introduction of an estimate of the basis set superposition error ͑BSSE͒ by means of the counterpoise method, 20 originates an energy barrier increase of about 2 kcal/mol for each basis set ͑Table II͒. Most of these values can be consistent with the experimental activation energy value ͑1.5Ϯ0.5 kcal/mol͒. 11 Nevertheless, a small down scaling will be necessary in order to get an analytical 3 AЈ PES suitable to be used in kinetic and dynamical studies. On the other hand, the use of a larger basis set, as the 6-311G(3d,2f ) one, implies a huge cost in disk space and computer time ͑each triatomic point takes about 23 h on the abovementioned workstation͒ to calculate dense grids of points on this surface. However, the main features of the 3 AЈ PES would probably not change significantly, with the single exception of the barrier height. The exoergicity only changes from Ϫ90.6 to Ϫ90.3 kcal/mol when going from the 6-311G(2d) to the 6-311G(3d,2f ) basis set, while the experimental value is Ϫ87.4 kcal/mol. 19

B. Analytical fit
A many-body expansion 21 has been used to obtain an analytical representation of the 3 AЈ potential energy surface, which can be written as: in which V ͑2͒ and V ͑3͒ are the two-body and three-body terms, respectively, and R 1 , R 2 , and R 3 are the CO, CS, and SO distances, respectively. In this equation, monoatomic terms V ͑1͒ have been omitted because all possible reaction channels for O͑ 3 P͒ϩCS͑X 1 ⌺ ϩ ͒ correlate with the atoms in their ground electronic states Extended Rydberg potentials up to the fifth order have been used to describe the two-body interactions ͑diatomic molecules͒ in Eq. ͑1͒, for the CO and CS diatomic curves: in which D e and R e are the dissociation energy and equilibrium bond length of the diatomic molecule, respectively, and is defined as being equal to RϪR e . The parameters a i have been determined by means of a nonlinear least squares procedure, 22 using the corresponding ab initio points ͑31 and 38 points for CO and CS, respectively͒. Figure 1 depicts both the ab initio data and fitted energy curves of the three diatomics. For the SO molecule the PUMP4 method is not suitable, as the wave function must be described with more than a single determinant, as occurs in the case of the isoelectronic O 2 molecule. Thus, in the SO case only a few points ͑8͒ around the equilibrium distance and the dissociation energy have been included in this fit, using an extended Rydberg potential only up to the third order. Nevertheless, as the aim of this work is the study of reaction ͑2a͒, the lower accuracy in the description of the SO diatomic does not represent a serious difficulty. Moreover, as the experimental endoergicity of reaction ͑2b͒ is 46.2 kcal/mol, 19 this second reaction channel will be closed under the usual experimental conditions. The root-mean-square ͑RMS͒ deviations for the diatomic energy curves of CO, CS, and SO are 0.081, 0.087, and 6.0ϫ10 Ϫ6 eV, respectively. The optimum extended Rydberg parameters of each molecule are given in Table III. The three-body term consists of a fourth order polynomial expressed in terms of three variables i ͓ i ϭR i ϪR i 0 , where the reference structure (R 1 0 ,R 2 0 ,R 3 0 ) has been taken as equal to the TS geometry͔, and a range function T(R 1 ,R 2 ,R 3 ) which cancels the three-body term as one of the three atoms is separated from the other ones where 1рiϩ jϩkр4 being i, j, and k positive integer numbers, and

͑6͒
The 38 three-body parameters (V 0 ,͕c i jk ͖,͕␥ i ͖) have been determined by a weighted nonlinear least squares procedure 23 using the following ab initio data: ͑a͒ energy, geometry, and harmonic force constants of the ab initio TS, and ͑b͒ a total of 739 ab initio points within the abovementioned triatomic geometrical intervals and with energies up to 2.0 eV over reactants. The inclusion of points with higher energies worsened the adjustment and, in fact, is not necessary for the study of reaction ͑2a͒ at the usual experimental conditions, thermal ͑300 K͒ or hot energies ͑relative translational energy of 0.16 eV͒. [13][14][15] Moreover, the energy barrier has been scaled to reproduce as closely as possible the experimental rate constant at 300 K and the activation energy, when these properties have been calculated using the variational transition state theory ͑see next section͒. About 150 analytical representations of the 3 AЈ PES have been carefully analyzed in order to obtain the best fitted surfaces. In a first stage of the fit, an initial guess for the ␥ i parameters has been used and the corresponding optimum polynomial parameters have been obtained by a linear least squares procedure. These parameters have later been used as input parameters in a nonlinear least squares procedure that leads to several analytical PES with smaller RMS deviations. In the last stage of the fitting procedure, a weighted nonlinear least squares procedure has been considered, taking as start-  Second partial derivatives at the TS geometry have also been introduced in the fit by using very small weights ͑0.01͒, because higher weight values greatly increased the RMS deviations of the ab initio points. During the fit, additional ab initio points have also been calculated to better describe the regions of the surface where spurious minima appeared ͑e.g., in the OSC and SOC collinear regions͒. For all analytical surfaces obtained, the spurious minima have been located and plots at different OCS angles have been carefully analyzed. Spurious minima with energies lower than 1 kcal/mol have been considered negligible. In fact, these shallow minima have also been found in the ab initio calculations, at the entrance and exit valleys, with energies lower than 0.5 kcal/mol. Three analytical 3 AЈ potential energy surfaces ͑PES 1, PES 2, and PES 3͒ have been selected from the fitting procedure. Their main difference is the barrier height for reaction ͑2a͒. Even though very different initial parameters have been considered in the fitting procedure, there has been a clear convergence to the three sets of optimum parameters here reported. The parameters for the three surfaces ͑Table IV͒ and their equipotential contour diagrams ͑Fig. 2͒ are very close to each other. The surface plots show a smooth fit of the ab initio data for each OCS angle without spurious artifacts. In Figs. 2 and 3 the energetically favored nonlinear approach of the attacking oxygen atom to the C end of the CS molecule can also be seen. Reaction channel ͑2b͒ is energetically forbidden at the usual experimental collision energies. Table V shows the differences between these 3 AЈ analytical PES and the quality of their fit. In general, the RMS deviations for each OCS angle in these surfaces are very close to each other. Moreover, when the RMS deviations for different intervals of energy between reactants and up to 2 eV over them are calculated, the values are almost constant for each energy range and surface. Thus, these surfaces show a homogeneous goodness in all regions which are fitted. The total RMS deviation is around 3 kcal/mol, a value that is quite good in comparison with the estimated accuracy of the ab initio data, evidenced by comparing the ab initio energy barrier and the exoergicity of reaction ͑2a͒ with the corresponding experimental data.
A comparison of the best analytical representations for the 3 AЈ PES obtained in this work with the earlier one 16 shows two important differences, which will probably modify to a large extent the calculated kinetics and dynamics of this reaction: ͑a͒ the new 3 AЈ PES initially favors a nonlinear O-C-S approach instead of the collinear one found for the earlier surface; ͑b͒ the new 3 AЈ surface does not present a nonlinear stable OCS minimum, which had a strong influence in the kinetic and dynamical properties of reaction ͑2a͒ on the earlier surface. Only a very shallow minimum with no TS in the exit channel has been found in the present ab initio calculations, and it has not been reproduced in the new analytical 3 AЈ surfaces.

III. CALCULATION OF THERMAL RATE CONSTANTS
The rate constants for reaction ͑2a͒ have been calculated, within the 150-1000 K temperature range, using the conventional transition state theory ͑TST͒ and variational transition state theory ͑VTST͒ with semiclassical tunneling, as implemented in the POLYRATE program. 24 VTST rate constants have been computed using the improved canonical variational theory ͑ICVT͒. Tunneling corrections have been taken into account by means of the microcanonical optimized multidimensional tunneling method ͑OMT͒, although differences with the small curvature tunneling method ͑SCT͒ are very small. Table VI shows the rate constants at 200 and 300 K computed using the three ab initio analytical 3 AЈ PES, and the previous MNDO/CI based analytical surface. 16 For each analytical PES, the rate constants calculated with the different levels of the transition state theory are quite close to each other, due to the shape of the PES and to the heavy atoms involved. The ICVT/OMT values are, however, in general slightly better when compared with the experimental values. The PES 3 surface produces the best rate constants in the 150-300 K interval. At 300 K only a factor of 2 below the experimental data occurs. This result is much better than the earlier one obtained from the MNDO/CI surface ͑where a factor of 50 was found 16 ͒. However, if it is assumed that only the lowest surface among the three which connect reactants with products ͑one 3 AЈ and two 3 AЉ͒ contributes to the reac- tivity, a statistical factor equal to 1/3 has to be included into the calculated rate constants and preexponential factors, such as the values are given in Table VI ͑see also Fig. 4͒.
For all PES and methods the Arrhenius plot of the rate constants shows a curvature between 150 and 1000 K ͑Fig. 4͒. Table VI also gives the Arrhenius parameters in two temperature ranges: 150-300 K and 300-1000 K. PES 1 is the surface that best reproduces the experimental activation energy in the interval 150-300 K. Nevertheless, it is necessary to take into account that only rate constants around 300 K have been measured. Only in one experimental work 12 were the rate constants at lower temperatures determined ͑Fig. 4͒ by means of a discharge flow system with mass spectrometer detection of CO product, and rather poor Arrhenius behavior was found. In fact, the Arrhenius parameters obtained from the experimental data shown in Fig. 4 ͑Aϭ9.82ϫ10 13 cm 3 mol Ϫ1 s Ϫ1 and E a ϭ1.35 kcal/mol͒ correspond to an r correlation coefficient of 0.944, indicating as well a small curvature which introduces an important error bar in the activation energy. In the same experimental work, 12 the 294 K rate constant was introduced twice. One value was measured monitoring the CS concentration and the other one by following the CO concentration. This fact modified the Arrhenius parameters appreciably ͓Aϭ͑15.7Ϯ2.4͒ϫ10 13 cm 3 mol Ϫ1 s Ϫ1 and E a ϭ1.5Ϯ0.3 kcal/mol͔, improving slightly the quality of the fit ͑rϭ0.948͒. Comparison of these parameters with the experimental recommended ones, 11 which included more rate constants at room temperature ͑Table VI͒, and with the calculated values for the three analytical surfaces derived from ab initio calculations ͑PES 1, PES 2, and PES 3͒ seems to indicate that even though PES 1 gives the activation energy closest to the recommended value, PES 3 could also furnish an activation energy consistent with the experimental data, once both the experimental and theoretical error margins are included.
A more accurate calculation of the rate constants would have to take into account the contribution of the two 3 AЉ surfaces which connect reactants with products. In a previous study, 17 the energy difference between the 3 AЈ and the lowest 3 AЉ surface was shown to be very small. As the lowest 3 AЈ and 3 AЉ surfaces present on the whole a similar shape ͑nonlinear OCS transition states with similar geometries͒, the differences in the energy barrier will be the main factor affecting the calculated rate constants. Thus, to evaluate the relative importance of the two lowest PES contributions to the thermal rate constant, the same ab initio and VTST method have been used for these surfaces.
First and second partial derivatives with respect to the interatomic distances are necessary to calculate ICVT rate constants. Even though the UMP4 or PUMP4 methods are more appropriate for calculating absolute rate constants, as was shown in Sec. II, for relative rate constants UMP2 calculations can furnish a good estimate. Thus, minimum energy paths for the two lowest surfaces have been obtained by means of the intrinsic reaction coordinate ͑IRC͒ method. 25 FIG. 3. Equipotential contour maps of the PES 3 for the O to CS approach with R CS ϭ1.5645 Å. The CS molecule lies along the X axis at Y ϭ0 with its center of mass at the origin. Each path begins at each transition state with a step of 0.015 bohr. A total amount of 83 UMP2/6-311G(2d) ab initio points have been calculated, 41 for the 3 AЈ PES and 42 for the lowest 3 AЉ PES, including first and second analytical derivatives ͑this represents around 500 h on the abovementioned workstation͒. Figure 5 represents the 3 AЈ/ 3 AЉ ratio of the ICVT/SCT rate constants at different temperatures. This ratio decreases from 46.3 at 150 K to 2.4 at 1000 K. The effect of this second surface ͑the lowest 3 AЉ PES͒ on the thermal rate constant calculation is shown in Fig. 4, where the third surface ͑the first excited 3 AЉ PES͒ is supposed not to be reactive. It was also assumed that the UMP2/ ICVT/SCT calculated rate constant ratios ͑Fig. 5͒ are valid for the calculation of PUMP4/ICVT/OMT rate constants. Thus, the effect of the lowest 3 AЉ PES is more important at higher temperatures. Adding this contribution to that of the 3 AЈ PES, both multiplied by a factor equal to 1/3, it is not possible to adequately reproduce the rate constant at room temperature ͑Fig. 4͒. Calculated rate constants for PES 3 at 300, 200, and 150 K, estimating the contribution of the corresponding lowest 3 AЉ PES as indicated above, are equal to 2.45ϫ10 12 , 1.48ϫ10 12 , and 1.03ϫ10 12 cm 3 mol Ϫ1 s Ϫ1 , respectively. At the lowest temperature considered, the agreement with the corresponding experimental value ͑1.02ϫ10 12 cm 3 mol Ϫ1 s Ϫ1 ͒ is excellent, but as temperature increases the FIG. 5. Relative thermal rate constants for O͑ 3 P͒ϩCS→COϩS͑ 3 P͒ on surfaces 3 AЈ and 3 AЉ calculated using UMP2 ab initio points and the ICVT/ SCT method. Solid line shows an exponential fit whose equation is indicated in the plot. accord progressively diminishes. From the UMP2 calculations it appears that the 3 AЈ/ 3 AЉ rate constant ratio decreases from 46.3 to 8.4 between 150 and 300 K. Hence, the 3 AЈ surface contribution clearly dominates over the 3 AЉ one at this temperature range.
VTST calculations indicate that the three optimal analytical 3 AЈ PES that have been determined in this work cannot simultaneously fit the experimental rate constant at room temperature ͑300 K͒ and the activation energy ͑150-300 K temperature range͒. Additional and more accurate ab initio calculations on lowest 3 AЉ surface would be necessary to shed more light on this problem. On the other hand, the role of the second 3 AЉ surface will probably be relevant at much higher temperatures than the ones here studied.

IV. CONCLUDING REMARKS
In this work we have presented three analytical representations of the 3 AЈ adiabatic potential energy surface ͑PES 1, PES 2, and PES 3͒, which fit 816 PUMP4/6-311G(2d) ab initio points and three different activation energies for the reaction O͑ 3 P͒ϩCS͑X 1 ⌺ ϩ ͒→CO͑X 1 ⌺ ϩ ͒ϩS͑ 3 P͒. These optimal analytical surfaces show small RMS deviations in all adjusted regions and a smooth shape for all OCS angles. Variational transition state rate constant values of the PES 3 surface with semiclassical tunneling correction at room and lower temperatures ͑300-150 K interval͒ are reasonably close to the experimental values, with excellent agreement at 150 K. After a great number of attempts, it has not been possible to obtain an analytical surface capable of simultaneously reproducing both the rate constant at 300 K and the activation energy for the 150-300 K range of temperatures. The effect of the lowest 3 AЉ surface on the thermal rate constant has to be considered at higher temperatures to compare in a meaningful way with the experimental data. The present results indicate that the PES 3 surface is a good candidate to be used for kinetic or dynamical studies on a single PES. A quasiclassical trajectory calculation on this surface is in progress, to study the CO vibrational distribution, several rovibrational distributions, and some vector correlations studied experimentally at 300 K and for hot ͑hyperthermal͒ O͑ 3 P͒ atoms.