A Theoretical Study on the Formation, Binding Energy and Monomer Dipole Moment of Small Water Cluster Systems

The formations, binding energies and monomer dipole moments of small water cluster systems (H2O)n with n = 1-12, have been investigated by the density functional calculations using B3LYP/6-311++(2d,2p) level theory. A new method based on reactivity indices from Fukui functions has been introduced to generate the initial structures. Constant adding one by one water molecule from monomer to the cluster systems (hydration reactions) have transformed the cluster shapes by following formation order: linear (n = 2), cyclic planar (n = 3-6) and 3-dimensional (n = 7-12) ones. The average binding energies of small water cluster systems have converged asymptotically to the intermolecular binding of bulk water, concerning the local binding energy fluctuation effects on the average binding energy trend. Based on total electronic energies, zero point energies and optimized transition structures energies (activation energies) analysis, we have predicted that the cyclic planar is the most stable hexamer formation to compete with the cage and the book ones. Considering the monomer dipole moment calculations with regard to the experimental data references, we have discovered new findings which have not clarified yet before, that the smallest piece of ice water is the cyclic tetramer (n =4) and the cage hexamer is the smallest stable liquid formation. [DOI: 10.1380/ejssnt.2009.871]


I. INTRODUCTION
Small water cluster systems have been largely studied for many years in the last decade since their potential applications are still growing up. Water clusters are involved in various areas of great scientific and technological importance, in atmospheric planetary, environment, biochemistry, and novel energy devices such as in the membrane design of Nafion-based Proton Exchange Membrane for Fuel Cell (PEMFC) [1] applications. Recently, the performance of PEMFC is highly expected to be affected by the clustering behavior of water molecules at the molecular level, which control the humidity conditions (wet/dry) inside the Nafion [2,3] membrane.
Wet and dry conditions are strongly related to the physical phase states of matters. In case of very small water cluster systems, it is still hard to identify the phase state of the clusters. A systematic scheme is needed to identify the structures (formation) stabilities, to analyze the intermolecular binding energy trend, and to evaluate the molecular dipole moment values. From our preliminary study, we have investigated the most stable formations and the average binding energies trend of small water clusters as the number of molecules are increased using the * Corresponding author: handoko@dyn.ap.ang.osaka-u.ac.jp density functional theory (DFT). The formation of water clusters have many possible conformations or local minima [4,5] due to their weak intermolecular H-bonding. Thus, finding the most stable formations as hydration products could be arduous, and they strongly depend on the initial orientations of the clusters. One of physical properties that make sense to distinguish between the liquid and solid phase in condensed matter is the stability of structures which are associated with the intermolecular binding energies. However, the property has no charge distribution information. Other identifying of the phase state which include the atomic charge distribution information had been made experimentally by measuring the monomer dipole moments (dipole moment is a vector) of bulk water system [6][7][8]. But it is still uneasy to obtain experimentally the properties on small scale of water cluster systems. Since the hydration reactions also have wide spread applications in chemistry, we have separated our investigation purposes into some objectives.
There are three objectives to be achieved for these studies; the first one is to find a reliable theoretical approach in order to design the initial orientation for intermolecular formation. This could be solved by determining the reactivity index to reduce search of the number of reactive sites between the host cluster and the guest single water molecule in the molecule reactions (hydrations). The second one is to analyze the local binding energy fluctuations in terms of the number of molecules and to observe the behavior trend of the average binding energies for different cluster sizes. The last one is to examine the water clusters physical state (phase) via their monomer dipole moments values with regard to the experimental values [6][7][8].

