Thermally Induced Hydrogen-Bond Rearrangements in Small Water Clusters and the Persistent Water Tetramer.

Small water clusters absorb heat and catalyze pivotal atmospheric reactions. Yet, experiments produced conflicting results on water cluster distribution under atmospheric conditions. Additionally, it is unclear which "phase transitions" such clusters exhibit, at what temperatures, and what are their underlying molecular mechanisms. We find that logarithmically small tails in the radial probability densities of (H2O) n clusters (n = 2 - 6) provide direct testimony for such transitions. Using the best available water potential (MB-pol), an advanced thermostating algorithm (g-BAOAB), and sufficiently long trajectories, we map the "bifurcation", "melting", and (hitherto unexplored) "vaporization" transitions, finding that both melting and vaporization proceed via a "monomer on a ring" conformer, exhibiting huge distance fluctuations at the vaporization temperatures (T v). T v may play a role in determining the atmospheric cluster size distribution such that the dimer and tetramer, with their exceptionally low/high T v values, are under/over-represented in these distributions, as indeed observed in nondestructive mass spectrometric measurements.

Introduction

Small water clusters, W ≡ (H2O), are supposedly prevalent in the atmosphere,[1,2] catalyzing chemical reactions[3,4] of profound environmental significance, such as acid rain formation,[5] and absorbing heat radiated from the earth’s surface, which contributes to global warming.[1] Theoretical assessments of water cluster densities in the atmosphere are commonly based on an equilibrium assumption that every cluster can convert to any other cluster on the equilibration timescale.[6−9] Under these conditions, the monomer dominates at all temperatures and the cluster population diminishes nearly exponentially with n at all altitudes.[8] This is seemingly in agreement with rotation and vibrational spectroscopy, which so far identified only the water dimer in water vapor under ambient conditions[10] and the trimer at 650 K.[11] In contrast, various types of “gentle” mass spectrometric methods (that avoid excessive cluster fragmentation) suggest that humid air contains a wide distribution of cluster sizes[12−16] from which the monomer and dimer are conspicuously absent while the tetramer probability is enhanced. There is no explanation to this observation (nor to the stark differences between methods).

Water clusters have multiple conformers. At low temperatures (T), the most stable n ≤ 5 conformers are cyclic, which is a manifestation of water’s cooperativity.[17,18] Their average hydrogen-bond (HB) strength increases monotonically up to n = 6, while the dissociation energy (for W → W + W1) reaches a maximum for n = 4.[18,19]

The hexamer represents the transition from monocyclic “two-dimensional” to multicyclic “three-dimensional” structures.[20] But what happens at higher T and are there stable water clusters in the tropospheric temperature range (ca. 220–300 K)? “Many insightful spectroscopic studies have been done at very low temperatures, ... conditions [that] are far from those found in the much warmer atmosphere of Earth”.[3] For example, vibration-rotation-tunneling (VRT) spectroscopy identified HB rearrangement events, such as the acceptor switch (AS), interchange (I), and bifurcation rearrangements (B).[21,23] The splittings in the VRT bands report on the tunneling barriers of these events, but otherwise, there is no direct time or temperature dependence in these data. Indeed, for neutral clusters, temperature control is difficult to achieve. Nevertheless, it recently became possible to measure shifts in the IR OH bands of “warm” neutral water clusters (up to rotational temperatures of 160 K).[24]

The large gap between the laboratory (∼160 K) and atmospheric conditions (220–300 K) suggests turning to theory for advice. In particular, molecular dynamics (MD) simulations can readily depict the temperature dependence. For example, MD simulations of water clusters observed a “melting transition” (at a melting temperature, Tm), for example, by using the empirical “Lindemann criterion” for the root-mean-squared displacement (RMSD) of the O atoms,[25−29] from the heat capacity,[29−34] or the Q4 order parameter.[35] Unfortunately, there is a huge variability with the water model,[34,36] which renders much of the current literature qualitative at best. In addition, the molecular mechanism behind cluster melting is rarely investigated,[27] and finally, there is not always a heat capacity maximum,[27,34] hence no evidence for any transition. Consequently, a more sensitive method for detecting transitions is required.

