Electronic structure, stability, and cooperativity of chalcogen bonding in sulfur dioxide and hydrated sulfur dioxide clusters: a DFT study and wave functional analysis

Density functional theory calculations and wave functional analysis are used to examine the (SO2)n and (SO2)n–H2O clusters with n = 1–7. The nature of interactions is explored by molecular electrostatic potentials, electron density distribution, atoms in molecules, noncovalent interaction, and energy decomposition analysis. The putative global minimum of SO2 molecules has a 3D growth pattern with tetrahedral. In the hydrated SO2 clusters, the pure hydrogen bond isomers are less stable than the O–S chalcogen bond isomers. The cluster absorption energy of SO2 on water increases with the size of sulfur dioxide, implying reactivity of sulfur dioxide with water increases with size. The presence of cooperativity was evident from the excellent linearity plot of binding energy/polarizability vs the number of SO2 molecules. Molecular electrostatic potential analysis elucidates the reason for the facile formation of S–O chalcogen than hydrogen bonding in hydrated sulfur dioxide. Atoms in molecule analysis characterize the bonds chalcogen and H bonds to be weak and electrostatic dominant. EDA analysis shows electrostatic interaction is dominated in complexes with more intermolecular chalcogen bonding and orbital interaction for systems with intermolecular H-bonding. The reduced density gradient (RDG) analysis of sulfur dioxide clusters has blue patches and green patches due to S–O chalcogen bonding O–O electrostatic interaction. The RDG analysis of hydrated sulfur dioxide clusters shows intensive blue patches and green patches for the existence of S–O chalcogen and hydrogen bonding respectively. Thus, the presence of strong electrostatic S–O chalcogen interaction and weak H bonds acts cooperatively and stabilize the hydrated sulfur dioxide clusters.


