Probing the network topology in network-forming materials: The case of water

Rings statistic has been widely used to investigate the network topology in numerically simulated network-forming materials in order to rationalize their physical and mechanical properties. However, different topologies arise depending on how rings are counted, leading to incomplete or even contrasting physical interpretations. Solving this critical ambiguity is of primary importance for the correct assessment of material properties. Here, we show how such differences emerge in water, a complex network-forming material endowed with polyamorphism and a directional network of hydrogen bonds whose topology is correlated with the anomalous behavior of water. We probe the network in the liquid state at several thermodynamic points under equilibrium conditions, as well as during the out-of-equilibrium first-order-like low density to high density amorphous transformation. We study three schemes for counting rings and show that each of them provides complementary insightful information about the network, suggesting that a single counting scheme may not be sufficient to properly describe network topologies and to assess material properties. Our results provide a molecular description of the rings in supercooled water and of the amorphous-to-amorphous transformation kinetics, hence shedding light on the complex nature of water. Nonetheless, our results expose how delicate the proper choice of method for counting rings is, an issue with important consequences for rationalizing the properties of network-forming materials at large.


I. INTRODUCTION
In his seminal paper, Zachariasen 1 rationalized the evidence that, under some thermodynamic conditions, the mechanical properties of glasses are comparable with those of crystals, suggesting that the way the atoms are connected with each other must play a crucial role in the determination of such properties. Nowadays, the connectivity in network-forming materials generated via numerical simulations is commonly probed via the rings statistic, i.e., by counting the number of rings (or closed loops) present in the network and formed following the links (bonds or interactions) between atoms. The rings statistic has been widely adopted in characterizing the network topology of several simulated network-forming materials, from amorphous systems 2,3 to clathrate hydrates 4 and chalgogenide glasses, 5 and is an essential tool to characterize continuous random networks. [6][7][8][9][10][11][12] Nonetheless, the rings statistic has been widely employed in water to characterize its different phases, [13][14][15][16][17] to investigate the origin of water anomalies, [18][19][20][21] and to inspect the dynamics of homogeneous nucleation. 20,22 Several definitions of rings and counting schemes have been reported in the literature, 2,3,23-27 which yield different and inconsistent results even for the simplest crystalline structures, not to mention more complicated networks, such as amorphous silicate structures 3,25,[28][29][30][31] or water. 17,18,21,22 In this work, we compare different ring counting schemes and show, for the first time, that each scheme carries complementary information leading to incomplete or even contrasting physical interpretations if considered individually. In particular, we probe the network topology of water, a complex 32,33 network-forming material with a directional network of hydrogen bonds (HBN) and endowed with polyamorphism, i.e., water can acquire more than one amorphous state. [34][35][36][37][38] The low density amorphous (LDA) ice is characterized by similar degrees of density and tetrahedrality to crystalline ice and can be produced, e.g., upon rapid quenching of liquid water down to low temperatures or upon the vapor deposition of water droplets onto cold surfaces, while the high density amorphous (HDA) ice is characterized by a density ∼20% higher than that of LDA and can be produced, e.g., upon isothermal compression of LDA. 14,39,40 Here, we study the HBN of liquid water under different thermodynamic conditions, as well as the kinetics of the HBN during the LDA-to-HDA phase transformation. Nonetheless, we also show that each counting scheme can provide information about the molecular nature of different rings, i.e., how many water molecules can donate-accept hydrogen bonds (HBs), the information that will be crucial for further investigations of homogeneous and heterogeneous nucleation and for a better understanding of water anomalous behaviors.
The definition of hydrogen bond (HB) follows Ref. 41. In this regard, any quantitative measure of HBs in liquid water is somewhat ambiguous since the notion of an HB itself is not uniquely defined. However, qualitative agreements between many proposed definitions have been deemed satisfactory over a wide range of thermodynamic conditions. 42,43 We construct rings by starting from a tagged water molecule and recursively traversing the HBN until the starting point is reached or the path exceeds the maximal ring size considered (10 water molecules in our case). We show that different counting procedures give different distributions, opening the doorway to incomplete, or even contrasting physical interpretations if such distributions are not carefully analyzed together. It is important to emphasize, at this point, that the three counting schemes here adopted are not exhaustive and other schemes can be implemented. On the other hand, the scope of this work is not to analyze all possible counting schemes, rather to show that using only one of them may lead to incomplete interpretations.
This article is organized as follows: In Sec. II, we briefly report the details of the numerical simulations. Section III is devoted to the presentation of the three counting schemes. In Sec. IV, we present the main findings of this work. Results for the case of supercooled water are reported in Sec. IV A and the results of the LDA-to-HDA phase transformation are reported in Sec. IV B. In Sec. V, we summarize our findings and discuss their implications.

