Isomerization, melting, and polarity of model water clusters: „H2O...6 and „H2O...8

Energetics, structural features, polarity, and melting transitions in water clusters containing up to eight molecules were studied using ab initio methods and empirical force field models. Our quantum approach was based on density functional theory performed at the generalized gradient approximation level. For the specific case of (H2O)6, we selected five conformers of similar energy with different geometries and dipolar moments. For these cases, the cyclic arrangement was found to be the only nonpolar aggregate. For (H2O)8, the most stable structures corresponded to nonpolar, cubic-like, D2d and S4 conformers. Higher energy aggregates exhibit a large spectrum in their polarities. The static polarizability was found to be proportional to the size of the aggregates and presents a weak dependence with the number of hydrogen bonds. In order to examine the influence of thermal fluctuations on the aggregates, we have performed a series of classical molecular dynamics experiments from low temperature up to the ...


I. INTRODUCTION
In recent years, the subject of water clusters with typical lengthscales in the nanometer range has conceited an increasing attention from the theoretical [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] and experimental [14][15][16][17][18][19][20][21][22][23][24][25] points of view.This interest is mainly due to the fact that investigations on small water clusters are a perfect means to characterize structural changes and bonding mechanisms in passing from isolated molecule to bulk states.From the theoretical point of view, a large body of computational work has been devoted to the study of structural features and binding energies of small water clusters at the Hartree-Fock, 1 MP2, [2][3][4]7,15 and DFT 8,10,11,14 levels. At pesent, there is little doubt that the minimum energy structures of the water trimer, tetramer, and pentamer correspond to quasiplanar, cyclic rings and that the heptamer, octamer, and larger clusters are expected to exhibit three-dimensional geometries.The specific case of the water hexamer minimum energy structure is more subtle since it appears to be in the crossover region between cyclic planar and tridimensional structures.
Results from molecular beam electric deflection experiments 18 performed more than 25 years ago seem to indicate that, under the particular experimental conditions, water clusters ranging from the trimer up to the hexamer exhibit nonpolar behavior.Furthermore, these findings have been interpreted as a sufficient proof to infer that the stable structures were cyclic. 10However, a recent terahertz laser vibration-rotation tunneling ͑VRT͒ spectroscopy study performed by Saykally et al. 16 at temperatures of the order of 5 K, provided a conclusive evidence that the most stable structure for ͑H 2 O͒ 6 corresponds to the so-called cage configuration ͓see Fig. 1, ͑H͔͒.With this state of the art spectroscopic technique, Saykally et al. were able to explore the manifold of different low lying energy minima produced by flipping the relative orientations of different hydrogen bonds.Considering the low temperature of these experiments, the relevant mechanism that drives these transitions is mainly quantum proton tunneling.These results have been confirmed by Quantum Monte Carlo simulations performed by Gregory and Clary, 9 as well.In addition, by inducing Stark splittings within VRT techniques, the Saykally group also succeeded in directly measuring the dipole moment of different water clusters; 25 for ͑H 2 O͒ 6 , they found that the magnitude of the dipolar moment was consistent with the cage structure derived from the previous VRT experiments.
At this point, there seems to be a paradox between the first deflection experiments and these new sets of spectroscopic results.In this article, we will try to provide one possible explanation to reconcile these apparent contradictory observations; more specifically, we will show that the simple consideration of thermal fluctuations may lead to a unified consistent picture.We have already mentioned that due to the low temperature involved, VRT experiments basically a͒ Also at Unidad Actividad Quı ´mica, Comisio ´n Nacional de Energı ´a Ato ´mica, Avenida Libertador 8250.͑1429͒ Capital Federal, Argentina.explore regions of the potential energy surface that correspond to minimum energy conformations.A completely different scenario may prevail in deflection experiments where the involved temperatures may be sufficiently high 26 to allow for a richer variety of aggregates including equilibrium between different conformers and even the possibility of solidliquid-like phase transitions.
Anyhow, the crux of the matter still resides in the viability of characterizing the structure of a given cluster using as a signature the magnitude of its dipole moment.If this is possible, one could conceive to use it as a tool to monitor for example, isomerization and phase transitions in nanoaggregates.This represents the main motivation of this study that includes a series of classical molecular dynamics ͑MD͒ experiments monitoring simultaneously different order parameters like the number of hydrogen bonds, the Lindemann's index 27 along with the corresponding dipole moments.Before doing so and as a preliminary test for our classical force fields, we have compared predictions from the MCY 28 and TIP4P 29 Hamiltonians with results from density functional theory ͑DFT͒ calculations for the interaction energies, bond lengths and dipole moments for selected stable structures of ͑H 2 O͒ n , with 1рnр8.
The organization of the article is as follows: in Sec.II we give technical details of both, electronic structure and molecular dynamics computations.In Sec.III we present the results for 0 K structures and simulation results covering a temperature interval from Tϭ10 K up to TϷ200 K.The concluding remarks are left for Sec.IV.

