Thermal stability of water polyhedra with square faces

The ability to form numerous crystalline modifications of ice and gas hydrate frameworks is a characteristic feature of water. In fact, this structural variety is much wider due to the proton disorder. Configurations with different arrangements of hydrogen atoms (protons) in hydrogen bonds are not equivalent in their properties. Polyhedral water clusters are convenient objects for studying the effect of proton disorder on the properties of ice-like systems. It was previously established that the stability of water polyhedra is determined by the competition of two factors. The geometric factor gives preference to tetrahedrally coordinated structures with a large number of pentagonal faces. The topological factor takes into account the number of energetically most favorable types of H-bonds. This number increases with the number of square faces. It was found that tetrahedrally coordinated structures are not the most stable. However, these calculations were performed without taking thermal effects into account (Kirov M. V., J Phys Chem A, 124:4463 − 4470, 2020). The purpose of the present article is to study the structural stability of various water polyhedra at different temperatures. In the course of modeling, using the Amoeba force field, the advantage of configurations with a large number of square faces is demonstrated. The structure and energetics of surface defects are studied. Several very stable structures of unusual shape were found, including polyhedra which contain 4-coordinated molecules and polyhedra whose O–H groups are directed to the cluster center. The comparative analysis of cluster stability includes the temperature intervals of melting-like transitions.


Introduction
Clathrate hydrates attract many researchers regarding the problems of global climate change and green energy. The large number of different gas hydrate frameworks is a continuation of the variety of ice forms. In the most common clathrate hydrates, many cavities have the shape of a pentagonal dodecahedron [1]. Most of the faces of the other cavities also form pentagons, which correspond well to the tetrahedral coordination of the hydrogen bond (H-bond) network. However, both ice and gas hydrate frameworks are also characterized by proton disorder [2]. The number of configurations that differ only in the location of protons in H-bonds (proton configurations) is huge [3,4]. An important circumstance is that these configurations differ from each other in energy and other characteristics [5]. Understanding the energetics of proton configurations of clathrate frameworks and ice can be of considerable scientific and practical interest.
It was found that the stability of water polyhedra is determined by the influence of two main competing factors [6,7]. The geometric factor considers the degree of deviation of the H-bond network from tetrahedral coordination. The topological factor considers the number of energetically more preferable fragments of the H-bond network and the specific location of protons in H-bonds. The fragment that determines the energy of water polyhedra is one of the configurations of an H-bonded pair of molecules. The number of types of such pairs (types of H-bonds) on the polyhedron surface is five. The discrete model of strong and weak effective H-bonds (the SWEB model), which takes into account the interaction of the first, second, and third neighbors in the network, is a physical justification for the advantage of one type of H-bond [8]. The competition between the two stability factors is that, on one hand, structures with tetrahedral geometry and with a large number of pentagonal faces should be more stable. On the other hand, structures with square faces have a greater number of more preferred H-bonds.
For the set of all configurations allowed by the Bernal-Fowler ice rules [3], the geometric stability factor is dominant. That is why pentagonal dodecahedra and pentagonal faces prevail in the structure of clathrate frameworks. However, if we compare only the most stable proton configurations, the picture is completely different. Among all water polyhedra with the number of molecules 20 and 24, the proton configurations of clusters with a large number of square faces are the most stable [7]. In such clusters, the maximum number of the most profitable types of H-bonds is greater. Calculations of the energy for water polyhedra were performed using the TIP4P [9] and AMOEBA [10] force fields. For the polarizable force field Amoeba, the determining role of square faces is manifested to a greater extent [7]. The conclusion about the special energy efficiency of the cluster with the number of molecules 24 containing no pentagonal faces was also made by Anick [11]. Note that clusters with a large number of square faces have another potential advantage. Square faces, in which opposite molecules are located in the planes of these faces, can take the shape of a diamond. Such spontaneous deformation gives an additional gain in energy. The number of faces with such an arrangement of molecules is an accompanying topological parameter, which also gives preference to polyhedra with a large number of square faces [7].
The motivation of this article is to extend the conclusion about the dominant role of the topological stability factor of water polyhedra to a wide temperature range. In this study, we will limit ourselves to seven polyhedral clusters with the number of molecules (vertices) N = 20 and 24. The structure of the article is as follows: The methodology section presents a method for obtaining the most stable proton configurations. This section also provides a brief description of molecular dynamic simulation methods and a description of methods for analyzing structural transformations. In the next result section, the temperature limits of structural stability, the classification of the most characteristic structural defects, and the description of the melting-like transitions are given for each of the selected configurations. The final section contains a discussion of the results obtained and conclusions.

