Atmospheric reaction of methyl mercaptan with hydroxyl radical as an acid rain primary agent

For the CH3SH + OH atmospheric reaction, we study the mechanism, potential energy surface, thermodynamic parameters of all stationary points, and rate of generation of the main product channels at high, low, and intermediate pressures. In this study, the UMP2, UM062X, UB3LYP, and CCSD(T) methods by Dunning and Pople basis sets are used and the results are compared with the experimental data. It is theoretically predicted that the reaction has fourteen possible pathways with eight different products in the gas phase. The thermodynamic results show that OH radical extracts predominantly the hydrogen of the SH functional group compared to the hydrogen of the CH3 group of CH3SH. Also, the rate constant calculations indicate that the extraction of the hydrogen atom of the SH group has a major role in 150–3000 K, while a good contribution is observed for the hydrogen of methyl group above 1200 K. Our results show that the used methods lead to good agreement with experiment. Finally, we demonstrated that why the main path is the main path.

Reactions of methanethiol with active atmospheric species such as CN, NO X (X = 1-3), and OH radicals have been noticed by experimental and theoretical researchers. Bao-En et al. have investigated the reaction mechanism of CH 3 SH + CN in the gas phase theoretically. This work showed that CN tends to react with the hydrogen of the SH group. Also, this channel was a barrier-less reaction 14 . Reactions of CH 3 SH plus NO and NO 2 could have a significant role in the atmosphere. The photolysis of CH 3 SH in the presence of NO radicals led to producing CH 3 SNO as an intermediate and (CH3S) 2 as a final adduc 15 . Also, The reaction of CH 3 SH + NO 2 led to generate CH 3 16 . A theoretical analysis of the kinetics and thermochemistry of the CH 3 SH + NO 2 reaction was presented by Brudnik et al. 17 . To generate CH 3 S + HNO 2 and CH 2 SH + HNO 2 products, they computed rate constants were k(T) = 7.9 × 10 -15 (T/300) 1.9 exp(− 8190/T) and k(T) = 7.9 × 10 -13 (T/300) 1.94 exp(− 16,290/T) cm 3 molecule −1 s −1 , respectively. The reaction of methyl mercaptan with NO 3 18 had less important in comparison with OH, due to the low rate with an order of 10. Also, NO 3 reacts predominantly with the hydrogen of the SH group.

SNO intermediate and OH or CH 3 S + HONO and/or CH 3 S(O)H + NO
The kinetics and mechanisms of the CH 3 SH + OH have been investigated extensively by experimental researchers [19][20][21][22][23][24][25][26][27] . The measured values for the rate constant at room temperature were in the range 2.1-9.04 in unit 10 -11 cm 3 molecule −1 s −1 . The reported rate constants will be discussed in more detail later. Masgrau et al. 28 have investigated theoretically the title reaction by MCCM-CCSD(T)-CO-2m/cc-pVDZ//MP2(full)/cc-pVDZ level. They have used VTST theory for the calculation of the rate constant and the obtained result at this level showed that the rate constants of H-abstraction from the SH and CH 3 groups at 298 K were 8.85 × 10 -12 and 2.95 × 10 -14 cm 3 molecule −1 s −1 , respectively. In addition to the main pathway, some of the addition-elimination reactions of the CH 3 SH + OH were proposed computationally using the G1 method by Muino 29 .
In this article, we argue the connection between excitation energies and the stability of prereactive complexes as a starter of the reaction in the presence of sunlight. Using these results along with kinetic calculations, we answer this question "which reaction has the most probability to occur in the atmospheric condition?'' . It is shown that the methodology used here has precise results in the kinetic calculations of the title reaction in comparison with the previous study. Using the validated meta hybrid M06-2X method, the effect of pressure is cleared that has been less attended in the mentioned theoretical and the most experimental studies until now. We demonstrate that the pressure and altitude have significant effects on the rate of the CH 3 SH + OH reaction. Moreover, similar to our previous work 30 on the reactions of methyl mercaptan, we pay attention (a) to complete investigation of the mechanism of CH 3 SH + OH reaction (b) to obtain more accurate theoretical intuition of the potential energy surface (c) to calculate pressure/temperature dependent rate constants with high accuracy, and (d) to get precise thermodynamic viewpoints on the title reaction. Therefore, the CH 3 SH + OH reaction pathways are studied using both wave function and density functional theories to obtain reliable potential energy surface, precise rate constants of main paths, and thermodynamic parameters of all components in the proposed PES. Therefore, the multi-well potential energy surface is calculated by higher-level calculations. The coupledcluster method including triple excitations such as the UCCSD(T)/aug-cc-PV(T + d)Z//UMP2/aug-cc-pVTZ and the UCCSD(T)/6-311 + + g (3df,3pd)//UMP26-311 + + g (3df,3pd) levels are used to get an accurate and complete treatment of the mentioned reaction pathways. The stability of all adducts is examined energetically using thermodynamic parameters at the UCCSD(T)/aug-cc-PV(T + d)Z(energies) + UMP2/aug-cc-pVTZ (thermodynamic corrections) level at room temperature. In addition, to get the nature of interactions of pre-reactive collision complexes, topological analysis of the electronic charge density is implemented. High-pressure limit rate constants are carried out by using transition state (TST) and VTST theories for the reaction paths with one transition state (second-order elementary reactions), and Rice-Ramsperger-Kassel-Marcus (RRKM) theory is used for the low-pressure limit rate constant and its behavior in the fall of regime region.