While the literature focused on cluster melting, there is hardly any study of the analogous “vaporization transition” (temperature Tv) where monomer dissociation begins (W → W + W1). Indeed, Tv > 220 K implies that the cluster is stable in some atmospheric layers. However, melting simulations often restrain the system in various ways to prevent dissociation;[34] hence, such vital information is presently unavailable. In particular, it is unknown whether vaporization has any characteristics of a phase transition.

In terms of methodology, we have discovered a new fundamental criterion for determining the water cluster transition temperatures from (previously unreported) logarithmically small tails in their OO radial distribution function (RDF), g(r). We use these, in conjunction with the most accurate water potential (see below), for charting the HB rearrangements as a function of T and determining TB, Tm, and Tv. Although the (average) HB strength increases up to the ring hexamer,[18] we find that the tetramer is thermally the most stable cluster. In particular, it is the only water cluster with Tv > 220 K, and therefore, it might play a major role in phenomena such as global warming and acid rain.

Methods

Accuracy Assessments

For the results to be relevant to atmospheric science, the transition temperatures must be determined within the atmospheric temperature range. We are particularly interested in an accurate determination of the vaporization temperature because a cluster is expected to be stable if T < Tv, where T is the outside air temperature. Because we are not aware of experimental or computational data of cluster vaporization, we test our methodology on the melting transition. As a reference, we consider Tm for bulk ice. Different methods for determining the water potential energy function (PEF) (e.g., empirical force fields vs density functional theory (DFT)) show a huge variability (>200 K) in Tm.[36] For example, the BLYP and PBE exchange functionals of DFT give Tm = 411 and 417 K, respectively.[37] Adding dispersion corrects only a part of this error. For several of the empirical PEFs, Tm is closer to 273 K (the experimental value), but mostly because they were adjusted to reproduce this value.[38] Therefore, it is not possible to determine how accurate the transition temperatures predicted by empirical force fields are going to be for clusters, where no experimental data exists. Further progress thus depends on an accurate first-principles PEF for water.

A breakthrough in accuracy was brought about with the development of many body (MB) PEFs,[38−40] such as the MB-pol of Paesani and co-workers.[41,42] “At present, MB-pol achieves unprecedented accuracy in describing water properties from the dimer to the condensed phase and is perhaps one of the all-around best water models to date”.[40] Our test, Table S1 in the Supporting Information, suggests that the MB-pol accuracy in predicting water cluster energies is even better than previously claimed (Figure in ref (41)). In particular, we have also shown[43] that with a suitable correction to nuclear quantum effects (NQE), MB-pol reproduces the intramolecular vibrational frequencies of small water clusters to within about 5 cm–1 of the experimental values (an order of magnitude better than previous results). Classical MB-pol trajectories predict that Tm = 263.5 K for bulk water,[41] which is 9 K too low.

Figure 3

HB rearrangement events seen in the water dimer. (A) Logarithmically scaled r2g(r) at TAS, TI, TB, Tm, and Tv (red lines). Blue lines show fits to eq with parameters in Table S3. (B) Schematic depiction of the AS, I, and B events (with their VRT activation barriers)[22] and the long excursions taken by water molecules near Tv.

There are several additional error sources affecting the accuracy of the simulations. First (Figures S1 and S2), classical MD at a constant particle number, volume, and temperature (NVT) requires a stable integrator to eliminate errors from the thermostating process. The recently developed g-BAOAB integrator was shown to be “an order of magnitude more accurate for configurational averages than existing alternatives”.[44] Because we also keep the timestep small (0.2 fs), the difference from conventional Langevin dynamics exists but is not that large (Figure S3). Another factor is the length of the trajectory. While it was previously observed that a 0.5 ns trajectory is sufficiently long for the task,[27] we have found that detecting rare events near transition temperatures requires trajectories of at least 3 ns (Figure S5). Undersampling will cause shorter trajectories to overestimate the transition temperatures.