Method
A number of polyhedra with the number of vertices N = 20 and 24 formed by square, pentagonal, and hexagonal faces are 23 and 59, respectively [12,13]. Only the so-called cubic polyhedra, which correspond to cubic graphs, are taken into account here [14]. Three edges converge at each vertex of such polyhedra. The most stable proton configurations were selected from a variety of proton configurations of different water polyhedra [7]. The structure of polyhedra forming the most preferred proton configurations is shown in Fig. 1.
Based on the discrete model of strong and weak effective H-bonds (the SWEB model), the configurations were preselected during combinatorial optimization [8]. Optimal configurations contain the maximum number of more preferred types of H-bonds (t1-bonds) [7]. The final selection of the most stable configurations was made during calculations with the TIP4P [9] and AMOEBA [10] force fields.
In this paper, we analyze the most stable proton configurations of three polyhedra (H 2 O) 20   The regular pentagonal dodecahedron {0,12,0} and the polyhedron {0,12,2} do not have square faces. They correspond to the small and large cavities of the hydrate structure I. Note that in cubic polyhedra, when the number of pentagonal faces decreases, the number of squares and hexagonal faces increases simultaneously. The number of pentagonal faces is always even, and the number of square faces varies from 0 to 6 [12].
To study the thermal stability of the selected configurations, the Tinker molecular modeling package [15] and the AMOEBA force field [10] were used. Molecular dynamic simulation was carried out using the Berendsen thermostat with the time damping constant maintained at 0.1 ps [16]. Starting from 0 K, the temperature gradually rose by 20 degrees. For each temperature, a simulation time was 100 ps with a standard time step of 1 fs. To reduce the impact of randomness and have a more meaningful comparison, we repeated all heating processes 10 times. MD trajectories including atomic coordinates are saved every 100-time steps for analysis.
The analysis of structural changes was performed using the adjacency matrix. The presence of H-bonds between the molecules of water polyhedra was checked by a geometric criterion that takes into account the values of the distances between the oxygen atoms of neighboring molecules (R OO < 3.5 Å) and the minimum distance between the oxygen atom of one molecule and the hydrogen atom of another molecule (R OH < 2.5 Å). Short-term structural changes were not taken into account. The temperature of loss of stability was considered to be the temperature at which the topology of the H-bond network changed irrevocably. The average value between the last temperature of unchanged shape T u and the destruction temperature T d was considered a temperature limit of stability T s = T d -10 (temperature increment is 20 K). To determine the nature of the defect, some sections of MD trajectories were visualized using the VMD software package [17].
We believe that the temperature stability limit T s is the lower boundary of the transition region between solidlike and liquid-like states of the clusters. To determine the boundaries of this region, caloric curves were calculated showing the change in the potential energy of the cluster as a function of temperature. In addition, calculations of the Lindemann index, which characterizes the fluctuation of interatomic distances (oxygen-oxygen), were carried out.

