Background

The ionic hydrogen bond, which has a stronger intermolecular interaction than the typical hydrogen bond, plays significant roles in various biological and chemical processes[1–4], such as protein folding[5] and molecular recognition[6, 7]. Meanwhile, the ionic hydrogen-bonded structures provide a powerful lens through which the important phenomena of ion hydration could be further understood[8–11]. Intensive studies[12–15], with experimental techniques[8, 16–22] and theoretical approaches[23–38], have been focused on the detailed structures of the hydrated halide ion clusters.

Experimentally, the formation enthalpies of small hydrated halide ion clusters in the gas phase have been studied using a pulsed electron beam high-pressure mass spectrometer[21]. Diffraction studies[22] have provided limited structural information on these clusters, for example, the average distance between fluoride/chloride and oxygen is 2.6 to approximately 2.9/3.1 to approximately 3.5 Å. Recently, vibrational spectroscopy using argon predissociation technique has provided more abundant structural information about the fundamentals and first two overtones for X−(H2O) n (X = F, Cl; n = 1–4) clusters[8, 16–18].

On the other hand, to understand ion hydration at the molecular level, intensive theoretical works have focused on the structures and binding energies of hydrated ion water clusters[23–26]. For example, Kim and coworkers reported various local minima and transition state structures for these clusters with n = 1 to 6 through high-level ab initio molecular orbital (MO) calculations[23–25]. Recently, Kamp et al. studied the solvation of fluoride and chloride ion water clusters with n ≤ 20 using the effective fragment potential (EFP) method and Monte Carlo simulations[26]. In addition, other theoretical efforts, such as quantum mechanical/molecular mechanical (QM/MM)[28] and molecular dynamics (MD) with potential including polarization term[29, 39], are also devoted to these clusters. The interaction potentials are essential to the theoretical studies and high-level ab initio calculations are considered to be the most reliable. However, due to the bottleneck of computational limitations, it is formidable to use such expensive potentials to deal with the energetics and geometries of these clusters with respect to a large number of water molecules.

As an alternative solution, semiempirical potentials have been successfully applied into various systems with extremely cheap computation cost[40, 41]. Recently, the PM6 semiempirical potential[42, 43] has been successfully applied into the quantum chemistry study of neutral hydrogen-bonded systems, such as glycine-water clusters[44] and hydrogen storage material MOF-5[45]. The PM6 results show reasonable accurate structures and stabilization energies compared to high-level ab initio calculations. More recently, new semiempirical potentials of PM6-DH[46], PM6-DH+[47], and PM6-DH2[48], which improve PM6 by including two important corrections of dispersion energy and hydrogen bonding, have shown promising results for extended hydrogen-bonded complexes. To the best of our knowledge, however, there is still no report on ionic hydrogen-bonded systems with PM6 or PM6-DH+ semiempirical potentials.

In this study, thus, we apply the semiempirical potentials to X−(H2O) n (X = F, Cl; n = 1–4) clusters for investigating the various possible geometries. On the other hand, ab initio MO and density functional theory (DFT) methods are also utilized to examine the performance of the semiempirical potentials, with respect to the dependence of the number of water molecules on the energetics and geometries of the ionic hydrogen-bonded clusters. Those are very important data for our near-future study of on-the-fly semiempirical MD or path integral MD simulation.

Methods

First, we have explored the optimized energies and geometries of X−(H2O) n (X = F, Cl; n = 1–4) clusters by semiempirical PM6. Then two modified semiempirical methods, PM6-DH+ and PM6-DH2, were also applied to these clusters. PM6-DH+ was applied to examine both the energies and geometries of the clusters. However, due to the limitation on geometry optimizations of PM6-DH2[48], we examined the interaction energies of PM6-DH2 on the basis of the optimized geometries from PM6. All the semiempirical calculations were performed using the MOPAC2009 (SCC, Colorado Springs, CO, USA) program package[49].

On the other hand, we have examined the above-obtained semiempirical results by applying ab initio MO and DFT calculations for such ionic hydrogen-bonded systems. The ab initio calculations were performed with Hartree-Fock (HF) and second-order Møller-Plesset perturbation theory (MP2)[50] including electron correlations. The DFT scheme utilized Becke's three parameter hybrid method using the Lee-Yang-Parr correlation function (B3LYP)[51, 52]. The aug-cc-pVDZ basis set was used for the comparison of all methods. Of the above-mentioned theoretical levels, MP2 is usually considered to be the most reliable method, while it requires the most computational cost. The computations were performed using the Gaussian 09 (Gaussian, Inc., Wallingford, CT, USA) program package[53].