Introduction
Sulfur dioxide (SO 2 ) is a primary gas molecule responsible for air pollutions. Its source of accumulation in the earth's atmosphere is being both anthropogenic and natural including volcanic eruptions [1,2]. Recently, IPCC assessment report has publicized that short-term radiative SO 2 exposure to the atmosphere is much more hazardous than CO 2 exposure [3]. The reason for this is linked high polarity of both SO 2 and water which leads to the substantial solubility of SO 2 leading to the formation of sulfuric acid by an oxidation process. Because of its prominence in atmospheric implication, the interaction between water and sulfur dioxide has been the subject of many experimental and theoretical studies [4][5][6]. It is generally accepted that during solvation, the SO 2 molecules to be existent in hydrated sulfur dioxide adducts [7].
Several studies were carried out to elucidate the structure and electronic properties of SO 2 monomers and dimers. The geometry of the SO 2 molecule determined its electron affinity, bond length, and angle, and vibrational frequency using photoelectron spectroscopy in solid Argon and Neon matrices was reported [8]. Very recently, using ab initio methods, the structure and electronic properties of neutral and anionic ground states of SO 2 molecule were reported [9]. Taleb-Bendiab et al., using microwave spectroscopy, reported the structure of SO2 dimer to have a low symmetry (C s isomer) with the first SO 2 unit defines a plane and the second to lie above the first [10]. The infrared spectroscopic studies on Ar matrices show a larger peak shift supporting the C s isomer [11]. On the contrary, in the N 2 matrix, a small shift was noticed supporting the C i isomer. Very recently, Ito and Hirabayashi reported SO 2 dimer to subsist in Cs isomer in Kr and Xe matrices in a matrix-isolation study [12]. Bernstein et al. investigated the (SO 2 ) n clusters by using single-photon ionization [13]. The distribution of SO 2 clusters decreases roughly exponentially with increasing cluster size.
To understand the atmospheric implication of SO 2 , the interaction of SO 2 with water molecules was studied in gas or aqueous aerosols by both experimental and theoretical methods. Infrared and Raman studies in both gas and solution phase support the existence of H 2 SO 3 to be a loosely aquated molecule [14,15]. To get molecular insight into the microscopic structure of SO 2 /H 2 O complexes, several simulations were carried out; however, there persist considerable ambiguity in the structure. Born-Oppenheimer molecular dynamics simulations to investigate the hydrated structures of SO 2 show that in the gas phase and at a temperature 300 K, the main interaction was electrostatic S-O interaction between SO 2 and water [16]. However, at the surface, hydrogen-bonding interaction was predominant. Microwave spectroscopy analysis suggests that H 2 O and SO 2 form a double-decker orientation, in which the oxygen atom of water is close to the sulfur atom of the SO 2 molecule [17]. Infrared spectra of (H 2 O)n-SO 2 complexes trapped in argon matrices show hydrogen bonds for 2:1 and 3:1 complexes [18].
The calculations on the structure and energetics of H 2 O/ SO 2 complexes were studied at the MP2 level by Cukras and Sadlej [19]. The global minima of the dimer have a double-decker orientation and the presence of strong electrostatic, dispersion, and induction interaction stabilizes the complex. The structure of sulfur dioxide with up to 8 water molecules was investigated by Steudel and Steudel at B3LYP/6-31G(2df,p) level of theory [20]. SO 2 resides at the surface of the water cluster, and the sulfur atom does not take part in the hydrogen-bonding network until the water molecules are less than six. The further increase in the water resulted in the isomerization of the complex to solvated furous acid. Besides the above studies, the interaction of sulfur dioxide clusters with water molecules is not known to the best of our knowledge, although the SO 2 complexation with CO 2 , formaldehyde and its derivate, thioformaldehyde, and cyclohexanol was reported in which intermolecular interactions also occur through hydrogen bonding [21][22][23][24].
Gathering fundamental knowledge on the intermolecular interactions and binding topologies of biological important small molecules helps us to acquire knowledge on their role in the in vivo process. Especially, the S-O interactions not only control molecular conformations but also affect ligand-receptor binding modes, which consecutively stabilizes the folded structures of protein and regulate the enzymatic function [25]. Among the various interactions know, the chalcogen atoms can form noncovalent complexes with electron donors [21][22][23]26]. These interactions mainly ascend near the chalcogen atom at a low electron density region that is perpendicular to the negative regions and is designated as π-hole and is tend to be highly directional [27]. Recent studies on systems with different types of interactions were found to cooperatively strengthen the each of interactions and make the complexes more stable [28][29][30].
The present work aims to shed light on the role of chalcogen bond cooperativity in the growth pattern, structure, and stability of sulfur dioxide and hydrated sulfur dioxide clusters. Herein, we report our systematic study on (SO 2 )n and (SO 2 )n-H 2 O (n = 1-7) clusters using density functional theory methods. The nature of interaction and stability of the cluster are addressed using energetic analysis, molecular electrostatic potential (MEP), electron density difference (EDD), noncovalent interaction-reduced density gradient (NCI-RDG) analysis, and atoms in molecule (AIM) analysis. The cooperativity effect and the role of hydrogen bonding in (SO 2 )n-H 2 O clusters in cooperativity are investigated in terms of geometry, pair-wise interaction energy, MEP and AIM analysis. Energy decomposition analysis (EDA) was also carried out for (SO 2 )n-H 2 O clusters to provide insight into the role of driving forces for the hydration of sulfur dioxide clusters.

Computational details
The initial guessed geometries for (SO 2 )n and (SO 2 ) n-H 2 O (n = 2-7) are generated using the ABCluster code, which searches the global as well as the local minima of atomic and molecular clusters [31]. For more details on the theoretical basis of the ABCluster code, the reader is referred to the work of Zhang and Dolg [32]. The coordinates obtained in the ABCluster codes are further optimized at the DFT level with the meta-hybrid GGA functional PW6B95-D3 and Def2TZVP basis set, without imposing any symmetric restrictions. Previous benchmark studies on (SO 2 ) 2 by Tasinato and Grimme show that PW6B95-D3 provides the best agreement for binding energy with a mean absolute deviation of 0.3 kcal mol −1 [33]. Recently, Malloum and Conradie found the smallest mean absolute deviation for binding energies for pronated acetonitrile clusters [34]. Harmonic frequency calculations were carried out to confirm the minimum energy structures are in real minima on the potential energy surfaces. The energies of the most stable structures were corrected for basis set superposition error using the counterpoise correction method [35]. The stability of the hydrated clusters was evaluated by calculating thermochemical quantities of model formation reaction between water as one fragment and sulfur dioxide cluster as the second fragment. The total binding energy and the incremental binding energy of the sulfur dioxide cluster are computed to assess their stability.
The total binding energy of water to the sulfur dioxide cluster is computed as where E((SO 2 ) 2 -H 2 O) is the total binding energy of the hydrated sulfur dioxide cluster, E(SO 2 ) 2 is the total energy of the sulfur dioxide cluster and E(H 2 O) is the total energy of the water molecule. All the reported DFT calculations were carried performed using the Gaussian 16 suite of programs [36].
The PW6B95-D3/Def2TZVP method was used for the construction of electrostatic potential surfaces, NCI analysis, and NBO analysis. The wave function analysis-surface analysis suite (WFA-SAS) program was used to calculate quantitative electrostatic potential at 0.001 electron/Bohr 2 isodensity and to visualize the 3D surface [37]. The noncovalent interaction analysis (NCI) calculations were carried out with the Multiwfn program using larger grids and the isosurface is visualized using the Chemcraft program [38,39]. Atoms in molecules (AIM) analyses were carried out using the AIMALL program package [40].