II. THEORETICAL AND COMPUTATIONAL METHOD
A systematic scheme has been used. We have employed a DFT approach using different basis set based on B3LYP (Becke-3 parameters-Lee-Yang-Parr) exchangecorrelation functional [9] and 6-311++G(2d,2p) [10] to optimize the cluster geometry. One of the reasons in utilizing this functional is that the three parameters of B3LYP had been fitted to gaseous water data.
The reactivity index from Fukui functions [11] have been used to predict the initial structures of the clusters. The basic idea of this prediction is that the chemical potential (µ) differences drive the electron transfer from one system to another. Because the physical meaning of Fukui functions f (r) is implied by its definition as functional derivative [δµ/δ(r)] N : it measures how sensitive a system's chemical potential is to external perturbation at a particular point, then the overlap of left f − (r) and right f + (r) derivative (reactivity indices) will give us prediction directions for the initial structures with more charge-transfer indicating a stable structures approximation. In these studies, the f − and f + reactivity indices have been approached using isosurfaces of the Highest Occupied Molecular Orbitals (HOMOs) and the Lowest Unoccupied Molecular Orbitals (LUMOs).
Calculated charge-transfer populations was estimated using natural bond orbital [12] (NBO) have been also considered to evaluate the trimer (n = 3) binding characters. A saddle point optimizations are also performed to get the activation energy of hexamers in order to complete their stabilities analysis. The average binding energies and local binding energies are determined by the following equation: , and n denote the average binding energy, the absolute (hydration) binding energy when a H 2 O approaches (H 2 O) n−1 cluster, the total energy of cluster system, the total energy of isolated water molecule, and number of molecules, respectively. Equation (1) pertains to the average binding energy needed to dissociate or hydrate all water molecules from/to their cluster. Equation (2) represent the local binding energy needed (in equilibrium) to add or remove only one water molecule to or from the cluster. Both quantities refer to the intermolecular bond strength of the clusters. However, the average binding energies could have significant roles as the references for the fluctuations of the first equation and also their trend will lead to the average intermolecular binding energy of ice water.
The overestimated bond strength of 1b 2 molecular orbital (MO) in the dimer (n = 2) can be reduced by adding more d functions to O atoms and more p functions to H atoms [13]. The obtained functions for O atoms and H atoms will strengthen the covalent characters and reduce the overestimated orbital bond strength indicated by means of the intermolecular binding. The B3LYP has been chosen in order to consider that both the non correlation model (Hartree Fock) [14] and the smaller correlation parameter model (BLYP) could not give higher intermolecular binding. In addition, we also compare the results from HF and BLYP models, and we obtain a lower binding energy using B3LYP in optimized dimer which shows a good agreement with experimental values [15].
The monomer dipole moments are calculated by point of charge approximation that taking account Atomic Polar Tensor (APT) population analysis (the APT charge was obtained from G03 software [16]) which is condensed in atomic charge. To be different from Mulliken population analysis, APT population is based on the electron distribution (electron density) in molecules that is not distorted by an arbitrary choice of the reference basis set as reported by Cioslowski [17] (one of Gaussian03 contributors). The grid of numerical integration is pruned (75,302) grid, having 75 radial shells and 302 angular points per shell. The monomer dipole moments using APT population can be represented by following equation: where µ, q, d, and n denote monomer dipole moment (in Debye), atomic charge (APT population) (in e), atomic position (inÅ), and number of molecules, respectively, and x, y, z are coordinates (projected positions).