Results and discussion

The schematic illustrations of optimized structures for X−(H2O)1–4 are shown in Figure 1. Here, X refers to F or Cl. These topological structures obtained from semiempirical calculations are consistent with those of the ab initio calculations. In addition, the global minimum structures of fluoride clusters, which are optimized by only semiempirical methods, are shown in Figure 2. The binding energies for fluoride and chloride anion water clusters are listed in Tables 1 and2, respectively. The binding energy is defined as follows:

ΔE = E x − + n E H 2 O − E x − H 2 O n ,
(1)

where EZ refers to the optimized energies of species Z = X−, H2O, and X−(H2O) n , respectively. The selected bond lengths and bond angles, which represent the ionic hydrogen bond strengths of X−(H2O)1–4 clusters, are shown in detail Tables S1 and S2 in Additional file1 for X = F and X = Cl, respectively. In addition, the values and directions of the imaginary frequencies of the corresponding transition state structures are shown in Table S3 and Figure S1, respectively, in Additional file1.

Figure 1
figure 1

Optimized geometries of fluoride/chloride anion water clusters at MP2/aug-cc-pVDZ level. The colors of green, red, and white represent fluoride/chloride, oxygen, and hydrogen atoms, respectively.

Full size image
Figure 2
figure 2

Optimized global minima of fluoride anion water clusters at PM6-DH+ level. The F…H and H…O distances in each cluster are shown in black and pink color, respectively. The optimized F…H and H…O distances in isolated HF and H2O molecules are 0.966 and 0.949 Å at the same calculation level, respectively.

Full size image
Table 1 Binding energies (kcal/mol) of fluoride anion water clusters
Full size table
Table 2 Binding energies (kcal/mol) of chloride anion water clusters
Full size table

We would focus on the performance of the semiempirical methods on the optimized geometries and binding energies of X−(H2O) n clusters with respect to the increasing water number through comparison with the ab initio and DFT calculations.

Fluoride anion water clusters

In the F−(H2O) cluster, the stable structure takes the Cs symmetry, and the transition state takes the C2v symmetry for all calculation methods. The potential barrier height between the transition state and stable structure is 7.7 kcal/mol at PM6 level with zero-point vibrational energy (ZPE) correction, which is exactly the same as that at MP2 level. We note here that PM6-DH+ and PM6-DH2 show the same stabilization energies as PM6. The ionic hydrogen-bonded structures of F−(H2O) n are exhibited through r(F…Hnear), r(F…O), and θ(F…H-O). In comparison with the C2v structure, the semiempirical results show that the r(F…Hnear) of Cs structure is shorter and θ(F…H-O) is much larger, which agree well with other calculation levels (see Table S1 in Additional file1 for details). The results indicate that the ionic hydrogen bond strength in the stable structure is stronger than that in the transition state structure. However, r(F…O) of the Cs structure is slightly longer for the semiempirical methods, which is just on the contrary at higher calculation levels. Here, we have to mention that the r (F…Hnear) predicted by PM6-DH+ (1.165 Å) is much shorter than that of MP2 (1.413 Å). The optimized F…H and H…O distances in isolated HF and H2O molecules are 0.966 and 0.949 Å at PM6-DH+ level, respectively. Thus, the global minimum would be close to the HF(OH−) structure (see Figure 2).

In the case of F−(H2O)2, ab initio calculations show that stable and transition state structures take the C 2 and C2v symmetries, respectively. The potential barrier height between these two structures is 0.2 kcal/mol at the MP2 level with ZPE correction. The energy barrier is quite small in comparison with the F−(H2O) cluster, which indicates the frequent changes of cluster geometry under certain thermal fluctuations. The semiempirical methods can also provide such small energy difference (not more than 0.1 kcal/mol) between C2 and C2v structures. However, the preferred global minimum structure by semiempirical methods takes the C1 symmetry, rather than the C2 symmetry by ab initio calculations. The results indicate that the water-water interaction in F−(H2O)2 is not well estimated by semiempirical methods. In addition, 2(C2v′) of the F−(H2O)2 structure is shown to be a saddle point with the smallest binding energies by semiempirical as well as ab initio and DFT calculations. The 2(C2h) structure was optimized as a saddle point by MP2/6-311++G** and MP2/TZ (2df, pd)++ in[24]; however, using a larger basis set of aug-cc-pVDZ in this work, it is shown to be a transition state at the MP2 level and to be a local minimum at the HF and B3LYP levels. The semiempirical methods support that it is a saddle point structure. On the other hand, the geometry differences of the C2, C2v, and C2h structures, such as r(F…Hnear), r(F…O), and θ(F…H-O), are very close for both semiempirical and ab initio calculations. The results imply that the ionic hydrogen bond strengths in these geometries are very close. Similar to the case of n = 1, short r(F…Hnear) is predicted by PM6-DH+ (1.087 Å), the C1 structure is probably HF(OH−)H2O.

