Impact of a single water molecule on the atmospheric oxidation of thiophene by hydroperoxyl radical

Water as an important assistant can alter the reactivity of atmospheric species. This project is designed to investigate the impact of a single water molecule on the atmospheric reactions of aromatic compounds that have not been attended to comprehensively. In the first part, the atmospheric oxidation mechanisms of thiophene initiated by hydroperoxyl radical through a multiwell-multichannel potential energy surface were studied to have useful information about the chemistry of the considered reaction. It was verified that for the thiophene plus HO2 reaction, the addition mechanism is dominant the same as other aromatic compounds. Due to the importance of the subject and the presence of water molecules in the atmosphere with a high concentration that we know as relative humidity, and also the lack of insight into the influence of water on the reactions of aromatic compounds with active atmospheric species, herein, the effect of a single water molecule on the addition pathways of the title reaction is evaluated. In another word, this research explores how water can change the occurrence of reactions of aromatic compounds in the atmosphere. For this, the presence of one water molecule is simulated by higher-level calculations (BD(T) method) through the main interactions with the stationary points of the most probable pathways. The results show that the mechanism of the reaction with water is more complicated than the bare reaction due to the formation of the ring-like structures. Also, water molecule decreases the relative energies of all addition pathways. Moreover, atoms in molecule theory (AIM) along with the kinetic study by the transition state (TST) and the Rice–Ramsperger–Kassel–Marcus (RRKM) theories demonstrate that the overall interactions of a path determine how the rate of that path changes. In this regard, our results establish that the interactions of water with HO2 (thiophene) in the initial complex 1WHA (1WTA or 1WTB) are stronger (weaker) than the sum of its interactions in transition states. Also, for the water-assisted pathways, the ratio of the partition function of the transition state to the partition functions of the reactants is similar to the respective bare reaction. Therefore, the reaction rates of the bare pathways are more than the water-assisted paths that include the 1WHA complex and are less than the paths that involve the 1WTA and 1WTB complexes.

that ignoring this percent makes large errors in atmospheric models 9 . The HO 2 …H 2 O complex has also been involved in the heterogeneous uptake of hydroperoxyl radical by water 10 . In the last decade, the effects of water vapor on the main reaction pathways of methane, formic acid, acetaldehyde, acetone, HOCl, glyoxal, dimethyl sulfide (DMS), propionaldehyde, sulfuric acid nitric acid, and HCl with hydroxyl radicals have been investigated by several researchers [11][12][13][14][15][16][17][18][19][20] . In the dimethyl sulfide, dimethyl sulfoxide, and dimethyl sulfone reactions with OH, the presence of the water molecule reduces the barrier height of the H-abstraction channel. Also, the results indicated that the transition states are more stabilized in the DMSO + OH…H 2 O reaction than the DMS + OH…H 2 O reaction compared to the respective naked reactions 17 . However, in some simulated atmospheric reactions, an increase in the calculated rate constant was obtained for the water-assisted reactions compared to corresponding bare reactions. For example, the rate of H 2 SO 4 + OH reaction in the presence of water is higher than the naked reaction (about 10 3 times) 19 . In contrast, water not only cannot accelerate the rate of the NCO plus HCHO reaction, but also it has a negative influence on that reaction 21 . The same trend is also reported for the HNO 3 + OH reaction rate. The main reason is a rearrangement in the geometry of the hydrogen-bonded pre-reactive complex, which causes the kinetics of the reaction to be more complicated, and ultimately entire reaction process is slow 12 .
The simple and volatile organic compounds containing a sulfur atom in the field of heterocyclic chemistry, thiophene, and its derivatives, are emitted into the atmosphere by numerous sources such as the combustion of fossil fuels, burning plants, and biomasses 22,23 . The detection of sulfur compounds like thiophenes over seaweed fields and oceans was reported previously. Also, oxidation of these compounds occurs (with a high likelihood) near the surface of oceans in the marine boundary layer. It should be noted that the water concentration in this region is too high. So, it may form hydrated complexes with those compounds 17,24 .
It is worth mentioning that the removal of aromatic compounds is the main atmospheric issue. So, the oxidation of these chemicals plays a pivotal role in both atmospheric chemistry and the combustion process. In this regard, the oxidations of thiophene with active species have been studied widely [25][26][27] .
As we know, atmospheric compounds react preferably with hydroxyl radicals. So, the gas phase oxidation of thiophene is widely considered by hydroxyl radical experimentally and theoretically [27][28][29][30][31][32][33][34][35] . Some research showed that other active radical species such as 3 3 29,35 , NO 3 41-43 , and Cl 44,45 contribute to removing the atmospheric thiophene. Computationally, in the thiophene plus 3 O, 3 O 2 , and NO 3 reactions, the most suitable paths are energetically related to the addition/elimination reactions because of having a small barrier than the H abstraction reactions 26,39,43 . So, theoretical results proved that the main reaction channel is for the addition/elimination pathways. It is better to say that hydroperoxyl radical has high reactivity in the gas phase after hydroxyl radical. So, HO 2 plays an important role in the degradation of atmospheric pollutants which has been considered by others [46][47][48] . Also, the atmospheric hydroperoxyl radical play a pivotal role in hydrocarbon combustions 49 . Therefore, HO 2 is a key atmospheric species.
In addition, in atmospheric chemistry, the well-known oxidants are HO X family (OH and HO 2 radicals). They play a crucial role in the degradation of atmospheric pollutants due to stratospheric ozone. Also, these species are important in tropospheric oxidation reactions 50,51 . Furthermore, it is demonstrated that oxidation reactions with HO 2 are the main way to lose many important molecules and radicals in the troposphere 52,53 .
On one hand, according to the experimental and theoretical results, up to 30% of HO 2 forms the hydrated complexes with water, and neglecting it causes significant errors in the atmospheric models 8,9 . On the other hand, as mentioned above, a large portion of sulfur-containing aromatic compounds are in the hydrated form. So, water will play a significant role in the reactions of thiophene-based compounds when the oxidant is HO 2 .
In this study, we will present an example to show how treatment of aromatic compound reactions changes in a wet gas phase media (a fact that goes beyond this particular reaction). A straightforward approach to accomplish this aim is using simple atmospheric models. But, as we know, the reactions of aromatic compounds are complex. Thus, we select the thiophene reaction with hydroperoxyl radical that has been studied by the same authors in the dry environment 54 . We showed that this reaction begins with the addition of HO 2 to thiophene and proceeds through different channels. Here, we consider the impact of water on addition reactions and epoxy formation due to OH elimination only. First, we will construct the PES of the title reaction and investigate the influence of water on the relative energies of all stationary points. The geometries of species involved in the PES will be analyzed to find the main reason for the stability of water-assisted reactions. It will be shown that strong interactions due to water presence play an important role in the formation of ring-like structures in the course of the simulated reaction. Also, NBO analysis will be used for elucidating water's role in the stability of stationary points. We will use the AIM theory for clarifying the nature and strength of the interactions in ring-like structures. Second, the kinetic of the designed channels will be investigated to elucidate how water position impacts the rate of reaction. Finally, the results of wet and naked reactions will be compared for judging the water's role in the reactions of aromatic compounds. This study will aid the understanding and rationalization of similar processes that may occur in the gas-phase oxidation of other heterocyclic and aromatic compounds in the atmosphere.