II. COMPUTATIONAL DETAILS
Our studies are based on classical molecular dynamics simulations of rigid water molecules described by the TIP4P/2005 interaction potential 44 in the isobaric (NpT) ensemble. We have performed simulations with the GROMACS 5.0.1 package 45 running on IBM POWER8 machines with NVIDIA Kepler K80 GPUs.

A. Liquid water under equilibrium conditions
For the system of liquid water under equilibrium conditions, the system is composed of 512 molecules and details of the numerical simulations can be found in Ref. 18. In this work, we analyze the trajectories reported in Ref. 18 at the pressures of p = 1 bar from ambient to T = 210 K, and p = 400 bars in the temperature window from T = 230 K to T = 195 K. The results shown here are the average at each state point of five independent trajectories. At each state point, we have computed and carefully monitored the decay of the self-part of the intermediate scattering function with time 46 to assure proper equilibration, and we have monitored the absence of crystallites using a sensitive order metric. 13

B. Amorphous states
We have prepared LDA by cooling a sample of liquid water composed by 50 000 molecules under equilibrium conditions from T = 300 K to T = 80 K at ambient pressure using a continuous cooling rate of qc = 1 K/ns. In order to prepare HDA, we have simulated the isothermal compression of the so formed LDA from ambient pressure to p = 200 MPa with a compression rate of qp = 1.0 MPa/ns. Further details of the simulations can be found in Ref. 47. As shown in Refs. 14, 40, and 47, the LDA-to-HDA transformation is remarkably sharp, as also observed in experiments.

III. COUNTING SCHEMES
In the three definitions for counting rings, we move from a highly constrained scheme (def1) to a less constrained scheme (def3). In def1 (upper panel of Fig. 1), we consider only the shortest primitive rings (rings that cannot be decomposed into smaller ones 23,27,48 ) generated by water molecule labeled as 1 in Fig. 1 donating an HB. In this example, molecule 1 donates one HB to molecule 2 and the overall network, emphasized by the arrows, forms an hexagonal ring. It is important to notice that, according to this definition, it is irrelevant whether other molecules involved in the ring accept or donate a bond, a further restriction that could be taken into account in another counting scheme. Therefore, this definition particularly emphasizes the intrinsic directional nature of the HBN in water that might be the kinetic driving force in processes like homogeneous and heterogeneous nucleation. In the second scheme, def2, lower panel of Fig. 1, we remove the constraint of directionality in the HBs of molecule 1 and focus on the shortest primitive rings regardless of molecule 1's donor/acceptor nature. According to this definition, the shortest ring obtained starting from molecule 1 is the pentagonal ring emphasized on the lower panel of Fig. 1. In this case, the counted ring is formed when molecule 1 accepts one HB (either from molecule 2 or molecule 5) instead of donating as in the previous definition. In the third scheme, def3, we remove the constraint of counting only the shortest rings. Within this definition, the water molecule from which rings are counted is typically part of several rings and, therefore, def3 can capture also the formation of longer rings that characterize the network of a disordered material. In the example reported in Fig. 1, according to def3, water 1 generates and belongs to both the hexagonal and the pentagonal rings. As such, def1, def2, and def3 represent decreasingly strict criteria for rings, but also probe distinct physical features of the HBN. As expected,

ARTICLE
scitation.org/journal/adv FIG. 1. Schematic representation of the network of HBs connecting water molecules. Red filled circles represent oxygen atoms, while red empty circles represent hydrogen atoms. Water molecule labeled as 1 is the starting molecule from which the rings are counted. Upper panel: according to def1, the shortest ring is the hexagonal ring whose network is emphasized by arrows. Lower panel: according to def2, the shortest ring is the pentagonal ring whose network is emphasized by arrows. According to def3, water molecule 1 generates both the hexagonal and the pentagonal rings.
but so far never discussed in the literature, these three definitions produce different distributions that can be interpreted in different ways.

