# Prediction of Mechanism and Thermochemical Properties of O3 + H2S Atmospheric Reaction

Authors: Morteza Vahedpour; Reza Baghary; Freshte Khalili
Journal: Journal of Chemistry
 (2013)

---

## Abstract

Ozone and hydrogen sulfide reaction mechanism including a complex was studied at the B3LYP/6-311++G(3df,3pd) and CCSD/6-311++G(3df,3pd)//B3LYP/6-311++G(3df,3pd) levels of computation. The interaction between sulfur atom of hydrogen sulfide and terminal oxygen atom of ozone produces a stable H2S-O3 complex with no barrier. With the decomposition of this complex, four possible product channels have been found. Intrinsic reaction coordinate, topological analyses of atom in molecule, and vibrational frequency calculation have been used to confirm the suggested mechanism. Thermodynamic data at T = 298.15 K and the atmospheric pressure have been calculated. The results show that the production of H2O + SO2 is the main reaction channel with ΔG° = −645.84 kJ/mol. Rate constants of H2S + O3 reaction show two product channels, SO2 + H2O and HSO + HOO, which compete with each other based on the temperature.


---

## Body

## 1. Introduction

Ozone and hydrogen sulfide are two important atmospheric pollutants. Ozone as a reactive species plays an important role in the upper atmosphere. In view of the relatively large ozone concentration in the middle atmosphere and the importance of ozone chemistry in the depletion of ozone layer, by a series of reactions with atmospheric species, its gas phase reactions are of interest. Recently, the reaction of ozone with atmospheric molecules and radicals has been studied extensively in the theory and experiment [1–6].Sulfur compounds are released into the atmosphere in various forms by means of biogenic (biomass rotting and volcanic eruptions) and anthropogenic (chemical industry) sources [7]. Hydrogen sulfide is one of the sulfur compounds that play an important role in the chemistry of the atmosphere [8]. Also, it is a large source of sulfur in the atmosphere. Hydrogen sulfide is highly toxic and is used in a number of chemical industries [9]. Therefore, its atmospheric reaction to turn other compounds with lower toxicity is a topical subject in a number of studies.The reaction of hydrogen sulfide with different molecules and radicals in the stratosphere was the subject of many experimental and theoretical studies [10–13]. In experimental studies, for many years, most of hydrogen sulfide is removed by the well-known Claus process [2, 12]. Also, H2S reacts with ozone to produce SO2 and H2O, which potentially is one of the sources for acid rains [2, 3]. In theoretical work, kinetics and thermodynamics of hydrogen sulfur atmospheric reaction with some radical species such as CH3, OH, and NO2 are studied [1].The reaction of ozone and H2S can aid to remove the toxic pollutant which is important in atmospheric chemistry. Moreover a slow gas-phase reaction of ozone with hydrogen sulfide may be significant if the reaction proceeds to produce (HSOH) and (HSOOOH) that may be an important source of OH, HS, and HOO (hydroperoxyl) radicals which are well known as oxidizing agents and take part in the processes of the destruction of ozone. Therefore, in this work, we have computationally investigated the details of the mechanism and thermodynamics of the hydrogen sulfide and ozone reaction. Also, the hydrogen thioperoxide molecule is indicated as the other new product of the reactions.Computational Methods. The geometries of the reactants, complexes (denoted as Cs), products, and transition states (denoted as TSs) involved in the reaction are optimized using density functional theory at the spin unrestricted B3LYP with the 6-311++G(3df,3pd) bases set. To obtain more reliable relative energy higher levels of electronic correlation method,  CCSD(T)/6-311++G(3df,3pd)//B3LYP/6-311++G(3df,3pd), is employed in the single point energy calculation to improve the accuracy of energetic information on minimum energy path. The molecular energy at the CCSD(T) level and the zero-point energy (ZPE) which is obtained from the B3LYP level are reported. Harmonic vibrational frequencies were obtained at the B3LYP level for verification of the optimized geometries. The nature of the stationary points and transition states are determined according to the number of negative eigenvalue of the Hessian matrix. Additionally, intrinsic reaction coordinate (IRC) was performed to examine connections between all the species involved in the reaction. The counterpoise procedure is used to correct the interaction energy for basis set superposition error, and CCSD(T) method is used to calculate the T1 diagnostic values for the structures.At the level of this method, it is possible to generate a wave function in a suitable form to execute a topological analysis of atoms in molecules [14] using the AIM2000 series programs [15]. Moreover, thermodynamic data have been calculated using the statistical mechanics. The rate constants of the reaction channels for thermodynamic and kinetic products have been calculated using TST and RRKM theories due to the B3LYP/6-311++G(3df,3pd) method. All calculations reported in this work were carried out with the Gaussian 03 program package [16].


## 2. Result and Discussion

