Barrier-free molecular reorientations in polyhedral water clusters

doi:10.21203/rs.3.rs-1421948/v1

Download PDF

Research Article

Barrier-free molecular reorientations in polyhedral water clusters

https://doi.org/10.21203/rs.3.rs-1421948/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

The structural stability of the weakest bonded configurations of water clusters in the form of gas hydrate cavities D and T, as well as in the form of Kelvin’s polyhedron, was tested using the local geometric optimization with the polarizable force field AMOEBA. A number of different barrier-free molecular rearrangements were observed. Most of the configurations changed their shapes, forming surface defects of the same type. But the main attention is paid to the coordinated molecular reorientations leading to the proton shift along hydrogen bonds. A number of statements about the topology of the hydrogen-bonded network are rigorously proved. It is concluded that the structural stability of configurations and the features of structural transformations are basically determined by the arrangement of the homodromic hydrogen-bonded water rings. The stability of the configurations, which retained their original shape during local geometric optimization, was studied by the molecular dynamics method with a gradual increase in temperature. Some interesting effects were found. This is the stabilizing locking of the most polarized configurations, as well as very intense amplitude oscillations of the free OH group at the end of two-bond flips.

Water clusters

Hydrogen bond reorientation

Proton transfer

For many years, the interest in studying the properties of water and ice has not weakened at all [1–7]. Polyhedral water clusters (PWCs) in the form of gas hydrate cavities are very convenient for the theoretical analysis of the properties of these systems. As in bulk ice, the arrangement of hydrogen atoms (protons) follows the well-known Bernal–Fowler ice rules [8]. The number of defect-free proton configurations increases exponentially with the number of molecules. The logarithm of this number determines the residual entropy [9]. PWCs are quite well studied using a variety of computer modeling methods. Probably the most well-studied PWC is the pentagonal dodecahedron (H₂O)₂₀, which coincides in shape with the small cavity D of hydrate structures I and II [10]. For this cluster, the total number of proton configurations [11] and the number of symmetrically nonequivalent configurations [12] are precisely calculated. For PWCs, mathematical approaches [12–14] and discrete physical models [15] have been developed, which make it possible to estimate the relative energy of a particular cluster configuration, taking into account the topology of the H-bonded network, and predict the set of the most stable configurations.

An important feature of water H-bond network is its variability. The extended jump model [16, 17] for liquid water, in contrast to the Debye model of rotational diffusion, implies rather strong changes in the orientations of molecules during breaks and the formation of new H-bonds. However, the molecular picture of the reorientation mechanisms, as well as proton transport, is far from complete. Quantum tunneling of protons through H-bonds [18, 19] eventually does not lead to their breaking, but only to a change in the direction of the bonds. The interest in PWCs research may be due to the fact that they can be used for detailed study of the spontaneous reorientation of molecules without changing the bond network as a whole. This kind of rearrangements was observed in the course of geometric optimization of some configurations of PWCs using quantum chemical calculation methods [20, 21].

Here we present a systematic study of the highest-energy proton configurations of clusters in the form of gas hydrate cavities D and T, as well as in the form of Kelvin’s polyhedron (cluster K). These configurations were found using the Strong–Weak–Effective–Bond (SWEB) discrete model, which takes into account the Coulomb interaction of the second and third neighbors in the H-bond network [15]. The stability of these configurations was studied using the AMOEBA polarizable atomic multipole water model [22]. A number of barrier-free and low-barrier structural transformations which result in a proton shift along H-bonded chains were found. In addition, proton configurations that are stable despite very high energies were found. The topological features of the H-bond network were established, which allowed us to predict the location and direction of structural changes.

This article is structured as follows. The methodological part includes a brief description of the discrete SWEB model underlying the classification of proton configurations by the stabilization energy, as well as a description of the computer modeling methods used. The main part of the article includes 1) verification of the stability of proton configurations, 2) proof of rigorous statements concerning the topology of the H-bonded network, and 3) interpretation of the structural changes. The final part is devoted to the discussion of the results obtained.

The exact statistics of proton disorder in PWCs include the total number of configurations that satisfy the Bernal–Fowler ice rules [8] and determine the residual entropy [9], as well as the number of symmetrically nonequivalent configurations. For a cluster D, these numbers are 3,600,000 and 30026, respectively [10, 11]. In spite of the complex nature of the intermolecular interactions, fairly accurate regression models linking the topology of the H-bond network and energy of different configurations have been developed for PWCs [12–14]. A discrete physical model that takes into account the interaction between the first, second, and third neighbors in the H-bond network has also been developed [15]. This approximate SWEB model agrees well with the results of ab initio calculations. It allows us to divide the entire set of proton configurations into groups by the value of a single structural parameter.

