The solvation of SOH group in hydrated HSO 4 − (H 2 O) n clusters

Abstract

The S=O stretching and SOH bending peaks in the vibrational spectra of HSO₄⁻(H₂O) _n , with n up to 6, are analyzed by both harmonic analysis and ab initio molecular dynamics simulation. The SOH bending mode is found to be much more sensitive to the extent of hydration and to the fluctuation of hydrogen bonds than the S=O stretching mode. The SOH donor hydrogen bond is gradually stabilized by n = 4, and further shortened up to n = 6, which is the key factor to understand the trend of evolution observed in the infrared multiple photon dissociation spectra.

Classification dependents by the distance of hydrogen bond.

Introduction

Bisulfate ion is one of the trace sulfur species in the atmosphere, in addition to H₂SO₄, SO₄²⁻, and their precursor SO₂ from both natural and anthropogenic sources.^1,2 Although not as abundant as its neutral cousins, HSO₄⁻ could play an important role in aerosol formation, and its interaction with water molecules is one of the important factors underlying the homogeneous nucleation processes in the atmosphere.^3–8 Hydrated bisulfate cluster ions, HSO₄⁻(H₂O) _n , are molecular models to probe such interactions and the subject of a number of experimental and theoretical studies. Theoretical calculations have identified many local energy minimum structures,^9–12 with varying harmonic frequencies and hydration energy. In the 600–1800 cm⁻¹ region, the vibrational profiles of HSO₄⁻(H₂O) _n have been measured by infrared multiple photon dissociation (IRMPD) spectroscopy, with n up to 16, and assigned by comparison with the harmonic spectra obtained for optimized structures. ¹³ More recently, the photoelectron spectra for HSO₄⁻(H₂O) _n , with n up to 2, have also been reported. ¹⁴

Bisulfate ion is also one of the conjugate bases, ubiquitous in aqueous solution, and often encountered in general chemistry lab. Micro-solvation details around such polyatomic ions are of fundamental interests, and recent spectroscopic studies have shed considerable lights into the geometry, electronic structures, and dynamics in such hydrated cluster ions.^15–23 Particularly relevant to HSO₄⁻(H₂O) _n are the recent studies on H₂PO₄⁻(H₂O) _n and HCO₃⁻(H₂O) _n .^24–28 All three singly charged anions have two kinds of groups, X=O and XOH, with the XO stretching (υ(X=O)) and the XOH bending (δ(XOH)) coupled together with considerable overlap. Furthermore, δ(XOH) is more sensitive than υ(X=O) to the fluctuation in hydrogen bond (HB) distances. As the number of water molecules increases and the XOH group becomes better solvated, it would produce a trend of changes in the pattern of υ(X=O) and δ(XOH) peaks, which could be explained by the dynamic of HB fluctuation.^26,28

Such dynamics is also expected to be important for HSO₄⁻(H₂O) _n . Moreover, with a pK_a value of 2.0, HSO₄⁻ is more acidic than H₂PO₄⁻and HCO₃⁻, and the acidic dissociation of HSO₄⁻, when the number of water molecules increases, has been reported by a recent research. ²⁹ In this paper, we report ab initio molecular dynamics (AIMD) simulation results on HSO₄⁻(H₂O) _n , with n up to 6, and show that the solvation of SOH is also the key to understand the trend of change in the υ(S=O) and δ(SOH) peaks with increasing cluster size, as observed in the IRMPD experiment. ¹³

Computational methods

The optimization of 0 K structures and harmonic analysis are performed by GAUSSIAN 16 package (Rev A.03) at the B3LYP/6-311+G(d,p) level.^30–32 Relative energies include zero-point vibrational energies. A scaling factor of 1.071 is used for low-frequency region below 2000 cm⁻¹. ³³ The resulting stick spectra are convoluted by a Gaussian line shape function with a full width at half-maximum line width of 15 cm⁻¹.

AIMD simulations at finite temperatures are performed by the CP2K package. ³⁴ A cluster under simulation is put at the center of a periodic cubic box with a length of more than twice the largest interatomic distance in the cluster and with the effect of the periodical charge density images eliminated by the decoupling technique developed by Martyna and Tuckerman. ³⁵ At every time step, the potential energy and atomic forces are calculated within the framework for density functional theory (DFT). The wave functions are represented by double-zeta Gaussians, and the electron density is represented using a mixed basis set of Gaussians and auxiliary plane waves with a cut-off of 320 Rydberg. The Goedecker–Teter–Hutter (GTH) typed pseudopotentials are employed for the atomic cores. For the exchange and correlation energy, BLYP functional are employed with dispersion corrected by Grimme’s D3 scheme. ³⁶