Theoretical methods
We have previously shown that, for the thiophene plus HO 2 reaction, the B3LYP 55 method has good accuracy for simulation results 54 . Consequently, the potential energy surfaces for the addition of HO 2 to thiophene in the naked and water-assisted reaction are examined by the DFT-B3LYP method in conjunction with the Pople type basis set 6-31 + G(d,p). However, full optimization of the probable structure of all stationary points with and without the presence of water accompanied by vibrational frequencies is carried out at the DFT-B3LYP/6-31 + G(d,p) level of theory. The obtained frequencies along with the intrinsic reaction coordinate 56 (IRC) calculations are employed to verify the validity of the computed saddle pints structures and their connectivity at the mentioned level of theory. The optimized structures of all saddle points contain just an imaginary frequency that shows an alteration along the reaction coordinate. Also, each transition state is connected to the respective minimum Determining reactive centers in thiophene and hydroperoxyl radical by the Fukui function. In DFT, the Fukui functional is a favorite concept for determining a more reactive site on a molecule 71,72 . This function for a compound containing N electron is defined by where ρ(r) is the electron density.N is the number of atoms in the studying molecule. υ(r) is a constant external potential. This function is not evaluated directly because of a limitation associated with discontinuity 73,74 . Thus, the modified equations based on the involved mechanism in the investigating reaction are used for calculating the Fukui function 75 .
(a) For the nucleophilic attack, the following equation is used where P X (N + 1) and P X (N) are the population of atom X in molecule Y that denote the electron states with − and 0 signs, respectively, obtained from the geometry of neutral species. (b) For the electrophilic attack, it can be written where P(N − 1) is for the state with + signs. (c) For the radical attack, we have The population of an atom is obtained by subtracting the atomic number ( A ) of the target atom from its atomic charge ( Z) We apply the Hirschfeld charges in the calculation of the population of atom X. It is worth mentioning that the charge transfer mechanism is studied by using these formulas. Therefore, the Fukui function is a useful parameter to identify the reactive sites on a molecule. Tables S14-S21 contain the Fukui analysis results for all stationary points of the title reaction.