Results and discussion
Thermal destruction of the initial structures Table 1 shows the temperature limits of structural stability for all selected configurations of clusters. As noted above, in cubic polyhedra, an increase in the number of square and hexagonal faces occurs simultaneously with a decrease in the number of pentagons. We will estimate the dependence of stability on the number of squares, the number of which in such polyhedra varies from 0 to 6 [7,12].
First, we note that for water polyhedra with the number of molecules N = 20, the temperature limits of structural stability are completely consistent with the stabilization energies at zero temperature [7]. The most stable configuration is still configuration no. 3 {444} with the largest number of square faces (Fig. 1). The configuration of regular dodecahedron no. 1, formed only by pentagonal faces, is the least stable of the three. For water polyhedra with the number of molecules N = 24, the agreement of the results is also good. The most stable configuration under the influence of temperature is still configuration no. 7 {608} with the maximum number of square faces. Unlike the results at 0 K, configuration no. 4 {0,12,2} is not the least stable. But the advantage of configuration no. 7 over all others is still very strong.
The following is a description of the main types of defects that destroy the initial polyhedral shape of the clusters. Although the clusters under study are empty, thermally Table 1 The temperature limits of structural stability (K) of water polyhedra nos. 1-7 for ten realizations of the heating process No  1  2  3  4  5  6  7   1  110  150  170  130  110  110  150  2  110  130  150  110  130  110  130  3  130  130  150  130  130  90  130  4  130  150  130  130  90  130  170  5  130  90  150  110  90  110  150  6  110  150  170  130  90  110  170  7  130  130  150  130  90  130  150  8  130  130  170  130  110  90  170  9  130  130  170  130  130  130  190  10  130  150  150  110  130  110  170  Average  124  134  156  124  110  112  158 induced rearrangements rarely immediately lead to a collapse of the initial structure with the appearance of distant jumper bonds and the movement of some molecules to the central part of the cluster. Structural changes begin with the appearance of surface defects. Thermal fluctuations often lead to temporary structural damage, after which the original structure is restored. We analyzed only fatal defects that lead to further destruction of the original structure. These defects can be divided into three main types.
The most common surface defect is the appearance of a bridge inside a pentagonal or hexagonal face (Fig. 2a). Adjacent square and triangle most often appear instead of a pentagon (sqtr-defect). Such a transformation is also observed in the simulation of an isolated cyclic pentamer [18]. In our case, it is common for a triangle to be open; one of the H-bonds is broken. At the site of such a defect, the surface of the polyhedron becomes concave. There are two main reasons for this transformation. This can be caused by intense thermal fluctuations of the free O-H group, especially if such groups are located nearby. Under the influence of temperature or the repulsive force from the neighboring free hydrogen atom, one of these atoms breaks the existing one and creates a new surface H-bond. In this case, the neighboring molecule forms a bridge bond inside the pentagonal face. In addition, the cause of such a bridge may be the deformation vibrations of the pentagonal face itself. This bond occurs when the pentagonal face becomes very narrow and elongated.
Here, it is important to emphasize that the formation of the sqtr-defect leads to the appearance of 4-coordinated molecules (Fig. 2a), although it causes strong deviations from the tetrahedral coordination of the H-bond network. The transformation of the pentagonal face into square one at the initial stage of the structure change is the most typical defect. In polyhedra with a large number of hexagons, a similar narrowing of the hexagonal ring to the pentagonal one often occurs.
The second type of fatal defect is the occurrence of an inverse O-H group, i.e., when the "dangling" hydrogen atom becomes directed towards the cluster center (Fig. 2b). This type of defect was observed quite often (in about three cases out of ten), but basically in clusters with 20 molecules. Such a rearrangement is accompanied by a local flattening of the polyhedral surface or a small deepening. Subsequently, this often leads to the appearance of a distant bridge and the destruction of the polyhedral shape. Stability over 1 ps for a structure with a reverse O-H group seems to be quite unexpected. This is primarily because the most stable proton configurations of water polyhedra were selected for the study.
Another, rarer type of surface defect is the topological transformation according to the scheme: pentagon and hexagon two squares and pentagon (Fig. 2c). Such a spontaneous structural transformation was observed, in particular, in cluster no. 3 {444}. In this case, a noticeable distortion of the tetrahedral coordination is compensated by the appearance of two 4-coordinated molecules. The resulting polyhedron {643} is no longer cubic, but its surface is convex. Having a fairly even polyhedral shape, it surpasses regular polyhedra in the number of H-bonds. It is not difficult to make sure that for an even number of molecules N in the initial cubic polyhedron, the number of emerging 4-coordinated molecules k is also even. This is because the doubled number of H-bonds is determined by the expression 2n = 4k + 3(N − k) = 3N + k . Each pair of 4-coordinated molecules increases the number of H-bonds by one (Fig. 1), since in cubic polyhedra the number of bonds n 0 is one and a half times greater than the number of vertices 2n 0 = 3N . All other types of fatal defects appeared extremely rarely. Such rare defects include an inversion with the simultaneous appearance of a distant jumper.
The spontaneous appearance of surface defects and their stability for some time means a small difference in the energies of these configurations from the initial ones. Recall that for the study, the lowest energy proton configurations were selected from a huge set of configurations that satisfy the Bernal-Fowler ice rules. Figure 3 shows configurations of the two lower energy levels of cluster no. 1 (regular pentagonal dodecahedron) following the discrete SWEB model (106 and 1541 configurations) [8]. The total number of defect-free configurations in this cluster is 30026 [5]. The calculations were performed using the Amoeba force field. The minimum energy E MIN = -41.303 kJ/mol. The energy of geometrically optimized structure containing the sqtrdefect E SQTR = -41.248 kJ/mol (Fig. 3, circle). This is the energy of the most stable sqtr-defect of those that spontaneously occurred during the simulation. It is noteworthy that it is only slightly inferior to the energy of the lowest energy defect-free configuration of regular water dodecahedron. The energy of the geometrically optimized structure of non-cubic polyhedron with two 4-coordinated molecules E NC = -41.180 (Fig. 3, rhombus), which is also quite surprising, since this is only one of the possible polyhedra. Special calculations were performed for the defect in the form of an inverse O-H group. All O-H groups were sequentially inverted for all 106 configurations of the lower energy level in the SWEB model (1060 structures; the number of dangling hydrogens is N/2 [5]). The minimum energy of the configuration with the inverse O-H group is noticeably lower: E OH = -41,008 ( Fig. 3, triangle). However, on the set of all defect-free configurations of a regular dodecahedron, such a configuration should still be considered among the most stable. Cartesian coordinates of all atoms of the three defect configurations can be found in the Supporting Information. Coordinates of all seven defect-free configurations of clusters nos. 1-7 can be found in ref. 7.