The atomic motion is treated by Newtonian mechanics, with the temperature controlled by a Nose–Hoover thermostat and with a time step of 0.5 fs. The temperature is first scaled at the desired value for 3 ps, followed by a data collection run in the NVE (micro-canonical ensemble) ensemble for 15 ps. The value of dipole moment is calculated and collected at each time step, so that a vibrational spectrum can be directly simulated by the Fourier transform of the dipole time-correlation function (DTCF)

$$\alpha \left(\omega \right)=\frac{2\pi \beta {\omega }^2}{3n\left(\omega \right)cV}\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}dt\overrightarrow{M}(t).\overrightarrow{M}(0)\exp (i\omega t)$$

where $\beta =1/kT$ , $n(\omega )$ is the refractive index, c is the speed of light in vacuum, and V is the volume. $\overrightarrow{M}$ is the total dipole moment of the system, calculated by the polarization including both ionic and electronic contributions. Dynamic and anharmonic effects are automatically taken into accounts in such a scheme, although quantum effects at low temperature are not included.

Results and discussion

HB fluctuation in HSO₄⁻(H₂O)

In Figure 1, three optimized structures are shown. The energy gap among them is less than 2 kJ mol⁻¹. 1-1 and 1-2 have been identified before, and in both structures, H₂O forms two donor HBs with the O atoms of two S=O bonds, which is usually labeled as a DD-H₂O (where D denotes the donor). In 1-1, H₂O further accepts an HB from the SOH group, making it a fully solvated ADD water (where A denotes the acceptor). But the energy of 1-1 is raised slightly relative to 1-2, due to the distortion in the two donor HBs. Such an ADD water is not observed before either in H₂PO₄⁻(H₂O) or HCO₃⁻(H₂O).^26,28 There are only two X=O groups in HCO₃⁻ or H₂PO₄⁻, and a DD water solvating these two groups is pointed away from the XOH groups, making it impossible for another acceptor HB. In contrast, there are three X=O groups in HSO₄⁻, arranged in a trigonal pyramid, and a DD water could be bent up slightly to accept an HB from the SOH group.

Figure 1.

Optimized structures of HSO₄⁻(H₂O). Relative energy in parenthesis is in kJ mol⁻¹.

In AIMD simulations, these two donor HBs in 1-1 or 1-2 are maintained most of the time during a 100-ps period at 150 K, as shown in Figure 2. This is an interesting contrast to the DD water in H₂PO₄⁻(H₂O) and HCO₃⁻(H₂O) during AIMD simulations,^25,28 where one of the donor HBs is easily broken. For these two molecules, the water molecule would then migrate to a position between X=O and XOH with X=O···HOH distance below 1.8 Å, resulting in a stronger typical HB. For HSO₄⁻(H₂O), a structure with the H₂O between SOH and S=O groups is also found as 1-3 in Figure 1. But the HB distance of S=O···HOH at 2.01 Å is of typical strength, due to the fact that the negative charge on HSO₄⁻ is distributed over three S=O groups, while on H₂PO₄⁻ and HCO₃⁻, it is distributed over only two X=O groups. AIMD simulations also show that 1-3 is thermally unstable and even at 20 K, it transforms into 1-1.

Figure 2.

Left panel: fluctuation in the dihedral-O₅S₁O₃H₆ during AIMD simulations for cluster 1-1. Right panel: fluctuation in the HB distances during AIMD simulations at 150 K for 1-1.

The dynamics of HSO₄⁻(H₂O) at finite temperature is dominated by the transformation between 1-1 and 1-2. At 20, 100, and 150 K, the cluster would end up with SOH···OH₂ (as in 1-1), even when an AIMD simulation starts without it (as in 1-2). But by 150 K, rotation around the S−OH bond is observed, which could occasionally break the HB of SOH···OH₂ (O3H6···O7), as shown in Figure 2. By 200 K, this HB cannot be maintained, as S−OH rotation becomes quite extensive, and the structure constantly changes back and forth between 1-1 and 1-2.

S=O stretching and SOH bending in HSO₄⁻(H₂O)