Structures of (SO 2 )n clusters (n = 1-7)
In this section, we will focus on the putative global minima (PGM) of (SO 2 )n clusters (n = 2-7), whose optimized geometries are shown in Fig. 1. The other minimum energy geometries up to 12 numbers are provided in the supporting information Figs. S1-S6. The PGM of SO 2 dimer has a low symmetry geometry with C s symmetry in which one of SO 2 molecule is perpendicular to the other with an intermolecular S-O distance of 3.331 Å. The average S-O bond length in the PGM was 1.429 Å, which is shorter than observed for the SO2 molecule. By gas-phase infrared spectra study and ab initio calculations on SO 2 clusters, Ito and Hirabayashi have the recommend that the PGM be stable geometry [12,41]. There exits two low lying geometry with C i symmetry 2S-1 and 2S-2 has a petite twist in the O-S-O-S dihedral angle and are 0.55 and 0.56 kcal mol −1 less stable than the PGM. The C 2 symmetry geometry 2S-3 is 0.76 kcal mol −1 less stable than PGM. The existence of several low-lying isomers within a few kcal mol −1 energy suggests that extreme low-temperature studies are needed for studying the SO 2 clusters.
For the trimer SO 2 cluster, a cyclic ring with C1 symmetric with two SO2 molecules binding one SO 2 molecule was found to be the PGM. The average S-O bond length .093 Å respectively. The above values are longer and shorter than the PGM dimer cluster. The low lying isomer 3S-1 was shown in Fig. S2, with a chair conformation is just 0.07 kcal mol −1 in energy. Besides these, next four isomers were also in cyclic configuration; however, they were found to have different orientations of SO 2 molecules. For example, the 3S2-2 isomer exists in chair form but the angle of orientation of SO 2 differs from one other, while the 3S-5 isomer exists in a cup shape configuration. It is interesting to observe that the configuration other than cyclic forms such as the bent configuration (3S8-3S-11) is less stable than the cyclic configurations. Of the entire configuration, the linear 3S-12 was found to be the least stable.
The PGM and the 12 low-lying geometries along with selected bond parameters for the SO 2 tetramer clusters are provided in Fig. 1c and Fig. S3 respectively. The PMG of the tetramer cluster has a tetrahedral geometry in which the SO 2 molecules occupy the corners. The average S-O bond length in these cluster remains at 1.430 Å, while there is an increase in average intermolecular S-O bond length to 3.138 Å compared to the trimer PGM cluster. The first five low-lying isomers of tetramer clusters shown in Fig. S3 (4S-1-4S-5) are found to exist in the tetrahedral geometry, with a small change in the angle of orientation of SO 2 molecules. It is interesting to find that these are just 0.04-0.74 kcal mol −1 energy differences than the PGM. The butterfly-like configurations shown in Fig. S3 (4S-7, 4S-8, 4S-10, 4S-11) have nearly twice less stable as the tetrahedral configurations. Moreover, the planar geometry shown in Fig. S3 4S-12 is the least stable conformers. This indicates that 3D growth is preferred in the SO 2 clusters. It is worth point out that our previous studies on DMSO clusters show them to have ouroboros structure with 2D growth [42,43].
With the increase in the number of SO 2 molecules and the existence of several closely positioning isomers, the complexity of finding the PGM for larger clusters becomes more complex and computationally expensive [44]. Hence, for the pentamer to heptamer cluster, the possible guessed structures are obtained by adding one, two, and three SO 2 molecules to the lowest energy tetramer and trimer clusters in all possible directions. The PGM for the pentamer is shown in Fig. 1d and the low-lying isomers are shown in Fig. S4. In the pentamer to heptamer clusters, we have noticed that those SO 2 molecules whose oxygen atoms both are in bonding with nearly by the sulfur atom of sulfur dioxide have a shorter S-O bond distance. The average S-O bond length in pentamer, hexamer, and heptamer clusters remains at 1.430 Å, while the intermolecular S-O bond distance in these clusters is 3.216, 3.257, and 3.187 Å, respectively. Furthermore, in these clusters, the SO 2 molecules exist in near tetrahedron geometry when one, two, and three sulfur dioxides are removed from the pentamer, hexamer, and heptamer clusters.