A. DFT calculations
The electronic structure calculations reported in this work have been carried out within the framework of the DFT formalism. 30The Kohn-Sham self-consistent procedure 31 was implemented to obtain the electronic density and energy through the determination of a set of one electron orbitals.Gaussian basis sets were used for both, the expansion of the one-electron orbitals and for the electronic density.Matrix elements of the exchange-correlation potential were evaluated by the numerical integration scheme proposed by Becke. 32During the self-consistent cycle, the integration was performed on a set of coarse, atom-centered, spherical grids with 25 and 20 radial shells for oxygen and hydrogen atoms, respectively, and 50 angular points per shell.At the end of the procedure, the exchange-correlation energy was evaluated using an augmented, finer grid with 35 and 30 radial shells for oxygen and hydrogen atoms, respectively, and 194 angular points per radial shell.Geometry optimizations have been performed without symmetry constraints by using a quasi-Newton minimization method in Cartesian coordinates with analytical energy gradients.The computations have been performed at the generalized gradient approximation level.The combination of gradient corrections for exchange of Perdew and Wang 33 and for correlations developed by Perdew, 34 has been used, along with the Dirac's exchange term and the parametrization of the correlation energy of the homogeneous electron gas due to Vosko. 35Basis sets with contraction patterns ͑73111/521/1͒ and ͑721/1͒ 36 augmented with two diffuse polarization functions complemented with an auxiliary bases sets with contraction patterns of ͑1111111/ 11/1͒ and ͑111111/1͒ have been employed for O and H, respectively. 11We have checked that these sets provide reasonably converged results for the computed dipole moments, geometries, vibrational frequencies and interaction energies of the water monomer and dimer. 11t was also interesting to compute static polarizabilities for water clusters.We did that using a finite-field scheme 37 that included an extra term into the electronic Hamiltonian depending on an external uniform electric field F. The resulting Kohn-Sham one-electron equations were of the type where in the last equation, the effective potential V eff (r) includes ionic, Coulombic, and exchange-correlation terms, 2ប is the Planck constant, e and m represent the electron charge and mass, respectively, and i (r) are one-electron orbitals.From the self-consistent solution of Eq. ͑1͒, the dipole moment and total energy can be determined as a function of the applied electric field.
The cluster polarizability ␣ is defined in terms of the trace of the polarizability tensor: where the tensor elements ␣ i j are defined as , i, jϭx,y,z.͑3͒ In the last equation, i represents the ith component of the total dipole moment of the cluster and E is the total energy; x, y, and z are arbitrary orthogonal directions.Field induced polarization functions 38 have been added to the basis sets. 36It is well known that the addition of these diffuse functions is required to compute polarizabilities accurately. 39In order to determine the stability of our calcu- lation with respect to the magnitude of the applied electric field, the influence of this parameter on the computed polarizability of the water monomer and dimer has been examined.The best agreement with linear response predictions was found for field satisfying 0.005 a.u.р͉F͉р 0.02 a.u.Within this range, polarizabilities computed from the first derivatives of the dipole moment and from the second derivative of the total energy agreed better than 1%; consequently the polarizability values reported in the following sections were computed for ͉F͉ϭ0.01a.u.Finally, note that nuclear positions have not been allowed to relax in the polarizability calculations.