Although SOH···OH₂ is a weak HB with a distance ~2.5 Å, it has a significant impact on the vibrational profile between 1000 and 1400 cm⁻¹. There are four peaks in this region: the symmetric stretching of the three S=O bonds, υ_S(S=O); the asymmetric stretching of the two S=O bonds, cis to the HSO group, υ_AS-C(S=O); the asymmetric stretching of the S=O bond, trans to the SOH group, υ_AS-T(S=O); and finally, the bending of S-O-H, δ(SOH). For the naked HSO₄⁻, the frequency ordering is

$${\upsilon }_{\mathrm{S}}\left(\mathrm{S}=\mathrm{O}\right)<\delta \left(\mathrm{SOH}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{C}}\left(\mathrm{S}=\mathrm{O}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{T}}\left(\mathrm{S}=\mathrm{O}\right)$$

and there are couplings between δ(SOH) and υ_AS-T(S=O). Upon solvation by one H₂O, when the SOH···OH₂ HB is not formed as in 1-2, there is little change in the ordering or the relative intensity of these four peaks. As shown in Figure 3, the harmonic spectrum of 1-2 is essentially the same as that of HSO₄⁻. In line with our previous studies on H₂PO₄⁻(H₂O) _n and HCO₃⁻(H₂O) _n ,^26,28 the X=O stretching is not very sensitive to the solvation of X=O bonds. However, when the SOH···OH₂ HB is formed as in 1-1, the change is quite significant. Even though at an equilibrium distance of 2.54 Å, it should be as weak HB. δ(SOH) is 1224 cm⁻¹ in 1-1, which is blue-shifted comparing to 1165 cm⁻¹ in 1-2 and 1166 cm⁻¹ in naked HSO⁴⁻. So that we now have υ_AS-C(S=O) < δ(SOH). υ_AS-T(S=O) is 1331 cm⁻¹ in 1-1, which is blue-shifted comparing to 1292 cm⁻¹ in 1-2, as it has a mixed component of δ(SOH). The resulted harmonic spectrum of 1-1 is quite different from that of 1-2 or of the naked HSO₄⁻ (Figure 3).

Figure 3.

Harmonic spectra for naked HSO₄⁻, clusters 1-1 and 1-2 with scaling factor of 1.071, and experimental IRMPD spectra for HSO₄⁻(H₂O). ¹³

The DTCF spectra for 1-1 produced by AIMD simulations at various temperatures are shown in Figure 4. The 20 K spectrum is basically the harmonic spectrum, when the cluster is frozen around the equilibrium structure. The strength of the SOH···OH₂ HB is overestimated by the BLYP functional employed in our AIMD simulations, which is common among pure functionals. As a result, both δ(SOH) and υ_AS-T(S=O) shift further to the blue. Up to 150 K, the δ(SOH) and the three υ(S=O) peaks are only slightly broadened. But at 200 K, the cluster is in a constant dynamic transformation between 1-1 and 1-2, and the peaks are mixed and much broadened, in agreement with the experimental IRMPD spectrum. It lends support to the previous suggestion that the HSO₄⁻(H₂O) clusters measured in IRMPD experiment contained the combined contribution of both 1-1 and 1-2.

Figure 4.

Harmonic and DTCF spectra for 1-1 at various temperatures comparing to experimental IRMPD spectrum. ¹³ A scaling factor of 1.071 is used in the harmonic spectra for this region. No scaling factor is applied to the DTCF spectrum.

The effects of SOH solvation on SOH bending and S=O stretching

The strength of SOH···OH₂ HB also provides the key to understand the evolution of spectral region between 1000 and 1400 cm⁻¹. In my work, I did some classification and quantification research by the distance of this HB: naked SOH with infinite HB distance; solvated SOH with HB distance above 2 Å; and solvated SOH with HB distance under 2 Å. Below n ⩽ 4, a number of clusters have been optimized with the SOH group being unsolvated, and harmonic spectra are shown in Figure 5. The S=O groups become more solvated with more water molecules, and meanwhile the inter-water HB network becomes more elaborate. But such changes do not have any bearings on the pattern of the δ(SOH) and υ(S=O) peaks. As long as SOH is unsolvated, the patterns of δ(SOH) and υ(S=O) peaks are all similar to those for unsolvated HSO₄⁻, with the same frequency of ordering

$${\upsilon }_{\mathrm{S}}\left(\mathrm{S}=\mathrm{O}\right)<\delta \left(\mathrm{SOH}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{C}}\left(\mathrm{S}=\mathrm{O}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{T}}\left(\mathrm{S}=\mathrm{O}\right)$$

Figure 5.

Optimized structures of HSO₄⁻(H₂O)_n, n = 0–4, with the SOH group being unsolvated. A scaling factor of 1.071 is used in the harmonic spectra for this region.