Results and discussion
Fraction of possible complexes formed between HO 2 and H 2 O, and thiophene and H 2 O molecules. In this section, we argue about the complexes that may take part in the title reaction. First, all possible complexes are found by closing water molecules from different sides to the hydroperoxyl radical and thiophene molecules separately. Then, using the kinetic theory of gases, the fraction percent of those complexes is calculated based on their stability. In this regard, the Boltzmann distribution is used as follows 76 : where N i is the number of molecules with energy E i and c is constant. k B is the Boltzmann constant. The fraction of molecules which has the energy E i is www.nature.com/scientificreports/ where N is the total number of molecules. This equation shows obviously that the most stable complexes will have a higher percentage. This is an important result since it helps us to find more abundance complexes and consider them for reaction design. Figure 1 shows all possible complexes between hydroperoxyl and water molecules. By the possible interactions that may form between HO 2 and H 2 O molecules, three stable H-bonded complexes are established and named 1WHA, 1WHB, and 1WHC. The most stable complex is 1WHA with a relative energy of − 9.99 kcal mol −1 and after that 1WHB (− 3.60 kcal mol −1 ). The origin of more stability for 1WA is the formed hydrogen bond ( ρ(lcp) = 0.0392 e bohr −3 and ∇ 2 ρ(lcp) 0.1153 e bohr −5 ) in which water is a hydrogen bond acceptor. In the 1WHB and 1WHC complexes, the water molecule is a hydrogen bond donor with the electronic charge density of ρ 1WHB (lcp) = 0.0192 and ρ 1WHC (lcp) = 0.0196 e bohr −3 , and its Laplacian ∇ 2 ρ 1WHB (lcp) = 0.0527 and ∇ 2 ρ 1WHC (lcp) = 0.0492 e bohr −5 . The percent of the 1WHA fraction is almost 100% at 300 K. This result proposed that all HO 2 molecules at room temperature are complexed with water in the 1WHA form. Also, Table 1 shows that this complex has a dominant fraction over the 300-3000 K temperature range.
As abovementioned, thiophene can be hydrated by water molecules. Thus, the probable complexes between thiophene and water are explored. The sketch of the complexes is depicted in Fig. 1. By closing a water molecule to thiophene, four different complexes are yielded named 1WTA, 1WTB, 1WTC, and 1WTD. These complexes have also near relative energies. The obtained relative energies forecast that the 1WHA and 1WHB (HO 2 … H 2 O) complexes are more stable than the thio… H 2 O complexes under the same conditions. This is related to the nature of interactions among the interacting molecules. In the 1WTA and 1WTB complexes, there is a van der Waals interaction between nonbonding electrons of water with α hydrogen ( ρ(lcp) = 0.0122 e bohr −3 and ∇ 2 ρ(lcp) = 0.0395 e bohr −5 ) and β hydrogen ( ρ(lcp) = 0.0109 e bohr −3 and ∇ 2 ρ(lcp) = 0.0361 e bohr −5 ) of thiophene. For the 1WTC complex, the interaction between the hydrogen of water with π electrons of thiophene is also van der Waals type ( ρ(lcp) = 0.0098 e bohr −3 and ∇ 2 ρ(lcp) = 0.0271 e bohr −5 ). Finally, there is a 1WTD Figure 1. Possible complexes established between hydroperoxyl radical and water molecules, and also thiophene and water molecules (unit of all distances is angstrom). , respectively. It should be noted that the 1WTC complex during the reaction when HO 2 is close to interacting with it will convert to the 1WHX complex (X = A, B, and C) due to the orientation of the complexed water molecule. Thus, we skip its assistance for reaction.