IV. RESULTS
In this section, we present the main findings for the test cases we have chosen: liquid water at 1 bar and at 400 bars at different temperatures under equilibrium conditions, and the out-of-equilibrium isothermal compression of LDA to generate HDA.

A. Equilibrium liquid water under ambient and supercooled conditions
In 1976, Speedy and Angell observed that, at variance with other liquids, the isothermal compressibility of supercooled water increases upon cooling suggesting the presence of a thermodynamic singularity at −45 ○ C. 49 In 1992, based on the results of classical molecular dynamics simulations, Poole et al. hypothesized that a liquid-liquid critical point (LLCP) located under deeply undercooled conditions and low/intermediate pressures could be the source of the observed thermodynamic anomalies. [50][51][52] In correspondence with the LLCP, two liquids [an high-density liquid (HDL) and a low-density liquid (LDL)] are in metastable equilibrium with each other. The existence of the LLCP has been proven via numerical simulations for a molecular model of water, 20 and several other computational [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] and experimental 37,70-83 investigations point toward the same confirmation. The outstanding importance of water anomalies is nowadays recognized to play a central role in processes in several fields, from physics and chemistry, to biology, material science, geology and climate modeling. 34,37,84,85 Recently, one of us has probed the HBN connecting LDLlike environments using def3, showing that the HBN within LDLlike environments is directly connected to-if not responsible forthe appearance of the thermodynamic anomalies of water. 18 Here, we report the distribution of rings obtained by probing the entire HBN, without distinctions between HDL-like and LDL-like environments using the three counting schemes above reported and for two pressures, namely 1 bar and 400 bars.
In each panel of Fig. 2, we report the probability P(n) of having a ring of length n, with n ∈ [3,10] for the cooling of liquid water from T = 290 K to T = 210 K at ambient pressure. All P(n)s have been normalized to unity and therefore do not reflect the total number of rings of a given size.
We can observe that def1 (upper panel) gives a distribution that, at high temperatures, is broad and slightly maximized in correspondence with the hexagonal rings. Short (n < 5) and long (n > 6) rings are present in non-negligible amounts because the thermal energy induces molecular diffusion and the network can explore a larger amount of configurations. Upon cooling, we observe a marked depletion of longer rings and squared rings. Correspondingly, we observe a slight increase in the pentagonal rings and a pronounced increase in the fraction of hexagonal rings that correspond