δ(SOH) is 1160–1166 cm⁻¹ and υ(S=O) is 1045–1047 cm⁻¹ for all these clusters with the unsolvated SOH group. It means the position of vibrational modes will not change with the size evolution.

Figure 6 shows three structures in which the SOH group is solvated by a typical donor HB with a distance above 2 Å. The acceptor H₂O, while forming two donor HB with S=O groups, is not solvated further by any other water molecule. Their spectra between 1000 and 1400 cm⁻¹ are similar, and the frequency ordering is switched to

$${\upsilon }_{\mathrm{S}}\left(\mathrm{S}=\mathrm{O}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{C}}\left(\mathrm{S}=\mathrm{O}\right)<\delta \left(\mathrm{SOH}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{T}}\left(\mathrm{S}=\mathrm{O}\right)$$

Figure 6.

Optimized structures of HSO₄⁻(H₂O)_n, n = 1–3, with the SOH group being solvated by an HB with a distance above 2 Å. A scaling factor of 1.071 is used in the harmonic spectra for this region.

δ(SOH) is 1224 cm⁻¹ for 1-1, 1236 cm⁻¹ for 2-2, and 1226 cm⁻¹ for 3-4, quite similar to each other. Also, the frequencies of vibrational modes are not sensitive to the cluster size.

Starting at n = 2, the water molecule accepting SOH in hydrogen bonding could form donator HB with other water molecules. For 2-1 and 3-1, shown in Figure 7, this water molecule forms one donor HB with another H₂O and another HB with one S=O group, and the HB distance of SOH···OH₂ is shortened to below 2 Å. As a result, their δ(SOH) is moved further to the blue, and the frequency ordering is

$${\upsilon }_{\mathrm{S}}\left(\mathrm{S}=\mathrm{O}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{C}}\left(\mathrm{S}=\mathrm{O}\right)<{\upsilon }_{\mathrm{AS}-\mathrm{T}}\left(\mathrm{S}=\mathrm{O}\right)<\delta \left(\mathrm{SOH}\right)$$

Figure 7.

Optimized structures of HSO₄⁻(H₂O) _n , n = 2–6, with the SOH group being solvated by an HB with a distance below 2 Å. A scaling factor of 1.071 is used in the harmonic spectra for this region.

δ(SOH) is 1406 cm⁻¹ for 3-1, 1447 cm⁻¹ for 4-1, 1454 cm⁻¹ for 5-1, and 1479 cm⁻¹ for 6-1. As the HBs shorten in this group of clusters, the position of δ(SOH) is more sensitive to size evolution comparing to naked SOH and HBs above 2 Å. However, it does not affect the overall pattern of υ_S(S=O) < υ_AS-C(S=O) < υ_AS-T(S=O) < δ(SOH).

Furthermore, δ(SOH) is decoupled from the υ(S=O) modes, in contrast to the previous two types of structures, for which there is considerable coupling between δ(SOH) and υ_AS-T(S=O). When the water molecule solvating an XOH group is connected to other water molecules by HBs, the impact of δ(SOH) is absorbed by surrounding water molecules and lose its effects on S=O stretching. Such decoupling between XOH bending and X=O stretching has been observed in the cases of H₂PO₄⁻(H₂O) _n and HCO₃⁻(H₂O) _n .^26,28

For 4-1, 5-1, and 6-1 (Figure 7), the water molecule accepting SOH in hydrogen bonding forms two donor HBs with two other water molecules, as the network of H₂O becomes more extensive with increasing cluster size. The three υ(S=O) peaks do not change much, as these modes involve the vibration of S and O atoms which are much heavier than H, while δ(SOH) is shifted further to the blue.

The evolution of SOH bending and S=O stretching with cluster size

With the effects of SOH solvation understood, it is more straightforward to assign the region between 1000 and 1400 cm⁻¹ for the larger clusters. The SOH group is well solvated for n ⩾ 5. As shown in Figure 8, peaks E, D, and C are due to the three S=O stretching modes, υ_S(S=O), υ_AS-C(S=O), and υ_AS-T(S=O), respectively, in agreement with the previous assignment. ¹³ Higher in frequency than these three peaks is peak B, due to δ(SOH), as the HB of SOH···OH₂ is strengthened by such cluster sizes. The three υ(S=O) peaks (C, E, and D) are less sensitive to n value increases. Since these modes do not directly involve the movement of H atoms. It is perfectly reproduced by both harmonic analysis and DTCF spectra obtained by AIMD simulations shown in Figure 8. Due to the overestimate of HB strength, the calculated gap between peak B/δ(SOH) and other υ(S=O) peaks in the DTCF spectrum is 1450 cm⁻¹ to 1046 cm⁻¹ for 6-1 in 100 K, larger than 1335 cm⁻¹ to 1040 cm⁻¹ in the experimental values for n = 6. For both n = 5 and n = 6, when temperature arises to 200 K, due to the intensity of vibration, the fluctuation of bond length will make the peaks broadening. That is why the peaks are not sharp in 200 K comparing to lower temperature.

