Introduction
Hydrated ion clusters widely exist in oceans and living organisms, especially hydrated sodium ion clusters, which are important in the control of blood pressure, cell permeability, neuronal activity, and other somatic functions (Jensen, 1992; Feller et al., 1994; Pohl et al., 2013). Understanding the behavior of hydrated sodium ion clusters is helpful to uncover the mechanism of some key biochemical reactions (Mano and Driscoll, 1999; Snyder, 2002; Dudev and Lim, 2010; Payandeh et al., 2011). A number of experimental (Dzidic and Kebarle, 1970; Tang and Castleman, 1972; Schulz et al., 1986, 1988; Blades et al., 1990; Hertel et al., 1991; Patwari and Lisy, 2003; Vaden et al., 2004; Mancinelli et al., 2007) and theoretical (Perez et al., 1983; Arbman et al., 1985; Lybrand and Kollman, 1985; Cieplak et al., 1987; Probst, 1987; Bauschlicher et al., 1991; Dang et al., 1991; Perera and Berkowitz, 1991; Hashimoto and Morokuma, 1994; Glendening and Feller, 1995; Kim et al., 1995; Ramaniah et al., 1998; Carrillo-Tripp et al., 2003; Lee et al., 2004; Rao et al., 2008; Neela et al., 2012; Biring et al., 2013; Dinh et al., 2014; Soniat et al., 2015; Fifen and Agmon, 2016) studies on hydrated sodium ion clusters have been reported, particularly on the global minima. A global minimum can be determined by obtaining the stabilization energies of isomers. In experiments, comparing enthalpies is the most direct method to obtain thermodynamic information to deduce stabilization energies. With a high-pressure mass spectrometer containing a thermionic alkali ion source, Dzidic and Kebarle reported the enthalpies and entropies of hydrated sodium ion clusters for n = 1-6 in gas phase (Dzidic and Kebarle, 1970). With the hydration number increasing, the binding energy per water molecule decreases. Glendening and Feller calculated the stabilization energies and stabilization enthalpies of Na+(H2O)n (n = 1-6) at various levels of theory (Glendening and Feller, 1995), in which the RHF and MP2 levels with the 6-31+G* basis set reproduced the experimental values obtained by Dzidic and Kebarle well (Dzidic and Kebarle, 1970).
For the structures of hydrated sodium ion clusters, the coordination number is of particular appeal to studies in different phases. In liquid water, the coordination number of a sodium ion is about 5.5 ± 0.5 based on molecular dynamics simulation (Mancinelli et al., 2007; Megyes et al., 2008; Bankura et al., 2013, 2014; Lev et al., 2013; Galib et al., 2017; Liu et al., 2019). However, in gas-phase clusters, the coordination number is 4 from ab initio calculations at 0 K (Kim et al., 1995; Neela et al., 2012; Soniat et al., 2015; Fifen and Agmon, 2016). The structures of Na+(H2O)n (n = 1-4), all the water molecules surrounding the sodium ion, were firstly reported by Bauschlicher et al. from ab initio calculation (Bauschlicher et al., 1991). For Na+(H2O)5, 4+1+0 (the structures of isomers are presented in the form of n1+n2+n3, where n1, n2, and n3 are the numbers of water molecules in the first, second, and third solvation shells, respectively) is supported as the global minimum by most ab initio calculations at 0 K (Hashimoto and Morokuma, 1994; Kim et al., 1995; Lee et al., 2004; Rao et al., 2008; Neela et al., 2012; Biring et al., 2013; Soniat et al., 2015; Fifen and Agmon, 2016), and 5+0+0 is deemed to be concomitant with 4+1+0 at 298 K (Kim et al., 1995; Fifen and Agmon, 2016). For n = 6, at 0 K, several recent ab initio calculations stated that 4+2+0 (with D2d symmertry) is the global minimum (Lee et al., 2004; Rao et al., 2008; Biring et al., 2013; Soniat et al., 2015; Fifen and Agmon, 2016), which was proposed by Lybrand and Kollman (1985) based on RWK2 potential (Reimers et al., 1982). However, Neela et al. sustained 5+1+0 to be the global minimum for n = 6 calculated at MP2/cc-pVTZ level of theory (Neela et al., 2012). At room temperature, Kim et al. found that 5+1+0 possesses better stability than 4+2+0 by HF/TZ2P (Kim et al., 1995). Differently, Fifen and Agmon indicated that 4+1+1 is dominant, calculated at MP2/6-31++G(d,p) (Fifen and Agmon, 2016).
The infrared (IR) spectra are available in distinguishing the cluster isomers. For hydrated sodium ion clusters, the feature peaks of O-H stretching mode could accurately provide the structure information of clusters (Huang and Miller, 1989; Patwari and Lisy, 2003; Vaden et al., 2004; Miller and Lisy, 2008a,b; Ke et al., 2015), in which the non-covalent interactions, including ion-water interaction and hydrogen bond, weaken the O-H bonds, causing redshifts for O-H stretching vibration modes and producing different feature peaks for different structures (Muller-Dethlefs and Hobza, 2000; Vaden et al., 2002, 2004, 2006; Kozmutza et al., 2003; Bush et al., 2008; Miller and Lisy, 2008a), Using a custom-built, triple-quadrupole mass spectrometer as well as ab initio calculations, Lisy et al. reported the IR spectra of Na+(H2O)n (n = 2-5) and Na+(H2O)nAr (n = 2-5) (Miller and Lisy, 2008a,b; Ke et al., 2015). For n = 4, 3+1+0 is the stable structure with bent hydrogen bonds (Miller and Lisy, 2008a). Recently, through straightforward IR spectra for n = 5, they speculated that 4+1+0 and 3+1+1 could be concomitant at 75 K (Ke et al., 2015).
In spite of many works having been conducted for Na+(H2O)n (n = 1-6), the global minima of n = 4-6 remains unclear. Moreover, the non-covalent interaction and electron transfer in hydrated sodium ion clusters have not been discussed in detail, which can elucidate the principle of the shifts of O-H stretching frequency. In this paper, the comprehensive genetic algorithm combined MP2 method is used to determine the global minima of Na+(H2O)n (n = 1-6) and simulate their anharmonic vibrational frequencies. Furthermore, charge transfer inside the clusters through NBO analysis and charge density difference are contained, aiming to reveal the principle of non-covalent interactions effecting on the O-H bonds.
Methods
In this work, some of the structures of the Na+(H2O)n (n = 1-6) clusters are adopted from previous literatures (Bauschlicher et al., 1991; Glendening and Feller, 1995; Ke et al., 2015; Fifen and Agmon, 2016). To obtain more isomers for n = 4-6, a global search with the comprehensive genetic algorithm (CGA, Zhao et al., 2016) combined with DMol3 program (Delley, 2000) based on DFT was executed. The CGA method is described in our previous review in detail (Zhao et al., 2016). For each cluster size with n ≥ 4, we took 10 independent global searches, and for each search, we maintained mating and mutation operations on a population of eight members of up to 3000 GA iterations. Since the Becke-Lee-Yang-Parr (BLYP, Becke, 1988; Lee et al., 1988) functional would provide similar relative energies to MP2 (Møller and Plesset, 1934) method (see Table S1), the generalized gradient approximation (GGA) with the BLYP functional and p- and d- polarization functions (DNP) basis sets were employed to optimize the clusters’ isomers in CGA search without symmetry constraint. Considering the calculation cost, zero-point energy (ZPE) was not contained in global search.
Our previous work proved that MP2 is a reasonable method to obtain the energies and properties of small hydrogen-bonded systems (Liu et al., 2013; Shi et al., 2017, 2018). Since the geometrical optimization at augmented correlation-consistent polarized valence double-zeta (aug-cc-pVDZ, Dunning, 1989; Kendall et al., 1992) and aug-cc-pVTZ provide almost the same structures (see Table S2), MP2/aug-cc-pVDZ method was utilized to optimize the isomer structures. Single-point energies of these clusters were computed at MP2/aug-cc-pVQZ and MP2/aug-cc-pVDZ levels.
Within harmonic approximation, MP2 calculation usually overestimates the frequencies relative to the experiment, especially for the high frequencies in IR spectra, and may leave out some peaks. Hence, we calculated the IR spectra with anharmonic correction at MP2/aug-cc-pVDZ level at 298 K via second-order vibrational perturbation theory (VPT2, Barone, 2005; Barone et al., 2010), as well as to obtain ZPE and thermal correction at 298 K.
For visualizing the bonding strength between the two atoms intuitively, NBO (Carpenter and Weinhold, 1988; Reed et al., 1988) was calculated at MP2/aug-cc-pVQZ level, as well as to obtain the Wiberg bond order (Wiberg, 1968). All the calculations aforementioned were performed in the Gaussian 09 package (Frisch et al., 2013).
Charge density differences of 1+0+0, 2+0+0, 3+0+0, 3+1+0, 4+0+0, and 4+1+0 were calculated using GGA and Perdew-Burke-Ernzerhof (PBE, Perdew et al., 1996) functional, the projector-augmented wave potentials (Blochl, 1994) with an energy cutoff of 500 eV, as implemented in Vienna Ab-initio Simulation Package (VASP, Kresse and Furthmuller, 1996). Only Γ point is k-point with a vacuum layer of over 15 Å was employed in our calculation. The charge density difference is given by:
where ρcluster, ρNa+, and ρH2O are the charge density of entire hydrate cluster, sodium ion and all the water molecules, respectively.
Results and Discussion Structures
We re-optimized all the isomer clusters obtained from CGA search using the MP2/aug-cc-pVDZ method. The optimized structures and symmetries of Na+(H2O)n (n = 1-6) are present in Figure 1. Due to the high computational cost, we used the total energies computed at the MP2/aug-cc-pVQZ//MP2/aug-cc-pVDZ+ZPE level to rank the energy order of all the isomers. Table 1 lists the relative energies at 0 K and 298 K calculated at MP2/aug-cc-pVQZ, MP2/aug-cc-pVDZ, and BLYP/DNP levels, respectively.
For n = 1-3, the global minima, i.e., 1+0+0, 2+0+0, and 3+0+0, all the water molecules surround the sodium ions and locate at equivalent positions without hydrogen bonds, which are similar to those in previous reports (Hashimoto and Morokuma, 1994; Kim et al., 1995; Rao et al., 2008; Neela et al., 2012; Soniat et al., 2015; Fifen and Agmon, 2016).
For n = 4, CGA has located 3+1+0 and 4+0+0 isomers. Among them, 3+1+0 was the better result in all ten CGA global searches, which possesses lower energy than 4+0+0 without ZPE correction (see Table S1). In most previous reports by ab initio calculations, 4+0+0 was argued to be the most stable of the structures (Hashimoto and Morokuma, 1994; Kim et al., 1995; Ramaniah et al., 1998; Lee et al., 2004; Rao et al., 2008; Kamarchik et al., 2011; Neela et al., 2012; Fifen and Agmon, 2016), whereas Miller and Lisy reported that the IR spectrum of Na+(H2O)4 is similar to that of 3+1+0 at 300 K (Miller and Lisy, 2008a). From Table 1, at 0 K, our MP2 calculation manifests that 4+0+0 is lower in energy by 1.36 kcal/mol, while 4+0+0 and 3+1+0 possess almost equal energies at 298 K. As shown in Figure 1, 4+0+0 has four equivalent water molecules surrounding the sodium ion, while 3+1+0 is evolved by 3+0+0 connecting a water molecule with two coordinated water molecules by hydrogen bonds.
For n = 5, six isomers are contained in our calculation: 4+1+0, 3+1+1, two 3+2+0 structures, and two 5+0+0 structures with different symmetries. 4+1+0 was the best structure in all ten CGA searches, and had the lowest energy at 0 K. 3+1+1, 3+2+0(1), and 5+0+0(1) possess the relative energies of 2.00, 2.09, and 2.77 kcal/mol, respectively. Besides, 3+2+0(2) and 5+0+0(2) have equal relative energies of 2.99 kcal/mol. At 298 K, 5+0+0(2) possesses a lower energy than 5+0+0(1), and the energetic order of the other isomers doesn’t change. At the BLYP/DNP level of theory, 3+2+0(1) becomes the second lowest energy structure, and 3+2+0(2) has lower energy than 5+0+0(1) and 5+0+0(2) at 0 K. In Figure 1, it is noteworthy that the global minimum structure 4+1+0 has an extra water molecule located at outer shell of 4+0+0. 3+1+1 structure has a water molecule connecting with the outer water molecule in 3+1+0 via a hydrogen bond. Similarly, 3+2+0(1), with Cs symmetry, has a water molecule connecting the isolated coordinated water molecule in 3+1+0, and 3+2+0(2) has a water molecule located beside a coordinated water molecule with a hydrogen bond of 3+1+0. 5+0+0(1) structure has three water molecules form a water cycle via hydrogen bonds, and the other two coordinated water molecules are isolated opposite the water cycle. Differently, 5+0+0(2) has only one isolated coordinated water molecule, with the other four water molecules constituting a quaternary water cycle.
For n = 6, we found eight isomers, 4+1+1, 6+0+0, two 5+1+0 structures, and four 4+2+0 structures with different symmetries: D2d, Cs, C2, and C1. In all ten CGA searches, 4+2+0(3) had the best solution with the lowest energy at MP2/aug-cc-pVQZ without ZPE (see Table S1). From Table 1, at 0 K, 4+2+0(1) has the lowest energy, while the relative energies of 4+1+1 and 4+2+0(2) are 1.05 and 1.25 kcal/mol, respectively. The other five isomers possess relative energies of over 2 kcal/mol. At 298 K, 4+2+0(3) becomes the second lowest energy structure rather than 4+1+1, with the relative energy of only 0.76 kcal/mol. Four isomers with four coordinated water molecules have lower stabilization energies both at 0 K and 298 K, indicating that four coordination is more favorable for n = 6. Compared to MP2/aug-cc-pVQZ, the calculations at BLYP/DNP shows that 4+2+0(4) has lower energy than 4+2+0(3) and 5+1+0(1). Besides, 6+0+0 possesses the highest relative energy of 6.56 kcal/mol, which is obviously higher than the total energy interval at MP2/aug-cc-pVQZ (4.32 kcal/mol). Combining with the relative energies of n = 5, BLYP/DNP gives the same global minima and similar energetic order to MP2/aug-cc-pVQZ. However, BLYP/DNP would overestimate the energies of the 5 and 6 coordinated structures, indicating that the CGA search tends to provide the isomers with 3 and 4 coordinated water molecules. Since the previous ab initio calculations show that the coordination number is 4 at 0 K (Kim et al., 1995; Neela et al., 2012; Soniat et al., 2015; Fifen and Agmon, 2016), the CGA search at BLYP/DNP could provide the global minima and other reliable isomers. As seen in Figure 1, 4+2+0(1) is a water molecule via hydrogen bonds connecting with the two coordinated water molecules in 4+1+0. Three coordinated water molecules in 4+2+0(2) connect the two outer water molecules via hydrogen bonds, and the oxygen atoms in these five water molecules locate in a flat with the sodium ion approximatively. 4+2+0(3) with C2 symmetry forms a water cycle via hydrogen bonds between two coordinated water molecules and the two outer water molecules. Besides, 4+2+0(4) with lowest symmetry has a coordinated water molecule without hydrogen bond. The 4+1+1 structure is a water molecule connecting with the outer water molecule in 4+1+0 via a hydrogen bond. In addition, 5+1+0(1) is an extra water molecule located at the outer shell of 5+0+0(1), connecting with the two isolated coordinated water molecules. 5+1+0(2) is just a water molecule connecting with the isolated coordinated water molecule in 5+0+0(2) via a hydrogen bond. 6+0+0 could transform from the perfect S6 symmetry to D3 symmetry, with two water cycles on two sides of the sodium ion, in accordance with previous calculations based on the polarizable electropole model (Perez et al., 1983).
The bond lengths present interesting variation trends as summarized in Table 2. The r¯(Na-O) increases strictly with the increasing of coordination number, indicating the decreasing of average ion-water interaction. For the structures with two water shells, if a coordinated water molecule is the proton-donor in a hydrogen bond system, the r(Na-O) should be shorter. In contrast, if the oxygen atom forms a hydrogen bond, the r(Na-O) should become longer. For the r(O–H)s, each average r(O–H) of water molecules in clusters is longer than the r(O–H) of free water molecules (0.966 Å), which stems from the non-covalent ion-water interaction. Meanwhile, the hydrogen bonds also stretch the O-H bonds and make the water molecules asymmetric.
Charge Analysis
For elucidating the non-covalent interactions in hydrated sodium ion clusters, Figure 2 and Table 3 show the NBO overlapping 3D schematic diagrams and electron transfers of 1+0+0. Figures 2A,B depict the 2s orbital of sodium ion overlaps the O-H anti-bonding orbitals of water molecule, resulting in electron transfer from the sodium ion to the water molecule. In contrast, Figures 2C,D depict the O-H bonding orbitals overlap to the empty orbital of sodium ion, resulting in electron transfer from the water molecule to the sodium ion. From Table 3, the amplitude of E(2) manifests that electron transfer from water molecules, including the electrons in O-H bonding orbitals and the oxygen atom’s lone pair electron orbitals, to sodium ions is larger than that from sodium ions to water molecules, which synergistically weakens and stretches the O-H bonds of 1+0+0.
For revealing the strength of O-H bonds intuitively, the Wiberg bond order in 1+0+0 (0.740) is smaller than that in a free water molecule (0.790), indicating that the sodium ion weakens the O-H bonds, in accordance with the results from NBO analysis.
To show the charge transfer of the whole clusters directly, the charge density difference of 1+0+0, 2+0+0, 3+0+0, 3+1+0, 4+0+0, and 4+1+0 is presented in Figure 3. The electrons from coordinated water molecules assemble at the location between sodium ions and water molecules near the side of oxygen atoms. The electron dissipation mostly happens near the hydrogen atoms, proving that the strengths of O-H bonds become weaker. In addition, the charge density difference of 3+1+0 and 4+1+0 show that the electrons also assemble at the location between the outer water molecules and the sodium ions, indicating that the ion-water interaction also reduces the strength of the O-H bonds in the outer water molecules.
Vibrational Spectra
Vibrational spectrum is an intuitionistic method, providing deeper insight into structure differences and non-covalent interactions (Fan et al., 2019), especially O-H stretching vibration modes for hydrated sodium ion clusters. Therefore, the experimental Na+(H2O)n isomers can be determined by comparing the simulated IR spectra and experimental spectra. Our discussion focuses on the high-frequency region (>3,200 cm−1) in IR spectra, which contains the O-H stretching vibration modes and can generally be measured in experiments (Ke et al., 2015).
At first, Figure 4 shows the IR spectra of 1+0+0, 2+0+0, 3+0+0, and a free water molecule with anharmonic correction. The two O-H stretching vibrational modes are asymmetric (the higher peaks near 3,700 cm−1) and symmetric (the lower peaks near 3,620 cm−1) modes for each structure, and the other peaks are caused by anharmonic correction. Compared to the free water molecule, the asymmetric vibration modes of the three clusters possess redshifts, stemming from the ion-water interactions. The redshifts become smaller with the increasing of coordination number and r¯(Na-O) in Table 2.
For n = 4, the spectra of 4+0+0 and 3+1+0, as well as the experimental spectrum of Na+(H2O)4 at 300 K (Miller and Lisy, 2008a) are given in Figure 5. The two modes of 4+0+0 reproduce the two outstanding peaks of experimental spectrum well. Moreover, the lower peak of the experimental spectrum confirms the small fraction of the existence of 3+1+0. Therefore, 4+0+0 dominates in the experiment at 300 K, and 3+1+0 is concomitant with 4+0+0, in accordance with the almost equal energies at 298 K in Table 1.
Figure 6 shows the IR spectra of six isomers for n = 5 with anharmonic correction, as well as the experimental spectrum (Miller and Lisy, 2008a). Apparently, no single structure could reproduce the experimental spectrum well. Among, the vibrational modes of 3+1+1 are able to correspond three peaks of the experimental spectrum (Miller and Lisy, 2008a), hence 3+1+1 possesses the most possibility of existing in experiment. However, no mode in 3+1+1 could reproduce the experimental peak near 3,560 cm−1, while all the other five structures have vibrational modes near 3,560 cm−1. Combining with the relative energies in Figure 1, 4+1+0, the global minimum, could be the main contributor to the peak at 3,560 cm−1, in accordance with the conclusion in previous reports (Ke et al., 2015; Fifen and Agmon, 2016). Therefore, 4+1+0 and 3+1+1 are concomitant in experiments, which are the two lowest energetic structures at 298 K in Table 1.
For n = 6, Figure 7 shows the IR spectra of all the eight isomers presented in Figure 1. Due to all the coordinated water molecules being equivalent to 4+2+0(1) in Figure 1, only two distinct peaks can be observed. Similar to 3+1+1, 4+1+1 possesses an obvious peak at 3340.5 cm−1 caused by O-H stretching of the water molecule in the second shell. 4+2+0(2) has no peak under 3,500 cm−1 because the hydrogen bonds are not strong enough to make the water molecules distinctly asymmetric, while the three feature peaks, 3371.5, 3395.4, and 3460.2 cm−1, of 4+2+0(3) are generated by the O-H stretching in the water cycle. Because of the hydrogen bonds between the outer molecule and coordinated water molecules in 5+1+0(1), the spectrum has two modes at 3513.3 and 3525.6 cm−1, which can’t be found in 5+0+0(1). Due to no equivalent water molecule in 4+2+0(4), the O-H stretching vibration modes with different frequencies make the spectrum more complex than the others. 6+0+0 has a spectrum similar to that of 5+0+0(2) in Figure 6, corresponding to the similar water cycles in both structures. Unlike 5+0+0(2), 5+1+0(2) has a significant peak at 3390.0 cm−1 which is the O-H stretching mode of the proton-donating coordinated water molecule, indicating that the hydrogen bond between the outer water molecule and the coordinated water molecule is strong in 5+1+0(2).
Conclusion
In this work, we investigate the geometries, energies, charges, and anharmonic vibrational properties of Na+(H2O)n (n = 1-6). The CGA search and geometrical optimization for the cluster isomers provide accurate stable structures of Na+(H2O)n (n = 1-6). At 0 K and 298 K, for n = 1-4, all the water molecules in global minima are coordination water molecules, surrounding the central sodium ions. Meanwhile, 4+1+0 and 4+2+0(1) are the global minima of n = 5 and 6, respectively. Thus, the coordination number of global minima of hydrated sodium ion clusters is 4 for n ≥ 4.
The non-covalent interactions, including ion-water interactions and hydrogen bonds, weaken the O-H bonds, resulting in longer bond lengths, lower bond orders, and redshifts of the O-H stretching mode in IR spectra. The simulated IR spectra with anharmonic correction can reproduce the experimental results well. The results show that 4+0+0 dominates for Na+(H2O)4, while 3+1+1 and 4+1+0 should be concomitant for Na+(H2O)5 in experiments. The present study executes a believable simulation of structures and vibrational spectra, and provides a comprehensive insight into the non-covalent interactions including ion-water interaction and hydrogen bonds of hydrated sodium ion clusters.
Data Availability
The datasets generated for this study are available on request to the corresponding author.
Author Contributions
PW and RS participated in the design and calculated the data of this study. PW, RS, and YS performed the statistical analysis. LT and XH improved the comprehensive genetic algorithm to adapt the system in this study. YS and JZ carried out the study and collected important background information. All authors have read and approved the content of the manuscript.
Funding
This work was supported by the National Magnetic Confinement Fusion Energy Research Project of China (2015GB118001), the National Natural Science Foundation of China (11674046, 11604039), the Science Challenge Project (TZ2016001), the Fundamental Research Funds for the Central Universities of China (DUT16LAB01, DUT16LK16, DUT17GF203, DUT18LK07), and the Supercomputing Center of Dalian University of Technology.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2019.00624/full#supplementary-material