Melting-like transitions
As the temperature increases, the distinction between the water clusters with the same number of molecules begins to disappear. This mainly occurs in the transition region from solid-like to liquid-like states of clusters. However, one important point needs to be made here. The states of the clusters at low temperatures are not equilibrium, because there are significantly more stable non-polyhedral structures. For example, the global minimum energy structure of a water cluster with 20 molecules represents three fused pentagonal prisms [19,20]. A spontaneous transition to this structure, as well as to other very stable structures (finite ice nanotubes), is extremely unlikely due to the high barriers between adjacent local minima on the potential energy surface. When water polyhedra are heated, these structures practically do not manifest themselves in any way, and the transition from solid-like to liquid-like state occurs without them. It is extremely difficult to calculate the actual thermodynamic properties of such clusters [21,22]. This is not the purpose of our work. We only evaluate the difference in stability of various polyhedral configurations at different temperatures including the region of melting-like transition. Figure 4 shows the change in the potential energy of clusters as a function of temperature for all ten realizations of the heating process. For each of the selected cluster configurations, when the heating process is repeated, the initial sections of the caloric curves exactly coincide (there is no spread of curves). The destruction of the initial configurations and the beginning of the broadening of the caloric curves corresponds to the lower boundary of the transition from solid-like to liquid-like state (see the red mark in Fig. 4). General bending in the caloric curves is noticeable in the figure.
The upper boundary of the transition region is better visible in Fig. 5 which shows the dependence of the Lindemann index on temperature. This indicator is equal to the relative mean-square-root fluctuation of the interatomic (O-O) distance: where r ij is the distance between the i-th and j-th atoms; N is the number of water molecules. Figure 5 also shows the temperature limit of structural stability (red mark). The upper boundary of the transition region for all clusters is approximately 200 K and Fig. 4 The caloric curves of water clusters nos. 1-7 for ten realizations of the heating process corresponds to δ ≈ 0.1. The obvious exception is cluster no. 7 {608}, for which both the lower and upper boundaries of the transition region are noticeably higher. This confirms the extremely high stability of this cluster. At the same time, an overestimated value of the upper bound in comparison with other polyhedra N = 24 means dependence on the initial structure and proves the non-equilibrium nature of the melting-like transition. With a further increase in temperature, the non-equilibrium decreases. It was established that the evaporation of molecules from all clusters occurs in the same range of 292-300 K (averaging over ten heating processes).

Conclusions
A molecular-dynamic study of water polyhedra allows us to conclude that on the set of the most stable proton configurations, the geometric factor of stability is inferior to the topological factor. The most tetrahedrally coordinated structures with the maximum number of pentagonal faces are not the most stable. As at zero temperature, structures with a noticeable deviation from tetrahedral coordination and with a large number of flat square and hexagonal faces are more stable. This provides an unusual picture of the stability of polyhedral water clusters and casts doubt on the idea that the regular pentagonal dodecahedron is the most stable structure for the cluster (H 2 O) 20 . At the same time, the polyhedra nos. 5-6 with N = 24 have average temperature limits inferior to polyhedron no. 4. This nonlinear dependence on the number of pentagons is a reflection of the competition between two stability factors. Another interesting outcome of our simulation is the discovery of very stable defect configurations of water polyhedra. They are structures with surface defects, convex noncubic water polyhedra and configurations with an inverse O-H group. Most of the defect structures are characterized by an increase in the number of square cycles and the appearance of 4-coordinated molecules. The deviation from the tetrahedral coordination is compensated by the appearance of additional H-bonds. Defect configurations do not stand out in energy among the many defect-free polyhedral configurations. Therefore, the idea of the presence of a small class of clearly the most preferred topologies (adjacent pentagonal prisms, finite nanotubes, regular polyhedron) is a strong simplification of the real picture [23,24]. The results obtained indicate that the morphological diversity of water clusters is more continuous, even if only the most energetically favorable structures are taken into account. This is consistent with the fact that the Atlas of the lowest energy configurations of water clusters covers a wide variety of topologies of the H-bond network [25].
The results obtained also suggest that the most stable water polyhedra can contain 4-coordinated molecules. After all, a cluster of this type that we accidentally discovered is only one of a huge set of such configurations, which differ from each other in the topology and the arrangement of protons in hydrogen bonds. At that, it is only slightly inferior to the most stable proton configuration.
Author contribution There were equal contributions of the authors to the completion of this work.

Funding
The present work was supported by the Basic Research Program of the Russian Academy of Sciences, Project No. 121041600040-3.

Data availability
The data that support the finding of this study are available within the article.
Code availability TINKER Software Tools for Molecular Design. Version 6.2, VMD 1.9.3.