Structure of (SO 2 )n-H 2 O (n = 1-7)
The sulfur dioxide molecule can interact with hydrogen atoms in water and make a hydrogen bond or the π-hole observed on the sulfur atom can interact with the electronrich oxygen atom of water to make a chalcogen bonding. The PGM of SO 2 -H 2 O hetrodimer is shown in Fig. 2a and the other three low-lying isomers (SH-1, SH-2, and SH-3) are shown in Fig. S7. The PGM of SO 2 -H 2 O was found to have the S-O chalcogen bonding with an intermolecular distance of 2.719 Å. The intramolecular S-O bond length has contracted from 1.440 in the free SO 2 molecule to 1.431 Å in the SO 2 -H 2 O complex. On the contrary, intramolecular O-H distance of water got stretched from 0.957 in free water to 0.959 Å. This may be envisaged due to the electron density transfer from the sulfur atom to the water molecule. We noticed three low-lying isomers (SH-1, SH-2, and SH-3) with H-bonding and intermolecular distances of 2.414, 2.430, 2.168 Å, with 1.72, 2.09, and 2.87 kcal mol −1 higher than the PGM, respectively. The lower stability of the H-bonded isomers shows that chalcogen bonds are stronger than the H-bonds (vide infra).
The (SO 2 ) 2 -H 2 O hetrotrimer can be in various isomeric forms. The PGM of the trimer cluster is shown in Fig. 2b and the other low-lying isomers are provided in Fig. S7. The trimer PGM has one S-O chalcogen with an intramolecular distance of 2.630 Å, which is shorter than observed in the SO 2 -H 2 O cluster. The average Ow-H-O distance between water and the SO 2 was 2.420 Å, while the average intramolecular S-O bond length was 1.432 Å which is slightly higher than observed in the SO 2 -H 2 O cluster. The low-lying isomer 2SH-1 has 0.78 kcal mol −1 higher in energy than PGM and exists in a cyclic chair conformation with two intermolecular chalcogen bond and one hydrogen bond. The next three low-lying isomers 2SH-2, 2SH-3, and 2SH-4 all exist in a cyclic conformation similar to the 2SH-1, but the orientation of the SO 2 molecule differ in these isomers. The linear isomers 2SH-5, 2SH-6, and 2SH-7 are less stable than the cyclic conformers. Comparison between the linear structures shows that 2SH-5 in which the two SO 2 molecules has chalcogen bond with water molecule is more stable than structures where there exists intermolecular chalcogen bonding between SO 2 molecules as in 2SH-6 and 2SH-7. From the above discussion, we can conclude that the structures which exist with chalcogen bonding alone and with only H-bonding are less stable and the existence of cooperativity between these two bonding's helps in stabilizing the PGM of SO 2 -water clusters.
The PGM of the hetrotetramer cluster occurs with two SO 2 molecules in hydrogen bonding with water and the other SO 2 to be in chalcogen S-O w bonding with water as shown in Fig. 2c. Interestingly, the molecules occupy the corners in a tetrahedral geometry. The SO 2 molecules are in intermolecular bonding with one another. The average intermolecular S-O bond length is shorter than the hetrotrimer but longer than the pristine SO 2 molecule. Furthermore, the intermolecular H-bonds between water and SO 2 got reduced to 2.161 Å compared to the 2.419 Å observed for the hetrotrimer cluster. The low-lying isomers of hetrotetramer are shown in Fig. S8. One may notice that all the low-lying hetrotetramer clusters have both hydrogen bonding and one or two chalcogen S-O w bondings with water. The PGM clusters of hetropentamer and hetrohexamer are shown in Fig. 2d, e have three intermolecular H-bonding and one S-O w bonding with water. The PGM clusters of hetroheptamer and hetrooctamer have three intermolecular H-bonding two chalcogen S-O w bonding with the water molecule.