It was previously established that there are only five types of H-bonds on the surface of PWCs [23]. First of all, the eclipsed H-bonds are divided into trans- and cis-configurations (Fig. 1, a, b). On the surface of PWCs, rotations of the H-bonded dimer by 120 and 240 degrees lead to the need to further divide these two configurations into six types, two of which are equivalent. The five types of H-bonds are shown in Fig. 1 (c). For t1-bond, the arrangement of electric charges is preferred not only in the pair of H-bonded molecules themselves. The Coulomb interaction in the entire fragment, including the nearest molecules, is also extremely advantageous. For cluster D, the classification of all proton configurations by the number of t1-bonds is presented in ref. [15]. 106 configurations of the lowest energy level in the SWEB model (t1 = 7) are usually of the greatest interest. According to the results of various calculations, the configurations of the lower level, in fact, turn out to be energetically most favorable. The separation of the lowest energy configurations from the second-level ones (t1 = 6) and from all other configurations can be easily seen, e.g., in figures of refs. [15, 24–26]. For this paper, we studied the properties of proton configurations that have no t1-bonds. According to the results of Kuo et al., the absence of t1-bonds in the cluster excludes self-dissociation of water molecules [27].

The Tinker molecular modeling package [28] and the non-additive potential AMOEBA [22], which takes into account the internal degrees of freedom of the molecule and atomic polarizability, were used. The stability of the highest-energy proton configurations without t1-bonds was checked during structural optimization. A number of barrier-free transitions to more stable configurations, also obeying the Bernal–Fowler rules, were discovered. The stability of the configurations that retained their original structure during optimization was further investigated using a stepwise temperature increase of 10-K interval. The average temperature was kept constant using a Berendsen's thermostat [29]. A molecular dynamics simulation at each temperature was performed at a duration of 100 ps and time step of 1 fs. Heating was performed to determine the temperature at which the structure began to change. To obtain more accurate estimates, five independent heating procedures were performed for each configuration. The temperature limit of stability of each configuration was determined by averaging the five values.

Note that the construction of the initial proton configurations is quite correct only for stable structures, when the regions of local minima are separated from each other by high barriers. For less stable configurations, neighboring minima are readily available. Therefore, the results appear to depend on the method of constructing the initial configurations. The geometric optimization was performed at an average distance of 2.84 Å between oxygen atoms in the initial structures. Surprisingly, for a distance of 2.75 Å, very similar results were obtained regarding the stability of the configurations. Yet, there was a dependence of our results on the initial coordinates of all atoms. Given this circumstance, as well as the fact that the potential used was approximate and introduced an additional error, we were primarily focused on identifying the main qualitative regularities.

Structural stability check for cluster D

For cluster D, the distribution of the energies and total dipole moments of proton configurations belonging to the highest energy level in the SWEB model (t1 = 0) is shown in Fig. 2. Initially, all 94 configurations were built according to the Bernal–Fowler ice rules and then geometrically optimized with AMOEBA potential. The blue color corresponds to complete geometric optimization, and the red color corresponds to the first few steps of optimization. It can be seen that after the initial optimization stage, the distribution had a typical, approximately parabolic shape, reflecting the influence of the dipole moment on the energy of the PWC. The distributions for the classes t1 = 7, 6, etc. had exactly the same parabolic form [15, 24–26]. The distribution of fully optimized configurations naturally moved down, but had a different form (Fig. 2). Configurations with average values of the total dipole moment shifted much more with a noticeable decrease in the dipole moment. It turned out that during the local optimization, all configurations, with the exception of the configurations with the minimum and, surprisingly, maximum dipole moments, significantly changed their structures. The fact that most of the weakest configurations of this cluster changed in the course of geometric optimization seemed quite natural. Interestingly, in some cases, the structural rearrangement exactly corresponded to the transfer of protons along the H-bonds and led to a new configuration, which also obeyed the Bernal–Fowler ice rules. It is also interesting that all transformations could be reduced to simple strict rules.

Homodromic faces

The analysis of structural transformations in the course of energy minimization allowed us to draw a conclusion about the completely new, key role of faces with homodromic (clockwise/counterclockwise) circular arrangement of H-bonds. It was not difficult to make sure that three homodromic faces could not be adjacent, so their number was bounded from above. There was also a nontrivial lower bound due to the fact that, in any PWC, the minimum number of the homodromic faces is equal to two [30]. The proof relies on Stallings' lemma [31]. According to this lemma, in any grid of one-way city roads and in the absence of intersections, at which all roads converge or diverge, there is always a block that can be bypassed in one direction. Any closed path along a unidirectional H-bonded chain divides the H-bond network of a PWC into two parts, for each of which the Stallings' lemma applies. It turned out that for the weakest configurations within the SWEB model (i.e., at t1 = 0) some simple mathematics made it possible to significantly simplify the understanding of the properties of these configurations. We emphasize that all of the following statements refer precisely to the class of configurations without a single t1-bond.