Methodology
In the reaction of CH 3 SH + OH, all stationary points on the ground-state of the multi-well PES were fully optimized using the UMP2 (wave function) 31 , UM06-2X 32,33 , and UB3LYP 34,35 (DFT) along with the Dunning's augmented correlation-consistent polarized valence triple zeta, aug-cc-pVTZ 36 , and split-valence augmented triple zeta, 6-311 + + G(3df,3pd), basis sets. The harmonic vibrational frequencies were also computed for all stationary points at the mentioned level of theories. The calculated frequencies were used to provide the zeropoint vibrational energy (ZPE), enthalpy, Gibbs free energy, and absolute entropy for each stationary point. Also, the vibrational frequencies were applied to confirm the nature of the corresponding stationary points as minimum structure or transition state structure in the PES. All thermodynamic parameters were computed by assuming each stationary point has an ideal gas behavior of the standard pressure of 1 atm temperature of 298.15 K. Moreover, the intrinsic reaction coordinate (IRC) 37,38 was implemented to ensure that the transition states have connected to the desired reactants, complex reactants, intermediate and complex products. IRC calculations were done at the UMP2/6-311 + + G(d,p) level of theory for all transition states. The energies, geometries, and the harmonic vibrational frequencies in different levels were obtained to examine the sensitivity of the stationary points to different methods and basis sets and to assess the reliability of the mentioned levels in the prediction of minimum structures and the transition state structures. Since, it is well-known that the single reference methods such as UMP2 and UM06-2X, UB3LYP without including some important excitations have some large error on the estimated energy (for example for H atom abstraction reactions the UB3LYP method underestimates the computed energy barriers). Therefore, we carried out (frozen core) single-point calculation using coupled-cluster (CC) theory. In CC calculation entire single and the double excitations accompanied by the perturbative treatment of the entire connected triple excitations were included in a single reference determinant of the UHF formalism,. The UCCSD(T) 39,40 calculations were carried out on the geometries optimized at the UMP2/aug-cc-pVTZ level. In the UCCSD(T) calculations the high exponent d functions for the sulfur atom [41][42][43][44] were included in combination with the aug-cc-pVTZ basis set to generate the aug-cc-pV(T + d)Z basis www.nature.com/scientificreports/ set. The topological analysis based on the theory of atoms in molecules (AIM) 45 was employed to determine the bonding features of in complexes using the analysis of the electronic charge density and its Laplacian at critical points. For this purpose, the Hassian Matrix as a first-order density matrix was computed at the UMP2/aug-cc-pVTZ level of theory. T1 diagnostic values play an important role in open-shell systems because of the need for multi-reference calculations if the threshold value is larger than 0.45. T1 diagnostic values of all stationary points were calculated. All T1 diagnostic values are below 0.037 [46][47][48][49] .
All of the high-pressure limit rate constants were calculated using Gpop program 50 based on the transition state (TST) and variational transition state (VTST) theories and the reaction rate in the low-pressure limit and the fall of the regime were computed using strong collision master equation/ RRKM theory by Ssumes program 51 . The Eckart tunneling correction 52 was used for the inclusion of quantum effects in rate constant calculation. The Gaussian suite program 53 was executed for the electronic structure calculation of the title reaction components.
When an unrestricted Hartree-Fock based method like the UMP2 method was applied to calculate the electronic structure of an open-shell system, the eigenvalues of the S 2 operator are key parameters to assess spin contamination. The quantity of the spin contamination is gained by the expectation value of the S 2 operator, <S 2 >. The expectation value of <S 2 > of the UHF base methods obtains from the spatial overlaps among all pairs of α and β spin-orbitals as follows: If all α and β orbitals are the same, there is no spin contamination in the considered system. If the amount of the expectation value <S 2 > is greater than 5% of the pure spin state, the used wave function is unacceptable for computing the energy 54 . Our results show that the spin is contaminated at the MP2/aug-cc-pVTZ level. Spin contamination makes the computed energies to be unacceptable at the mentioned level. Thus, we extract the energies at the PMP2/aug-cc-pVTZ level ( P refers to the projected wave function). All of the energies reported in the text are based on the UCCSD(T)/aug-cc-pV(T + d)Z level, unless we have mentioned the level of computation.