In the case of F−(H2O)3, ab initio calculations show that stable structures take the C3 and Cs symmetries, and the C3 structure is more stable. It is interesting from the ab initio results that when ZPE corrections are included, the energy of the transition state with the C3h symmetry could be even lower than that with the C3 symmetry. Unfortunately, this important indication could not be well supported by PM6 or PM6-DH+. Instead, the most stable C1 structure of F−(H2O)3 optimized by semiempirical methods takes a planar structure similar to C3h. What differs with the C3h symmetry is that the fluoride ions are located outside rather than inside the cluster. The semiempirical results show that the binding energy of the Cs structure is higher than that of the C3 structure, which may be the result of the overestimation of water-water interactions. On the other hand, in agreement with the ab initio results, semiempirical methods show that r(F…Hnear) of the C s structure is the shortest in the F−(H2O)3 clusters, which is followed by C3h and C3 symmetries. The largest r(F…O) also appears in the C3h structure at PM6 level, while other higher calculation levels favor the C3 structure. The largest θ(F…H-O) appears in the Cs structure at the PM6 level, while other higher calculation levels prefer the C3h symmetry. The C1 structure, which is optimized through the PM6-DH+ method, has a very short F…H distance (1.063 Å). The results indicate that the structure of HF(OH−)(H2O)2 is probably formed instead of F−(H2O)3. However, the ab initio calculations do not support these results. In addition, both the semiempirical and ab initio calculations show that due to the formation of the second hydration shell in the Cs structure, a shorter r(F…Hnear) and a larger θ(F…H-O) are obtained in comparison with the C3 symmetry. The results indicate that the ionic hydrogen bond strength in the structures with the second hydration shell is stronger than that in single hydrogen shell structures, when the water molecule number included in the cluster is the same. More interestingly, the ionic hydrogen bond strength in the Cs structure of F−(H2O)3 is stronger than that in the C2 structure of F−(H2O)2. The result indicates that the ionic hydrogen bond strength is strengthened when additional water molecules are connected via water-water interaction.

In the case of F−(H2O)4, PM6 and PM6-DH+ show that the most stable structure takes the C1 symmetry, and F− is located at the surface of the clusters. However, in the global minimum optimized by HF and B3LYP, F− is in the middle of the cluster with the C1 symmetry. Besides, MP2 favors the Cs symmetry with a 0.4-kcal/mol binding energy higher than the C1 symmetry. The binding energy of the C1 structure by PM6-DH+ is closer to MP2 in comparison with PM6. However, the binding energy differences between C1 and other symmetries by PM6 or PM6-DH+ are greater than 6 kcal/mol, as compared with the differences of less than 2 kcal/mol in the MP2 level. Since the energy differences between the stable configurations are quite small at the MP2 level, the cluster tends to be delocalized and has the coexistence of several configurations. However, the semiempirical results indicate a more localized structure of the C1 symmetry. In this structure, the F…H distance (1.040 Å) at the PM6-DH+ level is close to the isolated HF bond distance (0.966 Å). Thus, the structure of HF(OH−)(H2O)3 is probably formed. From the results of ab initio and DFT calculations, the longest r(F…Hnear) and r(F…O) appear in the transition state structure with the Ci symmetry. The largest and smallest θ(F…H-O) appear in the C4h symmetry and C1 structure, respectively. However, the semiempirical results show that the longest r(F…Hnear) and r(F…O) appear in the C1 and Cs structures, respectively, while the largest and smallest θ(F…H-O) both appear in the C2 structure. Similar to the case of n = 3, the ab initio and DFT calculations agree that the shortest r(F…Hnear) and r(F…O) appear in the Cs structure, which indicates that the strongest hydrogen bond exits in a second hydration shell structure.