Statement 1. All molecules of homodromic faces have the same number of dangling hydrogen atoms, either 0 or 1, i.e. all such faces are either absolutely smooth (H0-faces) or completely equipped with hydrogen atoms (H1-faces). Two homodromic faces are shown schematically in Fig. 3a, b. It is easy to make sure that with any choice of the direction of external bonds #1, the direction of bonds #2 (as well as of all other external H-bonds) should be exactly the same to avoid the appearance of t1-bonds (see Fig. 1). Recall that in Fig. 1, the arrows indicate the direction of the H-bond: from the donor to the proton acceptor.

Statement 2 (inverse). If the face is H0 or H1, then it is necessarily homodromic. To prove this, let us assume the opposite that the face is not homodromic. Then there are at least two vertices (molecules) with different numbers of dangling hydrogens (Fig. 3c), i.e. the face is neither of type H0 nor H1.

Statement 3. Two homodromic faces cannot be adjacent, because they form a t1-bond (see Fig. 1).

As noted, in any PWC the minimum number of homodromic faces is two. The question arises: what is the number of homodromic faces at t1 = 0? We already know that, in this case, homodromic faces cannot be adjacent. At the same time, all these faces are H0 or H1. In cluster D, the maximum number of non-adjacent faces is equal to three (Fig. 3d). If all these faces are homodromic, then two of them are inevitably H0- or H1-faces. Hence, two neighboring vertices (molecules) belonging to these two faces have the same number of dangling hydrogens. But this is impossible, because the intermediate H-bond excludes the similarity of these vertices. Therefore, in cluster D the number of homodromic faces is equal to two at t1 = 0, since it cannot be less than two. This proves the following statement, which is also true for t1 = 0.

Statement 4. In cluster D, the number of homodromic faces at t1 = 0 is exactly equal to two.

Statement 5. If the homodromic faces are connected by one intermediate bond (not located opposite to each other in the polyhedron), then one of them is H1- and the other is H0-face. The reason is that one of the molecules of this H-bond is an acceptor of one hydrogen atom and the second of two hydrogen atoms (protons).

Statement 6; If in cluster D the homodromic faces are located opposite each other and both are of the same type (H0 or H1), then the remaining molecules forming a single H-bonded cycle have the opposite number of free hydrogens (1 or 0). This follows from the fact that in any three-coordinated PWC, the number of vertices with and without dangling bonds is the same [12, 13].

Structural transformations in cluster D

In cluster D, the structure of the initial configurations with t1 = 0 was not too diverse. It was characterized by the arrangement of two homodromic faces. In 68 configurations, the homodromic faces were separated by one H-bond. In the remaining 26 configurations, the homodromic faces were located opposite each other. Eighteen of these 26 configurations were asymmetrical (C₁). A separate group consisted of 8 highly symmetric configurations with the C₅ axes. Let us first consider 68 configurations with non-opposite arrangement of homodromic faces. All of them changed during the local optimization. One of their faces was H1, i.e. dangling hydrogen atoms were arranged in a pentagonal ring. Of course, this was energetically disadvantageous. Note in advance that this hydrogen atom configuration was the main driver of all structural transformations.

There is a need for clarity here. C1-bonds with nearest dangling hydrogen atoms are considered as the main factor that reduces the binding energy of PWCs [12]. But this reflects the disadvantage of the interaction of only the nearest H-bonded molecules. In PWCs, the other two variants of cis-configuration are equally unprofitable (Fig. 1c). Electrostatic interactions are long-ranged and constitute the main contribution to the intermolecular interaction. Taking into account the interaction of the second and third neighbors, the energy of PWCs is determined basically by the number of t1-bonds [15]. The suitability of the criterion c1 ◊ min for structural optimization is due to a purely combinatorial equality c1 + t1 = N/2, where N is the number of molecules in the PWC [13]. The destruction of configuration is caused more by the presence of extremely unfavorably located molecules than by the low stabilization energy as a whole. There are configurations with significantly lower stabilization energy that are still structurally stable (see below).