The optimized geometries of reactants, complexes (as denoted C and a number), transition states (as denoted TS and a number), and products with the corresponding available experimental data for some species are represented in Figure1. The results of geometry optimization for all species that were obtained at the B3LYP/6-311++G(3df,3pd) level of theory in this work are in agreement with experimental values.Figure 1
Optimized geometries of reactants, products, complexes, and transition states of studied reaction in the B3LYP/6-311++G (3df,3pd) level (bond lengths is in Å and bond angels and dihedral angles are in degree).One reactant complex, H2S-O3 (denoted as C1), has been formed between the reactants on the singlet potential energy surface. In the complex, the terminal oxygen atoms of O3 are associated with sulfur atom of H2S and a four-member ring is formed. Molecular geometry indicates that the C1 complex has an asymmetric structure.In the C1 complex, the bonds length of newly formed S–O1 and S–O2 are 3.371 and 3.235 Å at B3LYP level, respectively, and the bond of O–O is about 0.004 Å longer than that of the corresponding bond in the parent ozone molecule. To identify newly bonds formation and existence of ring structure, we have used the topological analysis of electronic charge density. The values of bond critical point parameters and electronic density gradient for newly formed bonds and ring critical point for ring structure are schematized in Figure2. The Laplacians of electronic charge density value of S–O1 and S–O2 bonds are 0.0262 and 0.0265 au, respectively. Both low and positive values of Laplacian from AIM analysis for S–O bonds of C1 complex show that there is a weak van der Waals interaction between sulfur atom of H2S and terminal oxygen atoms of ozone. The electronic charge density corresponds with the newly bonds are 0.0095 and 0.0685 au, respectively. The calculated electronic charge densities confirm the weak interaction of ozone and hydrogen sulfide in C1 complex.Figure 2
The values of the AIM theory topological parameters for some selected bonds of species in reaction.The binding energy of C1 is 5.78 kJ/mol below the original reactants (H2S + O3) at B3LYP level. The base set superposition error (BSSE) value for the complex C1 energy is 0.01 kJ/mol. No transition state has been found for its formation. So, the formation of C1 is a barrierless process. Subsequently, via variety of C1, three kinds of products are obtained. In these paths, we have found five transition states, TS1, TS2, TS3, TS4, and TS5, on the potential energy surface of H2S + O3 reaction and another one, TS6, which is described as the transformation of SOO to SO2 and four other stationary points (C2, C3, HSOOOH, and SOO). As a matter of fact, the corresponding molecular graph (Figure 2) shows the existence of a BCP between bonded atoms. Also, ring structures are formed when ring critical points appear in all species. The total energies and relative energies of all species at the B3LYP/6-311++G(3df,3pd), under atmospheric pressure and 298.15 K are listed in Table 1 and CCSD(T) relative energies are shown in Figure 3. Higher levels of electronic correlation method CCSD(T)//B3LYP are employed in the single point energy calculation to improve the accuracy of energetic information on minimum energy path. Harmonic vibrational frequencies were obtained at the B3LYP level for verification of the optimized geometries. Also the results of the single point energy calculations at the CCSD/6-311++G(3df,3pd)//B3LYP/6-311++G(3df,3pd) level of computation are summarized in Table 1. All energies in Table 1 are corrected by BSSE.Table 1
Total energies (Hartree), relative energies (kJ/mol) (in parentheses) including ZPE corrections at the B3LYP/6-311++G(3df,3pd) and CCSD(T)//B3LYP methods, ZPE (Hartree), and T1 diagnostic for various species.
species
B3LYP/6-311++G(3df,3pd)
CCSD(T)/6-311++G(3df,3pd)
ZPE
T1
H2S + O3
−624.9053 (0)
−624.0099 (0)
0.0224
TS1
−624.8647 (106.53)
−623.9442 (172.63)
0.0237
0.099
TS2
−624.8809 (63.65)
−623.9678 (110.69)
0.0250
0.047
TS3
−624.8880 (45.33)
−623.9768 (86.92)
0.0235
0.043
TS4
−624.9069 (−4.41)
−624.0015 (22.14)
0.0262
0.041
TS5
−624.8634 (110.05)
−623.9495 (158.65)
0.0246
0.047
C1
−624.9063 (−2.74)
−624.0121 (−5.78)
0.0239
0.023
C2
−624.8890 (42.76)
−623.9869 (60.40)
0.0274
0.042
C3
−624.9690 (−167.30)
−624.0811 (−186.98)
0.0293
0.029
H2SO–O2
−624.9468 (−109.12)
−624.0741 (−168.55)
0.0270
0.014
HSOOOH
−624.9551 (−130.76)
−624.0678 (−151.96)
0.0290
0.021
TS6 + H2O
−624.8529 (137.47)
−623.9998 (26.46)
0.0248
TS6
−548.4097
−547.6707
0.0035
0.046
SOO
−548.5234
−547.7462
0.0056
0.027
SO2
−548.7090
−547.9449
0.0070
0.022
H2SO + 3O2
−625.0074 (−268.10)
−624.1236 (−298.39)
0.0262
H2SO + 1O2
−624.9461 (−107.25)
−624.0666 (−148.89)
0.0262
HSO + HOO
−624.9478 (−111.74)
−624.0812 (−187.24)
0.0245
H2O + SO2
−625.1521 (−648.19)
−624.2740 (−693.29)
0.0283
H2O + SOO
−624.9666 (−160.91)
−624.0753 (−171.62)
0.0269Profile of calculated potential energy surface of H2S + O3 reaction at CCSD(T)//B3LYP level of computation.
(a)
(b)The vibrational frequencies of all stationary points along with the available experimental values are listed in Table2. As can be seen, the calculated frequencies at B3LYP method are in acceptable agreement with available experimental data. All the reactants, products, and complexes have only real frequencies, while TS is identified with one and only one negative eigenvalue of the Hessian matrix.Table 2
Vibrational frequencies of all species at B3LYP/6-311++G(3df,3pd) level of theory.
Species
Frequency (cm−1)
H2S
1210 (1183), 2695 (2615), 2708 (2626)
TS1
326i, 102, 198, 265, 355, 602, 694, 911, 1234, 1432, 1927, 2680
TS2
750i, 196, 344, 358, 424, 481, 639, 685, 763, 1281, 2016, 3774
TS3
1111i, 190, 331, 398, 476, 644, 775, 1031, 1113, 1162, 1484, 2701
TS4
1432i, 358, 411, 469, 501, 566, 695, 880, 1094, 1125, 1589, 3785
TS5
500i, 70, 167, 213, 272, 476, 535, 752, 1208, 1328, 2684, 3107
TS6
918i, 378, 1162
C1
45, 52, 68, 127, 169, 218, 752, 1201, 1211, 1231, 2696, 2709
C2
227, 248, 287, 393, 538, 575, 904, 1040, 1221, 2665, 3766
C3
58, 124, 142, 196, 236, 419, 502, 826, 1093, 1633, 3754, 3876
H2SO–O2
32, 1014, 49, 65, 2616, 99, 1619, 129, 1190, 755, 3791, 472
HSOOOH
71, 264, 392, 435, 519, 585, 770, 968, 1028, 1427, 2593, 3684
O3
753 (716), 1213 (1089), 1264 (1135)
3O2
1646 (1556)
H2SO
487, 759, 1015, 1195, 2610, 3797
HSO
1013, 1091, 2426
HOO
1171, 1435, 3604Our calculation led to the identification of the reaction paths on singlet potential energy surfaces which can be summarized in the four following paths:P1(1): O3 + H2S→C1→TS1→C2→TS2→C3→SOO + H2O→TS6 + H2O→SO2 + H2OP1(2): O3 + H2S→C1→TS1→C2→TS2→C3→TS4→HSOOOH→HSO + HOOP2:   O3 + H2S→C1→TS3→HSOOOH→HSO + HOOP3:   O3 + H2S→C1→TS5→1O2 + HSOH→3O2 + HSOHwhere HSOOOH is a complex between the HOO and HSO molecules. The details of the reaction mechanism on the singlet potential energy surfaces will be discussed below.


### 2.1. Reaction Pathway Properties

We have only one collisional complex (C1) which was located on the potential energy surface for the H2S + O3 reaction which is formed when the sulfur atoms of hydrogen sulfide were attached to the terminal oxygen atoms of the ozone. For the H2S + O3 reaction, there are three possible products as SO2 + H2O, HSO + HOO, and HSOH + O2.


#### 2.1.1. SO2 + H2O Formation Channel

In path P1(1), C1 undergoes H–O1 and S–O2 bonds formation and S–O1, S–O3, and H–S bonds rupture process to form C2 via TS1 with the energy barriers of 178.41 kJ/mol. Then C2 transforms to C3 via TS2 with the energy barriers of 50.29 kJ/mol by O1–S and H–O3 bond formation and hydrogen transfers from sulfur atom to O1 atom.In the four-member ring structure transition region (TS2) of C2→C3 process, the electronic charge density and its Laplacian are 0.0356 and 0.1718 au, respectively.The energy height shows C1→C3 conversion is energetically feasible with two low energy transition states. The reaction path continues with the formation of water and SOO molecules from C3 without passing through any transition state. Finally, SOO undergoes O–S bond formation and O–O bond rupture process to the formation of SO2 molecule via TS6 with the energy barriers of 26.46 kJ/mol relative to original reactants and 198.08 kJ/mol relative to SOO molecule. The height energy of the final product (H2O + SO2) is 693.29 kJ/mol lower than the original reactants (H2S + O3). Because the formation of SO2 can be produced from the rearrangement of SOO and the number of their atoms is different from the others, we show this step of path P1(1) individually in Figure 3. The corresponding molecular graph of all species obtained from AIM calculations in path P1(1) in Figure 2 confirms that the formation and the rupture of the bonds are indicated by the appearance of BCP and the collapse of BCP, respectively. Also, the imaginary frequencies of TS1, TS2 and TS6, which denote the cleavage of some bonds and the formation of other bonds in the reaction coordinate, are 327i, 751i and 970i cm−1, respectively. In all the steps of path P1(1), the formation of products from reactants is confirmed using related IRC curves.


#### 2.1.2. HSO + HOO Formation Channels

For the HSO + HOO product, there are two pathways on the singlet potential energy surfaces denoted as P1(2) and P2 which are confirmed using vibrational frequencies of transition states, AIM analysis and IRC curves for transition states.In path P1(2), the formation of C3, a five-member ring structure with a low value of electron charge density, 0.0081 au, and Laplacian of charge density, 0.0366 au, are similar to path P1(1). C3 undergoes H–O1 bond rupture process and disappear of ring structure to form HSOOOH via TS4 with the energy barriers of 209.12 kJ/mol. Finally HSO + HOO molecules are formed by the breaking of O2–O3 bond without undergoing any transition state. The formation of products from C3 is confirmed using AIM results and IRC curve.In path P2, C1 transforms to the HSOOOH complex and its corresponding products (HSO and OOH) in the one step reaction via TS3 with the energy barriers of 92.70 kJ/mol.  Its transformation, accompanied with H2–O3 bond formation and S–O3 bond rupture process to form a five-membered ring transition state, TS3, molecule. Then, ring structure of TS3 disappears and form HSOOOH specie. Ring existing, H2–O3 bond, formation and S–O3 bond cleavage are confirmed using AIM calculation which are schematized in Figure2.


#### 2.1.3. HSOH + O2 Formation Channel