B. Molecular dynamics calculations
Molecular dynamics experiments were performed using two classical force fields for water: the MCY model of Matsuoka et al., 28 and the mean field TIP4P model developed by Jorgensen et al. 29 In both cases, H 2 O is modeled as a rigid, four-site molecule.The MCY pair potential 28 is the sum of exponential terms between oxygen and hydrogen sites plus Coulomb tails between partial charges located at the positions of the hydrogens and at an auxiliary site slightly shifted from the oxygen towards the hydrogens along the molecular C 2v axis, respectively.The TIP4P model, 29 possesses similar Coulomb tails but the exponential short-range terms are replaced by a single Lennard-Jones contribution between the oxygen sites.The molecular dynamics experiments consisted in rather long canonical runs 40 of about 3-5 ns at a series of different temperatures ranging from Tϭ10 K up to approximately Tϭ200 K; within this temperature range, we found liquid-solid phase transitions with a only a few episodes of evaporated molecules.In all experiments, the Verlet algorithm 41 was used to integrate the equations of motions using typical time steps of ␦tϭ0.5 fs; to handle intramolecular constraints in the molecules, the SHAKE algorithm was implemented. 42

A. Characterization of stable structures
The clusters considered in this work are sketched in Fig. 1.Binding energies, geometrical parameters and spectroscopic information of many of the structures ͑A, B, C, D, E, F, K, M, and N͒ have already been investigated in a previous DFT study. 11In Table I we compare binding energies for the different optimized DFT structures to those obtained using the TIP4P and MCY force fields.DFT results include basis set superposition errors evaluated within the counterpoise scheme. 43he true value of the water dimer binding energy still is not a completely settled issue.Assuming as valid references the experimental result 44 D e ϭϪ5.4Ϯ0.7 kcal/mol and the accurate ab initio prediction 45 D e ϭϪ4.8 kcal/mol, one could conclude that the present DFT scheme yields a slightly overestimated value D e ϭϪ5.8 kcal/mol. 46In a different context, it is well known that empirical force fields parameterized to reproduce experimental liquid state properties, normally overestimate the water dimer binding energy.The TIP4P model falls into this category; the particular case of the MCY potential is somewhat different since it has been deliberately fitted to reproduce ab initio configuration interaction calculations for the water dimer.It yields D e ϭ-5.88 kcal/mol, a value somewhat closer to the experimental results.
A relevant geometrical parameter to characterize optimized structures is the nearest Oxygen-Oxygen distance, d O-O .For the water dimer, the experimental estimate is d O-O ϭ2.98 Å. 44 Our DFT calculations predict d O-O ϭ2.907 Å which agrees reasonably well with the experimental result and with other previous ab initio calculations. 47On the other hand, most classical force field models yield values of d O-O somewhat smaller than the observed experimental value. 48he TIP4P and MCY estimates are: d O-O ϭ2.75 Å 29 and d O-O ϭ2.87 Å, 28 respectively.For larger aggregates, typical values of d O-O computed using DFT normally present a slight reduction as cluster size increases; this reduction is less marked in results obtained using classical nonpolarizable Hamiltonians and has been already observed in previous ab initio calculations. 1,49heoretical predictions for the structure of the ͑H 2 O͒ 6 may yield results that depend on the specific level of theory implemented.On one hand, several theoretical works using DFT, 8,10 as well as MP2 calculations 2 predict the global minimum to be a ring species.On the other hand, more recent theoretical calculations using MP2 techniques 14 would seem to indicate that the global minimum is the cage form.In this work, we focused our attention on five different water hexamers ͑referred to as ͑E͒, ͑F͒, ͑G͒, ͑H͒, and ͑I͒ in Fig. 1͒.Our DFT calculations show that the conformers ͑F͒ and ͑H͒ posses almost degenerate binding energies, slightly lower than the rest of the structures ͑E͒, ͑G͒, and ͑I͒.The energy landscape for the MCY and TIP4P mean-field potentials is similar to the DFT one; however, in the TIP4P calculation the energy gaps between the different isomers are somewhat larger than in the MCY case.
For the ͑H 2 O͒ 8 clusters, it is well established that cubiclike structures which maximize the number of hydrogen bonds are the most stable forms; [3][4][5] the rest of less coordinated isomers, lie at much higher energies.We have restricted our study to the following conformers: a ring form ͑J͒ with 8 hydrogen bonds, a trilobate form ͑K͒ with 10 hydrogen bonds, a propeller form ͑L͒ with 9 hydrogen bonds and two cubic-like structures with D 2d and S 4 symmetries, each one with 12 hydrogen bonds.As expected, both DFT and force field calculations predict that the almost degenerate cubic forms ͑M͒ and ͑N͒ are the most stable structures.
Information about the charge distribution within the aggregates as reflected by the magnitude of the cluster static dipolar moments is presented in Table II.DFT results for the isolated water molecule yields the value of ϭ1.822D, which is in very good agreement with the experimental result of ϭ1.847D. 50For the dimer, DFT predicts a value of ϭ 2.606D which again agrees reasonably well with the experimental result of 2.60D. 44The dipolar moment of the MCY and TIP4P monomers are larger than the experimental value since these potentials were optimized to fit bulk properties.The symmetric cyclic hexamer ͑E͒ is the only nonpolar, while the rest of the structures are considerably polar.Our predicted dipolar moments agree reasonably well with previous calculations. 10,25For the ͑H 2 O͒ 8 clusters, the nonpolar structures correspond to the symmetric cubic-like ͑M͒ and ͑N͒ forms and the ring-like structure ͑J͒.