In the course of geometric optimization of the 68 configurations with non-opposite arrangement of homodromic faces, the following was established. The H1-face molecule participating in the formation of an intermediate H-bond between the homodromic faces always started to turn first. This is due to the fact that a strong electric field was created at the place of this H-bond, the direction of which is shown by the arrow in Fig. 4a. This field was created by the molecules of the homodromic faces themselves, since the dipoles of all adjacent bonds were directed mainly downwards, and the sum of the dipoles inside the faces was zero. This field was further enhanced by the two molecules shown in Fig. 4 as an oval. Due to the Bernal–Fouler ice rules, two hydrogen atoms of these molecules were always directed downwards. The "electric breakdown" (concerted rotation of molecules) occurred, which could be carried out in different ways.

During the rotation, the free hydrogen atom of the H1-face pushed the hydrogen atom lying on the H-bond inside the polyhedron. Most often, the second atom then formed an additional H-bond inside the second homodromic face (Fig. 4b, central dotted line). In this case, the pentagonal H0-face turned into an adjacent square and triangle. This transformation is also typical for an isolated cyclic pentamer [32] and for other PWCs [33]. Following [32], this transformation will be denoted by sqtr-structure or sqtr-defect.

In the course of local optimization, a pure two-bond flip was observed in four cases with non-opposite arrangement of homodromic faces. The H-bond connecting two homodromic faces was always flipped first. Under the action of a strong electric field, the free hydrogen atom of the first molecule turned and pushed the hydrogen atom of the second molecule inside the polyhedron. The hydrogen atom of the second molecule, immersed inside the polyhedron, turned around and pushed out the hydrogen atom of the third molecule, which had previously formed the second H-bond. In this case, the original polyhedral shape of the cluster was preserved. In fact, this was a barrier-free proton transfer along two adjacent H-bonds without tunneling and without dissociation of the molecules. The molecules rotated concertedly, although the molecule of the H1-face started moving first. Figure 5a shows the time-dependence of the rotation angles of these three molecules. In reality, the maximum values of the rotation angles were slightly different from each other. The O-H bond also never strictly coincided with the O-O direction. For this reason, the values in Fig. 5a are normalized to the tetrahedral angle: arccos (-1/3) = 109.47⁰. This process can be seen in Movie S1.

When homodromic faces were placed opposite each other, the two-bond flips were observed more often. Of the 18 configurations with C1 symmetry, such flips were observed in 12. In three cases, a pure two-bond flip occurred. In 9 configurations, the two-bond flips were accompanied by additional rotation of some molecules. In the remaining configurations, surface sqtr-defects occurred. These defects were accompanied by the appearance of additional H-bonds as well as preferred t1-bonds, and some general deformation. The H1-face was still the driver of the structural rearrangement. But here, it was more difficult to predict in advance the reconstructed fragment of the cluster due to the symmetrical arrangement of the homodromic faces. In total, a pure two-bond flip was observed in seven configurations with both non-opposite (3 configurations) and opposite (4 configurations) arrangement of homodromic faces. This was accompanied by a significant increase in the stabilization energy, because t1-bonds were formed at the beginning and at the ends of the H-bonded chains (see below).

Hydrogen bond reorientations in clusters T and K.

This section presents the results of the study of barrier-free reorientations in two PWCs with the number of molecules N = 24. The first cluster Т has the form of a large cavity of the hydrate structures sI [10]. Along with 12 pentagonal faces, it forms two hexagonal faces, which are located opposite each other (Fig. 6a). The second cluster K has the shape of Kelvin's polyhedron (Fig. 6b). The peculiarity of the second cluster is that it does not contain a single pentagonal face. This topological feature determines its advantage over polyhedra D and T, because the maximum fraction of t1-bonds is noticeably higher (1/3) [13] and, as a consequence, the stability of the most energetically preferred proton configurations is higher [33]. We relied on the Statements presented above, of which only 4 and 6 relate exclusively to cluster D.

It was found that the reorientation in clusters T and K is also basically determined by the location of the homodromic faces. However, the analysis of these reorientations in these clusters becomes more complicated due to the fact that the number of homodromic faces may be greater. The total numbers of defect-free proton configurations in T and K clusters are 73041408 and 75400162, respectively [34]. At t1 = 0, the number of configurations is 65312 and 62144, of which 2739 and 1343 are symmetrically distinct (non-isomorphic). It was established that the numbers of non-isomorphic arrangements of homodromic faces are 20 and 16. In Fig. 7, faces H0 and H1 are shown in black and gray, respectively. Each coloring corresponds to a class of proton configurations with the same arrangement of homodromic faces, but differing in the direction of H-bonds in the cluster as a whole. The number of configurations in each class is shown at the diagrams. The number of configurations in classes II-5 and II-6 for cluster T (as well as some other classes of both clusters) is the same, because changing the direction of all H-bonds transforms configurations of one class to another, and vice versa.