Another simple step of the reaction is path P3 which is an elementary reaction. In path P3, C1 undergoes to the products (HSOH + O2) using oxygen abstraction of ozone by hydrogen sulfide and at the same time hydrogen shift from sulfur atom to newly oxygen atom via TS5 with the energy barriers of 164.43 kJ/mol and vibrational imaginary frequency of 500i cm-1. The oxygen abstraction and hydrogen shift lead to HSOH formation which are confirmed using intrinsic reaction coordinate calculation (IRC).In summary, the H2S + O3 reaction is energetically feasible to the formation of all the three products throughout corresponding paths which all steps of reaction confirm using IRC and AIM calculations.


### 2.2. Topological Analysis of Electronic Density of the HSOH, HSO, and SOO Molecules

AIM topological analysis of electronic charge density has been elaborated using the AIM2000 program package. AIM topological analysis has been obtained using the electronic charge density integrated over atomic basins (up to 0.001 e/bohr3 level) as well as in terms of electron density, ρ(Γc), density Laplacian, ∇²ρ(Γc), and bond ellipticity, ε=λ1/λ2-1, at bond critical points (BCP) where λ1 and λ2 are the eigenvalues of the Hessian matrix of the BCP electron density.According to the topological analysis of electronic charge density, in the theory of the atoms in molecules (AIM) [14], electron density ρ(Γc) and Laplacian of the electron density, ∇²ρ(Γc), are used to describe the strength and the characteristic of the bond, respectively. The Laplacian (∇²ρ(Γc)) is the sum of λ1, λ2, and λ3, where λi is the ith eigenvalue of Hessian matrix of the electronic density. If a critical point has two negative and one positive eigenvalues, it is called (3,−1) or the bond critical point (BCP). If a critical point has two positive and one negative eigenvalues it is called (3,+1) or the ring critical point (RCP), which indicates that a ring structure exists. According to the theory of AIM, the Laplacian of the electron density (∇²ρ(Γc)) describes the characteristic of the bond. In general, when ∇²ρ(Γc) < 0 the bond is covalent, but when ∇²ρ(Γc) > 0 the bond belongs to the electrostatic interaction.The topological analysis of the electronic density were performed for all bonds in the HSOH, HSO, and SOO species as relatively stable products of our studied reaction. The topological characteristics in the BCP’s of S–O, O–H, and H–S bonds are tabulated in Table3. Also, the ∇²ρ(Γc) values in Table 3 indicate that all the bonds in HSOH have covalent character (∇²ρ(Γc) < 0). Their values show the bond strength manner O–H > H–S > S–O is corresponding to the electronic charge densities 0.3781, 0.2225, and 0.1929 au, respectively. In HSO radical, the H–S bond has a covalent character while the S–O bond electrostatic interaction nature with electronic charge densities value 0.2174 and 0.2686 au, respectively. The electronic charge density of HS bond of HSOH in comparison with corresponding bond in HSO bond shows that the HSOH bond is stronger than the corresponding bond of HSO, while S–O bond in HSO has electron charge density higher than the corresponding bond in HSOH molecule which is similar to sulfur and oxygen double bond. In SOO species, the O–O bond has a covalent character while S–O bond has an electrostatic character with the electronic charge densities values of 0.4369 and 0.2048 au, respectively. The S–O bond strange in SOO molecule is weaker than the corresponding bond in HSOH molecule and stronger than the corresponding bond in HSO species. Its nature is the same with HSO molecule.Table 3
Topological analysis of HSOH, HSO, and SOO molecules in atomic units.
HSOH
HSO
SOO
H–S bond
S–O bond
H–O bond
H–S bond
S–O bond
O–O bond
S–O bond
λ
1
−0.47454
−0.28127
−1.9606
−0.46772
−0.45831
−1.1149
−0.24991
λ
2
−0.39663
−0.22767
−1.9184
−0.42171
−0.4015
−1.7519
−0.24141
λ
3
0.18834
0.390444
0.91642
0.225109
1.36699
1.9251
0.77949
ρ
0.22245
0.192901
0.37809
0.217379
0.26855
0.43688
0.20482
∇
2
ρ
−0.68282
−0.1185
−2.9627
−0.66432
0.50718
−0.26501
0.28817
ε
0.19643
0.235402
0.021986
0.10908
0.14150
0.036940
0.035514The ellipticity is a measure of the ratio of the rate of density decrease in the two directions perpendicular to the bond path at the bond critical point, the general shape of the bond, and the degree ofπ-character. A value of zero indicates a symmetrical distribution of density about the bond path, such as that found in standard single and triple bonds, while large values indicate a preference for density buildup in a particular orientation. The values for the ellipticity of various bonds for HSOH molecule are listed in Table 3. The high ellipticity (ε) for S–O bond in both HSOH and HSO molecules shows its high π bonding character, whereas low value belongs to the O–H and H–S bonds which are corresponding with the low π bond contributions. In SOO molecule, both bonds have almost the same ellipticity (0.0369 and 0.0352 au). This manner shows both bonds have almost the same symmetry in the distribution of π bonding character and its distribution is almost symmetrical between the two bonds. The comparison of the two bonds shows that the O–O bond has a higher π-character compared with the S–O bond which is confirmed by the experimental predictions.


### 2.3. Molecular Orbital Analysis

Electron affinity, ionization potential, chemical potential, chemical hardness, and some other parameters for chemical reaction are predicted by the HOMO and the LUMO orbitals analysis. The HOMO (highest occupied molecular orbital) energy characterizes the ability of electron giving; the LUMO (lowest unoccupied molecular orbital) energy characterizes the ability of electron accepting [17]. We can determine the way the molecule interacts with other species; therefore, they are called the frontier orbitals. The energy gap between HOMO and LUMO characterizes that the molecules are chemically stable or not, and it is a critical parameter in determining molecular electrical transport properties because it is a measure of electron conductivity. In addition, the energy of the HOMO is directly related to the ionization potential, while the amount of electron affinity is related to the energy of the LUMO. The energy gaps are largely responsible for the chemical and spectroscopic properties of the molecules [18].The features of the HOMO and LUMO are calculated by B3LYP/6-311++G(3df, 3pd) level in gaseous phase. The energy band gap (Δε) (transition from HOMO to LUMO), chemical potential (μ), chemical hardness (η), ionization potential, and electron affinity of the HSOH, HSO, and SOO species can be seen in Table 4.Table 4
Hardness (η), chemical potential (μ), HOMO energy (EHOMO), and LUMO energy (ELUMO) (in atomic unit) of HSOH and HSO molecules.
E
HOMO
E
LUMO
μ
η
Δ
ε
HSOH
−0.2561
−0.0509
0.1535
−0.1026
0.2053
HSO
−0.2651
−0.02933
0.1472
−0.1179
0.2358
SOO
−0.28904
−0.1953
0.2422
−0.04685
0.0937The HOMO-LUMO energy gap in Table4 predicts the HSO (0.2358 au) with a large energy gap (Δε) are more stable than the HSOH (0.2052 au) and SOO (0.0937 au) molecules.The frontier molecular orbitals of HSO and HSOH are depicted in Figure4. In HSO molecule, the LUMO is localized on the bond between oxygen-sulfur atoms. The HOMO is mainly composed on the p-type atomic orbitals of oxygen and sulfur atoms. In SOO molecule, the HOMO orbital is localized on the sulfur atom while the LUMO is composed in the π type orbitals of all atoms.Figure 4
Frontier molecular orbitals of HSO, HSOH, and SOO molecules.The chemical potential (μ) and chemical hardness (η) [18] have been determined using HOMO and LUMO energies which are calculated using the relations μ=-(I+A)/2 and η=(I-A)/2, where I is the ionization potential (εHOMO) and A is the electron affinity (εLUMO) of the molecular system. According to the maximum hardness principal, the most stable anions should have maximum hardness value which is a minimum energy structure at a constant chemical potential [18]. The data in Table 4 show that the maximum absolute value of chemical hardness, −0.1179 au, belongs to HSO molecule. Also, the minimum chemical potential, 0.1472 au, of HSO molecule is in agreement with the predicted chemical potential.Determination of molecular charge distribution is fundamental to the application of the molecular modeling method to a realistic chemical system. The variation in the molecular charge distribution with structure can be applied to explain to stability of species. Atomic charges found with our calculations inside atoms are shown in Figure5. All atomic charges for HSOH, HSO, and SOO molecules are obtained by the Mulliken method at B3LYP/6-311++G(3df, 3pd) level. Results show that sulfur atoms have positive charge in three species and the oxygen atoms take negative charges in both HSOH and HSO molecules and, in SOO, the terminal oxygen atom has negative charge while the other atoms have positive charge as shown in Figure 5. But absolute value of charge distribution of oxygen atom in SOO is smaller than the HSO molecule and its oxygen negative charges smaller than the HSOH molecule that corresponds with chemical potential and HOMO-LUMO molecular orbital analysis.Figure 5
HSO and HSOH calculated partial charge distribution with the Mulliken method.Dipole moment shows the molecular charge distribution and is given as a vector in three dimensions. Therefore, it can be used as a descriptor to depict the charge movement across the molecule. The direction of the dipole moment vector in a molecule depends on the centers of positive and negative charges. In HSO species, dipole moment vector (2.298 D) is from sulfur atom to oxygen atom parallel with the S–O bond length but in HSOH molecule this vector (1.594 D) is to the outside of the molecule and maybe can cause the low stability in the molecule. For SOO molecule, this vector (1.525 D) is from terminal oxygen atom to sulfur atom.Rotational constants of the HSOH molecule are 202.1, 15.2, and 14.7 GHz and for HSO are 298.3, 20.3, and 19.0 GHz and for SOO species are 95.2, 6.9, and 6.4 GHz. Because the experimental value of the rotational constants are not accessible, so our calculated results are reported in an absolute form.


