Structural features of high-local-density water molecules: Insights from structure indicators based on the translational order between the first two molecular shells

The two-liquids scenario for liquid water assumes the existence of two competing preferential local molecular structural states characterized by either low or high local density. While the former is expected to present good local order thus involving privileged structures, the latter is usually regarded as conforming a high-entropy unstructured state. A main difference in the local arrangement of such “classes” of water molecules can be inferred from the degree of translational order between the first and second molecular shells. This is so, since the low local density molecules present a clear gap between the first two shells while in the case of the high local density ones, one or more molecules from the second shell have collapsed towards the first one, thus populating the inter-shell region. Some structural indicators, like the widely employed local structure index, LSI, and the recently introduced ζ index, have been devised precisely on the basis of this observation, being successful in detecting well-structured low local density molecules. However, the nature of the high local density state has been mainly disregarded over the years. In this work we employ molecular dynamics simulations for two water models (SPC/E, extended simple point charge model, and TIP5P, fivesite model) at the liquid and supercooled regimes combined with the inherent dynamics approach (energy minimizations of the instantaneous configurations) in order both to rationalize the detailed structural and topological information that these indicators provide and, additionally, to advance in our understanding of the high density state.

Over the years, the two-liquids scenario proposition [15,16] has motivated a broad scrutiny for local order in water [21,22,[28][29][30][31][32][55][56][57][58][59][60][61]. But even when this picture, as already indicated, assumes the existence of domains of both structured (low local density) and "unstructured" (high local density) molecules, the emphasis has been overwhelmingly placed on the structured state. In this regard, the anomalies in water behavior have been recently reinterpreted in terms of locally favored structures [60,61]. Such privileged structures imply the existence of a structured state involving molecular arrangements with both low energy and degeneracy [60,61]. In contrast, the high density state has been expected to be constituted by thermally excited configurations characterized by high energy, high disorder and high degeneracy [60,61]. To shed light on this subject, different structural indicators have been devised in the last decades to discriminate between the two kinds of water "species" [21,22,[28][29][30][31][32][55][56][57][58][59][60][61]. As indicated, such indices have specifically focused on the detection of well-structured molecules, sweeping under the carpet the "unstructured" ones by classifying as "unstructured" everything that lacks a well-developed local tetrahedral coordination. This procedure precludes the comprehension of the nature of the high density state, usually simply regarded as a plethora of very different distorted and unrelated configurations. However, the possible existence of two extreme characteristic length scales of interaction in water, thus yielding two kinds of preferred local structures or arrangement topologies, is compelling, as suggested by the study of Walrafen-like pentamers [62,63]. Thus, an interesting question that emerges in this context is how unstructured is, actually, the widely overlooked high density state. This state is characterized by the intrusion of a second-shell neighbor or "interstitial" molecule towards the first coordination shell whose order is therefore altered. Thus, structured and unstructured states should differ in the degree of translational order between the first and second shells of the water molecules. This fact is at the heart both of the popular structural indicator known as the local structure index [21,22,[28][29][30][31][32] and of the recently introduced ζ index [60,61]. The availability of accurate and reliable structural indices is paramount for the context of supercooled water, in particular to test theories of glassy behavior and of ice formation and the two-liquids scenario. On one hand, proper deconvolution of the index distribution should be necessary for the determination of the fraction of structured and unstructured molecules. Moreover, in case that the index distribution yielded a bimodal behavior (which is not the usual situation for all of the currently employed indicators), it would be possible to further classify structured and unstructured molecules beyond determining their relative fractions. In particular, it is interesting to note that the ζ index and the LSI have been built upon idealized pictures where the fifth neighbor or first non-hydrogen bonded molecule is expected to play a central role (particularly its distance to the central molecule and its intrusion within the gap between the first and second coordination shells). However, more complicated local arrangements are usually involved in the distortion of the first two coordination shells, as we shall learn later on. Thus, in this work, we shall carefully study the instant configurations and the inherent structures (configurations from the real dynamics but subject to energy minimization) for SPC/E (extended simple point charge model) and TIP5P (five-site model) water models in order to reveal subtle aspects on the information provided by such indicators and also on the nature of the local arrangements encompassing the generally disregarded high density state. : CLASSIFYING MOLECULES BY  MEANS OF TWO STRUCTURAL  DESCRIPTORS, THE LOCAL STRUCTURE  INDEX AND THE ζ INDEX We conducted molecular dynamics simulations of both SPC/E and TIP5P water models by using the GRO-MACS package version 5.0.2. All the bonds were constrained using the LINCS algorithm and long range electrostatics were evaluated with the PME method. We used a modified Berendsen thermostat and a Parrinello-Rahman barostat at 1 bar as reference pressure. All the dynamics were run within cubic boxes of appropriate sizes with periodic boundary conditions and a cutoff of 1nm for the short range forces. The SPC/E system consisted of 216 water molecules while the TIP5P box contained 512 molecules. After equilibration for times much larger than the corresponding α relaxation time, production data runs were carried out for simulation times of at least 30 ns in each case. In this work we focus on two structural indicators that are sensitive to the translational order in the region between the first and second shell of the water molecules. Low density well structured molecules present a clear gap between the first and second shell, while the intrusion of one or more molecules from the second shell disturbs (in some cases even significantly) the tetrahedral order of the first shell of the high density molecules. The local structure index, LSI [21,22], has been devised to actually sense the presence or absence of such a gap between the first and second coordination shells. This index has recently regained considerable attention after we have applied it combined with local minimisations of potential energy, that is, the index is calculated not at the instantaneous configurations (real dynamics) but at the so-called inherent structures, IS. As long ago introduced by the seminal works of Stillinger and Weber [64,65] the potential energy surface (PES) of a many particle system can be partitioned into disjoint basins, where a basin is unambiguously defined as the set of points in configuration space connected to the same local minimum, called inherent structure (IS) via a minimization trajectory. Thus, an IS is the lowest energy configuration of a basin of attraction in the PES (the bottom of the basin) [17,[64][65][66]. In simple terms, a IS represents a structure of the system once the thermal energy distortions have been removed form the instantaneous configuration. In practical terms, for any given instantaneous configuration of the molecular dynamics we shall perform an energetic minimization by means of a steepest-descent algorithm (until energetic convergence is achieved) in order to get the corresponding IS. Within the IS scheme, the LSI is able to produce a clear bimodal distribution of structured and unstructured molecules both for the supercooled and normal liquid regimes [28] for the different water models where it has been applied, as for example, SPC/E [28,29], TIP4P-ice [67], TIP5P-Ew [32], TIP4P/2005 [68] and ab initio water model potentials within the framework of density functional theory [69].

II. METHODS
To compute the LSI index, which we denote as I(i, t), for a central molecule i at time t, one orders the rest of the molecules depending on the radial distance r j between the oxygen of the molecule i and the oxygen of molecule j : r 1 < r 2 < r j < r j+1 < · · · < r n(i,t) < r n(i,t)+1 , where n(i, t) is chosen so that r n(i,t) < 3.7Å < r n(i,t)+1 (that is, we consider the the n water molecules closest than 3.7 A to the central molecule). Then, I(i, t) is calculated as follows: where ∆(j; i, t) = r j+1 −r j and ∆(i, t) is the average over all molecules of ∆(j; i, t) (an average over all the n water molecules closest than 3.7Å to the central molecule). Thus, I(i, t) expresses the inhomogeneity in the radial distribution within the sphere of radius around 3.7Å. A high value of I(i, t) implies that molecule i at time t is characterised by a tetrahedral local order and a low local density, while on the contrary, values of I(i, t) ∼ 0 indicate a molecule with defective tetrahedral order and high local density. As already indicated, differently from Refs. [21,22], we calculate such index at the inherent structures, IS, by minimizing the potential energy of the corresponding instantaneous structure in order to get rid of the randomizing effect of the thermal vibrations, thus removing the fluctuations that prevent a proper identification of the local structure in the real trajectory. In turn, the ζ index [60,61] is computed by measuring the difference between the distance d j i of the closest neighbor molecule j not hydrogen-bonded to the central molecule i (for simplicity, from now on we shall call this first non-hydrogen bonded molecule as "1NHB"), and the distance d j i of the farthest neighbor molecule j that forms a hydrogen bond (HB) with molecule i (which we shall term as "LHB", last hydrogen bonded molecule): ζ(i) = d j i − d j i , where we shall consider that two water molecules form a HB when the O-O distance is below 3.5 A, and the O-H...O angle is larger than 140 • . This index has the advantage, from a conceptual point of view, of explicitly taking into account hydrogen bond coordination, a major contributor to water structure. It produced clear bimodal distributions for the TIP5P water model within the real dynamics, that is, the instant configurations (without the need to resort to the IS scheme) [61]. For TIP4P/2005, in turn, the distributions did not yield such a clear bimodality, but nevertheless could be easily decomposed in two gaussian functions [60]. This measure produces values slightly larger than 1Å for low density or structured molecules (for which the 1NHB molecule is located within the second coordination shell, which is separated from the first one by a clear gap), but significantly lowers for high density or unstructured molecules, where the 1NHB molecule is placed in the region between the first and the second shell (the index even reaches values close to zero or negative, when the 1NHB molecule intrudes within the region of the first coordination shell). . 1 shows the distribution of the ζ index for a range of temperatures for the two models studied, TIP5P and SPC/E (whose melting temperatures have been estimated to be respectively T m = 272 K and T m = 214 K [70]) both at the real (instantaneous) and inherent (inherent structures, IS) configurations. We note that this indicator has not been employed so far at the inherent structures level (by so doing, we shall now reveal relevant structural information concerning the high density state while also helping better rationalize the way in which the structural indicator is working). For TIP5P at the real dynamics we observe a nice bimodality with a right peak (whose population increases as T decreases) comprising the structured molecules and a left one (with an amplitude that decreases as T decreases) containing the unstructured ones. However, and at odds with the situation for the LSI index (whose distribution at the IS level is included for comparison), we note that with the minimization procedure (that is, the calculation of the ζ index at the IS level) it significantly looses bimodality, mainly since the left peak gets much broader and suffers a significant outwards displacement. When we turn to the SPC/E model (for which the ζ index has never been calculated before not even at the real dynamics) we find no bimodality (neither at the real dynamics nor at the IS level), with a simple peak that moves to higher values as T is decreased. Notably, the peaks for the different temperatures are much closer to each other at the inherent dynamics since, again, the peaks for the higher temperatures experiment clear displacements towards larger index values. This fact speaks of the existence of significant changes during the minimization procedure. More interestingly, the minimization seems to be affecting significantly the local arrangements of the unstructured molecules, while that of the structured ones might only suffer modest refinements.In Fig. 1 we also display the corresponding situation for the LSI index at the IS scheme for both water models. It is immediately evident that in both cases the index yields a marked bimodality. In turn, it is well known that the TIP5P model favors structuring by improving local tetrahedral hydrogen bond coordination while SPC/E represents an under structured model. To this end, it is interesting to consider Figs. 1 (b) and (d), that present the ζ index distributions for TIP5P and SPC/E at the IS scheme, and Figs. 1 (e) and (f), which display the corresponding cases for the LSI index. If we compare the peaks at low temperature in Figs. 1 (b) and (d) we learn that the ζ index for TIP5P implies very well structured molecules with values of the index larger than 1Å, while for the case of SPC/E this peak is well below such value. Additionally, the peak for structured molecules is also higher in the case of TIP5P. Finally, when we compare the structured peak in Figs. 1 (e) and (f) we can learn that the fraction of structured molecules (area under such peak) is much larger for TIP5P as compared to SPC/E at the lowest temperatures. However, it is also noticeable that TIP5P underestimates a bit the fraction of structured molecules at high T.

Fig
To investigate on the molecular basis of the local rearrangements responsible for the significant changes produced by the minimization process on the probability density distributions of Fig. 1 and, in turn, to shed light on the nature of the high density (unstructured) state, in Fig. 2 we show plots of the location (distance from the central molecule) of the two molecules involved in the ζ index calculation: the LHB (last hydrogen bonded) and the 1NHB (first non-hydrogen bonded) molecules. Fig. 2 (a) shows such distances when the molecules are classified as 1NHB and LHB at the real dynamics (and thier distances to the central molecule are also calculated the real dynamics), while the plot of Fig. 2 (b) corresponds to the situation where the molecules are classified at the inherent dynamics configurations (and distances are also measured at such IS level). However, since during the minimization procedure the molecules that embody both classes could change for a given central molecule, in Fig. 2 (c) we show the results when we classify molecules as 1NHB and LHB at the real dynamics and measure their distance to the central molecule at the corresponding IS configurations, that is, after minimization (regardless of the fact they remain or not being the 1NHB and LHB at the IS configurations).
At the real dynamics, Fig. 2 (a), we can note that the LHB molecules present a peak consistent with their engagement in a tetrahedral coordination with the central molecule (located within the first coordination shell and forming a good quality HB with the central molecule). However, the distribution also displays a tail to the right, presenting a non-negligible population at distances far-FIG. 2: Distances form the central molecule for the 1NHB and LHB molecules for SPC/E at T=240 K when the molecules are classified at: (a) the real dynamics (1NHB: blue curve with circles; LHB: red curve with squares); (b) at the inherent dynamics (1NHB: blue curve with stars; LHB: red curve with triangles). The case (c) corresponds to the situation when we classify the molecules at the real dynamics but when their distances to the central molecule are not calculated at the real dynamics as in cases (a) and (b) but at the corresponding inherent dynamics, that is, after the real configurations are subject to the minimization process (1NHB: blue curve with stars; LHB: red curve with triangles; we also include the curves of case (a) for comparison). Case (d) is idem to (c) but for the TIP5P water model at T=250 K. We note that the temperatures chosen imply that the SPC/E system is slightly above its melting temperature, T m , while the TIP5P system is slightly below it. We have deliberately elected them to illustrate different situations for different water models to show that they yield a qualitatively similar behavior. Similar results are also obtained for other temperatures close or below T m in both cases. In particular, for TIP5P we have chosen a temperature value below T m , since the ζ index displays a nice bimodality for the supercooled regime at the real dynamics (while for temperatures above T m such bimodality is practically lost).
ther than 3Å form the central molecule. In turn, the distribution of the 1NHB molecules exhibits a peak close to 3.4Å but it also "fattens" to the left, displaying significant population at distances lower than 3Å. However, from direct inspection of Fig. 2 (b), we learn that for both cases such region loses most of its population when the index is calculated at the inherent structure configurations: the distribution of the LHB molecules has mostly decayed at distances larger than 3.1 − 3.2Å while the 1NHB molecules do not display appreciable population at values shorter than that. Thus, the first lesson we can learn is that there are not many favored locations for these molecules at the inherent dynamics but only two different separation distances from the central molecule: the one corresponding to the first minimum in the radial distribution function (first-shell position) and a second location at roughly 3.5Å (characterized by a peak with a nice Gaussian shape in the plot of Fig. 2 (b)). Thermal fluctuations make the positions of the molecules deviate from these two locations at the real dynamics.
The significant differences between the plots of Fig. 2 at the real and the inherent dynamics imply, from one hand, that the LHB molecules move inwards upon minimization, improving their hydrogen bond with the central molecule, as expected. Besides, it could also be the case that the minimization procedure expelled certain 1NHB molecules (the interstitial molecules) from their disturbing positions (close to the first shell) towards roughly 3.5Å. While this happens in many cases, it is not always so, as we can learn from the fate of the molecules classified as 1NHB at the real dynamics (not at the IS) and analyzed at the IS configuration (Fig. 2 (c)). From such plot, we can see that the distribution for the 1NHB molecules splits into two clear peaks. In fact, certain 1NHB molecules move inwards to the position of a typical first neighbor (thus forming a good quality HB with the central molecule). This means that when we calculated the ζ index at the inherent dynamics for this case (that is, when we classified molecules as 1NHB and LHB directly at the IS; Fig. 2 (b)), the formerly 1NHB molecule (at the real dynamics) had lost its nature upon minimization transforming into a tetrahedrally hydrogen-bonded neighbor. This also means that another molecule had become the 1NHB molecule upon minimization. When we analyze the HB coordination of the central molecule for this situation, we learn on a clear prevalence of 3 HB at the real dynamics. Thus, the 1NHB at the real dynamics is not strictly speaking an interstitial molecule but represents the fourth neighbor that has lost its HB as a consequence of the distortion of the tetrahedron caused by the actual interstitial. It is noteworthy that the ζ index at the real dynamics (yielding small or even negative values) implies for this case the calculation of the difference between the distances of the forth and third neighbors, and not between the fifth and fourth ones as expected from the original proposition [60,61]. In turn, Fig. 2 (c) shows that other 1NHB molecules indeed move towards 3.5Å or more (second peak). Upon minimization, these fifth-neighbor (typically interstitial) molecules decrease or cease their disturbing effect on the first shell of the central molecule, which thus improves its local tetrahedral structure. For these cases, the central molecule shows a prevalence for 4 HB at the real dynamics and the local tetrahedron is improved by the minimization procedure. The ζ index calculation indeed involves the fourth and fifth neighbors in this case. These results are similar when we use TIP5P (Fig. 2 (d)) instead of SPC/E ( Fig. 2 (c)). While the distribution of the 1NHB molecules already displays certain bimodality at the real dynamics for TIP5P, the region around 3.0Å is again notably cleared by the minimization procedure. Such region is also depopulated in Fig. 2 (d) for the LHB molecules. This fact reveals the reluctance of the water molecules to stay in such region, since the minimization distributes them between a typical first shell position and an outer peak (which for SPC/E is centered at 3.5Å while for TIP5P it is displaced to a distance a bit larger). The clearing of such region upon minimization also explains why the LSI index produces bimodal distributions at the IS scheme. While the first shell is improved, the fifth neighbor is displaced outwards. Thus, if the minimization procedure places such molecule beyond 3.7Å, that is, at a distance commensurate with the second shell and that is used as a cutoff for the index definition, then this fifth neighbor is the last neighbor incorporated in the index calculation. And since its distance from the fourth neighbor is large, it contributes with a significant weight to the index calculation in Equation (1), thus producing a large index value characteristic of a well structured central molecule. But if the fifth neighbor falls below 3.7 A thus retaining an interstitial nature (the minimization procedure is not able to expel it to the second shell), the index calculation also includes the sixth neighbor thus incorporating an extra pair in the calculation of Equation (1). Additionally, since there are no large distances between any of the consecutive pairs considered, the index renders a small value consistent with a unstructured central molecule. Thus, while the LSI does a good job in unambiguously determining structured and unstructured molecules at the IS scheme for all water models studied so far, this classification evidently depends on the cutoff used [29]. Also, since for the TIP5P model the fifth neighbor is expelled to larger distances by the minimization procedure, the LSI index estimates a higher fraction of structured molecules as compared to the situation for the SPC/E model [32].
To further investigate on the reasons for the bimodality in the distribution of the 1NHB molecules at the IS level, Fig. 2 (c), we now study their orientation and possible interaction with their corresponding central molecule and their closest first-shell molecule, which from now on we shall call as the CFS molecule. The CFS molecule forms a HB with the central molecule, that is, it pertains to its first shell. Being the first-shell molecule that is closest to the 1NHB, it might form a hydrogen bond with it in certain configurations, as would be the case, for instance, in a well-structured ice Ih-like configuration (a well structured central molecule would also have its 1NHB molecule placed at a large distance so that it does not perturb the central molecule's tetrahedral first shell). When considering the orientation of the 1NHB molecule with respect to the central one we define the angle α as that formed by O i ...H...O 1N HB , where the subscript i represents the central molecule and the subscript 1NHB indicates that such oxygen atom corresponds to the 1NHB molecule. In turn, the H atom might belong either to the central molecule or to the 1NHB molecule. We must bear in mind that the formation of a good quality HB between the 1NHB and the central molecule (α > 140 • and an O-O distance lower than 3.5Å) is forbidden by definition. But in certain situations, a weak or deformed HB interaction cannot be ruled out. However, if the angle α is large (we consider an arbitrary threshold of α < 120 • , but the results do not change appreciably by small changes in such parameter, thus providing correct qualitative insights), we can safely regard such orientation as incompatible with any kind of HB interaction regardless of the distance between the two oxygens. We shall indicate the 1NHB molecules that belong to such class as C1 molecules. In Fig. 3 (a) and (b) we show the distributions of the distances to the central molecule of the C1 molecules and of the rest of the 1NHB molecules (all other angular orientations). We consider the cases when the molecules are classified as 1NHB in the real configurations and their distance to the central molecule is also calculated at the real dynamics ( Fig. 3 (a)) and, again, when they are classified at the real dynamics but their distances are evaluated at the corresponding inherent structures (that is, after the minimization procedure, Fig. 3 (b)). From Fig. 3 (a) and (b) we can learn that the C1 molecules (that is, 1NHB molecules with an orientation incompatible with even a weak HB interaction with the central one) display a simple peak at a distance close to 3.5Å at the real dynamics and that this peak is virtually unaltered by the minimization procedure. At variance, the rest of the 1NHB molecules display a bimodal distribution of distances to the central molecule in Fig. 3 (a) which is clearly improved by the minimization procedure ( Fig. 3 (b)). This result is compatible with the findings of Fig. 2. By analyzing the first coordination shell of the corresponding central molecules at the real dynamics, we find that the peak at small distance in the blue curve with circles of Fig. 3 (a) mainly corresponds to central molecules with only three good-quality HBs at the real dynamics (and thus are able to form a fourth HB with the 1NHB molecule upon minimization), while the peak at large distance mainly comprises central molecules already 4-fold HB-coordinated. In a similar manner we now analyze the 1NHB molecules in relation to the their interaction with their CFS molecule. Since this pair of molecules might form a good quality hydrogen bond, we simply classify 1NHB molecules in this situation as C2 molecules while the rest would be the ones that do not from such a good quality HB with the CFS molecule. Again, from Figs. 3 (c) and (d) we can learn that the distribution of distances of C2 molecules to the central one displays a peak slightly below 3.5Å at the real dynamics which is virtually unaltered by the minimization, while the rest of the molecules yield a broader peak that is split in two clear peaks by the minimization, in a way analogous to that of Fig. 2 (c).
In short, the results of Figs. 3 tell us that the 1NHB molecules that are hydrogen bonded to their CFS molecule and that are virtually unrelated to the central one are located at the larger distances from it and that, in turn, this location is not altered by the minimization procedure. However, 1NHB molecules that lack a good quality HB with their CFS molecule and that present an orientation compatible with a weak HB interaction with the central molecule, are usually found at closer locations from the latter (usually typical interstitial locations), while the minimization procedure might significantly alter their position. In fact, some of such 1NHB molecules restore a good quality HB interaction with the central molecule upon minimization.
In turn, to learn on the interactions between these three kinds of molecules (central, 1NHB and CFS molecules) on quantitative grounds, in Fig. 4 we now plot the interaction potential energy between the 1NHB molecule and the central one and also between the 1NHB molecule and its CFS molecule, as a function of the distance of the 1NHB molecule to the central one (at the real dynamics, including both coulombic and Lennard-Jones terms). We can observe that below 3.5Å the attraction of the 1NHB molecule with the CFS begins to decrease, while towards 3.1 or 3.2Å it starts to feel the attraction of the central molecule. In turn, the sum of these two curves displays a local maximum (a barrier) close to 3.1Å which is consistent with the above-described reluctance of the 1NHB molecules to reside at such distance. The interaction potential between the 1NHB and the central molecule at low distances corresponds to situations in which the latter has formed a good quality HB with the central molecule and, thus, it cannot be classified anymore as 1NHB molecule in such cases.
Additionally, Fig. 5 displays the probability distribution of the global stability of the 1NHB molecules (the sum of the interaction potential energy for the 1NHB; we sum all its pairwise potential interactions, both coulombic and Lennard-Jones terms, with the molecules located within a 8Å radius sphere) for two cases: Fig. 5 (a) corresponds to 1NHB molecules classified at the real dynamics at T=240K (for the SPC/E system) and evaluating their energy also at the real dynamics. Fig. 5 (b), in turn, corresponds to the same molecules (classified as 1NHB at the real dynamics) but evaluating their energy after the minimization process (at the corresponding IS configurations). We plot energy distributions as a function of the distance of the 1NHB molecule from the central one (since we are calculating probability values, these plots reflect free-energy contributions). The main feature immediately evident from Fig. 5 (a) is that there is only one minimum slightly below 3.5Å that corresponds to a region that covers both the interstitial position of the 1NHB molecule and a typical second shell location. This region comprises molecules compatible with the formation of a HB with a first neighbor of the central molecule (its CFS molecule), which begins to be lost as the 1NHB molecule is placed closer to the central one (energy rises). The right hand side of this minimum (distances much higher than 3.5Å) is consistent with 1NHB molecules of structured low local density central molecules (1NHB molecules that thus belong to typical second coordination shell positions) while the left hand side is compatible with interstitial molecules of high local density central molecules (1NHB located between the typical first and second shell positions). In turn, for the case when the 1NHB molecules are classified at real dynamics but evaluated at the IS (Fig. 5(b)), we can observe two energy minima. One is located at slightly above 3.5Å and comprises again both 1NHB molecules of well structured central molecules and interstitials (which are dragged by the minimization procedure from regions close to 3.0 or 3.2Å to a larger distance from the central molecule; as also evident from Fig. 2 (a) and Fig. 2  (b)). However, it also yields another minimum closer FIG. 3: Distance of the 1NHB molecules from the central one with respect to their mutual orientation at the real dynamics (a) and (b) after minimization (classified as 1NHB at the real dynamics but evaluating their location after minimization). With a red curve with squares we depict the C1 molecules (that is, orientation completely incompatible with even a deformed HB with the central molecule) while the rest of the molecules are represented by the blue curve with circles. We also show the distance of the 1NHB molecules from the central one depending on the orientation (good quality HB formation or not) of the 1NHB molecule with respect to its closest first neighbor of the central molecule (CFS molecule) at the real dynamics (c). The case (d) corresponds to the situation after minimization. The black curve with stars corresponds to the C2 molecules (1NHB molecules with a good quality HB with their CFS) while the green curve with triangles depicts the behavior of the rest. In all cases we employ the SPC/E water model. T=240K.
to the central molecule that exhibits significant population. These 1NHB molecules cannot form a HB with the CFS molecule and interact with the central one until beginning to form a good quality HB when below 3.0Å. Such molecules were the ones that were close to 3.0 or 3.2Å in Fig. 2 (a) (real dynamics) and then populated the left peak in Fig. 2 (c) (by means of the minimization procedure). At this point, these molecules have left being 1NHB at the IS and have completed the tetrahedron of the central molecule. It is interesting to note that the plots of Figs. 2 (a) and (b) include all the molecules, both high and low local density ones. The 1NHB molecules for both classes of water molecules fall within the same minimum of Fig. 5 (a). The only difference is that the low local density ones contribute to the right hand side of this minimum. This distance for low local density water molecules might seem very small as compared to low local density ice-like structure, but we must bear in mind that we are dealing with the closest of the second-shell water molecules (the 1NHB, that is, the first molecule not hydrogen bonded to the central one). If we include in the calculation all the second-shell molecules, we get the plot of Fig. 5 (c), with a more extended minimum that goes to larger distances form the central molecule.
To summarize, the above-displayed results imply that both for structured and unstructured molecules there exists only one minimum in the potential energy landscape at the real dynamics. This minimum, in turn, splits into two basins upon minimization: One consistent with a tetrahedral coordination with the central molecule (2.8 FIG. 4: (a) Interaction potential between the 1NHB molecule and the central molecule (blue curve with circles) and between the 1NHB and its closest first neighbor of the central one (CFS molecule, red curve with triangles). The green curve with stars adds up both contributions. We show the case for SPC/E. T=240K A) and an outer position located roughly at around 3.5 A or larger. From a topological point of view, upon minimization these molecules are either hydrogen bonded to the central molecule or to their CFS molecule respectively, arrangements that are, in turn, distorted by the thermal fluctuations at the real dynamics. This means that while low density well structured molecules conform to favored local arrangements, the high density "unstructured" molecules also display clear structural or topological preferences. While thermal energy at the real dynamics distorts the local molecular arrangements, these configurations belong to a definite basin of attraction in the potential energy landscape. To make this evident, we now study Walrafen-like pentamers [62,63]. We consider two pentamers comprising a central molecule and its 1NHB molecule: One pentamer is formed by the central water molecule acting as its center, that is, we consider the central molecule together with its four first neighbors (ideally H-bonded in a tetrahedral manner). The other pentamer is the closest pentamer where the 1NHB molecule (that by definition cannot be H-bonded to the central molecule) acts as one of its four corners (the fact that the centers of two first neighbor pentamers virtually cannot approach at 1NHB-central molecule distances [62,63] precludes the 1NHB molecule from acting as a center of the second pentamer). Two extreme possibilities (that thus distinguish two characteristic length scales of interaction) are the so-called "handshake" configuration (with two corners of one pentamer approaching two corners of the other pentamer, thus forming the corresponding HBs) and a "tango" configuration (with one pentamer rotated 90 • with respect to the other, thus breaking the two interpentamer HBs) [62,63]. In addition to possible handshake (two interpentamer HBs) and tango (no interpentamer HB) configurations, there could also be pentamers linked by just one HB, which we shall term as "half-tango". Fig. 6 shows schemes of these pentamer configurations.
To learn on the possible prevalence of some of these configurations for high and low local density molecules we study archetypal configurations by defining typically high density (structured) molecules and typically low density ones for the central molecules under consideration. For the case of the TIP5P water model, the highdensity molecules were taken as that with a ζ value within a 0.1Å neighborhood of the left peak of the ζ index distribution at the real dynamics, while the low density molecules were that within a 0.1Å neighborhood of the right peak. In turn, for the SPC/E model at the real dynamics (whose index distribution does not present the bimodality shown by the TIP5P model) we take as representative low-density molecules that with a ζ value between 0.7 and 0.9Å, while representative high-density ones would be those whose index value is between -0.15 and 0.05Å. In turn, for the inherent dynamics case, we consider 0.9 to 1.1Å and 0.15 to 0.35Å respectively. From Fig. 7 (a) we can learn that for the real dynamics scheme in TIP5P, the low density (structured) molecules considered as centers of pentamers present a preference for 1 inter-pentamer HB (1 HB between their pentamer and the pentamer with their 1NHB molecule at the corner), that is, half-tango configurations. However, a high population of 2 inter-pentamer HBs (handshake configuration) can also be noted. In turn, the high density molecules exhibit an equivalent preference for 0 interpentamer HB (tango) and 1 inter-pentamer HB (halftango). On the other hand, the minimization (IS scheme) makes the low density molecules to show a clear preference for 2 HBs between pentamers, while the distribution for the high density ones peaks at 1HB (since the interstitials might have recovered their HB with their corresponding CFS molecule). Thus, at the IS, the low density state favors hand-shake-like configurations while the high density state presents a prevalence of half-tango-like ones. That is, the IS scheme more clearly distinguishes the two characteristic length-scales of interaction of the water molecules. We also found similar results for the SPC/E model, as can be inferred form Fig. 7 (b). Again, not only the low density molecules conform to favored pentamer-pentamer configurations but also the high density ones display neat topological preferences. While the inter-pentamer study incorporates certain order beyond a single water molecule, it is still not directly comparable with models based on higher-order considerations, like the cage-like model [71]. This analytical model of water [71] not only includes model states with pairs of water molecules, being either pairwise hydrogen-bonded or non hydrogen-bonded (pairwise contact or pairwise noninteraction), but it also incorporates a higher-order structure model, a 12-water hexagonal unit cell or hydrogenbonded cage. We note that the structured state for a water molecule we characterized in this work is compatible with both the hydrogen-bond pairwise interaction state and with the higher-order cage state while, in turn, the unstructured molecular state resembles the local arrangement of the other (non hydrogen-bonded) pairwise models. Also, the low density hands-shake structure we presented with pentamers that are hydrogen-bonded to each other (favored at low temperatures) would be compatible with the 12-water cage, while the tango-like interpentamer arrangement represents the unstructured state.

IV. CONCLUSIONS
In this work we have analyzed the degree of translational order between the first and second molecular shells of high local density water molecules (for SPC/E and TIP5P molecular dynamical simulations), together with the performance of two structural indicators devised to be sensitive to this information: the widely employed local structure index and the recently introduced ζ index. We have calculated the latter for a water model where it has not been used before, SPC/E, and we have additionally studied it at the inherent dynamics level (both for SPC/E and TIP5P), a context where it has not been employed before. In particular, we show that energy minimizations (obtaining the corresponding inherent structure configurations) enable us to clearly rationalize the outcomes of both the LSI and the ζ index from structural and topological points of view. In fact, we demonstrate that the distributions at the inherent dynamics show strong differences from their corresponding results at the real dynamics. We also show that the fifth neighbor is not always the most relevant one in the calculation of both indices, as has been usually assumed. In fact, the ζ index sometimes implies the third and fourth neighbors while the LSI index incorporates the sixth neighbor in many cases. By focusing on the behavior of the first non-hydrogen bonded molecule to the central one (which we called 1NHB molecule) located at the interstitial region between the two first molecular shells, we find that while thermal fluctuations broadly populate this region at the real dynamics, the minimization procedure signif- icantly clears it by distributing these molecules between only two locations or basins of attraction of the potential energy. One of these regions implies typical first shell positions where the 1NHB molecule recovers a HB with the central molecule. This corresponds to very unstructured high density configurations when, by effect of thermal fluctuations, an interstitial molecule has spoiled the arrangement of the first coordination shell to a point that one of the four HBs with the central molecule has been broken. In turn, the minimization procedure recovers this HB with the central molecule and expels the other interstitial to distances of around 3.5Å or larger from the central molecule. The other basin of attraction of the potential energy represents precisely this latter shallow minimum that extends roughly beyond 3.5Å from the central molecule (the minimum of this basin in located at around 4.5Å). The left hand side of this minimum comprises typical interstitial molecules characteristic of high density configurations (for which the minimization procedure usually recovers or improves their H-bonding with a molecule belonging to the first shell of the central molecule, which we called CFS molecule) while the right hand side of this region is occupied by second neighbor molecules typical of low density configurations where they establish a good quality HB with their CFS molecule even at the real dynamics. Thus, the IS scheme shows that neighboring water molecules are reluctant to be located at the region around 3.1 or 3.2Å from the central molecule where they cannot exhibit good HB-coordination neither to the central nor to a first-shell (CFS) molecule. Hence, it is evident that not only low local density molecules display structural and topological preferences, but high local density ones also exhibit certain constraints in such regard. Furthermore, this behavior is also supported by the study of Walrafen-like pentamer configurations which exhibits that both classes of molecules display clear preferential topological arrangements. GAA SRA and JMMO acknowledge support form CONICET, UNS and ANPCyT(PICT2015/1893).