Hydrated Sodium Ion Clusters [Na+(H2O)n (n = 1–6)]: An ab initio Study on Structures and Non-covalent Interaction

Structural, thermodynamic, and vibrational characteristics of water clusters up to six water molecules incorporating a single sodium ion [Na+(H2O)n (n = 1–6)] are calculated using a comprehensive genetic algorithm combined with density functional theory on global search, followed by high-level ab initio calculation. For n ≥ 4, the coordinated water molecules number for the global minimum of clusters is 4 and the outer water molecules connecting with coordinated water molecules by hydrogen bonds. The charge analysis reveals the electron transfer between sodium ions and water molecules, providing an insight into the variations of properties of O–H bonds in clusters. Moreover, the simulated infrared (IR) spectra with anharmonic correction are in good agreement with the experimental results. The O–H stretching vibration frequencies show redshifts comparing with a free water molecule, which is attributed to the non-covalent interactions, including the ion–water interaction, and hydrogen bonds. Our results exhibit the comprehensive geometries, energies, charge, and anharmonic vibrational properties of Na+(H2O)n (n = 1–6), and reveal a deeper insight of non-covalent interactions.

In spite of many works having been conducted for Na + (H 2 O) 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 + (H 2 O) 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 + (H 2 O) 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 DMol 3 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., 2017Shi et al., , 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).

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 + (H 2 O) 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).
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: D 2d , C s , C 2 , and C 1 . 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/augcc-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 C 2 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 S 6 symmetry to D 3 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 protondonor 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.  The * in sodium ion orbitals represent that the orbitals are empty. The subscripts of the atoms correspond to the serial numbers in Figure 2.

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 + (H 2 O) 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 + (H 2 O) 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. FIGURE 5 | Anharmonic correctional IR spectra of 3+1+0, 4+0+0 calculated at MP2/aug-cc-pVDZ level of theory and experimental spectrum.
Frontiers in Chemistry | www.frontiersin.org 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 + (H 2 O) n (n = 1-6). The CGA search and geometrical optimization for the cluster isomers provide accurate stable structures of Na + (H 2 O) 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 + (H 2 O) 4 , while 3+1+1 and 4+1+0 should FIGURE 7 | Anharmonic correctional IR spectra of eight isomers for n = 6 shown in Figure 1 calculated at MP2/aug-cc-pVDZ level of theory.
be concomitant for Na + (H 2 O) 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.