### 2.4. The Thermodynamic Data in the O3 and H2S Reaction Process

The change of thermodynamic characteristics for each reaction channel is the difference between the corresponding thermodynamic properties of the products and the reactants. The thermodynamic data are calculated by G4 level of computation in gaseous phase and are corrected by ZPE for O3 + H2S reaction.The calculated relative internal energies, enthalpies, Gibbs free energies, and entropies of all products of the reactions in gas-phase at atmospheric pressure and temperature of 298.15 K, have been summarized in Table5.Table 5
Relative energiesΔE°, enthalpies ΔH°, free energies ΔG° and entropies TΔS° (kJ mol−1) for reaction of ozone and hydrogen sulfide with ZPE correction at B3LYP level of computation.
Reaction
Δ
E
°
Δ
H
°
T
Δ
S
°
Δ
G
°
H2S + O3→ HSOH + 1O2
−117.09
−117.09
−1.25
115.85
H2S + O3→ HSOH + 3O2
−235.23
−235.23
1.47
−236.70
H2S + O3→ HSO + HOO
−135.15
−135.16
7.98
−143.13
H2S + O3→ SO2 + H2O
−655.63
−655.63 (−660.72)
−1.93 (2.29)
−653.70 (−718.4)Table5 data show that the ΔH of all overall product channels are negative that means the reaction pathways are exothermic. This value for the SO2 + H2O path is greater than the others with high excess of heat. Also, Gibbs free energy values of all production channels are negative so, all product channels are spontaneous in gas phase at atmospheric pressure and 298.15 K. The SO2 + H2O path product channel thermodynamically is the most feasible in comparison with the other products by −653.70 kJ/mol in the Gibbs free energy. Our calculated data are compared with the experimental value for the SO2 + H2O reaction channel and the obtained results are acceptable in comparison with the experiment (9 percent relative error in Gibbs free energy).


### 2.5. Calculation of the Rate Constants

Rate constant is a key parameter in kinetic study of reaction that can be calculated from statistical mechanics and its equation is obtained from transition state theory (TST) based on statistical mechanics [17, 18]. For bimolecular reactions, one traditionally uses transition state theory, which requires knowledge of an accurate transition structure so the activation energy is known. The following is a familiar equation for calculating the rate constant for a bimolecular reaction, originally derived from Eyring [19]:
(1)kbim=LχQTSQAQBkBThe-E0/kBT,
where χ is the tunneling correction; L is the Avogadro constant; E0 is the energy barrier between TS and reactants including zero point energy corrections; QTS, QA, and QB are the partition functions of the transition structure and reactant species, respectively. The Q parameters for one species are the product of translational, rotational, vibrational, and electronic partition functions, kB is Boltzmann’s constant, and h is Planck’s constant.The translational partition function per unit volume of a molecule of massm in three dimensions is
(2)Qtr=(2πmkBTh2)3/2.The external rotational partition function for a nonlinear molecule has the following form:(3)Qrot=8π2(8π3IaIbIc)1/2(kBT)3/2σh3,
here, the I parameters are the moments of inertia and σ is a symmetry number. As two reactants approach each other, three (for nonlinear molecules) or two (for linear molecules) external rotations would transform to hindered rotations or harmonic vibrations as they approach the transition state.Vibrations in the transition state are usually treated as harmonic oscillators. The vibrational partition function for one degree of freedom is calculated as(4)Qvib=11-e-hυ/kBT,
where, υ is vibrational frequency. According to classical mechanics, a particle can cross a potential energy barrier only if its energy is higher than the barrier height, while quantum mechanics predicts that there is a probability for the crossing at lower particle energies. This transmission probability is very sensitive to the particle energy and mass, the barrier height, and the shape of the barrier. The tunneling correction factor is defined as the quotient of the quantum mechanical rate to the classical rate. Although the tunneling process is a multidimensional phenomenon, for simplicity this is treated as a one-dimensional process. Also, one possible cause for a deviation from Arrhenius behavior is quantum mechanical tunneling of reactants through classical barrier which was assumed conservation of vibrational frequencies and moments of inertia of the fragments during the course of reaction. In this study, we used the simplest approach to assess the role of quantum tunneling which is the Wigner correction to the reaction rate [20]. It is related to the magnitude νim of the imaginary frequency along the reaction coordinates at the transition state. The rate constant is enhanced by a factor of
(5)χ(T)=1-(124)(hνimckBT)2,
where νim is imaginary frequency with negative value that accounts for the vibrational motion along the reaction path, and c is the speed of light.The rate constant for unimolecular reaction channel is calculated by RRKM (Rice-Ramsperger-Kassel-Marcus) theory. In RRKM theory, the thermal rate coefficient is obtained by integrating overE from E0 to infinity of the following equation:
(6)k=∫E0∞k2(E)dk1(E)[M]k-1[M]+k2(E),
where ki, E0, E, and M are defined based on the Lindemann mechanism for the first-order chemical reactions. The various terms of the rate expression are now evaluated using statistical mechanics. In the high pressure limit, RRKM theory reduces to transition state theory and the results of the two theories coincide.The rate constants of the reaction have been calculated on the singlet potential energy surfaces for two pathways which are defined as thermodynamic path (P1(1)) and kinetic path (P2) by the RRKM [20] (multichannel reaction) and TST [19] (one step reaction) theories implemented in the Ssumes [21] and Gpop [22] programs.For multichannel reactions (path P1(1)), total rate constant are evaluated in the high pressure limit of the Lindemann mechanism. In the high pressure limit, RRKM theory reduces to transition state theory and the results of the two theories coincide. High pressure limit of total rate constant is defined askikjkk/k-ik-j based on the Lindemann mechanism, where ki, kj, and kk are the rate expressions for the forward steps of the reaction that and k-i and k-j are the reverse steps of reaction follow unimolecular reactions. For the channels with one transition state, P2, total rate constant is calculated using the transition state theory for bimolecular reactions (Eyring equation (1)).On the singlet potential energy surfaces, four reaction pathways, P1(1), P1(2), P2, and P3, are obtained. The rate constants have been calculated for two paths (P1(1) as thermodynamic path and P2 as kinetic path). The first pathways, P1(1), has the combination of two steps, so, it is a complex reaction. One is the formation of SOO and H2O as products and, the next, SOO molecule rearranged to SO2 molecule as a final product after passing TS6. In the rate constant calculation the final step (SOO→SO2) was ignored because it is a unimolecular reaction. So, it is very faster than the formation of SOO steps. Also, the number of the atoms in the transition state is different from the number of the atoms in the other TSs which is a limitation for the software. The HSO and HOO adducts produce using two channels. The first channel for production of HSO and HOO molecules is passing TS3 and the second is passing TS1, TS2, and TS4. Between them the mechanism that includes TS3 is more reliable than the other three-step channel and the rate constant is calculated for the reliable path (P2).For the multichannel reaction of H2S and O3, the calculated rate constant for two paths are denoted as k1 and k2 which correspond with P1(1) and P2 paths, respectively. Table 6 shows the individual rate constant of two paths at the temperature range of 200–2500 K at B3LYP/6-311++G(3df, 3pd) level. The rate constant values for the paths in Table 6 reveal that the path P2 is kinetically the main path in all the temperature ranges. Therefore the production of HSO and HOO molecules is the main product. Because of low level energy of transition state from TS3 channel, production of HSO + HOO is the most important adducts (path P2 and corresponding k2) in comparison with the other paths. The value of rate constant shows in competition between P1(1) and P2 paths, path P2 is more reliable than the path P1(1). On the other word, the reaction of H2S + O3 kinetically proceeds HSO + HOO as the main products on the singlet potential energy surface through path P2 via the rate constants of k2. In thermodynamic aspect, the H2O + SO2 adducts are the main products on the singlet potential energy surface through path P1(1) via the rate constants of k1. Figure 6 shows the Arrhenius plots for the k1 and k2 rate constants.Table 6
The overall rate constant and branching ratio of the reaction (in cm3 molecule−1 s−1), calculated at the temperature range of 200–2500 K at the B3LYP/6-311++G(3df,3pd).
T
k
1
k
2
k
tot
k
1/ktot
k
2/ktot
200
9.99E − 38
9.17E − 24
9.17E − 24
1.09E − 12
1.00E + 02
300
4.66E − 29
3.63E − 21
3.63E − 21
1.28E − 06
1.00E + 02
400
9.33E − 25
1.79E − 19
1.79E − 19
5.21E − 04
1.00E + 02
500
3.49E − 22
2.40E − 18
2.40E − 18
1.45E − 02
1.00E + 02
600
1.83E − 20
1.56E − 17
1.56E − 17
1.18E − 01
9.99E + 01
700
3.15E − 19
6.47E − 17
6.51E − 17
4.84E − 01
9.95E + 01
800
2.71E − 18
2.02E − 16
2.04E − 16
1.33E + 00
9.87E + 01
900
1.48E − 17
5.14E − 16
5.28E − 16
2.80E + 00
9.72E + 01
1000
5.82E − 17
1.13E − 15
1.19E − 15
4.90E + 00
9.51E + 01
1400
2.22E − 15
1.08E − 14
1.30E − 14
1.71E + 01
8.29E + 01
1900
2.94E − 14
6.31E − 14
9.25E − 14
3.18E + 01
6.82E + 01
2500
1.95E − 13
2.59E − 13
4.54E − 13
4.29E + 01
5.71E + 01Figure 6
Arrhenius plot for thek1 (–∎–) and k2 (−▲−) rate constants of H2S + O3 reaction.Arrhenius parameters for thek1 and k2 rate constants are as follows:
(7)k1=1.219×10-11exp(-100.184kJmol-1RT),(8)k2=2.468×10-13exp(-43.191kJmol-1RT).The overall rate constant for of H2S + O3 reaction is denoted as k which is the sum of rate constant of two paths (k1 + k2). The branching ratios for the paths P1(1) and P2 are k1/k and k2/k, respectively. Table 6 shows the individual rate constant of paths, overall rate constants, and the branching ratios at the temperature range of 200–2500 K. The rate constant and branching ration values for several paths in Table 6 reveals that the path P2 is kinetically the main path in temperature lower than 1000 K. Therefore, the production of HSO and HOO molecules is the main adducts. In temperature higher than 1000 K, the production of SO2 and H2O molecules from TS1 and TS2 channel (path P1 (1) and corresponding k1 and k1/k) is the most important channels can kinetically compete with the HSO and HOO products. However, the production of SO2 and H2O molecules is thermodynamically important production. Figure 6 shows the Arrhenius plots for the k1 and k2 rate constants.As shown in Table6, the rate constant values of H2S + O3 under atmospheric pressure increase with increasing the temperature. Therefore, the main products of reaction depend on the temperatures and the reaction mechanism, so, this reaction can play different roles in the formation of acid rain.