Energetics of (SO 2 )n and (SO 2 )n-H 2 O (n = 1-7)
To understand the stability of the SO 2 and hydrated SO 2 clusters, we have computed their total binding energy, the average binding energy per SO 2 molecule, and the incremental binding energy. The total binding energy is obtained as the difference in total energy of the cluster with the size of the cluster multiplied by the total energy of the optimized SO 2 molecule. The incremental binding energy is computed as shown in bellow where E(SO 2 )n is the total energy of the sulfur dioxide cluster, E(SO 2 )n-1 is the total energy of n-1 th cluster, and E(SO 2 ) is the energy of the sulfur dioxide monomer. The computed energy parameters for the sulfur dioxide clusters are provided in Table 1. The computed binding energy and binding energy per SO 2 molecule for the sulfur dioxide cluster increase monotonically with the increase in the cluster size. The incremental binding energy decreases with the increase in cluster size and nears the saturation limit.
For the hydrated sulfur dioxide clusters, the binding energy increases monotonically, while the binding energy per SO 2 shows no regular trend. The IBE of hydrated sulfur dioxide clusters shows the increase in the values with the cluster size except for the hetrohexamer with five SO 2 molecules. The cluster adsorption energy is a useful parameter to understand the adsorption of sulfur dioxide cluster onto the water molecule which is computed as where E((SO 2 )n-H 2 O) is the total energy of the hydrated sulfur dioxide cluster, E(H 2 O) is the total energy of the water molecule, and E((SO 2 )n) is the total energy of the sulfur dioxide cluster. The cluster adsorption energy increases with the sulfur dioxide cluster size and reaches a near saturation. This implies that small sulfur dioxide clusters get adsorbed The comparison of energetic parameters between SO 2 and hydrated SO 2 clusters would help in understanding their relative stability and reactivity. Compared to the binding energy of the pristine SO 2 cluster, the hydrated SO 2 cluster has 53-70% lower values. Similarly, the incremental binding energy of hydrated SO 2 clusters is nearly half the value observed for the pristine SO 2 clusters. This shows that hydrated SO 2 clusters are relatively less stable and are more reactive.
The plot of the total binding energy vs the number of SO 2 molecules for the pristine SO 2 and hydrated SO 2 clusters is shown in Fig. 3a, b, which have a perfect linear correlation with a coefficient of 0.999. Previous studies on acetonitrile and interhalogen derivatives account such behavior to the existence of cooperatively in those clusters [45,46]. Molecular systems are composed of various types of interactions including the weak noncovalent interactions. When a pair of noncovalent interactions strengthen each other, it is called cooperative while when they weaken each other they are defined as operating in an anticooperative manner. Cooperativity implies that the sum of at least two interactions is larger than the simple addition of the individual interactions. In the absence of the cooperativity, the binding energy variation would be minimal and will be closer to the binding energy of the dimer cluster. Recently, Bartkowiak et al. found that in formamide and HCN clusters, the computed mean polarizabilities show exceptional linearity with the cluster size [47,48]. Our previous studies on linear DMSO clusters also show the existence of cooperativity when the mean polarizability is plotted as a function of cluster size [41,42]. To test the above hypothesis, we plotted the polarizability as a function of cluster size for pristine and hydrated SO 2 clusters as shown in Fig. 3c, d, and observed excellent linearly with a correlation coefficient of 0.999, indicating the existence of cooperativity in these clusters.
One can also examine the presence of cooperativity in these clusters by computing their pairwise interaction energies [49]. The pairwise energies for the trimer to heptamer for SO 2 clusters and dimer to octamer in hydrated sulfur dioxide clusters and their results are presented in Tables 2 and 3 respectively, while their numbering pattern is provided in Figs. 1 and 2. In the sulfur dioxide trimer, the interaction between I and II is higher than the other two pairs. It is interesting to observe the SO 2 molecule I have the maximum bonding with the nearby cluster. A similar observation has been noticed in the tetramer cluster, wherein the interaction energy is higher from the (II, III) in with II has maximum bonding than the other molecules. In the pentamer cluster, the highest interaction is observed for the pair (I, II) and the least interaction occurs between the pair (II, V) are at the nonbonding distance. In the hexamer and heptamer clusters, pair (III, IV) and (V, VI) which are in close contact and with maximum bonding have the highest interaction energy. The least interaction energies for the above clusters are for the pairs which are at the nonbonding distances (Table 2).
In the hydrated (SO 2 ) 2 -H 2 O hetrodimer cluster, the highest interaction is observed for the pair (II, III), the interaction between the sulfur dioxides. The H-bonding interaction between the pair (I, II) has the least interaction energy, while the interaction between water and sulfur dioxide (I, III) has interaction energy higher than the H-bonding dimer but far less than the sulfur dioxide dimer interaction. A similar scenario has been seen in the other larger clusters, wherein the hydrogen bonding has lower interaction energy and the intermolecular sulfur dioxide have higher interaction energies. In the case of hetroheptamer clusters, we observe a destabilization interaction between the water and the SO 2 . The highest interaction energy is observed in the case of SO 2 which is in bonding with water than SO 2 which is nonbonding with water. It is worth to point out that the pair-wise interaction energy is higher for the H-bonding dimers in the DMSO cluster [42,43]. Moreover, the interaction energy decreases with the increase in distance between the water and the sulfur dioxide molecule. Despite the above, the total binding energy shows excellent linearly when plotted against the number of SO 2 molecules, which confirms the cooperativity effect which stabilized the SO 2 and hydrated SO 2 clusters.