Scientific Reports
The gas-phase reaction of C 4 H 4 S + HO 2 without and with water. All possible complexes between reactant species and water are displayed in Fig. 1. The main parameters of the optimized geometries for all reaction components computed by the B3LYP method are shown in Fig. 2. The potential energy surfaces (PESs) for the title reaction in the presence and absence of water are drawn in Fig. 3. the rate constants of all paths are depicted in Fig. 4 and are listed in Tables S9 and S10. The energetics of each stationary point is compared by considering the original reactants as the base of comparison and are listed in Tables 2 and 3. Also, for all interacted molecules, the extracted T1 diagnostic values from the CCSD(T) calculations are listed in Tables 2 and 3. The useful thermodynamic parameters, such as inner energies, enthalpies, Gibbs free energies, and entropies computed in the standard conditions are collected in Table 4 Table S1 of Supporting Information). Thus, the obtained energetics are acceptable. The optimized stationary points of the C 4 H 4 S + HO 2 …H 2 O reaction have ring-like structures due to intermolecular interactions, then the AIM analysis is utilized to extract beneficial information about the line and ring critical points. The results of AIM are reported in Table S2.
In the preceding study 54 , we offered all probable pathways for thiophene plus hydroperoxyl radical reaction by similar theoretical methods. It was proven that the reaction starts with the addition of HO 2 to the α (more preferable) or β carbon atoms of the thiophene ring. The electrophilic attack of oxygen on the S atom was investigated as well. The reaction is simulated in the presence of a single water molecule to have an insight into the influence of a wet environment on activation energies and so the final rate constants.
To understand the role of water vapor in the title oxidation reaction, all electronic structure calculations are carried out with and without water to construct the respective potential energy surfaces (PESs). Comparing the PESs clears how water changes the reactivity of thiophene and hydroperoxyl radical when reacts with each other. For this, we must have the same reaction pathways for bare and water-assisted reactions. According to our previous research 54 , the HO 2 -addition pathways include two pre-reactive complexes. After the formation of the van-der-Waals complexes and passing through desired transition states, the path is completed. But, the presence of water is accompanied by more complexity due to several pre-reactive complexes having different relative energies and structures. In the sequel, we focused on the wet reaction and compare its results with the bare reaction.
Addition to β-C in the presence of water. When 1WHA and thiophene, and 1WTA and HO 2 are gathered, two pre-reactive complexes CR-TS1-1WHA and CR-TS1-1WTB are formed, respectively, that are starting points of the β-C addition reaction. Since water can form ring-like structures, there is large stability in these complexes. The relative energies of CR-TS1-1WHA and CR-TS1-1WTB are − 16.24 and − 11.20 kcal mol −1 , respectively, which are − 10.32 and − 5.28 kcal mol −1 more stable than CR1(bare complex). This result is important because water position determines how much energy is released. Then, according to water position, two transition states TS1-1WHA and TS1-1WTB are found for the addition of HO 2 to β-C. The TS1-1WHA saddle point has a fivemembered ring structure which is stabilized by two hydrogen bonds. In these bonds, water and HO 2 play both hydrogen bond donor and hydrogen bond acceptor roles simultaneously. The relative energy of TS1-1WHA is 5.95 kcal mol −1 , which is 9.65 kcal mol −1 more stable than the naked reaction due to the ring structure containing hydrogen bonds. Also, for the TS1-1WTB, a six-membered ring structure is obtained including one hydrogen bond and two van der Waals interactions. So, the stability of this transition state is slightly less than the TS1-1WHA. After the mentioned transition states the product complexes CP-TS1-1WHA1 and CP-TS1- Table 1. The relative energy (kcal mol −1 ) and the percentage of fractions for the probable complexes formed between HO 2 and H 2 O, and thiophene and H 2 O computed at the BD(T) method.