ARTICLE
scitation.org/journal/adv to the typical network topology of both cubic (Ic) and hexagonal (Ih) ices. Therefore, def1 accounts for the longer rings and, upon cooling, describes a network that tends toward (but is far from) that of the stable crystal form. The ring distributions computed via def2 (middle panel) remarkably differ from the distributions computed via def1. Even at high temperatures, the distribution shows well-defined maxima in correspondence with the pentagonal rings, a lower contribution of the hexagonal rings, and a marginal contribution from the heptagonal rings with longer rings being completely absent. In light of the fact that the longer rings represent natural configurations to accommodate higher densities 21 (density fluctuations also play a crucial role in understanding water anomalies 20,50,67,86 ), the distribution in the middle panel of Fig. 2 suggests that important configurations are not taken into account in def2. Nonetheless, upon cooling the sample down to 200 K, the heptagonal rings gradually disappear in favor of the pentagonal rings and only marginally in favor of few hexagonal rings.
In order to rationalize the striking difference between the distributions obtained from def1 and the distributions obtained from def2, we observe that, accordingly with def2, the pentagonal rings present in the system favor configurations in which some water molecules are only hydrogen-acceptors. For a pentagonal ring, at most 2 water molecules can be solely hydrogenacceptors such that for each pentagonal ring, the number of hexagonal rings per water molecule that can be added to def1 is, at most, 2 ×n 6 wheren 6 is the average number of hexagonal rings per water molecule. Therefore, the differences between the distributions for def1 and def2 provide crucial information about the structure of the pentagonal rings, i.e., how many of the participating water molecules are only hydrogen-acceptors. This information could help in enlightening the role of pentagonal rings in the anomalies of water, 18,19 as well as in better understanding the role of pentagonal structures in frustrating against crystallization. 87,88 Moreover, the increasingly hexagonal-character of the HBN upon cooling shown in def1 suggests that the directionality of the HBs plays a fundamental role, possibly as kinetic driving force, in complex processes such as homogeneous and heterogeneous nucleation of water.
In the lower panel of Fig. 2, we report the ring distributions obtained by counting the rings with def3. As in the case of def1, the distributions are peaked over the hexagonal rings but very long rings are taken into account. In particular, at high temperatures, the network is almost equally described by 5-, 6-, and 7-folded rings, and the distribution is very broad accounting also for the longer rings. Such a distribution is conceivable with that of a network exposed to high thermal energy. Upon cooling, the hexagonal character of the network emerges and the longer rings are depleted until no 10-folded rings appear at T = 210 K. Nonetheless, even at T = 210 K, the network still hosts 8-and 9-folded rings, as one would expect in a diffusing medium. Such rings are absent in def1 and def2 because all water molecules are concurrently part of smaller rings.
In Fig. 3, we report the P(n) computed for liquid water at p = 400 bars in the temperature window T = 195-230 K, the temperature window at which liquid water crosses the lines of maxima of thermodynamic response functions at this pressure. 18 According to def1 (upper panel), at T = 230 K, the network of the sample is mostly arranged in 5-, 6-and 7-folded rings and does not contain longer rings. Upon cooling, we observe a gradual depletion of the heptagonal rings in favor of the hexagonal rings as expected since we are approaching thermodynamic conditions of relevance for Ih and Ic. A different scenario holds if we count the rings using def2 (middle panel), where we can appreciate how the distributions at different temperatures describe a network, where most molecules are part of at least one pentagonal or hexagonal ring. We can observe that at T = 230 K, the system is dominated by the pentagonal rings, while the heptagonal rings are already absent. Upon cooling, we observe a depletion of the pentagonal rings in favor of the hexagonal rings, with a switch to a distribution dominated by 6-folded rings at the lowest temperature inspected, T = 195 K. This observation indicates that the network becomes more hexagonal-like regardless of the acceptor/donor nature of the starting water molecule. This tendency contrasts with the behavior of def2 at p = 1 bar (see Fig. 2 middle panel), where, instead, we observe an increase in the pentagonal nature of the network, in agreement with Ref. 17 (different definitions of HB do not affect the overall results 42,43 ). Therefore, this effect can be purely ascribed to the increased pressure. The lower panel shows P(n) computed using def3. We can observe that, at the highest temperatures, the distribution is broad with an almost equal contribution of the hexagonal and heptagonal rings, and accounts also for the longer rings, as one would expect at higher pressures. Upon cooling the sample, we observe a depletion of longer rings, an almost constant profile for the heptagonal rings, and an increase in the hexagonal rings. Therefore, def3 is able to capture both the increase in the hexagonal character of the network upon cooling and the presence of longer rings that allows us to host higher densities.