Nature of intermolecular interaction
To understand the nature of intermolecular interactions between sulfur dioxides molecules in sulfur dioxide clusters and water and sulfur dioxide clusters in hydrated sulfur dioxides, we used quantitative molecular electrostatic potential analysis (MESP), electron density difference (EDD) atoms in molecule (AIM), noncovalent interaction-reduced density gradient (NCI-RDG), and energy decomposition analysis. The acquired results are presented below.  [50,51]. The V s,max associated with the π-hole has a magnitude of 33.98 kcal mol −1 The MEP of SO 2 dimer is illustrated in Fig. 4b, in which the sulfur atom which is at the free end has a magnitude of 36.83 kcal mol −1 a value closer to the pristine SO 2 molecule. The V s,max on the sulfur atom which gets adsorbed to the oxygen atoms of adjacent SO 2 molecule got reduced to 29.59 kcal mol −1 . A similar situation has been observed on all the studied SO 2 clusters, in which the sulfur atoms at the terminal nearly retain their V s,max values. The V s,min value on SO 2 was observed on the oxygen atoms with a value of − 19.85 kcal mol −1 . As the name π-hole was misleading, in recent year, π-hole interactions are redesigned as  [54]. The presence of V s,max and V s,min on the surface of the clusters helps for the additional SO 2 adsorption on the surface and growth of cluster. The MEP of hydrated sulfur dioxide clusters are shown in Fig. 4g-l. In the hydrated dimer cluster, the V s,max on bonding side of sulfur got vanished, while the opposite side has a magnitude of 24.91 kcal mol −1 much less than observed in the case of pristine SO 2 cluster. This implies that if SO 2 gets added to it will have lower stability. Moreover, the V s.max on hydrogen atom was 50.02 kcal mol −1 . Hence, the addition of SO 2 occurs at the site, leading to an H-bonding. In the case of the timer ((SO 2 ) 2 -H 2 O) cluster, V s,max on hydrogen atom was 37.39 kcal mol −1 which is larger than the π-hole value observed on the free SO 2 molecule. Thus, in the tetramer cluster ((SO 2 ) 3 -H 2 O), the SO 2 prefers to be an H-bonding state. In the larger clusters, the V s,max on hydrogen atom is far less than the values observed on the sulfur atom π-hole, and hence, bonding occurs at the site leading to the chalcogen bonding. The reduction of the V s,max on hydrogen atom upon complexation implies that electron density on oxygen atoms are polarized towards the positive hydrogen and tend to move from the oxygen atom on complexation. In the pentamer cluster, a V s,min value of 23.49 kcal mol −1 on the oxygen atom of water was observed. Henceforth in the octamer cluster, we observe the presence of two intermolecular chalcogen bonding between two sulfur dioxide and water molecule. Thus, the value of V s,max or V s,min dictates the direction of the interaction of new SO 2 molecule in larger clusters.

