Beyond Local Structures In Critical Supercooled Water Through Unsupervised Learning

Beyond Local Structures In Critical Supercooled Water Through Unsupervised Learning Edward Danquah Donkor,†,‡Adu Offei-Danso,†,‡Alex Rodriguez,∗,†,¶Francesco Sciortino,§and Ali Hassanali∗,† †The Abdus Salam International Center for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy. ‡Scuola Internazionale Superiore di Studi Avanzati (SISSA) – via Bonomea 265, 34136 Trieste, Italy. ¶Dipartimento di Matematica, Informatica e Geoscienze, Università degli studi di Trieste, via Valerio 12/1, 34127 Trieste, Italy §Dipartimento di Fisica, Sapienza Università di Roma, P.le Aldo Moro 5, 00185 Rome, Italy E-mail: alejandro.rodriguezga[email protected]; [email protected] Abstract The presence of a second critical point in water has been a topic of intense investigation for the last few decades. The molecular origins underlying this phenomenon are typically rationalized in terms of the competition between local high-density (HD) and low-density (LD) structures. Their identification often require designing parameters that are subject to human intervention. Herein, we use unsupervised learning to discover structures in atomistic simulations of water close to the Liquid-Liquid Critical point (LLCP). Encoding the information of the environment using local descriptors, we do not find evidence for two distinct thermody1 arXiv:2401.16245v2 [cond-mat.soft] 29 Mar 2024 namic structures. In contrast, when we deploy non-local descriptors that probe instead heterogeneities on the nanometer length scale, this leads to the emergence of LD and HD domains rationalizing the microscopic origins of the density fluctuations close to criticality. TOC image 2 The physics of the critical behavior of matter close to phase transitions remains one of the most cherished areas of study in both experimental and theoretical physics.1–4 One of the most lively areas of discussion in this regard, pertains to the microscopic origins of the complex phase diagram of water.5–8 Besides the rather well characterized liquid-gas critical point, a series of theoretical predictions over the last few decades have proposed the existence of another critical point - the Liquid Liquid Critical Point (LLCP) of water in the supercooled regime.9–11 The physics underlying this criticality is thought to be one of the essential ingredients for understanding the anomalies of water. Probing the molecular origins of this second critical point has been dominated by theoretical and numerical predictions due to the challenge of spontaneous nucleation of ice at supercooled conditions.12–14 Just over three decades ago, Poole and co-workers demonstrated using the ST2 water model,15 that deeply supercooled water showed the presence of two distinct liquid phases, with fluctuations between the two phases terminating at the LLCP.9Several groups have also shown in the last decade, using advanced sampling free energy calculations, that the ST2 model exhibits two distinct liquid phases.10,16–19 More recently, this has been bolstered by tour-de-force microsecond simulations of realistic classical models of liquid water,11 as well as ab initio neural network models of liquid water20 that give further evidence for a LLCP scenario at least, on the numerical front. On the experimental side, work on supercooled water under elevated pressures as well as pioneering sound velocity measurements appear to be breaking the boundary of the so-called no-mans land giving strong indications of the existence of a LLCP.21–23 One of the central holy-grails of understanding the possible polymorphic nature of liquid water, has been the use of locally-stable structures7,24–28 which are thought to be rooted in water’s unique hydrogen-bond network. Many of the anomalies in water have been rationalized in terms of a competition between two types of hydrogen-bonding structures.6,25 One is said to have a more ordered tetrahedral and therefore open structure, often referred to as a Low Density (LD) local configuration, and the other, disordered due to the presence of interstitial water molecules which 3 is referred to as a High Density (HD) local configuration. Numerous order parameters have been constructed in an attempt to identify and distinguish these two local environments.7,25,28–34 These order parameters have also been shown to be tightly coupled to the macroscopic density fluctuations that occur close to the critical point (CP35). Although they provide a manner in which to physically interpret the simulation data, these order parameters often require significant human intervention which necessarily involves chemical bias (and often some arbitrary cut-off in their definition). Furthermore, it is also not a priori clear whether the interpretations made through these parameters are transferable across different regions of the phase diagram. Recently, some of us proposed a protocol that streamlines an unsupervised learning procedure for liquids and applied it to study the structure of water at room temperature36 as well as the study of the excess proton in hydrochloric acid.37 In brief, the method involves a three-step process. We begin by encoding the information of local environments using local atomic descriptors computed from the Smooth Overlap of Atomic Positions (SOAP38), which preserve important symmetries when comparing different molecular structures.39 In the second step, these high-dimensional descriptors are subsequently processed through an algorithm that extracts the Intrinsic Dimension (ID)40 which is crucial to understand the embedding manifold of the data. In the final step, the ID is used to extract the high dimensional free energy of the system41 and identify the minima.42,43 For room temperature liquid water, we found a rather broad and rough landscape separated by small barriers on the order of thermal energy where the shallow minima arise from a continuum of local molecular structures that continuously connect the canonical low or high-density local environments.36 In this contribution, we apply this protocol to understand fluctuations in supercooled water. Specifically, we uncover the molecular origins of critical-like fluctuations using unsupervised learning, analyzing trajectories recently reported in reference [11] connected to the presence of a second-critical point in atomistic water models. The free energy landscape constructed using local SOAP descriptors results in a single minimum despite there being macroscopic fluctuations 4 of the global density. By systematically expanding the SOAP descriptor to include fluctuations on a length scale of up to 1 nanometer, we uncover non-local domains relevant to critical-like fluctuations in supercooled water. The free energy landscape close to the critical point evolves between the high and low-density macroscopic phases, through a complex topography which we link to collective fluctuations of chemical-based order parameters that include non-local information of the water network. Methods The trajectory used for our analysis obtained from reference [11] is a 40µslong NPT trajectory of 300 TIP4P/2005 water molecules produced using the Gromacs 5.1.4 software close to critical conditions (177K, 1751 bar). More information on the simulation conditions is detailed in the main and supplementary text of reference [11]. Additionally, we apply our analysis protocol to a large system (36424 TIP4P/2005 water molecules) in the NVT ensemble at a temperature and density of (180K, 1011.83 kg/m3) in order to explore the larger length-scale fluctuations in the density. Our unsupervised learning protocol developed in reference [36] involves encoding local environments of molecules in a local atomic descriptor, extracting the intrinsic dimension, and constructing a high dimensional point-dependent probability density function from which the thermodynamic information can be inferred. The details of this procedure are outlined in the paragraphs below. As indicated earlier, the first step in our analysis is to encode the water molecular environments in a local atomic descriptor. To this end, we use the Smooth Overlap of Atomic Positions (SOAP ) descriptor,38,39 which preserves rotational, translational, and permutational symmetries of our molecular environments. In brief, given an atomic environment χaround a central atom, one characterizes the local density as a sum of Gaussian functions with variance σ2centered on each of the neighbors of the central atom including the central atom itself: 5 ρχ(r) = ∑ j∈χ exp −r−rj 2 2σ2!(1) This atomic neighbor density can be expanded in terms of radial basis functions and spherical harmonics Ylm such that: ρχ(r)≈ nmax ∑ n=0 lmax ∑ l=0 l ∑ m=−l cnlmgn(r)Ylm(θ,φ)(2) where the cnlm are the expansion coefficients. The number of expansion coefficients one chooses to compute is bounded by the number of radial and angular basis functions (nmax,lmax). In practice, one defines a cut-off radius (rcut) for the atomic environment being considered. One can then define a rotationally invariant power spectrum (p), whose elements are: pnn′l=πr8 2l+1∑ m (cnlm)†cn′lm (3) Thus, the distance between two environments χand χ′is related to the SOAP kernel by the following expression: d(χ,χ′) = 1−KSOAP(p,p′)(4) where: KSOAP(p,p′) = p·p′ √p·p p′·p′(5) Using the Dscribe package,44 the local SOAP descriptor for a water molecule p(i)is formed by computing the power spectrum on only oxygen species within a cutoff radius (rcut = 3.7 Å) centered about each oxygen atom. The local SOAP descriptors encode fluctuations on the length scales around the first coordination shell. To explore non-local fluctuations, we form a glocal SOAP descriptor by taking an average of the SOAP descriptor for each molecule and its neighbors 6 within a distance rgloc, given as: prgloc (i) = 1 n n ∑ j=1 p(j)(6) Here, nis the number of neighbours within a distance rgloc of the water molecule i. We also consider pglobal(i)which is the SOAP descriptor obtained by averaging the descriptors of all molecules within a snapshot. Similar types of non-local descriptors have been previously used by Lechner and Dellago45 as a means to accurately include non-local structural information in crystalline solid-state systems. The quality, size and accuracy of the SOAP descriptors depend on the parameters that go into its definition. In particular, one needs to have a balance between the level of detail the descriptors encode and also the computational management of the datasets one uses. In this work, we compute the SOAP descriptors considering only oxygen species and with the following parameters: nmax = 8 , lmax =6 and σ=1.0 Å since it offers a good balance between the level of detail of the molecular environment encoded and the size of the descriptors. In section S1 of the Supporting Information we explain how the descriptors used in the ensuing analysis and the ones that are built to include hydrogens as well as the use of smaller σ, encode similar information. Besides the SOAP based local atomic descriptors, we are also interested in examining if and how well chemical based order parameters capture the relevant fluctuations in liquid water. Of the many order parameters, the ones of interest to us in this context were the qtet ,30,31 LSI,46–48 d5,49 ρvoro,50 ψ35 and ζ.6More details on these chemical order parameters are provided in section S6 of the Supporting Information as well as in the main texts of the referenced material. In data sets with numerous dimensions, the presence of correlations among variables describing each data point suggests that the system of interest likely lies on a manifold whose dimension (the Intrinsic Dimensionality of the data set) is much lower than the embedding dimension of the data. To illustrate this, consider a set of points in three dimensions – if distributed randomly, the Intrinsic Dimensionality (ID) would be three. However, correlations between coordinates could restrict data points to lie only on the surface of a sphere, resulting in an ID of 2. Computing the ID is closely tied to dimensionality reduction techniques,51–53 where the dataset 7 is projected into a lower-dimensional space for analysis, visualization, and interpretation. The ID denotes the minimum dimensionality in which the data can be projected by applying such techniques without significant information loss. A proper understanding of the ID guides the selection of the space to analyze system fluctuations. In our study, the ID is crucial for estimating a pointdependent density function, influencing the extraction of free energy, as elaborated later. In this work, we employed the Two-NN estimator,40 a recently developed technique estimating the ID based on information from the first and second nearest neighbors of data points. This method, successfully applied to various molecular systems,54–56 operates on the assumption that the density of a data point can be considered approximately uniform within the distance to the second nearest neighbor of a data point, demonstrating that the ratio of the second to the first nearest neighbor distances (µ=r2/r1) follows a specific distribution: P(µ) = d µd+1(7) Here, dis the ID. Assuming independence of sampled ratios µi, the ID can be estimated by maximum likelihood (other estimators are also possible) as: d=N ∑N i=1log(µi)(8) Where Nis the total number of samples in the dataset. Using SOAP distances, we estimated the ID of the water molecule environment. The ID represents the minimum number of independent order parameters needed to describe the environment, aiding in quantifying information gained or lost with different variables.57 The considerations of the previous chapter have a direct impact on the reconstruction of the free energy landscape of water. To this end, understanding relevant variables characterizing structural fluctuations is essential. A common strategy is to examine probability densities along chemicallyinspired variables like qtet, LSI, and d5.7,47,58 However, this assumes no information loss in the projection (something that cannot be strictly true if the number of variables employed is smaller 8 than the ID of the data) and that the variable correctly encodes the process of interest. Recent techniques automatically identify important degrees of freedom59,60 and construct free energies in high dimensions.41,61–64 For a detailed discussion, refer to a recent review.65 In this work, we employed the Point Adaptive k-nearest neighbor estimator (PAk),41 avoiding the need for projection and used successfully in studying complex molecular systems.55,56,66 The method uses the ID as a parameter to construct a point-dependent density (ρi). This density is computed by adding a linear correction to the standard k-nearest neighbor estimator, where the density is ρi=ki rd ki , and ki’s are computed for each data point as the larger neighborhood for which the density can be considered approximately constant. The rationale is that, at constant density, the variance of the density estimation scales with 1 √kiwhile the inclusion of regions with different densities introduces a bias term to the error, therefore the procedure controls the Bias-Variance trade-off. The point-dependent free energy is −Log(ρi). Previous work shows this method accurately estimates free energy errors up to dimensions as large as 8.41 With point-dependent free energies, independent minima in the free energy landscape (clusters) are determined using a modified density peak clustering algorithm (DPA),43 an extension of the original density peak clustering.42 In this procedure, cluster center candidates are chosen as those whose density is maximum within their kineighbors. Then, the saddle points between these free energy basins are computed and the clusters are considered as coming from statistical fluctuations (and therefore merged in one) if the free energy difference between the basin minima and the saddle point is lower than Ztimes the sum of the errors associated to these free energy estimates. The parameter Zis the only free parameter in DPA clustering and can be interpreted as a measure of the statistical confidence of the clustering partition. The higher its value, the more can be one sure that the clusters are not coming from statistical fluctuations, but, at the same time, the higher the probability of losing real clusters whose statistical confidence is low due to the limited number of data points. In this work, the choice of Zwas made by varying it in two independently generated datasets until the clusters were consistent. Finally, PAkand DPA results are visualized and interpreted using the uniform manifold ap9 descriptor. This is not surprising since the SOAP descriptor by nature is complete and contains information about the molecular orientations which the macroscopic density does not contain. All in all, the preceding results builds strong evidence to a picture where the LD-HD density fluctuations cannot be described in terms of local competing structures but instead, involves clusters of at least 100 water molecules. Thus, according to our unsupervised learning protocol, the density fluctuations underlying LD-HD transitions cannot be associated with properties assigned at the single molecule level. With this picture in mind, we can revisit the UMAP projections providing more chemical interpretability. In panels A and B of Figure 4, we show the UMAP projection of the SOAP data in 2 dimensions, now colored with the corresponding average d5(<d5>Global) and qtet (<qtet >Global) respectively. We confirm from this, that the two density peaks (or free energy minima) emerging from our clustering correspond to the HD and LD phases since one of the peaks overlaps with <d5>Global∼3.4 Å hence <qtet >Global∼0.8 (High Density), and the other peak overlaps with <d5>Global∼3.8 Å and thus <qtet >Global∼0.9 (Low Density). Figure 4: (A) 2D UMAP representation of the SOAP data manifold colored with the average d5. The peak on the left half of this panel corresponds to low average d5while the peak on the right half corresponds to relatively high average d5and (B) 2D UMAP representation of the SOAP data manifold colored with the average qtet. The peak on the left half of this panel corresponds to low average qtet while the peak on the right half corresponds to relatively high average qtet. One of the important signatures of critical behavior is the divergence in the structure factor in the low | k|limit which ultimately translates into an enhancement in long-range density fluctuations. To investigate this anomalous scattering behaviour of water close to the critical point, Debenedetti 16 and co-workers explored the properties of the static structure factor using a large system (36424 TIP4P/2005 water molecules at sub-critical conditions run in the NVT ensemble11). Their analysis indeed shows the signature of critical behavior. A question that remains however, is how exactly one rationalizes the relationship between the non-local structures that emerge from our preceding analysis and the long-range density fluctuations. The two clusters that have been automatically identified from the SOAP descriptors averaged on the nanometer lengthscale, (see Figure 1 earlier) provides a protocol for classifying water molecule environments in other contexts such as those used to construct the structure factors previously described. Using the k-nearest neighbour classifer (see details in section S7 of the Supplemental Information), we assign water molecules in the large box simulation to either LD or HD type depending on their respective similarities. Applying this procedure leads to the automatic identification of LD and HD domains. In the left-most panel of Figure 5, we show one snapshot with only oxygen atoms (for clarity) colored by the phase they have been assigned to - blue spheres represent molecules assigned as LD-like while red spheres are molecules assigned as HD-like. By visual inspection, one can see a tendency for the LD and HD water molecules to cluster together forming LD-like and HD-like domains. These domains extend over spatial distances of several nanometers and essentially percolate throughout the periodic box, and to the best of our knowledge, this is the first instance where the LD and HD domains have been identified in a completely unsupervised manner. If the density fluctuations close to the critical point are indeed creating LD and HD domains then this implies that there should be some signature of an interfacial region forming at the boundary of the domains. One signature of this would be that water molecules close to the boundary would not be classified as pure LD or HD environments. A manner in which this can be quantified is to measure the probability of identifying either an LD or HD environment and subsequently identifying pure LD environments as those with pLD >0.7 and the pure HD environments as those with pHD >0.7. Water molecules with 0.7>pLD >0.4 or 0.7>pHD >0.4 are then identified as those that are putatively assigned as boundary or interfacial points. 17 In the middle panel of Figure 5, we show the same simulation snapshot but now also coloring points that have been identified to be so-called boundary water molecules. From visual inspection we can see how the green molecules are typically located between LD and HD domains. These findings nicely demonstrate that our procedure of agnostically identifying environments with appropriately averaged SOAP descriptors on the nanometer lengthscale leads to the emergence of LD and HD domains which are identified at the same thermodynamic state point. In the right-most panel of Figure 5, we plot the PDF of the ρvoro order parameter constrained to the LD (blue full line), HD (red full line) and interfacial molecules (green dashed line). We note that the full distribution of ρvoro is unimodal and broad. However by restricting the distribution to the identified domains separately, we find that the peaks in the distributions are consistent with those associated with the LD and HD phases in the smaller box. It is also curious to observe interfacial molecules, which have ρvoro values peaked between the peaks of the LDL/HDL ρvoro distributions. Figure 5: (Left) Snapshot of the large system colored by which density phase it was assigned to; blue for LD and red for HD. We observe the LD and HD domains extend over 1nm spatial distance. (Middle) Same snapshot now coloring molecules that are found in the boundary between LD and HD domains in green. (Right) PDF of the Voronoi Density for all LD assigned molecules (Blue full line), PDF of the Voronoi Density for all HD assigned water molecules (Red full line) and the PDF of Voronoi Density values for all molecules assigned as interfacial molecules (Green dashed line). We note how the distributions are peaked towards low density, high density and intermediate density respectively, albeit with a huge overlap 18 Conclusion In this work, we have used unsupervised machine learning techniques to analyse data coming from molecular dynamics simulations of liquid water close to the second critical point, where one observes pronounced fluctuations of the global density between a High-Density (HD) and LowDensity (LD) liquid phase. We show that the free energy landscape in the space of local descriptors consists of one minimum despite the pronounced density fluctuations. This is rooted in the large heterogeneity of the local configurations sampled by the water molecules in both phases. However, by using descriptors that account for non-local information of the water network, bimodality in the free energy landscape emerges. We further confirm the importance of non-local information by deploying a statistical test which allows us to evaluate the strength of the mapping between different descriptors and the macroscopic density fluctuations. We find that the mapping is strongest when the descriptors are constructed to include information on approximately a nanometer length scale. Finally, armored with the non-local HD and LD structures that emerge from our analysis, we characterize the formation of HD and LD domains that is manifested in the anomalous scattering behaviour of water close to the second critical point. Our results bring forward important challenges in assigning and interpreting fluctuations of the hydrogen bond network in terms of single particle properties where longer range structural correlations are clearly more important. These findings should motivate more work32 in trying to understand the relationships between local-atomic descriptors and local-molecular chemically inspired parameters and how they change our understanding of fluctuations across the phase diagram of water. We believe our work provides a general framework for understanding water’s structural and dynamic properties in other scenarios where long range correlations may be important, such as at interfaces73–75 as well as under confinement.76,77 19 Supporting Information Information Imbalance between SOAP descriptors computed from different hyper-parameters, Clustering labels; from Local to Global descriptors, Intrinsic Dimension scaling; from local to global descriptors, Descriptor-Density coupling; comparison between sub-critical supercooled water and water at ambient conditions, SOAP descriptor - order parameter coupling; comparison between sub-critical supercooled water and water at ambient conditions, Description of chemical order parameters used in the main text, Test scores for the KNN model used in figure 5 of the main text. Acknowledgement AH thanks the European Commission for funding on the ERC Grant HyBOP 101043272. EDD thanks Alex Chen Yi Zhang for very fruitful discussions during this work. References (1) Corró, M.; El Hichou, A.; Cesari, E.; Kustov, S. Study of magnetic transitions in Dy by means of reversible Villari effect. Journal of Physics D: Applied Physics 2015,49, 015001. (2) Shamba, P.; Wang, J.; Debnath, J.; Kennedy, S. J.; Zeng, R.; Din, M. M.; Hong, F.; Cheng, Z.; Studer, A. J.; Dou, S. The magnetocaloric effect and critical behaviour of the Mn0. 94Ti0. 06CoGe alloy. Journal of Physics: Condensed Matter 2012,25, 056001. (3) Ekimov, E.; Sadykov, R.; Mel’nik, N.; Presz, A.; Tat’yanin, E.; Slesarev, V.; Kuzin, N. Diamond crystallization in the system B 4 C–C. Inorganic materials 2004,40, 932–936. (4) Noguera, C.; Mackrodt, W. A theoretical study of magnetic phase transitions in ultra-thin films: Application to NiO. Surface science 2006,600, 861–872. 20 (5) Shi, R.; Russo, J.; Tanaka, H. Common microscopic structural origin for water’s thermodynamic and dynamic anomalies. The Journal of chemical physics 2018,149, 224502. (6) Russo, J.; Tanaka, H. Understanding water’s anomalies with locally favoured structures. Nature communications 2014,5, 3556. (7) Cuthbertson, M. J.; Poole, P. H. Mixturelike behavior near a liquid-liquid phase transition in simulations of supercooled water. Physical review letters 2011,106, 115706. (8) Pipolo, S.; Salanne, M.; Ferlat, G.; Klotz, S.; Saitta, A. M.; Pietrucci, F. Navigating at will on the water phase diagram. Physical review letters 2017,119, 245701. (9) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Phase behaviour of metastable water. Nature 1992,360, 324–328. (10) Palmer, J. C.; Martelli, F.; Liu, Y.; Car, R.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Metastable liquid–liquid transition in a molecular model of water. Nature 2014,510, 385– 388. (11) Debenedetti, P. G.; Sciortino, F.; Zerze, G. H. Second critical point in two realistic models of water. Science 2020,369, 289–292. (12) Hou, Y.; Yu, M.; Shang, Y.; Zhou, P.; Song, R.; Xu, X.; Chen, X.; Wang, Z.; Yao, S. Suppressing ice nucleation of supercooled condensate with biphilic topography. Physical review letters 2018,120, 075902. (13) Murray, B.; O’sullivan, D.; Atkinson, J.; Webb, M. Ice nucleation by particles immersed in supercooled cloud droplets. Chemical Society Reviews 2012,41, 6519–6554. (14) Mishima, O. Liquid-liquid critical point in heavy water. Physical review letters 2000,85, 334. (15) Stillinger, F. H.; Rahman, A. Improved simulation of liquid water by molecular dynamics. The Journal of Chemical Physics 1974,60, 1545–1557. 21 (16) Liu, Y.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Low-temperature fluid-phase behavior of ST2 water. The Journal of Chemical Physics 2009,131. (17) Sciortino, F.; Saika-Voivod, I.; Poole, P. H. Study of the ST2 model of water close to the liquid–liquid critical point. Physical Chemistry Chemical Physics 2011,13, 19759–19764. (18) Kesselring, T. A.; Franzese, G.; Buldyrev, S. V.; Herrmann, H. J.; Stanley, H. E. Nanoscale dynamics of phase flipping in water near its hypothesized liquid-liquid critical point. Scientific reports 2012,2, 474. (19) Poole, P. H.; Bowles, R. K.; Saika-Voivod, I.; Sciortino, F. Free energy surface of ST2 water near the liquid-liquid phase transition. The Journal of chemical physics 2013,138. (20) Gartner III, T. E.; Piaggi, P. M.; Car, R.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Liquidliquid transition in water from first principles. Physical review letters 2022,129, 255702. (21) Amann-Winkel, K.; Gainaru, C.; Handle, P. H.; Seidl, M.; Nelson, H.; Böhmer, R.; Loerting, T. Water’s second glass transition. Proceedings of the National Academy of Sciences 2013,110, 17720–17725. (22) Kim, K. H.; Amann-Winkel, K.; Giovambattista, N.; Späh, A.; Perakis, F.; Pathak, H.; Parada, M. L.; Yang, C.; Mariedahl, D.; Eklund, T.; others Experimental observation of the liquid-liquid transition in bulk supercooled water under pressure. Science 2020,370, 978– 982. (23) Amann-Winkel, K.; Kim, K. H.; Giovambattista, N.; Ladd-Parada, M.; Späh, A.; Perakis, F.; Pathak, H.; Yang, C.; Eklund, T.; Lane, T. J.; others Liquid-liquid phase separation in supercooled water from ultrafast heating of low-density amorphous ice. Nature Communications 2023,14, 442. (24) Hamm, P. Markov state model of the two-state behaviour of water. The Journal of chemical physics 2016,145, 134501. 22 (25) Tanaka, H.; Tong, H.; Shi, R.; Russo, J. Revealing key structural features hidden in liquids and glasses. Nature Reviews Physics 2019,1, 333–348. (26) Tanaka, H. Liquid–liquid transition and polyamorphism. The Journal of Chemical Physics 2020,153. (27) Skarmoutsos, I.; Franzese, G.; Guardia, E. Using Car-Parrinello simulations and microscopic order descriptors to reveal two locally favored structures with distinct molecular dipole moments and dynamics in ambient liquid water. Journal of Molecular Liquids 2022,364, 119936. (28) Foffi, R.; Sciortino, F. Identification of local structures in water from supercooled to ambient conditions. The Journal of Chemical Physics 2024,160. (29) Errington, J. R.; Debenedetti, P. G. Relationship between structural order and the anomalies of liquid water. Nature 2001,409, 318–321. (30) Chau, P.-L.; Hardwick, A. A new order parameter for tetrahedral configurations. Molecular Physics 1998,93, 511–518. (31) Lynden-Bell, R.; Debenedetti, P. G. Computational investigation of order, structure, and dynamics in modified water models. The Journal of Physical Chemistry B 2005,109, 6527– 6534. (32) Martelli, F.; Leoni, F.; Sciortino, F.; Russo, J. Connection between liquid and non-crystalline solid phases in water. The Journal of Chemical Physics 2020,153. (33) Faccio, C.; Benzi, M.; Zanetti-Polzi, L.; Daidone, I. Lowand high-density forms of liquid water revealed by a new medium-range order descriptor. Journal of Molecular Liquids 2022, 355, 118922. (34) Montes de Oca, J. M.; Sciortino, F.; Appignanesi, G. A. A structural indicator for water built upon potential energy considerations. The Journal of Chemical Physics 2020,152, 244503. 23 (35) Foffi, R.; Sciortino, F. Correlated fluctuations of structural indicators close to the liquid–liquid transition in supercooled water. The Journal of Physical Chemistry B 2022,127, 378–386. (36) Offei-Danso, A.; Hassanali, A.; Rodriguez, A. High-dimensional fluctuations in liquid water: Combining chemical intuition with unsupervised learning. Journal of Chemical Theory and Computation 2022,18, 3136–3150. (37) Di Pino, S.; Donkor, E. D.; Sánchez, V. M.; Rodriguez, A.; Cassone, G.; Scherlis, D.; Hassanali, A. ZundEig: The Structure of the Proton in Liquid Water from Unsupervised Learning. The Journal of Physical Chemistry B 2023,127, 9822–9832. (38) Bartók, A. P.; Gillan, M. J.; Manby, F. R.; Csányi, G. Machine-learning approach for one-and two-body corrections to density functional theory: Applications to molecular and condensed water. Physical Review B 2013,88, 054104. (39) De, S.; Bartók, A. P.; Csányi, G.; Ceriotti, M. Comparing molecules and solids across structural and alchemical space. Physical Chemistry Chemical Physics 2016,18, 13754–13769. (40) Facco, E.; d’Errico, M.; Rodriguez, A.; Laio, A. Estimating the intrinsic dimension of datasets by a minimal neighborhood information. Scientific reports 2017,7, 12140. (41) Rodriguez, A.; d’Errico, M.; Facco, E.; Laio, A. Computing the free energy without collective variables. Journal of chemical theory and computation 2018,14, 1206–1215. (42) Rodriguez, A.; Laio, A. Clustering by fast search and find of density peaks. science 2014, 344, 1492–1496. (43) d’Errico, M.; Facco, E.; Laio, A.; Rodriguez, A. Automatic topography of high-dimensional data sets by non-parametric density peak clustering. Information Sciences 2021,560, 476– 492. (44) Himanen, L.; Jäger, M. O.; Morooka, E. V.; Canova, F. F.; Ranawat, Y. S.; Gao, D. Z.; 24 Rinke, P.; Foster, A. S. DScribe: Library of descriptors for machine learning in materials science. Computer Physics Communications 2020,247, 106949. (45) Lechner, W.; Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. The Journal of chemical physics 2008,129. (46) Shiratani, E.; Sasai, M. Growth and collapse of structural patterns in the hydrogen bond network in liquid water. The Journal of chemical physics 1996,104, 7671–7680. (47) Appignanesi, G. A.; Rodriguez Fris, J. A.; Sciortino, F. Evidence of a two-state picture for supercooled water and its connections with glassy dynamics. The European Physical Journal E2009,29, 305–310. (48) Malaspina, D.; Schulz, E.; Alarcón, L.; Frechero, M.; Appignanesi, G. Structural and dynamical aspects of water in contact with a hydrophobic surface. The European Physical Journal E2010,32, 35–42. (49) Saika-Voivod, I.; Sciortino, F.; Poole, P. H. Computer simulations of liquid silica: Equation of state and liquid–liquid phase transition. Physical Review E 2000,63, 011202. (50) Rycroft, C. H.; Grest, G. S.; Landry, J. W.; Bazant, M. Z. Analysis of granular flow in a pebble-bed nuclear reactor. Physical review E 2006,74, 021306. (51) Mika, S.; Schölkopf, B.; Smola, A.; Müller, K.-R.; Scholz, M.; Rätsch, G. Kernel PCA and de-noising in feature spaces. Advances in neural information processing systems 1998,11. (52) Jolliffe, I. T. Principal component analysis: a beginner’s guide—I. Introduction and application. Weather 1990,45, 375–382. (53) Kruskal, J. B.; Wish, M. Multidimensional scaling; Sage, 1978; Vol. 11. (54) Ansari, N.; Laio, A.; Hassanali, A. Spontaneously forming dendritic voids in liquid water can host small polymers. The journal of physical chemistry letters 2019,10, 5585–5591. 25 Figure S3: Information Imbalance between the different order parameters and the global density. Comparison between sub-critical supercooled water (left panel) and room temperature water (right panel). We show that the coupling between the order parameters and the global density is a feature of water in the supercooled regime as for water at ambient conditions we do not observe a decrease in the Information Imbalance with length scale of averaging. S5 Coupling Between different descriptors and the SOAP descriptor in sub-critical supercooled water in comparison to water at ambient conditions 32 Figure S4: Information Imbalance between the SOAP descriptors and the other order parameters, comparison between sub-critical supercooled water and room temperature water. We can see that the SOAP descriptors are tightly coupled to the other descriptors at larger length scales in the supercooled regime which is not the case at ambient conditions. S6 Chemically-Inspired Order parameters used in the study of the molecular structure of supercooled water •qtet measures the similarity between the first coordination layer and a tetrahedron.30,31 More precisely, qtet is defined by the following equation: qtet =1−3 8 3 ∑ i=1 4 ∑ j=i+1 (cos(φi j)+ 1 3)2(S1) where φi j is the angle formed by the lines joining the oxygen atom of the central water molecule to its four nearest neighbor oxygen atoms iand j. • The d5parameter is the distance between the fifth nearest neighbor to the central oxygen atom and reflects the extent of separation between the first and second solvation shells. A 33 larger value of d5is interpreted as being a more open and locally ordered structure49 and vice-versa, a more locally disordered and closed structure. • The LSI parameter was designed to distinguish environments with well-separated first and second coordination shells from those that are more disordered.46–48 Consider the distances between the oxygen atom of a central water molecule and the ith neighboring oxygen atom ordered in the following manner r1<r2< ... < ri<ri+1< ... < rn<3.7Å <rn+1. The LSI is then defined as: LSI =1 n n ∑ j=1 (∆(j)−¯ ∆)2(S2) where ∆(j) = rj−rj+1and ¯ ∆corresponds to the difference and average in consecutive distances respectively. Large values of LSI such as 0.3 correspond to structures with well separated first and second coordination shells while very low LSI values are consistent with interstitial waters between the two shells. • In order to measure local density variations, we computed the Voronoi density (ρvoro), which is the inverse of the Voronoi-volume associated with a water molecule. This volume is the sum of the volume of the oxygen and two hydrogen atoms.78,79 The Voronoi volume is found by performing a Voronoi tesselation on the water network. We carry out the Voronoi tessellation using the Voro++ code50 • The ψdescriptor, originally introduced in reference [35], is another descriptor which has been proposed to study topological arrangement of the water hydrogen bond network in supercooled conditions. More details on ψis presented in the original manuscript but in brief, ψmeasures the minimum physical distance between the reference oxygen atom of a water molecule and its neighbor located at a chemical distance D = 4. Where D is measured in units of number of hydrogen bonds. 34 • Finally, the ζdescriptor, introduced in reference [6] quantifies the separation between the first and second solvation shell around a water molecule by measuring the difference in the distance between the furthest hydrogen bonded molecule and the closest non-hydrogen bonded molecule around a central water molecule. Where two molecules are said to be hydrogen bonded using the geometric criteria by Luzar and co-workers.80,81 Several works have used this order parameter to classify water molecules into high ζ(Low Density) and low ζ(High Density) local environments.5,6,82,83 Classification of LD/HD environments from NVT simulations of sub-critical supercooled water Consider the data set consisting of N data points and D features. In our context each i∈{1,...,N} is the 10 Å glocal SOAP descriptor of a molecule sampled from the NPT trajectory that shows strong density fluctuations. From this data set our DPA clustering provides us with two clusters: an LD and an HD cluster. So for each data point, we have an associated label showing which phase it belongs to. Using these labels we perform a k-nearest neighbour classification (with k = 11) task on the 10 Å glocal SOAP descriptors computed from water environments sampled from the NVT trajectory with 36424 molecules. The aim is to identify the LD and HD domains. In figure S5 we show that the predictive power of the k-NN model does not depend significantly on the k used. However, we use a k = 11 to get good estimates of the probability (p) of being assigned as an LD or HD type water. After the classification task is carried out we obtain the probability (p) of being an LD or HD environment. The probability of being HD or LD are related (pLD +pHD =1). From these probabilities we classify the core LD environments as those with pLD >0.7 and the core HD environments as those with pHD >0.7. Then the interfacial or boundary molecules are labelled as those with 0.7>pLD >0.4 or 0.7>pHD >0.4. These cutoffs in the probabilities do not significantly affect the populations in the various categories of molecules. 35 Figure S5: Test Score as a function of the number of neighbors (k) used. Test scores are computed by splitting our 100000 10 Å glocal SOAP descriptors from the NPT trajectory into a training and testing set. The test score is thus evaluated on the test data for several ks. 36