B. LDA-to-HDA phase transformation
A remarkable evidence in favor of the LLCP is the polyamorphic character that water acquires at low temperatures and ARTICLE scitation.org/journal/adv low/intermediate pressures. The two amorphous states LDA and HDA are separated by a first-order phase boundary 14,40,70,[89][90][91][92] and correspond to the glassy states of LDL and HDL. A very high-density amorphous state has also been hypothesized. 93 Signatures of LDA and HDA have been found in the inherent potential energy surface (IPES, a collection of local potential energy minima resulting from systematic quenches along a liquid trajectory allowing to separate packing from thermal motion effects 94,95 ) of ab initio liquid water under ambient conditions, where the LDA-like and the HDA-like environments are characterized by having distinct HBNs. 21 The distinct HBNs in the IPES of ab initio liquid water indicate that there is an energy barrier associated with the different HBNs, and this barrier could be large enough at low temperatures and pressures to account for the first-order nature of the LDA-HDA phase boundary. Two distinct local structures have also been observed experimentally in liquid water from ambient to supercooled conditions, and interpreted as LDL-like and HDL-like environments 96 as well as under pressure. 97 In Fig. 4, we report the fraction of n-folded rings, n ∈ [5,8], as a function of the pressure when we simulate the isothermal compression of LDA to generate HDA. We can observe that all counting schemes capture the phase transition from LDA to HDA occurring in the interval p = 75-100 MPa, 14,36,40,98 but each definition describes a completely different kinetics (details of the HBN kinetics during the decompression of HDA and during compressiondecompression cycles as a function of different temperatures will be presented elsewhere). Recently, it has been observed that amorphous ices hide deep connections with high pressure crystalline ices and, in particular, with the metastable ice IV. 14,39 In order to unravel such connections, a detailed understanding of the kinetics of the HBN is necessary. Therefore, the issue raised in this work about the different possible pathways to probe the network topology is of remarkable importance. According to def1 (upper panel), the network of LDA resembles, as expected, the network of supercooled liquid water shown in the upper panel of Fig. 2 and is dominated by the hexagonal rings that, in correspondence with the phase transition, suddenly decrease from a fraction of ∼0.6 to ∼0.4 and keep decreasing in the network of HDA. On the other hand, the pentagonal rings retain an almost constant fraction of ∼0.3 independently of the pressure, while the heptagonal rings show an increment from a fraction of ∼0.1 to ∼0.3 with a slow but continuous increment with the pressure in HDA. The longer rings (n = 8) are almost absent in LDA and appear in the network only in correspondence with the phase transition, showing a continuous increase upon increasing the pressure. Since the longer rings account for the increase in the density characterizing the HDA phase, the network topology described by def1 is conceivable with the LDA-to-HDA phase transformation.
The middle panel of Fig. 4 shows P(n) as a function of the pressure for def2. At variance with def1, the network of LDA resembles that of supercooled liquid water shown in Fig. 2 middle panel, and is dominated by the pentagonal rings (∼0.6) that, in correspondence with the phase transition, show only a minor increase followed by a minor decrease upon further increase in the pressure. The fraction of hexagonal rings decreases from ∼0.4 to ∼0.2 in correspondence with the phase transition. Interestingly, the longer rings are almost completely absent upon increasing the pressure, a tendency incompatible with the (sudden) increment of the density in HDA.
The lower panel of Fig. 4 reports P(n) as a function of the pressure for def3. As for the supercooled liquid case reported in Fig. 2 lower panel, the fraction of each ring is more equally distributed along different lengths. The hexagonal rings account for the higher fraction (∼0.4) followed by the heptagonal (∼0.3), pentagonal (∼0.2), and octagonal (∼0.1) rings. Each fraction remains almost constant upon increasing the pressure and, in correspondence with the phase transition, we observe a sudden decrease of rings with the length n < 8, while the octagonal rings and longer rings (data not shown) increase to account for the increase in the pressure.