Reaction Coordinate
Thio+1WHA   A decrease in the value of f 0 X (X = O and C) from reactants to the TSs is related to the formation of a covalent bond. This means that after the formation of a covalent bond between the reacting centers, the reactivity of these centers is decarded.
In the sequel, the reaction can be continued by the TS5-1WHA and TS5-1WTB saddle points. Firstly, without any barrier, the CP-TS1-1WHA1 and CP-TS1-1WTB1 can rearrange to the intermediate complexes CP-TS1-1WHA2 and CP-TS1-1WTB2, respectively. After that, an epoxy structure can obtain if the oxygen-oxygen bond of HO 2 breaks via TS5-1WHA or TS5-1WTB. The ring structure of TS5-1WHA has five members containing two hydrogen bonds and one van-der-Waals interaction. The obtained values for the electronic charge density along with its Laplacian at the ring critical point are ρ(rcp) = 0.2663 e bohr −3 and ∇ 2 ρ(rcp) = 0.0569 e bohr −5 . Also, the existence of a six-membered ring with two hydrogen interactions and one van-der-Waals interaction in TS5-1WTB is confirmed by AIM analysis ( ρ(rcp) = 0.0878 e bohr −3 and ∇ 2 ρ(rcp) = 0.0371e bohr −5 ). The IRC path for the epoxy formation shows that the CP-TS5-1WHB and CP-TS5-1WTB are two final product complexes with some differences in their internal interactions. These differences are accompanied by relatively different relative energies for CP-TS5-1WHB (− 21.70 kcal mol −1 ) and CP-TS5-1WTB (− 20.32 kcal mol −1 ). It is better to note that CP-TS5-1WHB has a seven-membered ring structure ( ρ(rcp) = 0.0033 e bohr −3 and ∇ 2 ρ(rcp) = 0.0129 e bohr −5 ). AIM result demonstrated that in addition to two hydrogen bond interactions, this cyclic structure has an interaction between the O11 atom and the pi bond of C4 and C3 atoms. The CP-TS5-1WTB complex has also a seven-membered ring structure involving two hydrogen bonds. The computed AIM parameters for the CP-TS5-1WTB ring structure are ρ(rcp) = 0.0030 e bohr −3 and ∇ 2 ρ(rcp) = 0.0152 e bohr −5 . Also, for the discussed product complexes, the existence of an epoxy ring (a three-membered ring structure) between C1, O10, and C2 atoms is evaluated by AIM. The amounts of ρ(rcp) and ∇ 2 ρ(rcp) for the epoxy ring in CP-TS5-1WHB and CP5-1WTB are 0.2023 and 0.2003 e bohr −3 and 0. 3171 and 0.3172 e bohr −5 , respectively. The final products of these complexes are an oxirane-like bicyclic product plus hydroxyl radical (and water). From the thermodynamic point of view, the formation process of P5 (SCHOCHC 2 H 2 + OH) is exothermic by -5.22 kcal mol −1 in standard enthalpy and spontaneous by -2.72 kcal mol −1 in standard Gibbs free energy.
In the other pathway, the HO 2 addition to β-C is occurred by the OH part of HO 2 . This reaction starts by attacking the HO 2 moiety to the β-C of the IWTB complex or the addition of the 1WHA moiety to the β-C of    www.nature.com/scientificreports/ And there is only a strong interaction between shifting hydrogen and oxygen atom of water which is hydrogen type ( ρ(lcp) = 0.0376 e bohr −3 and ∇ 2 ρ(lcp) = 0.1028 e bohr −5 ). Also, the CP-TS4-1WTB has a seven-membered ring structure with one hydrogen bond (H14…O12) and one van der Waals interaction (O13…H7). The ring structure of the post-reactive complex CP-TS4-1WHA is similar to the CP-TS1-1WHA complex.
Addition to α-C in the presence of water. Before describing the reaction paths, it should be noted that the Fukui value for the α-C center ( f 0 α−C ) is − 0.143 e which is higher than f 0 β−C (see Tables S18 and S19). This result indicates that the α-C center is more reactive than the β-C center. The addition to α-C is considered by two paths. The orientation of the water molecule makes the distinction between these paths. The first starts with CR-TS1-1WHA and occurs through TS2-1WHA. TS2-1WHA has a ring-like structure with five members containing two hydrogen bonds and one van der Waals interaction. This ring was analyzed by AIM theory. For the ring critical Table 4. The calculated thermodynamic parameters (298.15 K) for all components of the water-assisted reactions using the BD(T) energies and the B3LYP corrections (unit of all numbers is kcal mol −1 ). www.nature.com/scientificreports/ point, AIM parameters are ρ(rcp) = 0.0154 e bohr −3 , and ∇ 2 ρ(rcp) = 0.0831e bohr −5 . For this reason, it is around 10.13 kcal mol −1 more stable than the TS2 (the corresponding transition state in the naked reaction). After TS2-1WHA, the post-reactive complex CP-TS2-1WHA is obtained in the IRC path with the relative energy of − 13.63 kcal mol −1 which is 9.50 kcal mol −1 more stable than the P2 (final) product. In the second path, CR-TS2-1WTA is the key per-reactive complex for the terminal oxygen (of HO 2 ) addition to the α carbon of thiophene. Geometrically, it has a ring-like structure with six members, including two covalent bonds, one hydrogen bond, and two van der Waals interactions. This route is completed after surmounting the saddle point TS2-1WTA and the formation of the post-reactive complex CP-TS2-1WTA in the exit channel. The saddle point TS2-1WTA is 3.99 kcal mol −1 above the initial reactants. Also, the product complex CP-TS2-1WTA is 10.08 kcal mol −1 more stable than the initial reactants. The six-membered ring structure is formed in CP-TS2-1WTA. In this ring, the O13 of water has a van der Waals interaction with the hydrogen of thiophene, and also the hydrogen of water acts as a hydrogen bond acceptor with the middle oxygen of HO 2 . The P2 and water are exothermic and nonspontaneous products by − 2.77 and 8.53 kcal mol −1 in enthalpy and Gibbs free energy, respectively. Finally, it should be noted that between the α-C and β-C additions, the addition process to α-C is favorable both thermodynamically and kinetically. It is better to say that a similar tendency was also established in the addition of the oxygen atoms of NO 3 and O 2 radicals to thiophene 26,43 .
It has been shown that in the bare reaction, the thiophene ring-opening and elimination of the OH group may simultaneously happen after the addition of HO 2 and generation of the P2 adduct, (a two-step reaction). In the presence of water, this process occurs through two pathways. In the first path, CP-TS2-1WHA rearranges to IN8a-1WHA and after passing through TS8a-1WHA converts to CP-TS8a-1WHA, and by dissociating the fragments of CP-TS8a-1WHA, the mentioned adduct is produced. The second path has a similar trend but the position of the water is different. respectively. The geometry of CR-TS2-1WTA has been discussed above. The CR-TS3-1WHA complex has an eight-membered ring-like structure in which the primary 1WHA complex is connected to thiophene by two van-der-Waals interactions. One of those interactions is located between the S atom of thiophene and the H14 atom of water, and the other is established between the O10 atom of HO 2 and the H2 atom of thiophene (see Fig. 2).
From an energetic point of view, there is high stabilization for transition states after the addition of water. The calculated relative energies for TS3-1WTA and TS3-1WHA are predicted to lie 10.72 and 9.36 kcal mol −1 below the energy of TS3, respectively. Besides, CP-TS3-1WTA and CP-TS3-1WHA are 14.41 and 14.86 kcal mol −1 below CP3, indicating more stability for post-reactive complexes than respective TSs. The main reason for the differences in the stabilization energies of the mentioned transition states goes back to their ring structures. The seven-membered ring ( ρ(rcp) = 0.0040 e bohr −3 and ∇ 2 ρ(rcp) = 0.0184 e bohr −5 ) in TS3-1WTA includes three van-der-Waals interactions and one hydrogen bond. The van der Waals interactions are between O13 (of water) and H6 (of thiophene), O11(of HO 2 ) and O10 (of HO 2 ), and O10 and S5 atoms. The hydrogen bond is between O11 and H14 atoms. But, in the TS3-1WHA structure, a six-membered ring is confirmed ( ρ(rcp) = 0.0044 e bohr −3 and ∇ 2 ρ(rcp) = 0.0195 e bohr −5 ) by AIM data with two van der Waals interactions and one hydrogen bond. The hydrogen atom (H14) of water also has a hydrogen bond with the O11 atom, and O11 has a van-der-Waals interaction with O10. In the TS3-1WHA structure, a six-membered ring is confirmed by AIM analysis with two van der Waals interactions and one hydrogen bond. The van der Waals interaction in TS3-1WHA is similar to TS3-1WTA, but the hydrogen bond is slightly different. The hydrogen atom of the water molecule acts as a donor and interacts with the O10 atom as an acceptor in the newly formed hydrogen bond. The characteristics of the new bond are ρ(lcp) = 0.0287 e bohr −3 and ∇ 2 ρ(lcp) = 0. 0827 e bohr −5 .
Similar trends are true about the energies of the product complexes CP-TS3-1WTA and CP-TS3-1WHA in comparison with CP3. The optimized structures of the CP-TS3-1WTA and CP-TS3-1WHA also indicate different ring structures in these complexes. Also, the calculated densities for these rings by AIM theory are ρ(r rcp ) = 0.0018 e bohr −3 ( ∇ 2 ρ(rcp) = 0.0087 e bohr −5 ) for CP-TS3-1WTA and ρ(rcp) = 0.0033 e bohr −3 ( ∇ 2 ρ(rcp) = 0.0131 e bohr −5 ) for CP-TS3-1WHA, respectively. The CP-TS3-1WTA has an eight-membered ring structure due to three established interactions. In this ring, the H12 atom of hydroxyl radical interacts with the connected oxygen (O10) to the sulfur atom. Also, the water molecule has one van der Waals interaction with the hydrogen of thiophene and one hydrogen bond with hydroxyl radical ( ρ(rcp) = 0.0309 e bohr − www.nature.com/scientificreports/ including two hydrogen bonds and one van der Waals interaction. The hydrogen atom of water has an interaction with the attached oxygen to the sulfur atom. Also, the hydrogen of the OH radical interacts with the oxygen of the water molecule. In addition, the radical site of OH interacts with the C2 atom of thiophene (the van der Waals interaction). Finally, these channels are ended by rupturing CP-TS3-1WTA and CP-TS3-1WHA to thiophene-S-oxide, hydroxyl radical, and water moieties. The computed standard thermodynamic parameters for this reaction show that the formation process of the final adduct (P3) is endothermic by 14.35 kcal mol −1 and nonspontaneous by 16.10 kcal mol −1 in Gibbs free energy in the gas phase.

Natural bond orbital (NBO) analysis.
To interpret the charge transfer processes 54 and characterize the inter-and intra-molecular interactions in molecular systems, a suitable quantum chemistry method is the NBO analysis. The energies of such interactions are expressed by the second-order perturbation theory between a donor NBO (i) and an acceptor NBO (j). The second-order stabilization energy E 2 for delocalization of electrons between i and j NBOs can be calculated as: Table 5. Analysis of the Fock matrix using the NBO calculations for some species in the thiophene plus hydroperoxyl radical reaction with and without a single water molecule. www.nature.com/scientificreports/ where q i is the orbital occupancy, ε i and ε j are diagonal elements of the Fock matrix, and F i,j is the off-diagonal elements 77 . The results of the Fock matrix obtained by the NBO calculations for some stationary points are listed in Table 5. By referring to part I 54 , we stated that two strong interactions exist in thiophene due to the n S(2) → π*(C 1 -C 2 ) (E 2 = 21.26 kcal mol −1 ) and n S(2) → π*(C 3 -C 4 ) (E 2 = 21.28 kcal mol −1 ) transitions. In encountering thiophene with HO 2 in wet (and dry) media, the energies of n S(2) → π*(C 1 -C 2 ) and n S(2) → π*(C 3 -C 4 ) transitions in CR-TS1-1WHA (and CR1) are 21.92 (and 21.50) and 21.40 (and 21.99), and in TS1-1WHA (TS1) are 5.68 (5.46) and 10.90 (10.88) kcal mol −1 respectively. The stabilization energies (E 2 ) of the same transitions in TS1-1WTB are 5.63 and 10.70 kcal mol −1 respectively. These results show that for β addition, water alters the n S(2) → π*(C 1 -C 2 ) and n S(2) → π*(C 3 -C 4 ) transitions. There are also transitions between the lone pairs of the 11O atom and the unoccupied π*(C 1 -C 2 ) and π*(C 3 -C 4 ) orbitals in TS1-1WHA (TS1) and TS2-1WHA (TS2), respectively. These transitions vary the interactions of the lone pair electrons of the S atom and the unoccupied π*(C 1 -C 2 ) and π*(C 3 -C 4 ) orbitals. Thus, the activation energies of that TSs are changed. The stabilization energies E2 value for the n O10(3) → π*(C 1 -C 2 ) transition in TS1-1WHA and TS1-1WTB (and TS1) are 47.94 and 44.88 (48.00) kcal mol −1 , respectively.
Overall, all the new interactions by a water molecule can influence directly and indirectly the transitions between the π*(C 1 -C 2 ) and π*(C 3 -C 4 ) orbitals and lone pair electrons that can change the stability of the thiophene ring. Finally, we can conclude that in the presence of a single water molecule the E 2 value is reduced for the discussed stationary points, relating to the involved intermolecular interactions (particularly hydrogen bonds by water molecules).
Kinetic rate constant calculation. To find the effect of an added single water molecule on the kinetics of the title reaction, all rate constants of the studied pathways are calculated in terms of the canonical transition state theory 78 . That paths start with the formation of a pre-reactive complex and then pass through just a transition state and are finally released to the final products can be considered by the following reaction mechanism where R 1 and R 2 are reactants. k 1 and k −1 are the rate constants of the association of reactants and dissociation of the pre-reactive collision complex, respectively. Also, k 2 is the rate constant for the formation of final products from CR. The steady-state approximation should be used to calculate the overall reaction rate. Thus, the following rate expression by considering the equilibrium between CR and reactants is obtained.
where the K eq is the equilibrium constant of the first step in the reaction of eq. 4. The following equations are applied to calculate the rate constant of two-step reactions where k B , T, h κ , and σ are the Boltzmann constant, temperature, the Planck constant, the tunneling correction, and the symmetry factor, respectively. Q R 1 Q R 2 ...H 2 O , Q CR , and Q TS represent the partition functions of the reactants (that is shown by the product of separated reactants), pre-reactive complex, and transition states, respectively. The E R 2 ...H 2 O + E R 1 , E CR , and E TS indicate the total energy of the reactants, the complexes, and transition, respectively.
The tunneling effect is a quantum mechanical phenomenon that is defined as the ratio of barrier crossing rates obtained by quantum mechanics and classical mechanics. In this study, the unsymmetrical Eckart potential www.nature.com/scientificreports/ energy barrier has been used to simulate approximately the potential energy curve with high accuracy during tunneling 79 . The Eckart formula which is a potential function uses an unsymmetrical barrier to estimate the tunneling factor for computing precise rate constants theoretically. The Eckart formalism uses the following integral over the threshold energy E to calculate the tunneling factor.
where E b is the barrier energy of the reaction path, and P T (E) is the tunneling transmission probability.
The tunneling transmission depends on three parameters that are as follows where ΔE is calculated from the difference of energy between the sum of products and the sum of reactants which is called the reaction energy. The definition of other parameters is found elsewhere 80 . Also, we use the statistical Rice-Ramsperger-Kassel-Marcus (RRKM) theory to compute the rate constants of two-step reactions, which are considered to be continued by unimolecular reaction channels 81 .
To gain more useful insight into the influence of a single water molecule on the title reaction, all rate constants for the naked and water-assisted reactions are calculated by the abovementioned theories. For the abovediscussed paths, the calculated rate constants are denoted by k X as follows:  Table S9 and S10. Also, for easier comparison, the ratios of rate constants in the presence and absence of water are collected in Table 6.
As seen in Fig. 4, for the calculated rate constants at the mentioned temperature range, the non-Arrhenius behavior is observed. Therefore, the three-parameter equation k = A( T 300 ) m exp(− E b RT ) is used to fit the calculated temperature dependence of rate constants. All the following rate expressions are bimolecular in the unit of L mol −1 s −1 . Table 6. The ratios of the computed rate constants at the BD(T) level for the water-assisted paths to the respective ones in the bare reaction. www.nature.com/scientificreports/ Also, the plots of these fitted equations are shown in Fig. 4. These results demonstrated that the water-assisted reaction has a similar trend the same as the naked reaction. Our computed rate constants show that the addition to the α carbon, α-C, is more favorable than β-C in the studied temperature range. For access to more detail about the naked reaction please see our previous paper 54 . It should be noted that according to the TST theory, in addition to the barrier energy of a path, the ratio of the total partition function of a saddle point to the product of the partition function of reactants is important. It shall be emphasized to this point that this ratio for the water-assisted reactions is smaller than the corresponding bare reactions.
As aforementioned, a single water molecule has different effects due to different positions. So, the rate constants k P1-1WHA and k P1-1WTB have different behavior by temperature variation. k P1-1WTB is more favorable than k P1-1WHA in a temperature range of 300-1600 K but above 1600 K, the rate constant k P1-1WHA is more. Also, the ratio of rate constants (k water /k nacked ) for k P1-1WHA /k P1 (and k P1-1WTB /k P1 ) is 3.58E−02 (1.01E+02), 4.97E−02 (2.92E+00), and 5.90E−02 (3.57E−01) at 300, 500, and 800 K, and is 6.23E−02 (1.88E−01), 6.69E−02 (7.41E−02), and 7.20E−02 (2.91E−02) and, at 1000, 1500, and 3000 K, respectively. These results show a meaningful difference in the values of computed rate constants for the P1 generation pathways in wet and dry media at high temperatures. Other paths to produce P1 include TS4, TS4-1WHA, and TS4-1WTB. The calculated rate constants for these paths (k P1′, k P1′-1WHA, and k P1′-1WTB ) are small due to having large energy barriers. It is better to say that the P1 formation through the mentioned paths in the presence and absence of water is not kinetically important up to 1500 K since the ratio of rates of k P1′-1WHA (k P1′-1WTB ) to k P1′ is small. This shows that this path even in the presence of water does not occur at low and moderate temperatures. The value of the rate constant for the P2 product, k P2-1WTA, is more than all pathways of the water-assisted reactions that coincide with the corresponding naked reaction because addition to the α carbon is the main reaction pathway among all addition reactions. Also, the computed rate constants for the HO 2 addition to the S atom (k P3-1WHA and k P3-1WTA ) in the water-assisted reaction indicate that the same as the respective naked reaction (k 3 ), this reaction does not occur at low temperatures. However, the computed rate constant for k P3-1WTA is higher than k P3-1WHA up to 800 K and after that the rate of P3 formation via the path that is procced by the 1WHA complex is high. The ratio k P3-1WHA /k P3 (k P3-1WTA /k 3 ) is 5.57E−02 (1.07E+03), 1.58E−01 (1.30E−01), 2.16E−01 (943E−03), 2.59E−01 (2.10E−03) at 300, 800, 1500, and 3000 K, respectively (see Table 6).
For two-step reactions, the rate constant of product formation is computed in both dry and wet media (see Tables S9 and S10). For P4 generation, the ratio of rate constant in the presence and absence of a water molecule k P4-1WHA /k 4 (k P4-1WTA /k 4 ) is 3.35E−02 (2.66E+01), 8.23E−03 (4.72E−02), and 8.68E−02 (3.02E−03) at 300, 1500, and 3000 K, respectively. These results demonstrate that the same as other paths the product formation through 1WHA formation has a lower rate constant at low temperatures. Regarding the ratio k P5-1WHA /k 5 , the results indicate that these ratios are smaller than others at low temperatures. The obtained value for k P5-1WHA /k 5 is 3.61E−03, 3.99E−03, 7.11E−02 and at 300, 400, and 500 K, respectively. But, for k P5-1WTB /k 5 the same values are 1.97E+02, 1.52E−01 and 3.03E+00, respectively. Finally, it can be concluded that there is a remarkable difference in the computed rate constants of a product formation in wet media due to the orientation of the water molecule. Also, there is a meaningful difference between the computed rate constants of the argued paths in wet and dry media because of a change in barrier energy over the 300-3000 K temperature range.

Conclusion
This work was accomplished to prove the crucial role of a water molecule in atmospheric reactions in which one component is an aromatic compound. To establish the impact of water, the addition mechanism in the HO 2 … H 2 O + thiophene system was selected as an example. The reaction begins with the association of the 1WHA and thiophene and 1WTA and thiophene, and then the formation of different pre-reactive complexes in the initial stage of addition reactions. Then, the reactions go on with adding HO 2 radical to the α and β carbons and also the sulfur atom. The obtained reaction pathways showed that the reaction mechanism with a single water molecule is more complicated than the naked reaction. Our results made clear that the presence of water molecules leads to the creation of ring-like structures. These ring structures were confirmed by AIM analysis and were observed due to hydrogen bonds and van der Waals interactions in all prereactive, post-reactive, and transition state structures. The main consequence of the strong interactions in ring structures was to lead to high stabilization of the stationary points. is similar for all paths, the amount of stabilization energy due to the presence of water in the 1WHA complex is higher than that of the water-assisted TSs. Thus, these paths have a smaller rate constant than the respective naked reactions. Despite this, the inverse treatment is seen for the paths that start with the 1WTA and 1WTB complexes. Therefore, the results of this study show that the presence of water accelerates or does not accelerate the rate constant of the addition pathways which is related to the reacting primary complex formed between water and one of the reactants.