In the course of the local geometric optimization with the AMOEBA potential, it was found that the reorientation of H-bonds without the polyhedral shape breaking occurs only in a small number of configurations. In clusters T and K, the flips were observed only in 82 and 57 configurations out of the total number of configurations 2739 and 1343, respectively (see above). The basic mechanism of concerted H-bond reorientation is the same as in cluster D. This is the pushing one of the OH-groups from H1-face and the sequential rotation of three molecules around H-bonds, which results in a change in the direction of two H-bonds. However, in clusters T and K, the distance between the homodromic faces is various. For the flips between the homodromic faces, we used a special designation, which, along with the number of reversed bonds, also takes into account the distance between these faces. So, in cluster D the flips between the homodromic faces with one intermediate H-bond are 2/1-flips. In clusters T and K, 2/2-flips are typical for configurations in which the distance between these faces is 2. In this case, only intermediate bonds are flipped, unlike the 2/1-flips. The designation 2/3 is used for two-bond flips located between faces H0 and H1 separated by three bonds.

Table 1

Characteristics of hydrogen bond flips in cluster T
Class of configurations	II-1	II-2	II-3	II-4	II-7	II-10
Number of Flips	7	16	33	20	3	3
Flip type	2/1	2/1	2/2	2/3	2/1	2/3
Standard flips, %	100	44	76	45	67	100
Complemented flips, %	100	81	100	85	100	100

A distribution of H-bond flips by the configuration classes, taking into account the homodromic faces arrangement, is presented in Tables 1 and 2. In cluster T, flips were found only in configurations with two homodromic faces. 2/1-flips were observed in configurations of classes II-1, II-2 and II-7. The third and fourth rows in Table 1 shows the fraction of single standard flips of the given type and the fraction of standard flips accompanied by any additional flip. The mechanism of the additional flips is the same. Only the pushing of a free OH-group does not occur from fully equipped H1-face, but from a face with some number of OH-groups. In rare cases, there were flips of H-bonds located not between homodromic faces. Naturally, any additional flip further stabilize the cluster.

It is easy to see that the flips in cluster T indeed due to the arrangement of homodromic faces. A smaller fraction of standard (predictable) flips in class II-4 is associated with a very distant and asymmetrical arrangement of homodromic faces. With a symmetrical, maximally distant location of homodromic faces in class II-10, all three 2/3-flips are oriented from H1- to H0-face. In class II-4, some unusual flips were observed. Pairs of simultaneous 2/3-flips were detected three times. The absolute exception is an extremely atypical 3/3-flip, in which three intermediate H-bonds between the homodromic faces are flipped. In this case, the rotations of the molecules were more complex.

A very similar picture was observed in cluster K (Table 2), although here flips were also observed in configurations with three homodromic faces. The maximum number of flips was observed in the largest class II-1 (see Fig. 7b). All of them were single standard 2/1-flips. In classes II-5 and II-6, most 2/2-flips were also standard, although sometimes they were accompanied by additional flips. Both reorientations in class II-2 are 2/2-flips, accompanied by additional flips of two bonds. In classes III-1, III-2 and III-3, standard 2/1-flips were also observed. In the first two of these classes, homodromic squares are arranged in a semicircle. Three of the six reorientations in class III-2 deserve special attention. These are converging 2/1-flips directed from opposite H1-bases of the polyhedron to the same H0-square located at the equator. In this case, the total dipole moment of the cluster has changed slightly and only in the perpendicular direction. The convergent direction of the flips towards each other is counterintuitive and is due to the local features of the electrostatic field (Movie S2). The only reorientation in class III-3 is interesting because it combines two standard flips: 2/1 and 2/2.

Table 2

Characteristics of hydrogen bond flips in cluster K
Class of configurations	II-1	II-2	II-5	II-6	III-1	III-2	III-3
Number of Flips	24	2	13	7	4	6	1
Flip type	2/1	2/4	2/2	2/2	2/1	2/1	2/1
Standard flips, %	100	0	85	29	100	100	100
Complemented flips, %	100	100	100	86	100	100	100

In most reorientations, the total dipole moment significantly decreased, and the number of t1-bonds increased, i.e., the stabilization energy of cluster increased. It is noteworthy that as a rule the number of t1-bonds increased to at least two. The appearance of only one t1-bond was noted only once in cluster D with a non-standard flip. In this case, two flipped bonds that are not located between homodromic faces turned over. For all optimized configurations after standard 2/1 and 2/2-flips, strict equality t1 = 2 is performed. In our opinion, this is an interesting manifestation of the topological properties of PWCs. Most often there were exactly two t1-bonds. With multiple flips, 3, 4 and even 5 t1-bonds appeared. Five bonds were found after triple two-bond flips. This was observed both in cluster T (three reorientations in classes II-3 and II-4) and in cluster K (both reorientations of class II-2).