As the water molecules in the F−(H2O) n clusters increase from 1 to 4, the semiempirical methods represent reasonable binding energies of the global minima. However, PM6-DH+ and PM6 probably favor an optimized global minimum structure of HF(OH−)(H2O)n−1. As n increases from 1 to 4, the F…H distance in the global minimum structure is getting closer to the HF bond distance (see Figure 2). Similar structures have been optimized by the EFP method[26] as local minima, while they are not global minima on the potential energy surface. However, the HF(OH−)(H2O)n−1 structures are not optimized by ab initio or DFT calculations. As shown in Figure 3, the ab initio and DFT calculations show that as n increases from 1 to 4, r(F…Hnear) and r(F…O) distances of the stable configurations with only one hydration shell increase, and the θ(F…H-O) angle decreases. The results imply that the ionic hydrogen bond strength decreases with the increase of cluster size. However, the opposite trend emerged for the semiempirical methods. The semiempirical methods are sufficient to show that the ionic hydrogen bond strength is enhanced when the second hydration shell forms, compared with the structure without the additional water molecules. However, the semiempirical methods are not accurate enough to study the structures of small fluoride ion water clusters, where strong ionic hydrogen bond exists.

Figure 3
figure 3

Comparison of semiempirical and ab initio /DFT calculations on ionic hydrogen-bonded structure of fluoride anion water clusters. The horizontal axis represents the water numbers and the vertical axis represents (a) F…H distance and (b) F…H…O angle.

Full size image

Chloride anion water clusters

In the Cl−(H2O) cluster, all calculations show that the most stable structure takes the Cs symmetry. The semiempirical methods also show the C2v symmetry to be a local minimum, though it is supposed to be a transition state by ab initio and DFT calculations. The energy difference between the Cs and C2v symmetries is 0.4 kcal/mol at the PM6-DH2 level with ZPE correction, which is slightly lower than 1.3 kcal/mol at the MP2 level. The small energy barrier indicates that the cluster geometry is very flexible. In comparison with the C2v structure, the Cs structure shows shorter r(Cl…Hnear) and r(Cl…O) and a much larger θ(Cl…H-O) by all calculation levels (see Table S2 in Additional file1 for details). The results indicate that the ionic hydrogen bond of Cl−(H2O) is stronger in the asymmetric Cs structure than that of the symmetric C2v structure.

In the case of Cl−(H2O)2, semiempirical methods show that the stable and transition state structures take the C1 and C2 symmetries, respectively. The results are compatible with the ab initio and DFT calculations. The potential barrier height is 1.7 kcal/mol at the PM6-DH2 level with ZPE correction, which is higher than 0.1 kcal/mol for the MP2 level. Such small energy difference indicates that the delocalized cluster geometry varies easily under certain thermal fluctuations. The 2(C2v) and 2(C2h) for Cl−(H2O)2 structures are not obtained by semiempirical calculations. In addition, the 2(C2v′) structure is optimized as a local minimum by semiempirical calculations, while it is a saddle point structure for ab initio and DFT methods. Despite the underestimation of θ(Cl…H-O) for C2 structure, the semiempirical methods could present the same trend as the ab initio and DFT methods, i.e., the bond lengths of r(Cl…Hnear) and r(Cl…O) in the C2 structure are close to the average of two asymmetric bond lengths in the C1 structure. The results imply that the ionic hydrogen bond strength of the transition state C2 structure is between the two ionic hydrogen bond strengths of the local minimum C1 structure.

In the case of Cl−(H2O)3, all calculations show that the most stable structure takes the C3 symmetry. By including ZPE correction, it is found that PM6-DH+ and PM6-DH2 could present binding energies very close to those of MP2 as well as experimental results. From the geometry point of view, the semiempirical PM6 and PM6-DH+ show that the shortest r(Cl…O) in the Cl−(H2O)3 clusters is in the Cs structure, which is consistent with the ab initio and DFT methods. The higher calculations also show that the shortest r(Cl…Hnear) and largest θ(Cl…H-O) are in the Cs structure, while it takes the C3 symmetry by semiempirical method. The higher calculation results indicate that the ionic bond strength is weakened in the sequence of Cs, C3h, and C3 structures. However, the semiempirical method PM6 or PM6-DH+ do not show a clear tendency.