V. CONCLUSIONS
In conclusion, we have inspected three different ring counting schemes adopting different constraints and we have investigated the network of HBs in the case of supercooled liquid water under equilibrium conditions and different state points as well as during the isothermal compression of LDA to generate HDA. We have adopted the geometrical definition of the HB reported in Ref. 41, whose validity has been assessed under the thermodynamic conditions explored here. 42,43 Each counting scheme provides important, complementary insights. In particular, def1 introduces an asymmetry into the measure based on acceptor/donor type of the starting molecule, which biases the analysis toward the intrinsic directionality of the HBN, while def2 and def3 gradually remove this constraint.
As a first test, we have probed the network of liquid water at p = 1 bar and at p = 400 bars and different temperatures. The scheme def1 applied to supercooled water shows that most of the hexagonal rings are generated by one water molecule donating one HB. Upon cooling liquid water, the network topology tends toward that of crystalline ice, as one would expect. The increasingly hexagonal nature of the HBN indicates that the directionality of the HBs should play ARTICLE scitation.org/journal/adv a fundamental role in driving the kinetics of homogeneous crystallization of water, a phenomenon whose kinetics is still the subject of intense investigation. We propose here that introducing further constraints in def1 such as donor-only or acceptor-only characters on the other molecules involved in the counting scheme could provide further, important insights in the kinetics of such processes. We also remark here that, although the network tends toward that of crystalline ice, all samples are liquids under equilibrium conditions 18 and we have verified the absence of crystallites with a sensitive order metric. 13 The counting scheme def2 does not distinguish between the acceptor/donor character of the starting molecule that generates the rings, and the distributions computed with this criterion describe a network in liquid water mostly characterized by the pentagonal rings. The central role played by the pentagonal rings in the anomalies of water has been discussed in the literature 18,19 as well as in the field of geometrical frustration, 87 and here we provide a molecular description of the composition of such pentagonal rings. Nonetheless, as shown in Fig. 3, middle panel, def2 reports a switch from a topology dominated by the pentagonal rings at higher temperatures, to a topology dominated by the hexagonal rings at lower temperatures. Considering the validity of the definition of HB under these conditions, 42,43 we infer that this observation has to be attributed to the balance between low temperature and low pressure that favor the formation of an hexagonal network. Therefore, taken together, def1 and def2 provide information about the nature of hexagonal and pentagonal rings and, in particular, of the number of participating water molecules that are only acceptors. Therefore, def2 can be considered as a useful tool to, e.g., design colloidal self-assembly. 99 Finally, def3 removes the constraint of counting only the shortest rings starting from a water molecule. This measure captures both the hexagonal character of the network and the presence of longer rings, expected in a diffusive disordered network, but the important information about the nature of the rings at the molecular level are washed out. None of the counting schemes here adopted capture any sensitive change in the topology of the HBN in correspondence with maximal thermodynamic response functions. In Ref. 18, one of us has shown that such correlations occur when using def3 to study the topology of the HBN connecting LDL-like environments only, while the network connection in HDL-like environments only does not show such correlation. Since the supercooled liquid samples always contain a majority of HDL-like environments, 18 we infer that the correlation with the maxima of thermodynamic response functions in the overall HBN here inspected is shrouded by the heavier contribution of the network connecting HDL-like environments. Nonetheless, other counting schemes might be able to capture such correlation between the overall HBN and maxima in the thermodynamic response functions.
Next, we have adopted the three counting schemes to study the kinetics of the LDA-to-HDA first-order-like phase transition at T = 80 K. According to def1, the transition to HDA is characterized by a significant drop in the hexagonal rings and an increase in the heptagonal rings. At the highest pressure here reported of 200 MPa, the network of HDA is characterized mostly by the pentagonal, hexagonal, and heptagonal rings. The presence of longer rings and the drop in hexagonal rings is in agreement with the evidence that the longer rings are required to host the increased density in HDA. On the other hand, the presence of only heptagonal rings (in a minor fraction compared with the hexagonal and pentagonal rings) among the longer ones is not enough to account for the high density of HDA. According to def2, the transition from LDA to HDA causes a drop in the fraction of hexagonal rings and an increment in the pentagonal rings. At the highest pressure, the network of HDA is mostly arranged in the pentagonal rings, with a minor fraction of hexagonal rings and longer rings accounting for only negligible fractions. Hence, as for the liquid case, the network described by def2 provides a description that, considered alone, does not agree with physical intuition. The most reliable kinetic profile of the LDA-to-HDA transformation is provided by def3, according to which the network of LDA is dominated by hexagonal rings followed by the heptagonal, pentagonal, and octagonal rings, respectively. In correspondence with the transition, we observe a sensitive drop in the hexagonal and heptagonal rings and an increase in the octagonal rings. At the highest pressure, the network is dominated by the octagonal rings, with an equal amount of shorter and longer rings. Further studies enlightening the kinetics during the compressiondecompression cycles at different temperatures will be presented elsewhere. 100 In conclusion, the three counting schemes presented here are not mutually exclusive but, rather, provide complementary information and hence should be investigated together, and other counting schemes might carry different important information. While our conclusions can be quantitatively applied only to the case of water, and other definitions might be more suitable to inspect the network topology in other materials, our results expose how delicate the proper choice of counting rings is, in general, an issue that has never been reported in the literature before and that is critical for the proper extrapolation of physical and mechanical properties in network-forming materials at large, such as molecular liquids, polymer melts, amorphous materials, and vitrimers. This issue becomes even more evident in the case of phase transitions, where pervasive network rearrangement occurs, 40 and unclear connections with different crystalline phases might be present. 14,39