Molecular electrostatic potential and electron density difference (EDD) analysis
To provide a visualization of the changes in the electron density on complexation, the electron density difference (EDD) maps [55] for the hydrated sulfur dioxide clusters have been plotted and depicted in supporting information Fig. S12. Red regions indicate regions of an increase of electron density, and blue regions are associated with regions of decreased electron density. These EDD maps are obtained by subtraction of electron densities of the complex and the corresponding monomers on their complex geometries. In the case of SO 2 clusters, the monomers are SO 2 molecule which binds with two oxygen as unit one and other as unit two. In the case of hydrated complexes, water is considered unit one and SO 2 cluster as unit two. In the pristine SO 2 clusters, the noticeable aspect is the build of blue regions encamping the regions where sulfur atom which is in S-O chalcogen bonding. In the larger clusters, these regions were found to merge generating a large surface. Also, notice that the oxygen adjacent to the sulfur atoms gains electron density envisaged by the red regions. Thus, a substantial electron density shift from the sulfur (donor) to the oxygen atom (acceptor) occurs during the cluster formation. This was also observed in the hydrated sulfur dioxide clusters. Besides the above, we also notice a presence of red regions fully surrounding the water molecules in hydrated clusters. This corresponds to the gain of electron density by water molecules similar to the previously reported hydrated complexes.

Atoms in molecules
The quantum theory of atoms in molecules provides detailed information about the existence, strength, and characteristics of different intermolecular interactions in various complexes [55][56][57][58]. The presence of a bond critical point (BCP) between sulfur dioxide and sulfur dioxide/water can be used to characterize the chalcogen bonds. Molecular graphs for all the pristine SO 2 cluster and hydrated SO 2 clusters are depicted in Fig. 5. The electron density properties at the BCPs for sulfur dioxide clusters are provided in Table S1 and for hydrated clusters in Table S2. The electron density ρ(r) at the intermolecular S-O in sulfoxide clusters is in the range of 0.005-0.0133 au which fall in the accepted range of chalcogen bonding [59,60]. However, these values are much smaller than observed in hydrogen bonds (0.002-0.035 a.u.), which indicates that these interactions are weak. The Laplacian of electron density were all positive values in the range 0.019-0.052 au and are closed-shell weak interactions. The ratio of kinetic and potential energy density at BCP, -G(r) /V(r) is greater than unity at all the BCPs revealing the presence of noncovalent interactions. Ellipticity is a measure  Table 4 The summation of QTAIM parameters corresponding to intermolecular bonding between the SO2 molecules and water and SO2 molecules   [63,64]. The computed summation of electron density properties at the BCPs is provided in Table 4. The results indicate that the summation of (r) and ∇ 2 (r) values of clusters increases with the size of the cluster and is always higher than that of dimers. Recently, Esrafili and Mohammadian-Sabet observed that the average (r) and ∇ 2 (r) values of larger clusters of OCS and OCSe are greater than those of dimers and associated it with the existence of chalcogen bond cooperativity due to the redistribution of electronic density in the clusters [54]. In our present study, we observe a very (r) and ∇ 2 (r) value for the BCPs with H-bonding, while the other sulfur dioxide part of hydrated clusters were found to have higher low (r) and ∇ 2 (r) compared to the pristine sulfur dioxide clusters. Thus, the presence of H-bond and chalcogen bond act cooperatively which resulted in the redistribution of electronic density in the hydrated clusters.

EDA analysis
Decomposition of the total interaction energy into repulsive and attractive terms opens a window to know the nature of the interaction. In Table 5, we provide the energy decomposition analysis results for the hydrated system obtained using the Morokuma's bond energy decomposition method [65]. In the Morokuma's method, the attractive terms are partitioned into electrostatic, orbital, dispersion energies. It can be seen from Table 5, the destabilization Pauli's terms increase with the cluster size and are the highest for the pentamer cluster (SO 2 ) 4 -H 2 O, presumable due to the addition of one more