Results and discussion
All of the computed relative energies (in kcal/mol) at the different methods are improved by ZPE correction. The corrections on the absolute energies calculated at the MP2, M06-2X, and B3LYP methods are done by the zero-point energy that is computed at the same level. Also, the ZPEs at the MP2 level are implemented for the correction of the UCCSD(T) energies for kinetic calculations. The relative energies of different species at the mentioned levels are tabulated in Table 2. The potential energy surface for all reactions is plotted by relative energies at the UCCSD(T)/aug-cc-pV(T + d)Z level and depicted in Fig. 4. The geometrical parameters of the optimized structures such as the reactants, the primary complexes, the transition states, the post reactive complexes, and the products are represented in Fig. 3 at the UMP2/aug-cc-pVTZ level of theory. For all stationary points, the extracted T1 diagnostic values and the expectation values of the S 2 operator are listed in Table S4 at the UCCSD(T)/aug-cc-pV(T + d)Z and the UMP2/aug-cc-pVTZ levels, respectively. All thermodynamic parameters are calculated at the UMP2/aug-cc-pVTZ level and collected in Table 3. Calculated equilibrium constants and rate constant values using the mentioned theories are collected in Table 1. The Arrhenius plot of the P1 and P2 products are depicted in Fig. 1. Throughout the paper, all thermodynamic values are reported at temperature 298.15 K and pressure 1 atm.
Enthalpy and free energy formation, relative stability, and nature of interactions in the prereactive complexes. The enthalpy formation of each pre-reactive complexes are calculated by the following equation where K i is the equilibrium constant of i-th complex (i = 1, 2 and 3 in respect of three pre-reactive complexes), X • 298 (S) is the absolute enthalpy or free energy of S species. All equilibrium constants computed at the UCCSD(T)/ aug-cc-pV(T + d)Z + UMP2/aug-cc-pVTZ level are sketched in Fig. 2 and tabulated in Table 1 and all of the reported X 0 298.15 (Y ) values (Y is reaction components) in this study are based on the results of the UCCSD(T)/ aug-cc-pV(T + d)Z + UMP2/aug-cc-pVTZ level (the total energies from UCCSD(T) method and thermodynamic corrections from the UMP2 method) and are collected in Table 3.
In Fig. 3, It is shown that three collisions pre-reactive complexes CR1, CR2, and CR3 are formed resulting from the collision of two reactants. The CR1 complex has relative stability of 4.47 kcal/mol in comparison with the original reactants. This value is 1.16 kcal/mol lower than the energy of the corresponding structure in the Muino study. In this complex, the atom S interacts with the atom O, and the bond distance S…O is 2.012 Å. The formation of CR1 complexes is exothermic and nonspontaneous with standard enthalpy and Gibbs free energy changes of − 2.27 and 6.97 kcal/mol, respectively. The CR2 complex is more stable than the CR3 complex. In the structure of the CR2 complex, the atom 7H of the radical hydroxyl is close to the sulfur atom of the methyl sulfide. The distance 7H-5S is 2.421 Å. The CR2 complex is formed without any barrier and has − 4.48 kcal/mol relative energy. The computed thermodynamic parameters, ∆H 0 and ∆G 0 , for the CR2 complex are − 3.36 and 3.26 kcal/ mol, respectively, that show the formation of CR2 is also an exothermic and nonspontaneous process. In the third initial complex, CR3, the atom O of hydroxyl radical interacts with the atom C of methyl sulfide, in which the distance O…C is 3.22 Å and has the relative energy of − 1.05 kcal/mol (see Table 2). Also, the formation of www.nature.com/scientificreports/ the CR3 is an exothermic and nonspontaneous process with ∆H 0 = − 0.17 and ∆G 0 = 3.52 kcal/mol. On the bases of the sketched PES, the reaction of methyl sulfide with hydroxyl radical is started with the three pre-reactive complexes on the doublet potential energy surface, and finally, 13 different products are predicted. The PES profile (relative energy versus reaction coordinate) is plotted in Fig. 4 at the ground state of all species based on the calculated energies at the UCCSD(T) level. The details of the reaction pathways are described in the next section. The equilibrium constants of the CR1, CR2, and CR3 complexes are calculated in the range of 230-3000 K. The equilibrium constant of the CR1 is decreased until 900 K, and it remains roughly constant after this temperature. www.nature.com/scientificreports/ Also, the same behavior was observed for the CR2. At the mentioned temperature range, the equilibrium constant of the CR3 is increased. The computed equilibrium constants for the CR1, CR2, and CR3 complexes at room temperature are 7.86 × 10 -23 , 2.06 × 10 -21 , and 2.35 × 10 -22 cm −3 molecule −1 , respectively. On the basis of the calculated equilibrium constants at the UCSD(T)/aug-cc-pV(T + d)Z level, we conclude that the CR2 and CR1 are more stable than the CR3 in low temperatures, but at the moderate and high temperatures, the CR3 equilibrium constant is more than the others. The topological analysis (AIM calculation) of the electron charge density in the mentioned complexes reveals the existence of bond critical points among the abovementioned interactions. In the CR1complex, the bond critical point between the S…O bond shows a noncovalent interaction with the charge density of ρ = − 0.1080 e bohr −3 and the charge density Laplacian of ∇ 2 ρ = 0.1259 e bohr −5 . The bond critical point in CR2 located between the atoms 7H and S (ρ = 0.01847 e bohr −3 and ∇ 2 ρ = 0.0437 e bohr −5 ) indicates a van der Waals interaction. www.nature.com/scientificreports/ www.nature.com/scientificreports/ The interaction between atoms O and C of the third complex is also a van der Waals (ρ = 0.0040 e bohr −3 and ∇ 2 ρ = 0.0204 e bohr −5 ), but it is weaker than the interaction of the CR2 complex.