Thermally induced structural transformations of stable configurations

Eight highly symmetric configurations of cluster D with opposite arrangement of the homodromic faces formed a completely separate group. Four of them had C₅ symmetry and the other four had S₁₀ symmetry. In this case, the homodromic faces could be different or the same (H0/H1). These eight configurations can be divided into three groups according to the types of faces (Fig. 8). Within these groups, the configurations differed in the direction of H-bonded cycles: clockwise or counterclockwise. Different directions of H-bonds within the three cycles (upper and lower pentagonal faces as well as an equatorial line) are shown at the bottom. The third group had only two options, since H-bonds along the equator were alternately oriented. All of these eight configurations retained their structure during local optimization. Further, their stability was tested with a gradual increase in temperature (by 10 K) in the course of molecular dynamics simulation. The simulation was run for 100 ps at each temperature, with an integration step size of 1 fs. To reduce random effects, the heating of each configuration was repeated five times.

Table 3

The temperature limits of structural stability (°K) for eight configurations of cluster D with opposite homodromic faces (Fig. 8) for five realizations of the heating process.
No.	1	2	3	4	5	6	7	8
Average	32	40	36	30	32	22	20	10

Table 3 shows the average temperature at which the original configuration changed. Configurations no. 1–3 with opposite H0-faces were the most stable. Configurations no. 4–6 with opposite H1-faces were slightly less stable. Configurations 7 and 8 with different homodromic faces were the least stable. For the last two configurations, even the mechanical stability at T = 0 was quite unexpected, given the highest values of the energy and dipole moment (Fig. 2). The stability was largely due to the symmetry of these configurations. In addition, in this case, the hydrogen atoms of the H1-face were locked by the hydrogen atoms located just below. This caused the self-locking effect. Three two-bond flips occurred during the heating of configuration no. 8. Two of these flips were closely related. Figure 6b shows the rotation angles of the three molecules of the first two-bond flip. Due to the natural spread, the maximum rotation angles were not equal to each other. The first two molecules started moving almost simultaneously, because in this case the H1-face molecules were locked and could not start independently. In Fig. 5b, one can see that the third molecule started moving later. The Berendsen thermostat was used with a coupling constant of 0.1 ps. During the rapid structural changes, the temperature did not have time to stabilize. Therefore, when the potential energy was released, there was a natural noticeable increase in temperature (kinetic energy). In the final stage of the two-bond flip, the kinetic energy was largely localized in the last hydrogen atom. This created a whip-like effect in the H-bonded chain. The hydrogen atom of the last molecule moved more quickly. These motions can be seen in Movie S3. According to recent results obtained using ultrafast spectroscopy, molecules on the surface of water lose energy in a similar way [35].

Thermally induced structural transformations in configurations no. 1–3 with two oppositely located H0-faces (Fig. 8a) almost always resulted in the formation of one or two sqtr-defects. Sometimes, a pair of sqtr-defects appeared simultaneously. A very similar pattern was observed when heating configurations with homodromic H1-faces.

In both clusters T and K, locally stable configurations were also found. Of eleven stable configurations in cluster T, three configurations of classes II-9 and II-11 differ in the direction of H-bonds in homodromic bases and along the equator. There was one maximally polarized configuration of class II-10 with the effect of self-locking. In addition, four configurations of class II-6 turned out to be stable. In cluster K, there were eight stable configurations with only opposite homodromic squares. Half of them belong to class II-3, and the other half belongs to class II-4. Similar configurations of classes II-9 and II-10, as well as maximally polarized configurations of classes II-2 and II-7 with different opposite faces, turned out to be unstable. In clusters T and K, the general picture of the loss of stability during heating, including the spontaneous occurrence of H-bond flips and the manifestation of whip effects, was exactly the same as in cluster D.

As a rule, in computer simulation of molecular systems, only the most stable structures are of interest. At the same time, the least stable configurations of ice-like systems are of particular interest because of the possibility of observing consistent water molecule reorientations. Spontaneous reorientations in PWCs lead to the displacement of the proton along the H-bond chains, as in quantum tunneling. The molecular dynamic trajectory of the system during structural rearrangement is very similar to the dominant tunneling path (instanton) [19]. Both of these trajectories provide a better understanding of the coherent collective nature of the process. Of no less interest are the thermo-induced low-barrier transformations, which also lead to the shift of protons along the H-bond chains. Note that the method of obtaining the proton configurations 7 and 8 in cluster D with the maximum dipole moment is probably the simplest. In our study, we did not take into account possible ionization of molecules. Although the above noted whip-like movement of H-bonded molecules perhaps, deserves a more accurate analysis, taking into account formation of t1-bonds, which can themselves lead to self-dissociation of water molecule [27].

