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 Hbonds 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 Hbonds. 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 oneway 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 Hbonded chain divides the Hbond 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 t1bond.
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 (H0faces) or completely equipped with hydrogen atoms (H1faces). 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 Hbonds) should be exactly the same to avoid the appearance of t1bonds (see Fig. 1). Recall that in Fig. 1, the arrows indicate the direction of the Hbond: 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 t1bond (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 nonadjacent faces is equal to three (Fig. 3d). If all these faces are homodromic, then two of them are inevitably H0 or H1faces. Hence, two neighboring vertices (molecules) belonging to these two faces have the same number of dangling hydrogens. But this is impossible, because the intermediate Hbond 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 H0face. The reason is that one of the molecules of this Hbond 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 Hbonded cycle have the opposite number of free hydrogens (1 or 0). This follows from the fact that in any threecoordinated 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 Hbond. In the remaining 26 configurations, the homodromic faces were located opposite each other. Eighteen of these 26 configurations were asymmetrical (C1). A separate group consisted of 8 highly symmetric configurations with the C5 axes. Let us first consider 68 configurations with nonopposite 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. C1bonds 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 Hbonded molecules. In PWCs, the other two variants of cisconfiguration are equally unprofitable (Fig. 1c). Electrostatic interactions are longranged 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 t1bonds [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 nonopposite arrangement of homodromic faces, the following was established. The H1face molecule participating in the formation of an intermediate Hbond 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 Hbond, 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 H1face pushed the hydrogen atom lying on the Hbond inside the polyhedron. Most often, the second atom then formed an additional Hbond inside the second homodromic face (Fig. 4b, central dotted line). In this case, the pentagonal H0face 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 sqtrstructure or sqtrdefect.
In the course of local optimization, a pure twobond flip was observed in four cases with nonopposite arrangement of homodromic faces. The Hbond 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 Hbond. In this case, the original polyhedral shape of the cluster was preserved. In fact, this was a barrierfree proton transfer along two adjacent Hbonds without tunneling and without dissociation of the molecules. The molecules rotated concertedly, although the molecule of the H1face started moving first. Figure 5a shows the timedependence of the rotation angles of these three molecules. In reality, the maximum values of the rotation angles were slightly different from each other. The OH bond also never strictly coincided with the OO 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 twobond flips were observed more often. Of the 18 configurations with C1 symmetry, such flips were observed in 12. In three cases, a pure twobond flip occurred. In 9 configurations, the twobond flips were accompanied by additional rotation of some molecules. In the remaining configurations, surface sqtrdefects occurred. These defects were accompanied by the appearance of additional Hbonds as well as preferred t1bonds, and some general deformation. The H1face 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 twobond flip was observed in seven configurations with both nonopposite (3 configurations) and opposite (4 configurations) arrangement of homodromic faces. This was accompanied by a significant increase in the stabilization energy, because t1bonds were formed at the beginning and at the ends of the Hbonded chains (see below).
Hydrogen bond reorientations in clusters T and K.
This section presents the results of the study of barrierfree 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 t1bonds 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 defectfree 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 (nonisomorphic). It was established that the numbers of nonisomorphic 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 Hbonds in the cluster as a whole. The number of configurations in each class is shown at the diagrams. The number of configurations in classes II5 and II6 for cluster T (as well as some other classes of both clusters) is the same, because changing the direction of all Hbonds 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 Hbonds 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 Hbond reorientation is the same as in cluster D. This is the pushing one of the OHgroups from H1face and the sequential rotation of three molecules around Hbonds, which results in a change in the direction of two Hbonds. 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 Hbond are 2/1flips. In clusters T and K, 2/2flips are typical for configurations in which the distance between these faces is 2. In this case, only intermediate bonds are flipped, unlike the 2/1flips. The designation 2/3 is used for twobond flips located between faces H0 and H1 separated by three bonds.
Table 1
Characteristics of hydrogen bond flips in cluster T
Class of configurations

II1

II2

II3

II4

II7

II10

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 Hbond 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/1flips were observed in configurations of classes II1, II2 and II7. 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 OHgroup does not occur from fully equipped H1face, but from a face with some number of OHgroups. In rare cases, there were flips of Hbonds 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 II4 is associated with a very distant and asymmetrical arrangement of homodromic faces. With a symmetrical, maximally distant location of homodromic faces in class II10, all three 2/3flips are oriented from H1 to H0face. In class II4, some unusual flips were observed. Pairs of simultaneous 2/3flips were detected three times. The absolute exception is an extremely atypical 3/3flip, in which three intermediate Hbonds 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 II1 (see Fig. 7b). All of them were single standard 2/1flips. In classes II5 and II6, most 2/2flips were also standard, although sometimes they were accompanied by additional flips. Both reorientations in class II2 are 2/2flips, accompanied by additional flips of two bonds. In classes III1, III2 and III3, standard 2/1flips were also observed. In the first two of these classes, homodromic squares are arranged in a semicircle. Three of the six reorientations in class III2 deserve special attention. These are converging 2/1flips directed from opposite H1bases of the polyhedron to the same H0square 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 III3 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

II1

II2

II5

II6

III1

III2

III3

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 t1bonds increased, i.e., the stabilization energy of cluster increased. It is noteworthy that as a rule the number of t1bonds increased to at least two. The appearance of only one t1bond was noted only once in cluster D with a nonstandard 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/2flips, 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 t1bonds. With multiple flips, 3, 4 and even 5 t1bonds appeared. Five bonds were found after triple twobond flips. This was observed both in cluster T (three reorientations in classes II3 and II4) and in cluster K (both reorientations of class II2).
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 C5 symmetry and the other four had S10 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 Hbonded cycles: clockwise or counterclockwise. Different directions of Hbonds 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 Hbonds 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 H0faces were the most stable. Configurations no. 4–6 with opposite H1faces 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 H1face were locked by the hydrogen atoms located just below. This caused the selflocking effect. Three twobond 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 twobond 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 H1face 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 twobond flip, the kinetic energy was largely localized in the last hydrogen atom. This created a whiplike effect in the Hbonded 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 H0faces (Fig. 8a) almost always resulted in the formation of one or two sqtrdefects. Sometimes, a pair of sqtrdefects appeared simultaneously. A very similar pattern was observed when heating configurations with homodromic H1faces.
In both clusters T and K, locally stable configurations were also found. Of eleven stable configurations in cluster T, three configurations of classes II9 and II11 differ in the direction of Hbonds in homodromic bases and along the equator. There was one maximally polarized configuration of class II10 with the effect of selflocking. In addition, four configurations of class II6 turned out to be stable. In cluster K, there were eight stable configurations with only opposite homodromic squares. Half of them belong to class II3, and the other half belongs to class II4. Similar configurations of classes II9 and II10, as well as maximally polarized configurations of classes II2 and II7 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 Hbond flips and the manifestation of whip effects, was exactly the same as in cluster D.