Pathways description. P1 (H 2 O + CH 3 S) and P2 (H 2 O + CH 2 SH) production pathways. The R1 and R2
pathways are H abstraction paths from SH and CH 3 groups of CH 3 SH to generate the P1 and P2 products, respectively. The possible paths are as follows: The CR1 complex is transformed via TS1 to the CP1 product complex by passing through the barrier-less process with the relative energy − 1.74 kcal/mol. Also, this path was considered by Muino. The computed relative energy for a similar structure was − 2.029 kcal/mol at the G1 level. IRC calculation shows that the atom 8H is shifting from the atom S of CH 3 SH to the atom O of OH fragment. In the structure of the TS1, the H-S bond is breaking up between the atoms 8H and 5S with a length of 1.398 Å and, the O-H bond is forming between the atoms 6O and 8H with a length of 1.454 Å. The imaginary frequency of this saddle point in the mentioned conversion is 1051 cm −1 at the MP2/aug-cc-pVTZ level of theory. The CP1 complex without passing any transition state converts to the P1 as a final product. The reverse reaction has less important due to the need for a high energy barrier (35.48 kcal/mol). It may relate to the low electronegativity of the sulfur atom in comparison with the atom O. So, the bond between the atoms O and H (123.58 kcal/mol) in isolated water is stronger than the bond between the atoms S (91.00 kcal/mol) and H in methanethiol.
Regarding the above-mentioned path, the CR2 complex via TS2 by surmounting on the barrier of 6.63 kcal/ mol converts to the CP2 product complex. The transition state TS2 is 2.15 kcal/mol higher than CH 3 SH + OH on the proposed PES. Then, CP2 yields the P2 product by breaking up the bonds between CH 2 SH + H 2 O moieties. The transition state 2 involves cleavages the bond between the atoms 3C and 5H with a length of 1.725 Å and also at the same time, it includes the formation bond between the atoms 1O and 5H with a length of 1.364 Å. The reverse reaction of this path same as the previous path requires high energy to occur (27.27 kcal/mol). Here the reason is different than the first path. As we know, the electronegativity of the carbon atom is lower than the oxygen atom, and also the stability of methyl radical is lower than OH radical. Because the high electronegativity of OH radical cause the radical electron roughly stable in comparison with that electron when it locates on the carbon atom.
On the basis of the collected data in Table 3, we can conclude that the P1 and P2 formation processes are thermodynamically favorable with the standard enthalpy of − 31.12 and − 22.09 kcal/mol and Gibbs free energy of − 32.03 and − 23.90 kcal/mol. This is an exothermic and spontaneous process. Our computed enthalpy is www.nature.com/scientificreports/ comparable with that values found by Muino ( H 0 = − 33.90 kcal/mol), but the Gibbs free energies have a large difference ( G 0 = − 34.01 kcal/mol).