Finally, for classical MD simulations with light H atoms, it is important to estimate the NQE on water cluster melting and vaporization. Due to zero point energy and tunneling, the NQE can be large or negligible, depending on the property studied and conditions.[45] It is common to assume that, experimentally, the isotope effect on water is a measure of its NQE. Table shows the experimental melting and boiling temperatures of bulk water isotopes under standard conditions (http://www1.lsbu.ac.uk/water/water_properties.html). As the differences reduce sharply with mass, T2O is already in the classical limit,[46] suggesting a 4.5 K decrease in Tm and 1.5 K decrease in Tv due to NQE.

Table 1

Water Isotope Melting and Boiling Temperatures at 101.325 kPa (in °C)

isotope	melting	boiling
H2O	0.00	100.0
D2O	3.82	101.40
T2O	4.49	101.54

Previously, a large NQE for bulk water melting (27 K) was calculated from quantum path integral simulations on empirical, rigid water PEFs.[47] This outcome is apparently due to inadequacies in the PEF because when a flexible water model was parametrized on the basis of quantum path integral simulations (q-TIP4P/F), the NQE reduced to only 8 K.[47] Extending these calculations to isotope effects gave ΔTm ≡ Tm(T2O) – Tm(H2O) = 8.2 K,[46] supporting the conjecture that the melting NQE ≈ ΔTm (experimentally 4.5 K).

For MB-pol, quantal simulations of bulk water melting have not yet been performed, but its improved accuracy suggests NQE in the range of 4.5 – 8 K (6 K on average). Classical MB-pol simulations gave Tm that is already 9 K too low,[41] and the NQE would further reduce it, resulting in an underestimation of Tm by 9 + 6 = 15 K, which we term the “residual MB-pol error”.

Proceeding to the clusters, a quantum statistical simulation of the water hexamer[35] found about a 10 K increase in Tm from (H2O)6 to (D2O)6. Based on Table , we add an estimated 20% for (T2O)6, obtaining a 12 K NQE. By a direct comparison of classical and quantal MB-pol simulations for (H2O)6, we find (below) NQE of ≈10 K. Thus, we estimate that for correcting both the residual MB-pol error (+15 K) and the lack of NQE in classical cluster simulations (−10 K), one should add 5 K to our classical Tm values.

Performing quantum statistical simulations for equilibrium conformations over the whole temperature range might be doable but will not produce the time dependence, which is needed for determining the molecular mechanisms underlying the transitions. While methods such as centroid molecular dynamics (CMD) will generate the dynamics of quantal nuclei, below 240 K they become highly inaccurate.[48,49] Thus, quantal MD is presently not feasible. It is also not needed at this stage as the NQE is small and we have established a procedure for correcting it. Hence, we proceed with classical MD.

Computational Method Summary

We ran classical NVT MB-pol trajectories for the W clusters (n = 2 – 6) using the g-BAOAB integrator for 3 ns or more (see trajectory length column of Table S2). Excluding a brief equilibration period, we utilize the oxygen atom coordinates to obtain the OO RDF, g(r), using the gofr function of VMD (ver. 1.9).[50] For a liquid, g(r) → 1 at large r (uniform distribution). For a finite cluster with no boundary conditions, g(r) → 0, and this allows to observe its logarithmically small tails that would otherwise be submerged in the continuum.

We are interested in the radial probability density (RPD) P(r) ≡ r2g(r) that, up to normalization by Z = ∫0∞r2g(r)-dr, is the probability density for observing O atoms in a spherical shell around a chosen O atom. P(r) for W2 and W3 has a single peak so that ≡ ∫0∞rP(r)-dr/Z is the average HB length. For W4 and W5, there are two well-separated peaks corresponding to the first- and second-shell neighbors, while most W6 isomers have three solvation shells. In these cases, the HB length is obtained by restricting the integrals to the first solvation shell.

Often, g(r) or P(r) are fitted to a Gaussian function, which is inadequate for markedly anharmonic HBs. Anharmonicity leads to asymmetric distributions, which were shown[51] to be adequately described by the skewed normal distribution (SND)where x ≡ (r – μ)/σ and erf(x) is the error function (note the missing in ref (51)). The four parameters, A, μ, σ, and ξ, obtained from fitting eq to the simulated P(r), are presented in Tables S3–S8. A is for normalization, μ is the peak location, σ is its width, and ξ is its asymmetry. Due to it, the average r is

The SND is fitted to each RPD peak separately. For the first peak, ξ > 0 increases with T; μ decreases but always increases with increasing T.[51] These smooth variations represent thermal expansion in which HBs continuously lengthen, whereas HB rearrangements manifest themselves as upward deviations from the SND with abrupt onsets (“transitions”). Thus, the SND serves as a reference whose subtraction from P(r) defines the RPD tails.

Results and Discussion

Inadequacy of Current Methodologies

As discussed in the Introduction, temperature-dependent structural transitions in water clusters were usually elucidated using either caloric curves or the Lindemann criterion. Caloric curves depict versus T, and their slopes, CV = (∂/∂T)V, are the constant volume heat capacities. A maximum in CV(T) characterizes a phase transition.[29−34]Figure shows that the caloric curves for n = 2 – 6 are actually linear, and hence, CV is constant (which depends on n, see inset). Note that the temperature range extends to above Tv, and yet, there is no indication for any transition. Sometimes, the averaged potential energy is used,[29,30,32] and Figure S4 shows a linear behavior in this case as well. Considering the Lindemann criterion, we have calculated both and as a function of T (where r is the OO distance), finding that these RPD moments can show transitions under certain conditions (see Computational Methods in the Supporting Information). Nevertheless, this evidence alone is not always sufficient (Figure S6).

Figure 1

Temperature dependence of the average total energy from classical g-BAOAB trajectories using the MB-pol PEF for W2 (black bars), W3 (red triangles), W4 (green squares), W5 (blue circles), and W6 (cyan asterisks). The best-fit linear functions are shown as lines, with their equations written near each line. The black asterisks in the inset are the slopes of the five lines, d/dT kJ/(mol K), seen to be a linear function of n. An analogous plot for the average potential energy is given in Figure S4.

RDF Tails

Since the RPD moments do not always show the transitions, we have considered the full distribution. The top panels in Figure show the computed OO RPD for W4 at two values of T. As expected, there are two peaks corresponding to nearest and second-nearest neighbors. While on a linear scale there seems to be nothing unexpected with P(r), on a logarithmic scale (bottom panels) we observe low amplitude tails that extend to large distances. We shall now show that the emergence of these tails marks the onset of new HB rearrangement events, particularly those connected with cluster melting and vaporization.

Figure 2

(Top) linear vs (bottom) logarithmic plots of the RPD for the water tetramer (red line) at two different temperatures (calculated from 3 ns MB-pol trajectories using the g-BAOAB integrator). The tails extend beyond SND2 (blue), namely, the SND fit to the second RPD peak. Note that even the relatively heavy tail at 239 K is easily missed in the linear plot.

Water Dimer

We first exemplify the analysis for W2. Figure A shows P(r) at the onset temperatures for the AS, I, and B events (i.e., the lowest T where these are observed in our 3 ns trajectories, see Movies S1–S3), which are denoted by TAS, TI and TB, respectively. Figure B shows their atomistic depiction and estimated barriers according to VRT spectroscopy.[22] AS corresponds to the donor H atom switching between the two lone-pairs of the acceptor. An I event switches the roles of the donor and acceptor by concerted rotation, while B exchanges the two donor H atoms via a bifurcated HB. Figure A also shows P(r) at the melting (Tm) and vaporization (Tv) temperatures, with vaporization dynamics depicted in Movie S4.

At 20 K (not shown) P(r) coincides with the SND, but at 30 K there are slight upward deviations due to the onset of AS (Movie S1). This disturbs the HB structure very mildly, with the donor H atom sliding along the electron density connecting the two lone pairs.[52] The I and B onsets occur at 40 and 50 K, respectively, as verified from Movies S2 and S3. For W2, B is actually composed of two consecutive I events, written schematically asWe term this an “I-assisted” B-event. Here the initial geometry is as in Figure S1 of ref (43), H2′ being the free H atom of the HB donor. It interchanges with H2 via two concerted rotations.

While these onset temperatures may decrease when NQEs are included, their classical values are interesting because they reflect the barrier heights. Indeed, we find that TI/TAS = 40/30, identical to the VRT barrier height ratio, 207/157 = 1.32. However, TB/TI = 50/40 < 394/207 = 1.9, hence, I assistance results in a smaller barrier as compared to pure B events. The general agreement with the VRT HB rearrangements at low T and their order of appearance raise hopes that the experimentally inaccessible melting and vaporization transitions are also properly described by the present theory.

The first moment of P(r), (the HB length), is shown in Figure (W2). Up to TB, increases linearly with T, which represents thermal HB expansion without the HB cleavage (i.e., the effect of AS and I on this curve is negligible). Near TB, we observe a small increase in slope (dashed ochre line). At higher T, two new HB rearrangement events induce larger slope changes. We attribute the one near Tm = 110 K to “melting”. Here, we observe more rapid HB rearrangements so that the large r tail extends beyond the SND into distances characterizing second-shell water molecules (Figure A, 110 K).

Figure 4

Average OO distance, (), in W2 and W3 (W4 and W5) (red symbols). Calculated from ∫∞rP(r)dr/∫∞P(r)dr, where rmin = 0 for and 3.26 Å (3.40 Å) for of W4 (W5) (see Figure S6). Lines are piecewise linear fits (Figure S7) with parameters presented in Table S9: up to TB (bold ochre), between TB and Tm (dashed ochre), between Tm and Tv (dashed blue), and above Tv (dashed red). Gray circles depict eq with parameters obtained by fitting eq to the first (second) peak of P(r) for W2/W3 (W4/W5).

Of particular interest is the vaporization transition at Tv ≈ 140 K. It is the most conspicuous, with both an abrupt increase in slope and a jump in (Figure ). The large r tail extends beyond 10 Å (Figure A), even when the pair does not eventually separate to infinity. As Movie S4 shows, the pair may return to HB distances so that the “volume” of the cluster undergoes huge fluctuations near Tv, which is reminiscent of the volume fluctuations at the critical point of bulk liquids.[53] Above Tv, the cluster lifetimes diminish rapidly with increasing T (Table S2). For W2, Tv is exceptionally low, which might explain its absence from the mass spectra of water vapor in which care was taken to eliminate cluster fragmentation.[12−16]

The Cyclic Clusters, Trimer to Pentamer

The lowest energy isomers for n = 3 – 5 are cyclic, with monotonically increasing HB strengths.[18] Thus, one might expect the transition temperatures to increase monotonically with n. Figure shows P(r) for these three clusters from which the onset of bifurcation, melting, and vaporization can be gleaned as upward deviating tails in comparison to the SND. Contrary to the above expectation, all transition temperatures are higher for W4 than for the other clusters, with W5 behaving similarly to W3. This is corroborated by the changes of slope in Figure . In particular, only W4 has a Tv within the tropospheric temperature range.

Figure 5

Logarithmically scaled P(r) (red lines) from MB-pol simulations for (A–D) W3, (E–H) W4, and (I–L) W5 near TB (A, E, I), near Tm (B, F, J), just below Tv (C, G, K), and just above Tv (D, H, L). Blue lines show fits to eq with parameters in Tables S4–S8.

While in W3 all three O atoms are first-shell neighbors, in W4/W5, there are second-shell neighbors as manifested by the doubly peaked P(r) in Figure . Because the tails extend here beyond the second shell, we show in the bottom panels of Figure the average r for the second shell, . It is indeed more sensitive to the transitions than the overall (see Figure S6). Two exceptional behaviors are noted in Figure . First, for W4, there is a negligible change in slope of at TB (b for “solid1” vs “solid2” in Table S9). Rather, B is noted in filling up the trough between the two peaks in Figure E. Experimentally, all clusters showed the B rearrangement except W4.[21] The simulations do show B events also for W4 (Movie S5). Possibly, TB (hence the B barrier) is too high to manifest itself in the VRT lines. Secondly, for W5, decreases with T, in contrast to W4 (and to ). This may be due to the deviation of the W5 ring from planarity (ca. 50) that may increase with T, thus reducing .

The melting transition (Movies S6–S8) is manifested by the low amplitude tail in P(r) extending to OO distances characteristic of second- (W3) and third-shell (W4/W5) water neighbors (Figure B, F, J). As T increases further, a plateau develops due to open chains having no clear minima (W3) and/or populating higher energy isomers with third-shell neighbors (W5). However, up to Tv, the tails stop abruptly around 7–8 Å (Figure C, G, K), as if there is a small residual barrier preventing complete dissociation.

At Tv, new OO tails appear, extending to notably larger distances (Figure D, H, L), resulting from the huge fluctuations in cluster volume (Movies S9 and S10), and reminiscent of critical behavior in bulk water.[53] The long-range intermolecular forces between water molecules in the clusters extend to 10–12 Å. Thus, the spectacle of dissociation reversal at larger distances arises exclusively from the thermal “bath” in the NVT ensemble. The thermostat mimics an inert gas with a Boltzmann distribution of velocities and, at the transition, the dissociating pair recedes at sufficiently low velocities that can be reversed by the thermal collisions.

This can be further clarified by a depiction of a real time trajectory just above Tv (Figure ). This trajectory first undergoes dissociation reversal (near 0.3 ns). However, the ultimate fate at T > Tv is the evaporation of a monomer, leaving behind the tetrameric fragment, and this occurs near 6.26 ns. Shortly thereafter, their separation is nearly reversed, but a final collision sends the fragments off to distances from which they never recombine. The raggedness of the receding pair path is due to the thermal collisions, and indeed, at constant energy (NVE ensemble), the separation proceeds at constant relative velocity (not shown).

Figure 6

An exemplary trajectory of W5 dissociation at 219 K, W5 → W4 + W1, just above its vaporization temperature (210 K), showing the time dependence of the distance between the centers of mass of the W4 and W1 fragments. (a) Short time segment showing dissociation reversal from over 40 Å, leading to the long-distance tail observed in P(r). The inset shows a magnification thereof. (b) Behavior at longer times shows dissociation with no subsequent reversals.

Clearly, above Tv, the cluster has a finite lifetime (here, 6.26 ns). As Tv is approached from above, the lifetime becomes longer until it diverges at Tv. The trajectory length in Table S2 is an integer number below Tv meaning that it has completed its preset lifetime without dissociating. At Tv, the signature for vaporization is attempted dissociation with concomitant reversal, manifested statistically as long-distance tails of P(r) (Figure L). Because the reversal is a random event, such clusters are bound to dissociate at sufficiently long times.

Transition Temperatures

For AS, I, and B, the transition temperatures were deduced by visualizing the low-temperature trajectories, searching for the lowest T at which each event appears. For melting and vaporization, the transition temperatures were determined from the first appearance of the RDF tails (Figure ), agreeing semiquantitatively with the temperatures at which a change of slope occurs in Figure . This is summarized in Table .

Table 2

Transition Temperatures (in K) for HB Rearrangement Events in Small Water Clusters, Using ≥3 ns MB-pol Trajectories

cluster/event	AS	B	meltinga	vaporizationa
dimer	30b	50c	110	140
trimer	20	75	140	213
tetramer	50	100	200	240
pentamerd	10	85	125	210
hexamer			115	205

Add 5 K to correct for residual MB-pol error and NQE.

I events occur as of 40 K for the dimer, but are not observed for the cyclic clusters.

I-assisted B event.

Sole previous MD study reporting multiple transitions found for W5: TAS = 35 K, TB = 120 K, and Tm = 150 K.[27]

Bifurcation, melting, and vaporization involve momentary or more permanent HB cleavage events, and thus, TB, Tm, and Tv show a similar n dependence (Table ), with the transition temperatures notably higher for the tetramer than for the other clusters. For Tv, this could have been anticipated from the qualitative trend in monomer dissociation energies,[18,19] but quantitative results for high T are difficult to obtain from 0 K quantum chemistry calculations. To check this, consider the high-accuracy CCSD(T)/CBS calculations of Temelso et al.[19] for the free energy change, ΔG(T), for W → W + W1 at 1 atm pressure, using scaled harmonic frequencies and assuming that the reaction occurs from the average W isomer. From the equilibrium condition ΔG(Tv) = 0, we obtain, by interpolating in Table 5 of ref (19), Tv = 190, 250, 316, 291, and 226 K for n = 2,...6, respectively. While the trend is the same as in Table , these temperatures are markedly higher. Of course, the two calculations are not precisely the same, as we only simulate the dissociation reaction in vacuum (and not the recombination), albeit with the accurate anharmonicities embedded in MB-pol. However, it is possible that the largest discrepancy arises from Boltzmann averaging ΔG over all the 0 K isomers,[19] whereas the reactions in Table occur from excited conformations having lower dissociation barriers (see the vaporization mechanism below). This two-step pathway is preferable when the excitation energy is lower than the ground state barrier, leading to the lower Tv values (and showing that for quantitative results, there is no alternative to the dynamical approach).

Underlying Molecular Mechanisms for Melting and Vaporization

Figure shows a sequence of HB rearrangements for W4 near Tm (Movie S7). First, a ring HB temporarily cleaves, generating a chain with two free H atoms at one end. One of these attacks the next nearest water, producing a cyclic trimer with one external water molecule donating an HB to it (the “3 + 1 isomer” in Figure B). This singly HBed molecule undergoes facile rotation to become an HB acceptor (panel C). Now, it has two free H atoms, one of which attacks the next nearest neighbor regenerating a W4 ring in which two water molecules have been exchanged (Figure D). Thus, the tetramer is the smallest entity in which water exchange occurs, a predecessor of water self-diffusion in the bulk.

Figure 7

Possible sequence of HB rearrangements leading to water monomer exclusion from the tetrameric ring (A–C) followed by (D) exchange (Tm and above) or (E) vaporization (Tv and above). Double-headed arrows depict HB cleavage whereas single headed arrows depict HB formation.

Near Tv, the HB of the external water molecule elongates, leading to dissociation by a single HB cleavage. If vaporization was a one-step reaction involving double HB cleavage, one would expect Tv ≈ 2Tm. As can be seen in Table , this is not the case. For W4, Tv is just 20% higher than Tm, which is explained by the two-step mechanism, each step governed by a single HB cleavage.

The Hexamer

This cluster has a plethora of low lying isomers (see Figure A in ref (35)). It is the smallest water cluster whose lowest energy isomer (the cage isomer) has a 3D structure, rather than being a planar ring. The tetramer stability is manifested in the hexamer low-lying isomers (cage, prism, and book), which are composed from several tetrameric rings. Next, in order of stability, is a group of isomers constructed of a pentameric ring with a water monomer forming two HBs to it. As of 90 K, we observe structures from this group of isomers. At 110 – 120 K, the sixth water molecule cleaves one of its HBs, remaining bound through a single HB (“5 + 1” structure, see Movie S11). This is likely the origin of the tail seen, for example, in Figure A at 119 K. Quantum statistical MB-pol simulations of (H2O)6 found a melting transition between 90 and 120 K (average 105 K).[35] Our classical MD value is about 10 K higher, in excellent agreement with the 12 K NQE deduced (Methods) from the isotope effect on Tm. It is also in excellent agreement with Tm= 116 K obtained using the TIP4P-ice water model,[34] which (contrary to MB-pol) has been parametrized to reproduce the experimental melting point of bulk ice.

Figure 8

Logarithmically scaled P(r) from MB-pol simulations of the water hexamer (A) near the melting temperature and (B) near the vaporization temperature. The dynamics giving rise to the pronounced tail (cyan line) just above Tv are depicted by Movie S12. The starting geometry in these simulations is that of the cage isomer.

The hexamer vaporization temperature is around 205 K (see Figure B). In spite of its 3D structure, we find a similar dissociation mechanism: W6 dissociates from a “5 + 1” structure, similar to W5 and W4 dissociating from “4 + 1” and “3 + 1” structures, respectively. Movie S12 shows the dynamics at 214 K (just above Tv), with the water monomer making long excursions (up to 25 Å) before returning to a book structure. In this example, dissociation did not occur up to t = 3 ns (the preset trajectory length), but we know it must occur at some later time, when thermal collisions do not succeed in reversing the dissociation attempt.

Conclusions

We have presented a first detailed dynamical description of thermal structural rearrangements of small water clusters (W) with a reliable PEF, from the low-temperature AS, I, and B events and up to vaporization, when the cluster disintegrates. We have introduced a new operational measure for the melting and vaporization transitions, manifested by logarithmically small tails of the RDF. For n = 4 – 6, these tails are attributed to the formation of an (n – 1) ring structure to which a monomer is attached. At Tm, such structures first appear, whereas at Tv, the single HB connecting the monomer to W cleaves. For W4, this mechanism (“3 + 1”) is the least favorable because it produces a strained W3 ring. As a result, both Tm and Tv are notably higher for W4 than for the other clusters.

The first appearance of the RDF tails is a particularly accurate monitor of vaporization, offering a first quantitative determination of Tv for water clusters. Tv marks the threshold of cluster stability toward unimolecular dissociation and thus has characteristics of a phase transition, as manifested also by the abrupt change in vs T and the critical-like distance fluctuations. Interestingly, all investigated clusters have their Tv below the atmospheric temperature window except the tetramer for which, after corrections (Methods), Tv ≈ 245 K.

The Persistent Tetramer

Theoretical assessments of water cluster densities in the atmosphere are commonly based on an equilibrium assumption that every cluster can convert to any other cluster on the equilibration timescale,[6−9]Our results cast doubts on the equilibrium assumption in the sense that some of the rate constants in this scheme vanish. First, Table S2 shows that the dimer at 180 K is already very short-lived, so by 220 K, it will not live long enough for undergoing bimolecular collisions. Generally speaking, a typical molecule experiences about 1010 collisions per second at standard pressure and 25 °C,[54] so that the typical time between collisions is larger than 100 ps. Hence, the dimer at 220 K will not form trimers and will dissociate faster than two water monomers collide. Second, for T ≤ 245 K, the reaction W4 → W3 + W1 would cease to occur. Thus, the cluster distribution is expected to be depleted in dimers and enriched in tetramers, as indeed observed in the “gentle” mass spectrometric measurements of water cluster distributions in humid air.[12−16]

More quantitative work for unraveling water cluster distributions in the atmosphere could involve the simulation of a large number of MB-pol water molecules in the gas phase for various densities and temperatures. Such simulations could test the hypothesis that, under appropriate conditions, the tetramer is an important atmospheric “workhorse”, absorbing IR radiation in its expanded ring system, energy which becomes available for carrying out chemical reactions.