## 2.1. Reaction Pathway Properties

We have only one collisional complex (C1) which was located on the potential energy surface for the H2S + O3 reaction which is formed when the sulfur atoms of hydrogen sulfide were attached to the terminal oxygen atoms of the ozone. For the H2S + O3 reaction, there are three possible products as SO2 + H2O, HSO + HOO, and HSOH + O2.


### 2.1.1. SO2 + H2O Formation Channel

In path P1(1), C1 undergoes H–O1 and S–O2 bonds formation and S–O1, S–O3, and H–S bonds rupture process to form C2 via TS1 with the energy barriers of 178.41 kJ/mol. Then C2 transforms to C3 via TS2 with the energy barriers of 50.29 kJ/mol by O1–S and H–O3 bond formation and hydrogen transfers from sulfur atom to O1 atom.In the four-member ring structure transition region (TS2) of C2→C3 process, the electronic charge density and its Laplacian are 0.0356 and 0.1718 au, respectively.The energy height shows C1→C3 conversion is energetically feasible with two low energy transition states. The reaction path continues with the formation of water and SOO molecules from C3 without passing through any transition state. Finally, SOO undergoes O–S bond formation and O–O bond rupture process to the formation of SO2 molecule via TS6 with the energy barriers of 26.46 kJ/mol relative to original reactants and 198.08 kJ/mol relative to SOO molecule. The height energy of the final product (H2O + SO2) is 693.29 kJ/mol lower than the original reactants (H2S + O3). Because the formation of SO2 can be produced from the rearrangement of SOO and the number of their atoms is different from the others, we show this step of path P1(1) individually in Figure 3. The corresponding molecular graph of all species obtained from AIM calculations in path P1(1) in Figure 2 confirms that the formation and the rupture of the bonds are indicated by the appearance of BCP and the collapse of BCP, respectively. Also, the imaginary frequencies of TS1, TS2 and TS6, which denote the cleavage of some bonds and the formation of other bonds in the reaction coordinate, are 327i, 751i and 970i cm−1, respectively. In all the steps of path P1(1), the formation of products from reactants is confirmed using related IRC curves.


### 2.1.2. HSO + HOO Formation Channels

For the HSO + HOO product, there are two pathways on the singlet potential energy surfaces denoted as P1(2) and P2 which are confirmed using vibrational frequencies of transition states, AIM analysis and IRC curves for transition states.In path P1(2), the formation of C3, a five-member ring structure with a low value of electron charge density, 0.0081 au, and Laplacian of charge density, 0.0366 au, are similar to path P1(1). C3 undergoes H–O1 bond rupture process and disappear of ring structure to form HSOOOH via TS4 with the energy barriers of 209.12 kJ/mol. Finally HSO + HOO molecules are formed by the breaking of O2–O3 bond without undergoing any transition state. The formation of products from C3 is confirmed using AIM results and IRC curve.In path P2, C1 transforms to the HSOOOH complex and its corresponding products (HSO and OOH) in the one step reaction via TS3 with the energy barriers of 92.70 kJ/mol.  Its transformation, accompanied with H2–O3 bond formation and S–O3 bond rupture process to form a five-membered ring transition state, TS3, molecule. Then, ring structure of TS3 disappears and form HSOOOH specie. Ring existing, H2–O3 bond, formation and S–O3 bond cleavage are confirmed using AIM calculation which are schematized in Figure2.


### 2.1.3. HSOH + O2 Formation Channel

Another simple step of the reaction is path P3 which is an elementary reaction. In path P3, C1 undergoes to the products (HSOH + O2) using oxygen abstraction of ozone by hydrogen sulfide and at the same time hydrogen shift from sulfur atom to newly oxygen atom via TS5 with the energy barriers of 164.43 kJ/mol and vibrational imaginary frequency of 500i cm-1. The oxygen abstraction and hydrogen shift lead to HSOH formation which are confirmed using intrinsic reaction coordinate calculation (IRC).In summary, the H2S + O3 reaction is energetically feasible to the formation of all the three products throughout corresponding paths which all steps of reaction confirm using IRC and AIM calculations.


## 2.1.1. SO2 + H2O Formation Channel