P3 (CH 3 + HSOH), P4 (H + CH 3 SOH) and P5 (SH + CH 3 OH) production pathways.
The following schemes show the production pathways of P3, P4, and P5 adducts: P3 adduct is produced by two paths. The R3a and R3b paths occur through only one transition state. The barrier energy of path R3a and path R3b are 11.84 and 20.46 kcal/mol, respectively. The CR1complex after passing the transition states TS3 and TS5 transforms into the CP5 complex. The mentioned transition states include methyl elimination from CH 3 SH and OH addition to the SH group. Firstly, the OH group contacts to the S atom. Secondly, the methyl group separates from CH 3 SH. The differences between the mentioned saddle points are related to the orientation of the methyl group and the distance of the atoms C and S. In the saddle point TS3, the angle C-S-O, and the distance C-O are 163.25° and 2.92 Å, but the corresponding angle in TS5 is 83.25° and 2.178 Å.
The CR1 is transformed into the CP4 by passing through TS4 with an energy height of 21.18 kcal/mol at the UCCSD(T) method. Then, without entrancing any barrier energy, the CP4product complex converts to the P4 as a final product. This path is the same as the P3 production paths, but the hydrogen atom is separated in this reaction. In other words, firstly, OH radical is closed to the atom S and crates a covalent bond and after that, the hydrogen atom of the SH functional group is separated. The final result of the R4 path is the replacement of the OH fragment with the H atom. This path is also investigated by Muino. He calculated relative energy for the same saddle point is 15.42 kcal/mol, which our result is 5.76 kcal/mol higher than that value.
From a geometrical point of view, the structure of TS4 involves the bond formation among the O and S atoms with a length of 1.690 Å and breaking up the bond between the atoms 1H and S atoms with a length of 2.053 Å. The imaginary frequency of TS4 at the UMP2 method is 668 cm −1 in the reaction coordinate. According to Fig. 4, between the separation of a methyl group and a hydrogen atom through P3-P5 paths, our results show that methyl elimination is energetically favorable than hydrogen elimination. These observations are related to the stability of the separated part. It is well known that the hydrogen atom is very unstable than many active atmospheric species. So, the CP4 product complex and also the P4 product are energetically unstable than other minimum structures. Therefore, it may back to reactants. The product complex CP3 and P3 product are also unstable than the original reactants, the reaction continuation from CP3 is favorable than the reverse reaction.
The product complex CP5 obtains via TS6 with the energy barrier of 27.19 kcal/mol relative to the CR3 complex. In this pathway, S N 2 reaction has occurred, in which the OH functional group comes from one side and the SH functional group goes out from the other side. In the structure of TS6, the bond O-C with a length of 1.826 Å is weakening and the bond S-C with a length of 2.10 Å is forming. This reaction has a large energy barrier among the other discussed first step bimolecular reactions. So, in the kinetic point of view, it is slower than the other, but thermodynamically it is favorable than the products of the paths R3a-R4. Therefore the back reaction is unfavorable kinetically. The reason may relate to the stability of radical electron in the SH group, corresponding to the electron cloud distribution. Other paths for the generation of CP5 are the paths R5b and R5c. Until CP3 formation is discussed above. This complex by surmounting on the energy barrier of 21.42 kcal/mol through TS8 converts to CP5. Eventually, in the last step of the R3a-R5c paths, the P3, P4, and P5 products are generated without passing through an energy barrier from the corresponding product complexes.
The production processes of P3 and P4 in standard conditions are endothermic according to the computed enthalpies 5.22 and 3.84 kcal/mol and nonspontaneous concerning to the calculated free energies 15.74 and 18.30 kcal/mol, respectively. The large unstability of P4 may relate to the production of atomic hydrogen. But, the production process of P5 is exothermic with H 0 = − 6.98 kcal/mol and spontaneous with G 0 = − 16.93 kcal/mol.