In the case of Cl−(H2O)4, the semiempirical methods show that the most stable structure takes the C4 structure, which is consistent with the higher calculation results without ZPE corrections. When considering ZPE corrections, B3LYP shows that the Cs structure becomes the most stable, while HF and MP2 still favor the C4 structure. The PM6 and PM6-DH+ methods show that the shortest r(Cl…Hnear) and r(Cl…O) and largest θ(Cl…H-O) come up in the second hydration shell's Cs′ structure. The results are consistent with the other calculations and indicate that the strongest ionic hydrogen bond exits in the second hydration shell structure of Cl−(H2O)4. The longest r(Cl…Hnear) and r(Cl…O) appear in the Cs structure for semiempirical methods and in the C2 structure for higher calculations, respectively. The smallest θ(Cl…H-O) is obtained in the stable structure with C2 symmetry for all calculations. We note that the Cs′ structure of Cl−(H2O)4 is reported here for the first time. It is a different configuration from the Cs structure found previously. In the Cs′ structure, it is interesting that the strength of the single ionic hydrogen bond is stronger than that of the other two symmetric ionic hydrogen bonds, while it is contrary to the case of the Cs structure. This new local minimum obtained from PM6 and PM6-DH+ methods could be verified through ab initio and DFT results.

The binding energies and topological characteristics of the global minimum structures presented by PM6 and PM6-DH+ are in good agreement with those of other ab initio and DFT results. Thus, it is reasonable to apply semiempirical methods on Cl−(H2O) n clusters for on-the-fly MD studies from both energetic and geometric considerations. Due to the complex coupling of ion-water and water-water interactions, the ionic hydrogen bond strength is strengthened by the formation of the second hydrogen shell. As the number of water molecules in the Cl−(H2O) n cluster increases from 1 to 4, the tendency of weakening the ionic hydrogen bond strength with the increasing number of water molecules could be obtained both from semiempirical and ab initio and DFT calculations (see Figure 4). We note that the mean values of r(Cl…H) and θ(Cl…H-O) are adopted for cases of structures with more than one ionic hydrogen bonds. Especially for Cl…H distances, PM6-DH+ shows better results than PM6 as n increases in comparison with MP2.

Figure 4
figure 4

Comparison of semiempirical and ab initio /DFT calculations on ionic hydrogen-bonded structure of chloride anion water clusters. The horizontal axis represents the water numbers and the vertical axis represents (a) Cl…H distance and (b) Cl…H-O angle.

Full size image

Comparison of fluoride and chloride clusters

For both fluoride and chloride clusters, PM6-DH+ could substantially improve the binding energies of PM6 in comparison with experiments as well as MP2 results. The semiempirical methods are sufficient to show the well-known results that the binding energy of fluoride cluster is larger than that of chloride with similar geometries.

From the geometric point of view, the optimized geometries by PM6 and by PM6-DH+ show little distinction for both X = F and X = Cl. The results show that the HF(OH)−(H2O)n−1 structures are probable to be global minima for X = F. However, the ab initio and DFT results all show that F prefers to be located inside the small F−(H2O) n clusters. The strong ion-water interactions are dominated in these clusters, while they are not well estimated by PM6 or PM6-DH+. On the other hand, the semiempirical results show that it is quite ubiquitous to find the cooperation of ion-water and water-water interactions in Cl−(H2O) n clusters. The results are compatible with those of the ab initio and DFT calculations and that the global minima take the same topological characteristics.

Conclusions

We have performed semiempirical calculations to analyze the stabilization energies and ionic hydrogen-bonded structures of X−(H2O) n (X = F, Cl; n = 1–4). The results were also compared with the high-level ab initio and DFT calculations.

For fluoride and chloride clusters, the stabilization energies by PM-DH+ are better than those by PM6 in comparison with those by MP2 as well as of the experimental results. On the other hand, PM6-DH+ presents tiny differences in optimized geometries from PM6. In the case of X = F, the optimized global minima are more probable to take the HF(OH)−(H2O)n−1 structures, which are not supported by the ab initio or DFT calculations. The results may be due to the fact that the strong ion-water interactions are not well estimated by PM6 or PM6-DH+ methods. However, the optimized global minima for X = Cl present the same topological characteristics as MP2, which take Cs, C1, C3, and C4 symmetries for n = 1 to 4, respectively.

As the number of water molecules increases, the ionic hydrogen bond strength becomes weaker for X−(H2O) n clusters. The ionic hydrogen bond strength of fluoride cluster is stronger than that of the chloride cluster with similar geometry. It is interesting to find that the emergence of the second hydration shell enhances the ionic hydrogen bond strengths, as compared with the structure without the additional water molecular. Meanwhile, we have reported a new stable structure of Cl−(H2O)4 with Cs′ symmetry with the second hydration shell. The aforementioned conclusions obtained from semiempirical methods have been verified by ab initio and DFT calculations. The study of on-the-fly semiempirical MD or path integral MD simulation with PM6 or PM6-DH+ potentials is in progress.