B. Molecular dynamics simulations
One of the most common means to monitor solid-liquid transitions in nanoclusters is through the examination of relative fluctuations in the interparticle distances as they are reflected by the Lindemann parameter ␦, defined as where r i j represents the distance between the ith and jth oxygen atoms, ͗ . . ..͘ represents a canonical ensemble aver- age, ␦xϭxϪ͗x͘ and the sum includes the N constituents of the aggregate.Usually, it is assumed that melting occurs when ␦ attains Ϸ0.1.In Fig. 2 we present results of ␦ for ͑H 2 O͒ 6 ; for both, the MCY and TIP4P Hamiltonians, the melting temperature can be located at T m Ϸ50 K, where ␦ shows a sharp rise of almost one order of magnitude in a narrow temperature interval of about 20 K.The analysis of the number of hydrogen bonds, N HB , computed as the number of O-H•••O moieties with a d O-O smaller than 3.1 Å, also shown in Fig. 2, is consistent with the previous results of the Lindemann parameter.In this case, as the cluster melts at TϷ50 K, N HB presents a steady drop from its low temperature plateau value of N HB Ϸ8 down to N HB Ϸ6.5 for temperatures of the order TϷ200 K when spontaneous evaporation of some of the constituents begin to appear.Although our previous description based on the Lindemann parameter and the N HB is fully valid, one can still wonder whether the melting process, or equivalently, any solid-solid isomerization process in clusters, could be monitored by another order parameter of more direct experimental access.Perhaps the first and most natural choice would be to consider a quantity related to the electrical response of the aggregates.At least two types of cluster beam experiments provide information about the electric moments of the aggregates: ͑i͒ those in which a transverse inhomogeneous electrostatic field exerts a deflecting force on polarized clusters; 51,52 ͑ii͒ the second type of experiments involve an homogeneous electrostatic field used together with a radiation field provoking a transition between two Stark spectral splitted states. 25e start by analyzing the temperature dependence of the time averages of the dipolar structure of the different clusters.The computation of a suitable scalar characterizing the polarity of a given aggregate deserves a brief comment.For solid-like structures, i.e., those in which the molecules remain almost immobile maintaining practically the same intermolecular distances and mutual orientations, details about the internal polarization within the clusters can be obtained by simply computing the average value of modulus of the dipole moment ͉͉͗͘.However, for liquid-like entities in which there can be considerable diffusion and/or important changes in the interparticle orientations, the single consideration of ͉͉͗͘ may not longer be appropriate.To rationalize this assertion, imagine a liquid-like aggregate whose molecular spatial distribution looks elongated along a given z axis whose direction coincides with that of the instantaneous largest moment of inertia of the cluster.Furthermore, consider that at a given time t 0 , the cluster presents an instantaneous charge distribution that yields a net dipole vector (t 0 ) pointing towards the positive z direction.Imagine now an spontaneous fluctuation at time t 1 , that involves almost no change in the spatial arrangement of the cluster molecules except a global flip of all the individual molecular dipoles into opposite directions, i.e., (t 1 )ϭϪ(t 0 ).It is clear that such two configurations would contribute to a nonzero value of ͉͉͗͘, while from a statistical point of view, their contri- butions should actually cancel out.The previous considerations led us to consider one possible alternative scheme to characterize polarity in liquid-like environments: at each time step of the simulation experiment, we considered an instantaneous local coordinate system characterized by versors ͕e ˆa ,e ˆb ,e ˆc͖ coinciding with the directions of the princi- pal moments of inertia of the cluster 53 and computed the corresponding projections: ϭ a e ˆaϩ b e ˆbϩ c e ˆc , where i ϭ•e ˆi , iϭa,b,c.A reasonable estimate for the cluster polarity may be obtained from an average value of the instantaneous projections ¯defined as