P6 (CH 4 + OSH), P7 (SH 2 + CH 2 OH) and P8 (H 2 O + H + CH 2 S) production pathways.
The following paths are designated to create the products P6, P7, and P8 www.nature.com/scientificreports/  www.nature.com/scientificreports/ The R6a and R6b paths are two-step reactions and so are occurring among two transition states (TS3 and TS7 or TS5 and TS7). Until CP3 generation was discussed above. The low-level transition state, TS7, with 6.39 kcal/ mol height is a barrier between the CP5 and CP7. Also, the IRC calculation shows the CP6 complex product is created from the related transition state TS7. Abstracted hydrogen in the structure of the saddle point 7 has a weak bond with O atom with a length of 1.80 Å and a bond with C atom with a length of 1.305 Å.
There are two paths with three steps for the production of P7 (R7a and R7b) as shown above. The reaction mechanism until CP5 generation is investigated in the R5 pathway. The reaction is continued by surmounting on TS9 with an energy barrier of 12.23 kcal/mol. The IRC calculation shows that TS9 leads to the formation of CP7. From a geometrical point of view, to continue the reaction by TS9, the HS fragment in the CP5 is closed to the CH 3 OH fragment and gets a hydrogen atom from the methyl group of methanol. AIM calculations indicate the presence of a bond critical point between the atoms C and 4H ( ρ = 0.1196 e bohr −3 and ∇ 2 ρ = − 0.1558 e bohr −5 ) and the atoms S and 4H ( ρ = 0.1541 e bohr −3 and ∇ 2 ρ = − 0.3084 e bohr −5 ) in the TS9 structure, which shows the C-H covalent bond is weaker than that of in the CP5 ( ρ = 0.3028 e bohr −3 and ∇ 2 ρ = − 1.288 e bohr −5 ) and arising a weak covalent bond between the atoms S and 4H.
In the exit channels of the R6a-R7b paths, the CP6 and CP7 product complexes are turned to relevant the P6 and P7 products without passing any energy barrier. The relative energies of these complexes are − 25.51 and − 13.04 kcal/mol, respectively.
The final path (R8) leads to produce the H 2 O + H + CH 2 S species. This path is a bimolecular reaction with two steps. The first step is CP1 formation by path R1, and the second step is the generation of P8 by path R8. The CP1 complex may keep on the reaction by surmounting on TS10 and directly without entrancing complex product converts to P8. The barrier energy of TS10 is 58.96 kcal/mol.
According to the information of Table 3, the P6 and P7 formation processes with the standard enthalpies of − 23.09 and − 11.50 kcal/mol and Gibbs free energies of − 24.24 and − 12.91 kcal/mol are exothermic and spontaneous. These parameters for P8 are H 0 = 20.35 kcal/mol and G 0 = 12.21 kcal/mol which indicate an endothermic and nonspontaneous process. Table 3. Thermodynamic parameters of all species (in kcal/mol) in the CH 3 SH reaction computed at the UCCSD(T)/aug-cc-pV(T + d)Z (total energies) + UMP2/aug-cc-pVTZ (thermodynamic corrections) level.   13 at 299-426 K (k(T) = 8.89 × 10 -12 exp(790 ± 300)/RT), and Tyndall and Ravishankara 21 in the 244-430 K range with the value of (3.3 ± 0.40) × 10 -11 cm 3 molecule −1 s −1 . As mentioned above, the rate expression measured by the Atkinson group implies negative activation energy. This result is also observed in our computed rate constant for the R1 pathway (barrier less reaction). Also, our computed rate constant, 3.6 × 10 -  www.nature.com/scientificreports/ 11 cm 3 molecule −1 s −1 by Cox and Sheppard 15 . It is 6.71 times larger than the predicted value by the VTST rate (1.44 × 10 -11 cm 3 molecule −1 s −1 ). Also, good agreement is observed between the kinetic data of VTST theory, and  [17][18][19] , that have very excellent agreement with the computed value of rate in this temperature. Using electron pulsed laser photolysis-pulsed laser-induced fluorescence technique, Hynes and Wine have investigated the kinetics of OH reaction with CH 3 SH at the temperature of 300 K and a pressure of 700 Torr. They measured value is (3.27 ± 0.36) × 10 -11 cm 3 molecule −1 s −1 which has good agreement with the calculated rate constant at the mentioned temperature with the value of (1.44 ± 0.36) × 10 −11 cm 3 molecule −1 s −120 .
Masgrau et al. have calculated theoretically the rate constants of the R1 and R2 paths at the MCCM-CCSD(T) level 22 . The obtained results at this level show the rate constants of mentioned paths are 8.85 × 10 -12 and 2.95 × 10 -14 cm 3 molecule −1 s −1 at 298 K and are 6.965 × 10 -12 and 1.75 × 10 -13 m 3 molecule −1 s −1 at 500 K and are 1.03 × 10 -11 and 1.04 × 10 -12 cm 3 molecule −1 s −1 at 1000 K, respectively. These values indicate that our computed rate constants agree well with the experimental results.
In summary, a comparison between our calculated rates by TST and VTST theories with the measured rates using several experimental techniques shows that our applied computational level, the UM06-2X/Aug-cc-pVTZ, is adequately precise in describing the kinetic of the CH 3 SH plus OH reaction.
Low-pressure limit rate constant and its behavior in the falloff regime. To study the pressure-dependent behavior of the rate constant for the CH 3 SH + OH reaction in the falloff range and low-pressure limit, the Rice-Ramsper- Table 6. Excited state parameters of all pre-reactive complexes in the CH 3 SH + OH reaction computed at the TD-PBE0/aug-cc-pvtz level of theory. The units of vertical excitation energies (E v ) and wavelength (λ) are eV and nm, respectively. www.nature.com/scientificreports/ ger-Kassel-Marcus (RRKM) theory was applied with a weak collision approach. In the calculation of pressuredependent rate constant the chemical activation mechanism was defined as follows: here M is (the third body) N 2 molecule. Applying the steady-state approximation to the concentration of CR1*, the pressure-dependent rate of the CR1 → CH 3 S + H 2 O conversion is At high and low-pressure limits where [M] → ∞ and [M] → 0, k(T,p) is the first-order rate constant and the second-order rate constant, respectively. The expression k i (T, p) at the high and low-pressure limit is as follows: Finally, to compute the pressure-dependent rate constant, the following equation is used: where κ and K(T) are tunneling correction and temperature dependent equilibrium constant, respectively.  www.nature.com/scientificreports/ The Eckart tunneling correction was used to achieve a reliable pressure-dependent rate constant. Nitrogen molecule is used as a third body because, in the atmosphere, the nitrogen molecule is the most abundant species. The average amount of energy transferred per collision E , and Lenard Jones parameters including collision diameter, σ (Å), and energy parameter, ε k B (K) are necessary for the calculation of the effect of pressure on the rate constant of reactions by the proposed theories. The amount of E during both up and down energy transfer collision for N2 bath gas is taken as 74 cm −1 . Lenard-Jones parameters for CH 3 SH, OH, and N 2 are 3.900 Å and 350.00 K for CH 3 SH, 2.750 Å, and 80.00 K for OH 55 . In the low-pressure limit where P → 0, k(T,p)/[N 2 ] is called k 0 (T) that is a termolecular rate constant with the units of cm 6 molecule −2 s −1 . Also, k 0 (T) is called the pseudo-third-order rate constant. Our calculations show that k 0 (T) is 4.45 × 10 -27 , 4.44 × 10 -30 , 2.80 × 10 -31 , and 1.71 × 10 -31 cm −6 molecule −2 s −1 at 150, 298, 500, and 700 K, respectively. Literature reviews show there are no experimental results for k 0 (T) of the CH 3 S + H 2 O production pathway.
Our computed rate constants in the temperature range of 150-800 K at the pressure range of 10 -4 -10 +4 bar in the falloff regime are depicted in Fig. 5 and listed in Table S7. In the calculation of k(T,p), Eq. (8) shows that the ratio of k ∞ /k 0 is important. Therefore, in the investigated pathway, with increasing pressure, the rate constant increases at a definite temperature in the mentioned temperature range (see Fig. 5). The k ∞ /k 0 ratio is increasing with increasing temperature. For example, at the temperatures of 150, 298.15, 500, and 700 K, this ratio is 7.59 × 10 2 , 4.01 × 10 4 , 4.29 × 10 5 , and 1.30 × 10 6 , respectively. If the ratio of k(T,p)/k(T, 1 bar) is defined as the reduced rate constant, we can obtain better insight about the effect of pressure on this reaction at a constant temperature. The reduced values for P = 10 -3 , 0.1, 10 and 10 3 bar are 2.09 × 10 -2 , 3.10 × 10 -1 , 257, and 7.32 at 150 K, and 6.79 × 10 -3 , 1.152 × 10 -1 , 5.66, and 75.20 at 298.15 K, respectively. These results showed that the rate constant increases with increasing pressure on the investigated temperature range.
Fate of methanethiol in the atmosphere. As argued above, the pathways leading to form the CH 3 S + H 2 O and CH 2 SH + H 2 O products play key roles and act as the main sink for methanethiol. If the lifetime of methanethiol in the investigated elementary bimolecular title reaction is calculated, it can be predicted the fate of the CH 3 SH compound in good fashion. The computed lifetimes for methanethiol are listed in Table 5. The lifetime values are computed at 0-50 km height from the surface of the earth and in an environment of OH radicals. On one hand, as our results show, the rate constant of the reaction decreases with decreasing of pressure that is companied by altitude increasing. On the other hand, as Table 5 shows the concentration of hydroxyl radicals is independent of height. Therefore, the calculated lifetime in the atmosphere varies from 2.08 × 10 +6 s −1 at 0 km to 2.25 × 10 +7 s −1 at 50 km. The obtained lifetime is short. Thus, the hydrogen abstraction path from the SH group is a dominant mechanism in the atmosphere for the degradation of CH 3 SH species. The computed lifetime by the high-pressure limit rate constant is 2.26 × 10 +4 s −1 , which the lifetime obtained by the high-pressure limit rate is 100 times higher than the pressure-dependent rate constant at 0 km altitude.
Which path could be the main route? In this section, the photolysis of all complexes is argued, based on the TDDFT calculations in the gas phase. It is proved that the PBE0 57 method has the most accurate results for the TDDFT calculations 58 . The geometries obtained at the MP2/aug-cc-pVTZ level are used for the TD-PBE0/aug-cc-pVTZ calculations. The obtained results are listed in Table 6. The ground state and excited state orbitals for prereacctive complexes are depicted in Figs. 6, 7 and 8 and for post-reactive complexes are shown in Figs. S2-S8. Also, we have performed similar calculations based on the TD-M06-2X/aug-cc-pVTZ//M06-2X/ aug-cc-pVTZ level (see Table S8). The vertical excitation energy (E v ), oscillator strength (f), and wavelength (λ) related to the first six excited states of all complexes are listed in Table 6. On the basis of Table 6, among the prereactive complexes, CR2 and CR3 are photolyzed more than CR1 due to comparing the first three excited state energies. Therefore, CR1 is more stable than CR2 and CR3 as a beginner of reactions in the atmosphere. It may relate to unpaired electrons of the sulfur atom and the radical electron of the oxygen atom. In comparison with S atom, the unpaired electrons (but radical electron) of OH radical due to the large electronegativity of the oxygen atom need for higher energy to excite. So, we only focus on the unpaired electrons of the S atom and radical electron. In the CR1 complex, some of the unpaired electrons of S atom are shared with the radical electron of OH moiety and they need to get large energy for excitation, but in other complexes, they are free for excitation. This result agrees with the computed relative energies of the mentioned complexes at the PMP2, Ub3lyp, UM06-2X, and UCCSD(T) methods in conjunction with augmented triple zeta basis sets. Because there is a bond between the oxygen atom of OH moiety and the atom S of CH 3 SH moiety with good strength and interaction (− 0.1080 e bohr −3 and 0.1259 e bohr −5 ). As we know, when excitation occurs, the electrons that are shared in a bond or an interaction can move to the upper layer, so the bond or interaction is broken. There is also good agreement with the kinetic results obtained from the experiments. The experimental results indicate that under the atmospheric conditions, especially under sunlight, the R1 pathway is the most probable path to happen. Therefore, we can conclude that the CR1 complex has the most possibility for the continuation of the reaction. Therefore, in the experiment the only production of the CH 3 S + H 2 O adduct is feasible.
This does not mean that other reactions do not occur. On the other hand, hydrogen abstraction from the methyl group is somewhat kinetically feasible, under certain conditions, such as collisions with suitable energy and orientation, they may occur. For the second step reactions occurring from CP1 and CP3, CP1 is unstable than CP3, relating to the calculated vertical excitation energies. This unstability is related to both the radical electron and unpaired electrons of the atom S in CH 3 S. They require smaller energy for excitation (especially radical electron on S atom). So, the reaction between HSOH and CH 3 is feasible than the CH 3 S + H 2 O reaction. Another reason for the possibility of the HSOH + CH 3 reaction is the energy barrier of this reaction. For the reaction between CH 3 OH and HS (production of H 2 S + CH 2 OH, which starts from CP5), the Cp5 complex is also an unstable complex under the sunlight of the atmosphere. In the presence of sun rays, the electrons of the atom S Scientific Reports | (2020) 10:18081 | https://doi.org/10.1038/s41598-020-74767-6 www.nature.com/scientificreports/ in HS become excited, so the HS…CH 3 OH complex is splintered. Finally, for a bimolecular reaction with a high energy barrier and for multi-step reactions in which prereactive, post reactive, and intermediate complexes have low stability due to excitation, the sun's solar radiation not only has no important role in the supplying activation energy of that reaction but also prevent continuing the reaction. But, in the reactions with stable wells, sunlight acts as a significant source for activation energy. For first step reactions with unstable prereactive complex, a collision with suitable energy and orientation may influence the reaction occurring, but for multi-step reactions with the same feature, the strong collisions with suitable orientations may lead to the reaction occurring. From the above statements, we concluded that the photo-oxidation of sulfur-containing compounds is important to the sulfur chemistry of the atmosphere. And we also concluded that many of the above-discussed paths, but the P1 production pathway, may occur under a condition that they complex is stable or be under certain collisions.