III. RESULT AND DISCUSSION
A good initial prediction of the nano structures chemistry based on reactivity index from Fukui functions can be found in the stereo-selection of water dimer structures. In order to form a stable dimer structure, there are two logical approaches that a water molecule could take; one approach is where the planes of the water molecule are perpendicular to each other; and the other, the planes are in parallel. In both cases, the O atom of the approach- I: Average charge-transfer from monomer to others (CTAV), average charge-transfer within monomer/isolated (NCTAV) and average binding energies (∆Ē b ) in dimer and cyclic trimer within the cyclization reaction from dimer to trimer of water clusters. Schematics of the initial and final states of the reaction are also shown. "stg" (surrounded by dashed curves) and "ecl" (surrounded by real curve) denote staggered binding character and eclipsed binding character, respectively. f − and f + denote reactivity indices from Fukui functions.  ing molecule connects with H atom of the other molecule. For this reason, the former (eclipsed formation) will be more stable because H atoms are closer to the f − Fukui functions (reactivity toward an electrophile reagent) than the latter. The stability of the eclipsed dimer has been also examined with another dimer conformer by scanning the energy profile as the H atoms are rotated for about 180 degrees with O atom as the centre of rotation. The other conformer (staggered formation) is found on unstable valley curve (Fig. 1). The maximum energy difference from this rotation is about 2 kcal/mol. It indicates that H-bonding energy has a relative wide-quality based on its structures. Here, we have found that the staggered struc-ture is unstable and the eclipsed is the most stable configuration because the staggered structure has higher repulsive interaction between H atoms than the eclipsed one. Based from these structures, we state the binding tendencies of the molecules as staggered and eclipsed characters, as commonly used terminologies in chemistry. The cyclization is essential in the small water cluster reactions, since the experimental results showed that clusters from 3 to 6 molecules are stable in the cyclic formations [18]. However, using structure selections from the Fukui functions, the initial structure of cyclic trimer can be prepared well as shown at the dimer row in Table I. The figure at trimer row inside Table I   structure of water trimer after optimization. The preparation initial structures have shortened the optimization process, which means that the iterative steps are lesser and this will reduce the computational cost. Our calculations show that staggered binding character has a weaker H-bond and longer O-O distance as compared with the eclipsed ones, as summarized in Table I. These results are also supported by the charge-transfer analysis. We have found that majority of the donor orbitals are from the lone pair orbitals or the HOMOs of each water molecule unit and majority of the acceptor orbitals are from the LUMOs of the other neighbor unit molecule. With the similar ways, we also have prepared the initial structures for water tetramer (n = 4, as shown in Fig. 2) and the others higher number of molecules based on the reactivity index from Fukui functions. The summary steps of the initial structure preparation based on the reactivity index are: (a) Calculating the f − and f + isosurfaces of the host and guest system, (b) determining position of the biggest isosurfaces of the f − and f + , (c) moving the guest molecule close to the host system with respect to the great overlapping between the f − and f + , and (d) if possible, avoiding both the f − and f − and also f + and f + isosurfaces meeting. Figure 3 shows the optimized cluster structures from monomer (n = 1) to dodecamer (n = 12) and Fig. 4 manifests the average binding energies trend and local binding energies fluctuations. For higher number of clusters, it can be observed that the former planar formations have been transformed to 3-dimensional formations. The shape transformation results also agreed with other theoretical studies that was shown based on stability energy, the cyclic hexamer was changed to prism34 [19]. However, the number of water molecules in that ref. [17] (n = 1-10) is inadequate to predict the average binding energy trend with regard to the local fluctuation effects. For cyclic formations with n = 3 and n = 4, the additions of eclipsed and staggered characters stabilize significantly their binding energies. However, as the number of molecules increases (n = 5 and n = 6), the diameters of planar forma- tions lengthen and as a consequence the inter-molecular separation is also increased; hence the values of average binding energies is decreased.
The water hexamer (n = 6) has two interesting points to be considered. First, by monomer dipole moment probing wherein current studies show that the cluster is supposed to be in the liquid phase [20], a transition from gas to solid phase. Second, the most stable formation of hexamer is still under debate since the current studies have found that based from the theoretical binding energies [21] the book/cyclic formation is more stable as compared with the cage formation. This finding is ambiguous with the fact that very low temperature experimental [22] investigation could probe the cyclic formation beside the previous cage structures finding on hexamer cluster but not for the book formation. Our calculations also show that the cyclic formation is the most stable structure followed by the cage and then the book formation. Zero point energies calculations (cyclic = −287835.44 kcal/mol, book = −287833.90 kcal/mol, and cage = −287832.55 kcal/mol) have changed the order of stabilities as cyclic>book>cage. Nevertheless, after considering the electronic energy, the zero point energies and the activation energies (an energy barrier) through the optimized transitional states as shown in Fig. 5, we believe that the cyclic formation is more favored and stable than the book and the cage ones. In this case, our calculation results are in agreement with the experimental findings which were also reported that the cyclic hexamer is formed by cyclization process and that rapid quenching inhibits (it suppose to have an energy barrier) its rearrangement to the cage structures.
The symmetry (well-ordered) of the cluster is broken in heptamer (n = 7) structures and the cyclic formation from hexamer does not expand, but instead to re-arrange into 3-dimensional formation. The 3-D formation is continued in the optimized octamer (n = 8) formation but its local binding energy has significantly decreased as shown in Fig. 4. It happens because besides having high symmetrical (ordered) structures, the cubic-like formation of the octamer could reduce the intermolecular separations.
In case of the nonamer (n = 9), the addition of one water molecule forms two H-bonds but breaks one of the other H-bond. It has responsibility for increasing the local binding energy. The decamer (n = 10) has two most stable structures: prism and fish formations. Both formations have high ordered structures and tend to lower (stabilize) the local binding energy values. The prism-like formation of dodecamer (n = 12) looks like a two-layer cyclic hexamer structures but the two-cubic-like formation has the most stable one, as opposed to the previous hexamer stability examination. Our explanation is that the existences of additional four weak H-bonds at the middle of the two-cubic-like formation prohibit the two-layer book formation. Therefore, the higher symmetrical two-cubiclike formation is favored. From the observed fluctuations, it can be deduced that odd number of clusters have lower absolute binding energies than the even number ones. The fluctuation energies on higher number of molecules cannot be easily estimated from the local binding energy graph. However, we can predict the average binding energy trend (Fig. 4) that for n ≥12 tend to asymptotically converged at intermolecular average binding value of ice cluster. The reason is that, by constant increasing one molecule at a time would give smaller contribution at higher number of molecules. Then, the fluctuation effects on the average binding energies are reduced.
The other essential property of small water cluster systems is the monomer dipole moment. Experimentally [20], one way to distinguish whether water in gas, liquid or ice phase, is measured by the monomer dipole moment property. The dashed line in Fig. 6 shows the monomer dipole moments level for gas phase water is around 1.855 Debye, for the liquid phase is around 2.45 Debye and for solid phase is around 2.6 Debye. The monomer dipole moments have been calculated using Eq. (3) and fitted by dipole moment value of gaseous water. Based on these properties, the tetramer, pentamer, and hexamer are close to the solid phase. However, the cage hexamer is supposed to be the smallest stable liquid formation. The heptamer is close to liquid phase and others higher number of molecules are in the transition between the liquid and solid phase. For dimer and trimer, their phases are in the transition between gas phase and liquid phase. Experimentally [20], the monomer dipole moment can be approximated by calculating second order perturbation energy of Stark effect. These studies implement an electric field to the cluster system, which it will polarizes each cluster monomer. Then after the electric field was removed, the value of monomer dipole moment should be proportional to the second-order Stark energy.
The discussion on the water phases with regard to the monomer dipole moment properties is still hard to elucidate the phase of small water cluster system. However, the data from Table II strongly support that the cyclic hexamer is more rigid than the other hexamers. The cyclic formation has the closest H-bonding distance and the strongest H-bonding energy. As opposite, the cage one has the longest H-bonding distance and the lowest H-bonding energy. Based on the H-bonding distances, Hbonding strength and the monomer dipole moment calculations, we predict that the cyclic hexamer has solid phase state and the cage one has liquid phase formation.