͑5͒
In the middle panel of Fig. 2 we present results for ¯for ͑H 2 O͒ 6 ; note that the size of the spontaneous thermal fluctuations at Tϭ10 K are sufficient to make ¯decrease from its 0 K value of 1.9D down to 1.45D; basically, this change can be ascribed to small variations in the intermolecular orientations.Moreover, an even more dramatic drop occurs during the melting transition, where ¯reduces to 0.2D, vanishing for all practical purposes beyond TϷ120 K.The simultaneous analysis of the temperature behavior of the three parameters ␦, N HB and ¯led us to formulate the following picture for the thermal dependence of the water hexamer spatial and dipolar moment: at low temperatures, water hexamers are likely to exhibit caged-solid-like structures characterized by roughly eight hydrogen bonds and a dipole moment of about ¯ϭ 1.5-2.0D.As the cluster attains the melting transition, the orientational distribution of the aggregates with respect to the local instantaneous system ͕e ˆa ,e ˆb ,e ˆc͖ becomes more isotropic leading to smaller values of ¯.We tend to believe that this observation may provide a plausible explanation for the apparent paradox mentioned in previous paragraphs: the VRT experiments are performed at temperatures sufficiently low that only cage hexamers are present while in the deflection beam experiments, the temperatures may be high enough to lead to a preponderance of liquid-like entities.Consequently, it would be incorrect to interpret the lack of polarity in terms of the prevalence of ring nonpolar solid like structures.
Interestingly, we were also able to tune the temperature of the aggregates so as to detect a spontaneous premelting interconversion equilibrium between two essentially solidlike structures corresponding to a cage ͑H͒ and a book-like ͑I͒ structures at Tϭ48 K.For this particular case, we focused our attention on the free energy profile W(I M ) associated with the equilibrium To distinguish between the two conformers and to monitor the dynamics of the interconversion, we found convenient to adopt the magnitude of the cluster largest moment of inertia I M as the appropriate order parameter.More specifically, we computed the probability distribution P(I M ) for the cluster largest moment of inertia:

͑7͒
where ␤ Ϫ1 ϭk B T is Boltzmann's constant times the temperature.In the top panel of Fig. 3 we present the time evolution of I M during a 5 ns trajectory.Note that for most of the time, the magnitude of the largest moment of inertia presents small fluctuations around 4.9 amu nm 2 ; that corresponds to hexamers exhibiting cage-like structures.From time to time, this regularity is interrupted by a few spikes in which I M rises up to values of the order of 7 amu nm 2 ; an analysis of several snapshots shows that these episodes can be ascribed to interconversions to book-like structures.The free energy profile is shown in the bottom panel of Fig. 3; it was constructed from a histogram of I M .The profile shows that cage and book structures are separated by a free energy barrier of Ϸ0.5 kcal/mol measured from the most stable cage conformer.Our estimation for the free energy difference between the two stable structures is ⌬W cage-book ϭ0.3 kcal/mol; this value should be compared to the corresponding energy difference ⌬E cage-book ϭ0.6 kcal/mol ͑see Table I͒ leading to a positive value of ⌬S cage-book which is a reasonable result considering that the less coordinated forms are expected to be entropically favored. 5he characteristics of the caloric curves for ͑H 2 O͒ 8 aggregates contrast sharply with those of the hexamer cases since they are rather featureless.Results for the Lindemann's parameter and the N HB are shown in Fig. 4; MCY and TIP4P models yield similar thermal behavior and predict a melting transition at around the same temperature T m Ϸ160 K. Low temperature, solid-like clusters correspond to cubic aggre-gates characterized by almost no dipolar moment and melt into liquid-like structures with negligible polarity.N HB decreases monotonically with temperature for both, TIP4P and MCY simulations; in the MCY case, the drop looks somewhat more pronounced, and this may be related to the smaller energy gap that exists between the cubic-like and more open structures, like the book-like ͑L͒ or ring-like ͑J͒ isomers.
Similarly, as it was found for the hexamer case, we could also detect a solid-solid premelting transition between the ͑M͒ and ͑N͒ cubic-like forms that present D 2d and S 4 symmetries, respectively.Note that the structures present two kinds of three-coordinated molecules: ͑i͒ double proton donor-single proton acceptor ͑DDA͒ and ͑ii͒ single proton donor-double proton acceptor ͑DAA͒.DAA molecules are connected exclusively to DDA, and vice versa.These two structures have been identified in recent infrared experiments in isolated clusters 24 and in experiments involving benzene attached to ͑H 2 O͒ 8 . 15From a theoretical point of view, coexistence of these aggregates has also been discussed in the context of MC simulations performed by Tsai and Jordan. 54he spontaneous isomerization between these two structures is not observed at low temperatures, say, below Tϭ100 K probably due to the nonergodicity inevitable when transitions require an escape through a free energy barrier much larger than k B T. However, at TϷ 160 K, the magnitudes of the thermal fluctuations are sufficiently large so that the interconversion between the two forms starts to occur in a time scale comparable to the length of the simulation experiments.
As in the case of the water hexamer, we adopted the largest moment of inertia as a valid tool to monitor the interconversion between the two cubic forms.The top panel of 5 shows the time evolution of I M along a 5 ns trajectory at Tϭ160 K; this time, the spikes correspond to a collection of liquid-like, open structures, with N HB close to 10-11.Yet, a description based solely on the magnitude of I M is still incomplete since this parameter does not reflect the differences in symmetry of the stable structures.To be able to distinguish between these structures, we found convenient to define a second parameter defined as The calculation of cos( i ) required first, the identification of the four versors v ˆj taken along the directions of the hydrogen bonds where DAA molecules behave as donors.Secondly, we grouped the DAA molecules in N pairs; two DAA mol- ecules belong to the same pair if both are recipients of hydrogen bonds from the same DDA molecule.Finally, we computed cos( i )ϭv ˆj•v ˆk , where j and k belong to the same ith pair.We remark that for Tϭ0 K, D 2d structures are characterized by Ӎ1 and Nϭ2, while S 4 clusters are characterized by Ӎ0 and Nϭ4.
In the bottom panel of Fig. 5, we present results for the free energy associated with the interconversion between the D 2d and S 4 structures.The procedure through which we generated this plot is the following: we started by building up two histograms for I M : one for configurations with 0.5р р1 (D 2d structures͒ and a second one for those with 0р р0.5 (S 4 structures͒.Next, in order to match these two curves, we had to characterize a moment of inertia (I M † ) associated with the transition state.The analysis of the dynamical trajectory provided us with some helpful guidelines to make a reasonable choice.In particular, we noticed that isomerization would not take place unless I M did not surpass a threshold value of about 15 amu nm 2 ; without being too rigorous, we decided to locate the transition state, I M † , at that value. 55The magnitude of the resulting barrier is close to 1.3 kcal/mol.We also note that the stable states do not present the same free energy, namely the D 2d conformer lies 0.1 kcal/mol below the S 4 structures.If we consider the energy difference between the structures given by the MCY mean field model potential ͑see Table I͒ ⌬EϷ0.3Ϫ0.5 kcal/mol, the difference in free energy can be interpreted in terms of the more entropic ͑less symmetric͒ S 4 structure.In closing this analysis, we also remark the presence of two secondary minima located at I M Ϸ10 amu nm 2 ; although we tend to believe that they are not the product of an artifact of the simulations, we were unable to make a proper identification of them with any of the rest of the different minimum energy structures shown in Fig. 1.
Of course, the description of the polarization structures is not limited to the sole analysis of the permanent dipoles; in fact, deflection beam experiments can also be used to determine higher order terms of the electric response, such as, for example, the static polarizability. 51,52For example, an analysis of the polarizability of clusters composed by semiconductors 51 shows strong oscillations in the electrical response as a function of the number of atoms, a fact that would indicate an interplay between the cluster structure and its electrical properties.In Table II we have included DFT results for ␣ corresponding to the minimum energy structures ͑A͒-͑N͒.The computed value for the water monomer is in very good agreement with experimental 56 and previous DFT results; 38 to the best of our knowledge, no experimental or computational results are available for larger aggregates.The most important observation in this respect is the almost insensitivity of ͗␣͘/n with the cluster size.A weak depen- dence of the computed value on the N HB was also noted: the larger the number of hydrogen bonds, the smaller the computed polarizability.Consequently, we tend to believe that, at least for the case of small water clusters, the static polarizability does not seem to yield any relevant structural information beyond the size of the aggregates.

IV. CONCLUDING REMARKS
In this article we have presented a series of DFT and molecular dynamics results that provide a comprehensive picture of the relationship that exists between spatial and polarization fluctuations occurring in small water aggregates.In particular, we focused our attention on two cluster sizes, namely ͑H 2 O͒ 6 and ͑H 2 O͒ 8 .For the smallest aggregate, we have constructed the caloric curve for two widely used pseudopotentials: the MCY and TIP4P to locate the melting transition.These model Hamiltonians have the peculiarity that predict the correct energy landscape for the different ͑H 2 O͒ 6 conformers at 0 K.
The simultaneous analysis of the thermal dependence of spatial and dipolar fluctuations has revealed that the simple consideration of the net dipolar moment in the aggregate was sufficient to provide an unambiguous signature not only of solid-liquid like phase transitions, but also of interconversions between solid-like isomers.We also analyzed the case of the water octamer.For this cluster size, the only solid-like aggregates exhibiting a net polarity lie several k B T above the cubic-like minimum energy conformers; consequently, this cluster melts into liquid-like forms, which are isotropic and nonpolar in average, before higher energy conformers become favored.The fact that liquid-like water hexamers and octamers exhibit negligible polarity is not at all an obvious result.One could speculate about a possible coupling between spatial anisotropy due to the presence of a free surface and the special characteristics of the water intermolecular potential that would result in a net dipolar pattern even in aggregates with non-negligible intercluster diffusion.This phenomenon has been found, for example, in the case of clusters composed by Stockmayer model molecules. 57till, there are many important aspects-such as for example, the explicit consideration of local polarization effects-that have not been included in the models presented in this study.Incidentally, we remark that we have also performed several test runs using more sophisticated polarizable model Hamiltonians, such as the fluctuating charge model proposed by Rick et al. 58 However, for this force field, the minimum energy configuration for the water hexamer corresponds to the cyclic structure ͓see Fig. 1͑E͔͒.These drawbacks are not totally unexpected considering that this more elaborate model was parametrized to reproduce thermodynamic information for bulk water at normal conditions; consequently, we refrained from proceeding any further with this model.Anyhow, although the study of polarization effects will surely deserve further investigation, we still believe that the results obtained with the TIP4P and MCY models are able to capture the main features that govern the interplay between phase transitions and global polarization in two water aggregates that have been extensively studied in recent years.

FIG. 1 .
FIG. 1. Water clusters studied in the present work.

FIG. 3 .
FIG. 3. Top panel: Largest moment of inertia for ͑H 2 O͒ 6 at Tϭ48 K as a function of time.Lower panel: free energy profile as a function of the largest moment of inertia for the same cluster.

TABLE I .
Binding energies of water clusters ͑in kcal/mol͒.The number of hydrogen bonds is given in parentheses.

TABLE II .
Dipole moments and average polarizability per molecule of ͑H2O͒ n .