In path P1(1), C1 undergoes H–O1 and S–O2 bonds formation and S–O1, S–O3, and H–S bonds rupture process to form C2 via TS1 with the energy barriers of 178.41 kJ/mol. Then C2 transforms to C3 via TS2 with the energy barriers of 50.29 kJ/mol by O1–S and H–O3 bond formation and hydrogen transfers from sulfur atom to O1 atom.In the four-member ring structure transition region (TS2) of C2→C3 process, the electronic charge density and its Laplacian are 0.0356 and 0.1718 au, respectively.The energy height shows C1→C3 conversion is energetically feasible with two low energy transition states. The reaction path continues with the formation of water and SOO molecules from C3 without passing through any transition state. Finally, SOO undergoes O–S bond formation and O–O bond rupture process to the formation of SO2 molecule via TS6 with the energy barriers of 26.46 kJ/mol relative to original reactants and 198.08 kJ/mol relative to SOO molecule. The height energy of the final product (H2O + SO2) is 693.29 kJ/mol lower than the original reactants (H2S + O3). Because the formation of SO2 can be produced from the rearrangement of SOO and the number of their atoms is different from the others, we show this step of path P1(1) individually in Figure 3. The corresponding molecular graph of all species obtained from AIM calculations in path P1(1) in Figure 2 confirms that the formation and the rupture of the bonds are indicated by the appearance of BCP and the collapse of BCP, respectively. Also, the imaginary frequencies of TS1, TS2 and TS6, which denote the cleavage of some bonds and the formation of other bonds in the reaction coordinate, are 327i, 751i and 970i cm−1, respectively. In all the steps of path P1(1), the formation of products from reactants is confirmed using related IRC curves.


## 2.1.2. HSO + HOO Formation Channels

For the HSO + HOO product, there are two pathways on the singlet potential energy surfaces denoted as P1(2) and P2 which are confirmed using vibrational frequencies of transition states, AIM analysis and IRC curves for transition states.In path P1(2), the formation of C3, a five-member ring structure with a low value of electron charge density, 0.0081 au, and Laplacian of charge density, 0.0366 au, are similar to path P1(1). C3 undergoes H–O1 bond rupture process and disappear of ring structure to form HSOOOH via TS4 with the energy barriers of 209.12 kJ/mol. Finally HSO + HOO molecules are formed by the breaking of O2–O3 bond without undergoing any transition state. The formation of products from C3 is confirmed using AIM results and IRC curve.In path P2, C1 transforms to the HSOOOH complex and its corresponding products (HSO and OOH) in the one step reaction via TS3 with the energy barriers of 92.70 kJ/mol.  Its transformation, accompanied with H2–O3 bond formation and S–O3 bond rupture process to form a five-membered ring transition state, TS3, molecule. Then, ring structure of TS3 disappears and form HSOOOH specie. Ring existing, H2–O3 bond, formation and S–O3 bond cleavage are confirmed using AIM calculation which are schematized in Figure2.


## 2.1.3. HSOH + O2 Formation Channel

Another simple step of the reaction is path P3 which is an elementary reaction. In path P3, C1 undergoes to the products (HSOH + O2) using oxygen abstraction of ozone by hydrogen sulfide and at the same time hydrogen shift from sulfur atom to newly oxygen atom via TS5 with the energy barriers of 164.43 kJ/mol and vibrational imaginary frequency of 500i cm-1. The oxygen abstraction and hydrogen shift lead to HSOH formation which are confirmed using intrinsic reaction coordinate calculation (IRC).In summary, the H2S + O3 reaction is energetically feasible to the formation of all the three products throughout corresponding paths which all steps of reaction confirm using IRC and AIM calculations.


## 2.2. Topological Analysis of Electronic Density of the HSOH, HSO, and SOO Molecules

AIM topological analysis of electronic charge density has been elaborated using the AIM2000 program package. AIM topological analysis has been obtained using the electronic charge density integrated over atomic basins (up to 0.001 e/bohr3 level) as well as in terms of electron density, ρ(Γc), density Laplacian, ∇²ρ(Γc), and bond ellipticity, ε=λ1/λ2-1, at bond critical points (BCP) where λ1 and λ2 are the eigenvalues of the Hessian matrix of the BCP electron density.According to the topological analysis of electronic charge density, in the theory of the atoms in molecules (AIM) [14], electron density ρ(Γc) and Laplacian of the electron density, ∇²ρ(Γc), are used to describe the strength and the characteristic of the bond, respectively. The Laplacian (∇²ρ(Γc)) is the sum of λ1, λ2, and λ3, where λi is the ith eigenvalue of Hessian matrix of the electronic density. If a critical point has two negative and one positive eigenvalues, it is called (3,−1) or the bond critical point (BCP). If a critical point has two positive and one negative eigenvalues it is called (3,+1) or the ring critical point (RCP), which indicates that a ring structure exists. According to the theory of AIM, the Laplacian of the electron density (∇²ρ(Γc)) describes the characteristic of the bond. In general, when ∇²ρ(Γc) < 0 the bond is covalent, but when ∇²ρ(Γc) > 0 the bond belongs to the electrostatic interaction.The topological analysis of the electronic density were performed for all bonds in the HSOH, HSO, and SOO species as relatively stable products of our studied reaction. The topological characteristics in the BCP’s of S–O, O–H, and H–S bonds are tabulated in Table3. Also, the ∇²ρ(Γc) values in Table 3 indicate that all the bonds in HSOH have covalent character (∇²ρ(Γc) < 0). Their values show the bond strength manner O–H > H–S > S–O is corresponding to the electronic charge densities 0.3781, 0.2225, and 0.1929 au, respectively. In HSO radical, the H–S bond has a covalent character while the S–O bond electrostatic interaction nature with electronic charge densities value 0.2174 and 0.2686 au, respectively. The electronic charge density of HS bond of HSOH in comparison with corresponding bond in HSO bond shows that the HSOH bond is stronger than the corresponding bond of HSO, while S–O bond in HSO has electron charge density higher than the corresponding bond in HSOH molecule which is similar to sulfur and oxygen double bond. In SOO species, the O–O bond has a covalent character while S–O bond has an electrostatic character with the electronic charge densities values of 0.4369 and 0.2048 au, respectively. The S–O bond strange in SOO molecule is weaker than the corresponding bond in HSOH molecule and stronger than the corresponding bond in HSO species. Its nature is the same with HSO molecule.Table 3
Topological analysis of HSOH, HSO, and SOO molecules in atomic units.
HSOH
HSO
SOO
H–S bond
S–O bond
H–O bond
H–S bond
S–O bond
O–O bond
S–O bond
λ
1
−0.47454
−0.28127
−1.9606
−0.46772
−0.45831
−1.1149
−0.24991
λ
2
−0.39663
−0.22767
−1.9184
−0.42171
−0.4015
−1.7519
−0.24141
λ
3
0.18834
0.390444
0.91642
0.225109
1.36699
1.9251
0.77949
ρ
0.22245
0.192901
0.37809
0.217379
0.26855
0.43688
0.20482
∇
2
ρ
−0.68282
−0.1185
−2.9627
−0.66432
0.50718
−0.26501
0.28817
ε
0.19643
0.235402
0.021986
0.10908
0.14150
0.036940
0.035514The ellipticity is a measure of the ratio of the rate of density decrease in the two directions perpendicular to the bond path at the bond critical point, the general shape of the bond, and the degree ofπ-character. A value of zero indicates a symmetrical distribution of density about the bond path, such as that found in standard single and triple bonds, while large values indicate a preference for density buildup in a particular orientation. The values for the ellipticity of various bonds for HSOH molecule are listed in Table 3. The high ellipticity (ε) for S–O bond in both HSOH and HSO molecules shows its high π bonding character, whereas low value belongs to the O–H and H–S bonds which are corresponding with the low π bond contributions. In SOO molecule, both bonds have almost the same ellipticity (0.0369 and 0.0352 au). This manner shows both bonds have almost the same symmetry in the distribution of π bonding character and its distribution is almost symmetrical between the two bonds. The comparison of the two bonds shows that the O–O bond has a higher π-character compared with the S–O bond which is confirmed by the experimental predictions.


## 2.3. Molecular Orbital Analysis