Non-covalent interaction-reduced density gradient analysis
To visualize the noncovalent interactions, including weak van der Waals interaction, hydrogen bonding, and steric repulsion, NCI-RDG analysis was carried out on the sulfur dioxide and hydrated sulfur dioxide clusters [66][67][68].
The 2D NCI graphs with inserted RDG isosurfaces for the sulfur dioxide and hydrated sulfur dioxide clusters are provided in Fig. 6 and supporting information in Fig. S13. In the NCI graphs, low density and low gradient spikes that appear on the negative λ 2 (r) values are for attractive interactions (specifically electrostatic/hydrogen bonding) and positive values are for the repulsive interactions. For the SO 2 molecule, the NCI graphs is shown in Fig. 6a, spikes are not visible due to the absence of attractive and repulsive forces. In the dimer (SO 2 ) 2 clusters, there were spikes at sign ( λ 2 ) = -0.006X and 0.006X au appeared due to the attractive chalcogen bonding and sulfur-sulfur repulsive forces respectively. These interactions are evident from the RDG isosurfaces shown at the right top of Fig. 6b with green and yellowish-brown patches. In the trimer (SO 2 ) 3 cluster, we also notice the appearance of a pale green patch surrounded by a yellowish-brown patch between the oxygen-oxygen atoms. Their corresponding spikes were also noticed in the NCI graph. Interestingly, the presence of such interaction is not visible in the AIM analysis. The number of S-O bonds increases with the increases in cluster size; thus, the number of spikes in the attractive regions also got increased with the cluster size. The negative λ 2 (r) values are very close and in the range -0.006 to -0.015 au. In the cluster sizes above n > 5, the spikes got merged and hence not visible as distinctly for each of the chalcogen bonding.
In the hydrated sulfur dioxide cluster SO 2 -H 2 O, we notice the presence of a spike at −0.02 au and a broad spike in the region 0.018 au. The corresponding RDG isosurfaces show the presence of a blue patch indicating the presence of a strong electrostatic component between the water and SO 2 molecule. Besides, a red patch at the edges of the blue patch is also noticed which is due to the repulsive forces between them. In the trimer cluster (SO 2 ) 2 -H 2 O, the S-O chalcogen bond between a water molecule and SO 2 molecules appears in the more negative region at −0.028 au and green spikes due to intermolecular chalcogen bonding between the SO 2 appears between 0.010 and 0.015 au. It is worth point out the intermolecular chalcogen regions are observed in the range -0.006 to -0.015 a.u for the pristine SO 2 clusters. In the tetramer and pentamer clusters, the blue spikes are shifted to more negative regions along with the increase in the number of green spikes. In the larger clusters, the spike due to the intermolecular chalcogen bonding between water and SO 2 molecules gets splattered, indicating the intermolecular SO 2 chalcogen bonding's cooperatively act to enhances the increased interaction between the water and SO 2 molecule. These results are in good agreement with the geometry, pairwise energetic, MESP, AIM, and EDA analysis.

Conclusion
In conclusion, the pristine SO 2 and binary intermolecular clusters formed by SO 2 molecules and water was investigated using density functional theory and wave functional methods. The putative global minimum of SO 2 molecules has 3D growth pattern with tetrahedral configurations and several closely positioning isomers were observed. In the hydrated SO 2 clusters, the structures with chalcogen bonding alone or with only hydrogen bonding are less stable than structures with both bonding. The total binding energy of the SO 2 and hydrated SO 2 clusters increases with the increase in the number of SO 2 molecules. For the pristine SO 2 clusters, incremental binding energy decreases with the increase in cluster size and nears the saturation limit. The cluster adsorption energy for SO 2 clusters on water increase with the sulfur dioxide cluster size and reaches a near saturation. Relatively less binding energy and incremental binding energy of hydrated SO 2 clusters intend that they are less stable/highly reactive compared to pristine SO 2 clusters. The pristine SO 2 and the hydrated SO 2 cluster show cooperativity effects which are evident from the excellent linearity plot of binding energy/polarizability vs the number of SO 2 molecules which was supported by the pairwise energies computed on the PGM clusters.
Molecular electrostatic potential (MEP) analysis shows the oxygen atoms which are at the terminal have enhanced values generating an intermolecular S-O chalcogen bond between the SO 2 molecules helps for the growth of clusters. MEP also explains the reason for the facile formation of S-O chalcogen than hydrogen bonding in hydrated sulfur dioxide. Atoms in molecule (AIM) analysis characterize the bonds to be weak and electrostatic dominant. The observed ρ(r) values at BCPs with hydrogen bonds are lower than the values observed for the S-O chalcogen. EDA analysis shows electrostatic interaction is dominated in complexes with more intermolecular chalcogen bonding and orbital interaction for systems with intermolecular H-bonding. The reduced density gradient (RDG) analysis of sulfur dioxide clusters has blue patches and green patches due to S-O chalcogen bonding O-O electrostatic interaction. The RDG analysis of hydrated sulfur dioxide clusters shows intensive blue patches and green patches for the existence of S-O chalcogen and hydrogen bonding respectively. Thus, the presence of strong electrostatic S-O chalcogen interaction and weak H bonds acts cooperatively and stabilize the hydrated sulfur dioxide clusters.