Conclusion
The PES of the CH 3 SH + OH reaction was described on the doublet state in more detail using the single point calculations at the UCCSD(T)/aug-cc-pV(T + d)Z level, and the geometries, thermodynamic corrections, and zero points corrections at the MP2/aug-cc-pVTZ level. Therefore, the mechanisms and kinetics of the reaction were cleared using the mentioned levels. Our results showed that conventional hydrogen abstraction reaction from the SH functional group has more contribution than the others in the atmospheric degradation of methanethiol. Thus, on one hand, its pathway kinetically plays an important role. In another word, this path is a kinetic pathway. On the other hand, the obtained product from this path has more negative standard Gibbs free energy ( G 0 = − 32.03 kcal/mol) than the others.
The electronic structure calculations along with the kinetic calculation showed that another favorable path after R1 was R2 production, which was involved hydrogen abstraction from the methyl functional group of methanethiol compound by hydroxyl radical.
On the basis of the computed PES, it was concluded that the main path of the reaction was barrier-less, so, this system is ideal to investigate pressure effect. Up to now pressure effect on the title reaction had remained obscure, so our results clear it. The chemical activation mechanism for the R1 pathway was confined to the P1(CH 3 S + H 2 O) reaction product. It was confirmed that the overall rate of title reaction depends on not only the temperature but also the pressure. Therefore, the importance of variation of pressure on the gas-phase reaction of CH 3 SH + OH was proved in the atmospheric condition. The negative temperature dependent rate constant was observed for barrierless pathways and positive for another, while the pressure-dependent rate constant for the main pathway was positive. A comparison of the computed rate constants with experimental results indicated that the UM06-2X method with the aug-cc-pVTZ basis set was a good case to obtain reliable pressure and temperature dependent rate constants.
The fate methyl mercaptan in the atmosphere was determined by using the computed pressure-dependent rate constants in different altitudes and also using the computed high-pressure limit rate constant.
Finally, It was shown that the excitations of complexes in the beginning and in the course of reactions could play a very important role in the atmospheric reactions.

Author contributions
H.D. and S.M. performed calculations; H.D. wrote the paper, analysed data, and revised the manuscript; M.V. corrected the manuscript, reviewed the revised manuscript, and supervised the whole study.