Electron affinity, ionization potential, chemical potential, chemical hardness, and some other parameters for chemical reaction are predicted by the HOMO and the LUMO orbitals analysis. The HOMO (highest occupied molecular orbital) energy characterizes the ability of electron giving; the LUMO (lowest unoccupied molecular orbital) energy characterizes the ability of electron accepting [17]. We can determine the way the molecule interacts with other species; therefore, they are called the frontier orbitals. The energy gap between HOMO and LUMO characterizes that the molecules are chemically stable or not, and it is a critical parameter in determining molecular electrical transport properties because it is a measure of electron conductivity. In addition, the energy of the HOMO is directly related to the ionization potential, while the amount of electron affinity is related to the energy of the LUMO. The energy gaps are largely responsible for the chemical and spectroscopic properties of the molecules [18].The features of the HOMO and LUMO are calculated by B3LYP/6-311++G(3df, 3pd) level in gaseous phase. The energy band gap (Δε) (transition from HOMO to LUMO), chemical potential (μ), chemical hardness (η), ionization potential, and electron affinity of the HSOH, HSO, and SOO species can be seen in Table 4.Table 4
Hardness (η), chemical potential (μ), HOMO energy (EHOMO), and LUMO energy (ELUMO) (in atomic unit) of HSOH and HSO molecules.
E
HOMO
E
LUMO
μ
η
Δ
ε
HSOH
−0.2561
−0.0509
0.1535
−0.1026
0.2053
HSO
−0.2651
−0.02933
0.1472
−0.1179
0.2358
SOO
−0.28904
−0.1953
0.2422
−0.04685
0.0937The HOMO-LUMO energy gap in Table4 predicts the HSO (0.2358 au) with a large energy gap (Δε) are more stable than the HSOH (0.2052 au) and SOO (0.0937 au) molecules.The frontier molecular orbitals of HSO and HSOH are depicted in Figure4. In HSO molecule, the LUMO is localized on the bond between oxygen-sulfur atoms. The HOMO is mainly composed on the p-type atomic orbitals of oxygen and sulfur atoms. In SOO molecule, the HOMO orbital is localized on the sulfur atom while the LUMO is composed in the π type orbitals of all atoms.Figure 4
Frontier molecular orbitals of HSO, HSOH, and SOO molecules.The chemical potential (μ) and chemical hardness (η) [18] have been determined using HOMO and LUMO energies which are calculated using the relations μ=-(I+A)/2 and η=(I-A)/2, where I is the ionization potential (εHOMO) and A is the electron affinity (εLUMO) of the molecular system. According to the maximum hardness principal, the most stable anions should have maximum hardness value which is a minimum energy structure at a constant chemical potential [18]. The data in Table 4 show that the maximum absolute value of chemical hardness, −0.1179 au, belongs to HSO molecule. Also, the minimum chemical potential, 0.1472 au, of HSO molecule is in agreement with the predicted chemical potential.Determination of molecular charge distribution is fundamental to the application of the molecular modeling method to a realistic chemical system. The variation in the molecular charge distribution with structure can be applied to explain to stability of species. Atomic charges found with our calculations inside atoms are shown in Figure5. All atomic charges for HSOH, HSO, and SOO molecules are obtained by the Mulliken method at B3LYP/6-311++G(3df, 3pd) level. Results show that sulfur atoms have positive charge in three species and the oxygen atoms take negative charges in both HSOH and HSO molecules and, in SOO, the terminal oxygen atom has negative charge while the other atoms have positive charge as shown in Figure 5. But absolute value of charge distribution of oxygen atom in SOO is smaller than the HSO molecule and its oxygen negative charges smaller than the HSOH molecule that corresponds with chemical potential and HOMO-LUMO molecular orbital analysis.Figure 5
HSO and HSOH calculated partial charge distribution with the Mulliken method.Dipole moment shows the molecular charge distribution and is given as a vector in three dimensions. Therefore, it can be used as a descriptor to depict the charge movement across the molecule. The direction of the dipole moment vector in a molecule depends on the centers of positive and negative charges. In HSO species, dipole moment vector (2.298 D) is from sulfur atom to oxygen atom parallel with the S–O bond length but in HSOH molecule this vector (1.594 D) is to the outside of the molecule and maybe can cause the low stability in the molecule. For SOO molecule, this vector (1.525 D) is from terminal oxygen atom to sulfur atom.Rotational constants of the HSOH molecule are 202.1, 15.2, and 14.7 GHz and for HSO are 298.3, 20.3, and 19.0 GHz and for SOO species are 95.2, 6.9, and 6.4 GHz. Because the experimental value of the rotational constants are not accessible, so our calculated results are reported in an absolute form.


## 2.4. The Thermodynamic Data in the O3 and H2S Reaction Process

The change of thermodynamic characteristics for each reaction channel is the difference between the corresponding thermodynamic properties of the products and the reactants. The thermodynamic data are calculated by G4 level of computation in gaseous phase and are corrected by ZPE for O3 + H2S reaction.The calculated relative internal energies, enthalpies, Gibbs free energies, and entropies of all products of the reactions in gas-phase at atmospheric pressure and temperature of 298.15 K, have been summarized in Table5.Table 5
Relative energiesΔE°, enthalpies ΔH°, free energies ΔG° and entropies TΔS° (kJ mol−1) for reaction of ozone and hydrogen sulfide with ZPE correction at B3LYP level of computation.
Reaction
Δ
E
°
Δ
H
°
T
Δ
S
°
Δ
G
°
H2S + O3→ HSOH + 1O2
−117.09
−117.09
−1.25
115.85
H2S + O3→ HSOH + 3O2
−235.23
−235.23
1.47
−236.70
H2S + O3→ HSO + HOO
−135.15
−135.16
7.98
−143.13
H2S + O3→ SO2 + H2O
−655.63
−655.63 (−660.72)
−1.93 (2.29)
−653.70 (−718.4)Table5 data show that the ΔH of all overall product channels are negative that means the reaction pathways are exothermic. This value for the SO2 + H2O path is greater than the others with high excess of heat. Also, Gibbs free energy values of all production channels are negative so, all product channels are spontaneous in gas phase at atmospheric pressure and 298.15 K. The SO2 + H2O path product channel thermodynamically is the most feasible in comparison with the other products by −653.70 kJ/mol in the Gibbs free energy. Our calculated data are compared with the experimental value for the SO2 + H2O reaction channel and the obtained results are acceptable in comparison with the experiment (9 percent relative error in Gibbs free energy).


## 2.5. Calculation of the Rate Constants