Figure 8.

DTCF spectra for 5-2 and 6-1 comparing to experimental spectra for HSO₄⁻(H₂O) _n , n = 5–6. ¹³

The overall pattern for these four peaks changes little as n increases from 9 to 11, except that the height of peak B/δ(SOH) drops appreciably. The flattening of the δ(XOH) with increasing n peak has been observed in the cases of H₂PO₄⁻(H₂O) _n and HCO₃⁻(H₂O) _n .^26,28 It can be again attributed to the strengthening of the XOH···OH₂. As shown in Figure 9, 3-1 shows the changes in values of δ(SOH), υ_AS-T(S=O), and υ_AS-C(S=O), when the HB distance of SOH···OH₂ is constrained to a specific value. While υ_AS-T(S=O) and υ_AS-C(S=O) change insignificantly, δ(SOH) is sensitive to the fluctuation in HB distance, especially when its value is below 2 Å.

Figure 9.

Simulated harmonic frequency of HB SO-H6···O7 for 3-1 with constrain bond distance, 1.63–2.23 Å with scaling factor of 1.071.

The spectral assignment for the smaller cluster with n = 2−4 is more complicated. For example, the SOH···OH₂ distance for the optimized 3-3 shown in Figure 10 is below 2 Å, and both its harmonic and 100 K DTCF spectra are similar to the spectra shown in Figure 8 for larger clusters. However, when the simulation temperature is raised to 200 K, the HB of SOH···OH₂ cannot be maintained, while the position of peak B is shifted to the red and its intensity is raised relative to peak C and D, in agreement with the experimental result. ¹³ Another example is 4-2, shown in Figure 10, in which there is a four water ring and the SOH group does not form a donor HB. During AIMD simulation, 4-2 is maintained even at 200 K. In both cases, peak B is caused by υ_AS-T(S=O), similar to spectra shown in Figures 6 and 7. Overall, the height of peak B can be taken as an indication for the solvation of SOH group by a donor HB. For n ⩽ 4, this HB is either weak or broken, producing a relative strong peak B which is relatively υ_AS-T(S=O). For n ⩾ 5, this HB is strong and the acceptor H₂O is connected to other water molecules by HB network, producing a lowered peak B which is caused by δ(SOH). Similar to n = 5 and n = 6, in 200 K, peaks will also broaden for n = 3 and n = 4 due to the fluctuation of bonds, but it is still easy to clarify the peak assignments.

Figure 10.

Harmonic, DTCF, and experimental IRMPD spectra for 3-3 and 4-2. ¹³ Scaling factor of 1.071 is used in the harmonic spectra. No scaling factor is applied to the DTCF spectra.

The evolution of SOH in larger cluster size

In our previous study, ²⁹ we found that although bisulfate anion is a weak acid with a pKa of 2.0, HSO₄⁻(H₂O) _n is completely dissociated at a moderate size of n = 16 and also dissociates in the liquid environment. For such charged clusters, the size can be precisely selected, making it possible to probe the molecular details about the acidic dissociation from its onset around n = 12 to its completion around n = 16. However, in the small cluster size of n = 1–6, the SOH group remains without dissociation.

Conclusion

In summary, the evolution of the IRMPD spectra for HSO₄⁻(H₂O) _n versus the cluster size can be attributed to the solvation of SOH group. For n ⩽ 4, SOH is not well solvated, which is indicated by a relative strong peak B, due to υ_AS-T(S=O). There is a significant coupling between υ_AS(S=O) and δ(SOH) peaks. For 5 ⩽ n ⩽ 6, SOH is well solvated, with a donor HB shorter than 2 Å. The spectral pattern is υ_S(S=O) < υ_AS-C(S=O) < υ_AS-T(S=O) < δ(SOH), with δ(SOH) decoupled from υ_AS(S=O) and broadened. It can be concluded that the SOH donor HB is the key factor to understand the trend of evolution observed in IRMPD spectra.