IV. SUMMARY AND CONCLUSION
The following items give the summary of the cluster formation stabilizations, the average binding energies and local binding energies of small water cluster systems: • By preparing reasonable initial structures using both reactivity index from Fukui functions and the symmetry considerations, we could find a more efficient way in obtaining the most stable formations of small water cluster systems.
• The significantly decreased binding energies on dimerization and cyclization reactions around −2.49 kcal/mol and −2.63 kcal/mol respectively, are strongly caused by the weak H-bond additions to the host water cluster system, leading to the linear to planar transformation.
• For the cyclic planar formations, the smaller binding energy fluctuations are caused by removing and adding the staggered and eclipsed character and also changing the average of intermolecular separation by the diameter of cyclic ring.
• From hexamer to higher number of water molecules, the planar formation transforms to 3-dimensional structures. The binding energy fluctuations follow the even and odd number of molecules. Water clusters which have the odd molecule numbers tend to have low symmetrical structures and weaker binding energies as compared with the even number of molecule ones, which have high symmetrical structures and stronger binding energies.
• The average binding energies trend for n = 1-12 tend to asymptotically converged at the average intermolecular binding value of bulk water.
• By monomer dipole moments calculations, we expect tentatively that water clusters with number of molecules 4-6 are stable in solid phase and n = 7 stable in liquid phase. For the higher number of molecules n > 7, they are in the transition between liquid and solid phase state. And for the clusters with number of molecules n < 4, they are in the transition between gas phase and liquid phase.
• The smallest piece of ice water is cyclic tetramer and the cage hexamer is the smallest stable liquid formation.