According to the discrete SWEB model, the highest energy configurations have very specific topological properties. This greatly facilitated the interpretation of the results obtained. We emphasize that the established regularities can be easily transferred to other water polyhedra. So, in a triangular water prism as well as in a cubic water cluster, there are two highest-energy configurations with t1 = 0. It is easy to make sure that they are characterized by the opposite arrangement of homodromic faces [36]. One of them is a H0-face and another is a H1-face because they share H-bonds. Calculations with the TIP4P and AMOEBA potentials showed that these configurations are stable.

Supplementary Information The online version contains supplementary material available at https://doi.org/xxxxxxxxxxxxx. Three movie files animating trajectories for the barrier-free and low-barrier concerted flips of the hydrogen bonds.

Acknowledgments This research was carried out according to state assignment No. 121041600040-3 and was supported by the Russian Science Foundation [grant number 22-23-00092].

Availability of data and material The data that support the finding of this study are available within the article.

Author contribution There were equal contributions of the authors to the completion of this work.

Code availability TINKER Software Tools for Molecular Design. Version 6.2, VMD 1.9.3.

Competing interests The authors declare no competing interests.

Debenedetti PG (2003) Supercooled and glassy water. J Phys: Condens Matter 15:R1669–R1726. https://doi.org/10.1063/1.1595053.
Smith JD, Cappa CD, Wilson KR, Messer BM, Cohen RC, Saykally RC (2004) Energetics of hydrogen bond network rearrangements in liquid water. Science 306:851–853.https://doi.org/10.1126/science.1102560.
Malenkov G (2009) Liquid water and ices: understanding the structure and physical properties. J Phys: Condens Matter 21:283101. https://doi.org/10.1088/0953-8984/21/28/283101.
Gallo P, Arnann-Winkel K, Angell CA et al (2016) Water: A Tale of Two Liquids. Chem Rev 116:7463–7500. https://doi.org/10.1021/acs.jpca.1c08020.
Shephard JJ, Slater B, Harvey P, Hart M, Bull CL, Bramwell ST, Salzmann CG (2018) Doping-induced disappearance of ice II from water’s phase diagram. Nature Physics 14:569–572. https://doi.org/10.1038/s41567-018-0094-z.
Gohlke M, Moessner R, Pollmann F (2019) Polarization plateaus in hexagonal water ice Ih. Phys Rev B 100:014206. https://doi.org/10.1073/pnas.2018837118.
Drechsel-Grau C, Marx D (2014) Quantum Simulation of Collective Proton Tunneling in Hexagonal Ice Crystals. Phys Rev Lett 112:148302. https://doi.org/10.1103/PhysRevLett.112.148302.
Bernal JD, Fowler RH (1933) A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions. J Chem Phys 1:515–548. https://doi.org/10.1063/1.1749327.
Pauling L (1935) The structure and entropy of ice and of other crystals with some randomness of atomic arrangement. J Am Chem Soc 57:2680–2684. https://doi.org/10.1021/ja01315a102.
Sloan ED Jr (1998) Clathrate Hydrates of Natural Gases; Marcel Dekker, New York.
Kirov MV (1994) Residual entropy of polyhedral water clusters. Exact relations. J Struct Chem 35:126–128. https://doi.org/10.1007/BF02578511.
McDonald S, Ojamae L, Singer SJ (1998) Graph theoretical generation and analysis of hydrogen-bonded structures with applications to the neutral and protonated water cube and dodecahedral clusters. J Phys Chem 102:2824–2832. https://doi.org/10.1021/jp9803539.
Anick DJ (2002) Polyhedral water clusters, I: Formal consequences of the ice rules. THEOCHEM 587:87–96. https://doi.org/10.1016/S0166-1280(02)00101-X.
Anick DJ (2002) Polyhedral water clusters, II: correlations of connectivity parameters with electronic energy and hydrogen bond lengths. THEOCHEM 587:97–110. https://doi.org/10.1016/S0166-1280(02)00100-8.
Kirov MV, Fanourgakis GS, Xantheas SS (2008) Identifying the most stable networks in polyhedral water clusters. Chem Phys Lett 461:180–188. https://doi.org/10.1016/j.cplett.2008.04.079.
Laage D, Hynes JT (2006) A molecular jump mechanism of water reorientation Science 311:832–835. https://doi.org/10.1126/science.1122154.
Leone SR, Cremer PS, Groves JT, Johnson MA (2011) Reorientation and Allied Dynamics in Water and Aqueous Solutions. Annu Rev Phys Chem 62:395–416. https://doi.org/10.1146/annurev.physchem.012809.103503.
Hassanali A, Giberti F, Cuny J, Kühne TD (2013) Parrinello, M. Proton transfer through the water gossamer. Proc Natl Acad Sci 110:13723–13728. https://doi.org/10.1073/pnas.1306642110.
Vaillant CL, Wales DJ, Althorpe SC (2019) Tunneling Splittings in Water Clusters from Path Integral Molecular Dynamics. J Phys Chem Lett 10:7300–7304. https://doi.org/10.1021/acs.jpclett.9b02951.
Anick DJ (2003) Application of database methods to the prediction of B3LYP-optimized polyhedral water cluster geometries and electronic energies. J Chem Phys 119:12442. https://doi.org/10.1063/1.1625631.
Chihaia V, Adams S, Kuhs WF (2004) Influence of water molecules arrangement on structure and stability of 5¹² and 5¹²6² Buckyball Water Clusters. A theoretical study. Chem Phys 297:271–287. https://doi.org/10.1016/j.chemphys.2003.10.032.
Ren P, Ponder JW (2003) Polarizable atomic multipole water model for molecular mechanics simulation. J Phys Chem B 107:5933–5947. https://doi.org/10.1021/jp027815+.
Kirov MV (1996) Conformational combinatorial analysis of polyhedral water clusters J Struct Chem 37:84–91. https://doi.org/10.1007/BF02578574.
Anick DJ (2010) Topology-energy relationships and lowest energy configurations for pentagonal dodecahedral (H₂O)₂₀X clusters, X = empty, H2O, NH3, H3O+: The importance of O-topology. J Chem Phys 132:164311. https://doi.org/10.1063/1.3397812.
Parkkinen P, Riikonen S, Halonen L (2013) (H₂O)₂₀ water clusters at finite temperatures. J Phys Chem A 117:9985–9998. https://doi.org/10.1021/jp4003092.
Iwata S, Akase D, Aida M, Xantheas SS (2016) Electronic origin of the dependence of hydrogen bond strengths on nearest-neighbor and next-nearest-neighbor hydrogen bonds in polyhedral water clusters (H₂O)_n, n = 8, 20 and 24. Phys Chem Chem Phys 18:19746–19756. https://doi.org/10.1039/C6CP02487D.
Kuo JL, Ciobanu CV, Ojamae L, Shavitt I, Singer SJ (2003) Short H-bonds and spontaneous self-dissociation in (H₂O)₂₀: Effects of H-bond topology. J Chem Phys 118:3583–3588. https://doi.org/10.1063/1.1538240.
Ponder JW (2012) TINKER Software Tools for Molecular Design. Version 6.2. Washington University School of Medicine, Saint Louis, MO. http://dasher.wustl.edu/tinker. Accessed 6 September 2021
Berendsen HJC, Postma JPM,Van Gunsteren WF, DiNola A, Haak JR (1984) Molecular dynamics with coupling to an external bath. J Chem Phys 81:3684–3690. https://doi.org/10.1063/1.448118.
Kirov MV (2002) Sequential Orientation of Cycles in Gas Hydrate Frameworks. J Struct Chem 43:274–283. https://doi.org/10.1023/a:1019600524366.
Stallings JR (1987) In: Combinational Group Theory and Topology, ed. S. M. Gersten, J. R. Stallings. Princeton Univ. Press, Princeton.
Wales DJ, Walsh TR (1996) Theoretical study of the water pentamer. J Chem Phys 105:6957–6971. https://doi.org/10.1063/1.471987.
Gudkovskikh SV, Kirov MV (2021) Thermal stability of water polyhedra with square faces. J Mol Model 27:366. https://doi.org/10.1007/s00894-021-04996-7.
Kirov MV (2020) Energetics of Water Polyhedra with Square Faces. J Phys Chem A 124:4463–4470. https://doi.org/10.1021/acs.jpca.0c02835.
Inoue KI, Ahmed M, Nihonyanagi S et al (2020) Reorientation-induced relaxation of free OH at the air/water interface revealed by ultrafast heterodyne-detected nonlinear spectroscopy. Nat Commun 11:5344. https://doi.org/10.1038/s41467-020-19143-8.
Kirov MV (2016) Classification of hydrogen bond flips in small water polyhedra applied to concerted proton tunneling. Phys Chem Chem Phys 18:27351–27357. https://doi.org/10.1039/C6CP04960E.

Download PDF

Version 1

posted

You are reading this latest preprint version

Barrier-free molecular reorientations in polyhedral water clusters

Status:

Version 1

Abstract

Figures

Introduction

Method

Results And Discussion

Structural stability check for cluster D

Homodromic faces

Structural transformations in cluster D

Conclusions

Declarations

References

Supplementary Files

Status:

Version 1