Rate constant is a key parameter in kinetic study of reaction that can be calculated from statistical mechanics and its equation is obtained from transition state theory (TST) based on statistical mechanics [17, 18]. For bimolecular reactions, one traditionally uses transition state theory, which requires knowledge of an accurate transition structure so the activation energy is known. The following is a familiar equation for calculating the rate constant for a bimolecular reaction, originally derived from Eyring [19]:
(1)kbim=LχQTSQAQBkBThe-E0/kBT,
where χ is the tunneling correction; L is the Avogadro constant; E0 is the energy barrier between TS and reactants including zero point energy corrections; QTS, QA, and QB are the partition functions of the transition structure and reactant species, respectively. The Q parameters for one species are the product of translational, rotational, vibrational, and electronic partition functions, kB is Boltzmann’s constant, and h is Planck’s constant.The translational partition function per unit volume of a molecule of massm in three dimensions is
(2)Qtr=(2πmkBTh2)3/2.The external rotational partition function for a nonlinear molecule has the following form:(3)Qrot=8π2(8π3IaIbIc)1/2(kBT)3/2σh3,
here, the I parameters are the moments of inertia and σ is a symmetry number. As two reactants approach each other, three (for nonlinear molecules) or two (for linear molecules) external rotations would transform to hindered rotations or harmonic vibrations as they approach the transition state.Vibrations in the transition state are usually treated as harmonic oscillators. The vibrational partition function for one degree of freedom is calculated as(4)Qvib=11-e-hυ/kBT,
where, υ is vibrational frequency. According to classical mechanics, a particle can cross a potential energy barrier only if its energy is higher than the barrier height, while quantum mechanics predicts that there is a probability for the crossing at lower particle energies. This transmission probability is very sensitive to the particle energy and mass, the barrier height, and the shape of the barrier. The tunneling correction factor is defined as the quotient of the quantum mechanical rate to the classical rate. Although the tunneling process is a multidimensional phenomenon, for simplicity this is treated as a one-dimensional process. Also, one possible cause for a deviation from Arrhenius behavior is quantum mechanical tunneling of reactants through classical barrier which was assumed conservation of vibrational frequencies and moments of inertia of the fragments during the course of reaction. In this study, we used the simplest approach to assess the role of quantum tunneling which is the Wigner correction to the reaction rate [20]. It is related to the magnitude νim of the imaginary frequency along the reaction coordinates at the transition state. The rate constant is enhanced by a factor of
(5)χ(T)=1-(124)(hνimckBT)2,
where νim is imaginary frequency with negative value that accounts for the vibrational motion along the reaction path, and c is the speed of light.The rate constant for unimolecular reaction channel is calculated by RRKM (Rice-Ramsperger-Kassel-Marcus) theory. In RRKM theory, the thermal rate coefficient is obtained by integrating overE from E0 to infinity of the following equation:
(6)k=∫E0∞k2(E)dk1(E)[M]k-1[M]+k2(E),
where ki, E0, E, and M are defined based on the Lindemann mechanism for the first-order chemical reactions. The various terms of the rate expression are now evaluated using statistical mechanics. In the high pressure limit, RRKM theory reduces to transition state theory and the results of the two theories coincide.The rate constants of the reaction have been calculated on the singlet potential energy surfaces for two pathways which are defined as thermodynamic path (P1(1)) and kinetic path (P2) by the RRKM [20] (multichannel reaction) and TST [19] (one step reaction) theories implemented in the Ssumes [21] and Gpop [22] programs.For multichannel reactions (path P1(1)), total rate constant are evaluated in the high pressure limit of the Lindemann mechanism. In the high pressure limit, RRKM theory reduces to transition state theory and the results of the two theories coincide. High pressure limit of total rate constant is defined askikjkk/k-ik-j based on the Lindemann mechanism, where ki, kj, and kk are the rate expressions for the forward steps of the reaction that and k-i and k-j are the reverse steps of reaction follow unimolecular reactions. For the channels with one transition state, P2, total rate constant is calculated using the transition state theory for bimolecular reactions (Eyring equation (1)).On the singlet potential energy surfaces, four reaction pathways, P1(1), P1(2), P2, and P3, are obtained. The rate constants have been calculated for two paths (P1(1) as thermodynamic path and P2 as kinetic path). The first pathways, P1(1), has the combination of two steps, so, it is a complex reaction. One is the formation of SOO and H2O as products and, the next, SOO molecule rearranged to SO2 molecule as a final product after passing TS6. In the rate constant calculation the final step (SOO→SO2) was ignored because it is a unimolecular reaction. So, it is very faster than the formation of SOO steps. Also, the number of the atoms in the transition state is different from the number of the atoms in the other TSs which is a limitation for the software. The HSO and HOO adducts produce using two channels. The first channel for production of HSO and HOO molecules is passing TS3 and the second is passing TS1, TS2, and TS4. Between them the mechanism that includes TS3 is more reliable than the other three-step channel and the rate constant is calculated for the reliable path (P2).For the multichannel reaction of H2S and O3, the calculated rate constant for two paths are denoted as k1 and k2 which correspond with P1(1) and P2 paths, respectively. Table 6 shows the individual rate constant of two paths at the temperature range of 200–2500 K at B3LYP/6-311++G(3df, 3pd) level. The rate constant values for the paths in Table 6 reveal that the path P2 is kinetically the main path in all the temperature ranges. Therefore the production of HSO and HOO molecules is the main product. Because of low level energy of transition state from TS3 channel, production of HSO + HOO is the most important adducts (path P2 and corresponding k2) in comparison with the other paths. The value of rate constant shows in competition between P1(1) and P2 paths, path P2 is more reliable than the path P1(1). On the other word, the reaction of H2S + O3 kinetically proceeds HSO + HOO as the main products on the singlet potential energy surface through path P2 via the rate constants of k2. In thermodynamic aspect, the H2O + SO2 adducts are the main products on the singlet potential energy surface through path P1(1) via the rate constants of k1. Figure 6 shows the Arrhenius plots for the k1 and k2 rate constants.Table 6
The overall rate constant and branching ratio of the reaction (in cm3 molecule−1 s−1), calculated at the temperature range of 200–2500 K at the B3LYP/6-311++G(3df,3pd).
T
k
1
k
2
k
tot
k
1/ktot
k
2/ktot
200
9.99E − 38
9.17E − 24
9.17E − 24
1.09E − 12
1.00E + 02
300
4.66E − 29
3.63E − 21
3.63E − 21
1.28E − 06
1.00E + 02
400
9.33E − 25
1.79E − 19
1.79E − 19
5.21E − 04
1.00E + 02
500
3.49E − 22
2.40E − 18
2.40E − 18
1.45E − 02
1.00E + 02
600
1.83E − 20
1.56E − 17
1.56E − 17
1.18E − 01
9.99E + 01
700
3.15E − 19
6.47E − 17
6.51E − 17
4.84E − 01
9.95E + 01
800
2.71E − 18
2.02E − 16
2.04E − 16
1.33E + 00
9.87E + 01
900
1.48E − 17
5.14E − 16
5.28E − 16
2.80E + 00
9.72E + 01
1000
5.82E − 17
1.13E − 15
1.19E − 15
4.90E + 00
9.51E + 01
1400
2.22E − 15
1.08E − 14
1.30E − 14
1.71E + 01
8.29E + 01
1900
2.94E − 14
6.31E − 14
9.25E − 14
3.18E + 01
6.82E + 01
2500
1.95E − 13
2.59E − 13
4.54E − 13
4.29E + 01
5.71E + 01Figure 6
Arrhenius plot for thek1 (–∎–) and k2 (−▲−) rate constants of H2S + O3 reaction.Arrhenius parameters for thek1 and k2 rate constants are as follows:
(7)k1=1.219×10-11exp(-100.184kJmol-1RT),(8)k2=2.468×10-13exp(-43.191kJmol-1RT).The overall rate constant for of H2S + O3 reaction is denoted as k which is the sum of rate constant of two paths (k1 + k2). The branching ratios for the paths P1(1) and P2 are k1/k and k2/k, respectively. Table 6 shows the individual rate constant of paths, overall rate constants, and the branching ratios at the temperature range of 200–2500 K. The rate constant and branching ration values for several paths in Table 6 reveals that the path P2 is kinetically the main path in temperature lower than 1000 K. Therefore, the production of HSO and HOO molecules is the main adducts. In temperature higher than 1000 K, the production of SO2 and H2O molecules from TS1 and TS2 channel (path P1 (1) and corresponding k1 and k1/k) is the most important channels can kinetically compete with the HSO and HOO products. However, the production of SO2 and H2O molecules is thermodynamically important production. Figure 6 shows the Arrhenius plots for the k1 and k2 rate constants.As shown in Table6, the rate constant values of H2S + O3 under atmospheric pressure increase with increasing the temperature. Therefore, the main products of reaction depend on the temperatures and the reaction mechanism, so, this reaction can play different roles in the formation of acid rain.


## 3. Conclusion

The reaction of H2S and O3 has been investigated at the B3LYP/6-311++G(3df,3pd), CCSD/6-311++G(3df,3pd)//B3LYP/6-311++G(3df,3pd), and G4 levels of computation. By dissociation of H2S-O3 complex, four product channels have been found, two channels for production of HSO and HOO products, one channel for production of SO2 and H2O, and another one for the production of HSOH and O2 molecules via corresponding transition states and complexes. AIM analyses were used to elucidate the reaction mechanism by comparing ρ(r) of BCP in complexes and transition states. Furthermore, AIM indicates ring structure for some transition structures and complexes. Additionally intrinsic reaction coordinate (IRC) calculation was performed to examine the correspondence of the calculated TS to the reactants and products. The energy profile of reaction mechanism is obtained. The thermodynamic results show that the SO2 + H2O product is the main product of H2S + O3 reaction with the most negative value of ΔG. Also the rate constant and branching ratio are calculated and discussed. Results show two channels are competing with each other based on the temperature. The HSO + HOO adducts are the main kinetically products while the SO2 + H2O molecules are the main products thermodynamically.



---

*Source: 659682-2012-12-30.xml*