Size-dependent effect of nano-confinement of water in an ionic liquid matrix at low temperature

Abstract

One of the leading hypotheses explaining water's anomalies is a metastable liquid–liquid phase transition (LLPT) at high pressure and low temperatures, which remains experimentally elusive due to homogeneous nucleation. Infrared spectroscopy experiments have shown that adding hydrazinium trifluoroacetate to water induces a sharp, reversible LLPT at ambient pressure, potentially originating from the same underlying mechanism as in pure water. In a previous work, we demonstrated that this transition can be attributed to the behavior of pure water only when nanosegregation of the aqueous component is brought into play. Here, by means of molecular dynamics simulations and the structural order parameter ζ, we explicitly analyze the effect of the ionic compound on the structure of liquid water at low temperature, both in a mixed solution and nanoconfined in spherical clusters of varying size. Our findings indicate that the ions surrounding the water induce structural perturbations that disrupt the water hydrogen-bond network up to a depth of approximately 0.70–0.75 nm from the surface toward the center of the sphere. This suggests that, in order to preserve a low-density liquid state within this ionic matrix, and more in general highly ionic matrices, water must be confined within pockets with radii greater than approximately 0.70–0.75 nm.

---

Introduction

Water is a key focus in many areas of research, including physics, biology, and chemistry. However, even under normal conditions, its microscopic structure remains a topic of intense debate.¹ At room temperature and pressure, water exhibits different local molecular arrangements, along with fluctuations in its hydrogen bond (HB) network.² This complexity increases when water exists as a supercooled liquid below its melting point, but still in a metastable state. Supercooling amplifies waters well-known thermodynamic and dynamic anomalies, driving significant research into their origins. Among various proposed theories the liquid–liquid critical point (LLCP) hypothesis is the most prominent explanation for these anomalies.^1,3

The LLCP scenario points to the presence of a first-order liquid–liquid phase transition (LLPT) between two phases: the low-density liquid (LDL) and the high-density liquid (HDL).³ In addition to their density differences, the LDL and HDL states also display distinct local structures. The LDL structure closely resembles that of ice (Ih), with four nearest neighbors arranged in a regular tetrahedral configuration. In contrast, the HDL structure presents a distorted HB pattern, featuring a fifth nearest neighbor molecule, often referred to as the interstitial water molecule. According to the LLCP scenario, the significant increase in thermodynamic response functions during supercooling is linked to the presence of the LLCP in the so-called “no man's land”, a region of the phase diagram in which ice nucleation occurs very rapidly, thus preventing the observation of the LLPT with conventional measurement techniques. Although there are several computational studies supporting this hypothesis, there is still no definitive experimental evidence for the existence of the LLPT.^4–7

Experimentally, a commonly adopted strategy to prevent ice crystallization consists in the use of binary aqueous solutions to shift the freezing point out of the “crystallization curtain”.^8–10 The addition of solutes, on the other hand, may lead to the suppression of neat water anomalies. However, in recent experiments it was shown that the use of hydrazinium trifluoroacetate (HYD-TFA 15.6% mol) aqueous solutions can effectively suppress crystallization, while keeping at the same time the heat capacity and density anomalies as for neat water.^8,9 By means of computed infrared (IR) spectra and comparison with available experimental data, our previous work¹¹ provided evidence that the HDL–LDL transition observed in these experiments is coupled to a mixing-demixing process. The hypothesis was that as the system is cooled, ions are progressively expelled from the LDL-like phase of water, leading to nanoscale separation into ion-poor/LDL-like and ion-rich/HDL-like regions. In a subsequent study combining molecular dynamics simulations with calorimetric measurements and X-ray diffraction,¹² we showed that even at 300 K, where the HYD-TFA solution is considered “fully mixed”, one third of the water molecules were found in dynamic clusters comprising, on average, 25 molecules, extending beyond a single hydration shell. “Soft” confinement of water (on the sub-nanometer scale) in supercooled aqueous sorbitol solutions has also been reported.¹⁰

Regarding aqueous solutions in the presence of other ionic liquids, as for example 1-butyl-3-methylimidazolium nitrate ([BMIM][NO₃]) and 1-butyl-3-methylimidazolium dicyanamide ([BMIM][DCA]), experimental studies have shown that water tends to cluster into structures commonly referred to as “water pockets”.^13–16 The size of these water pockets typically ranges on the order of a few nanometers, although the exact dimensions are strongly influenced by the properties of the ionic liquid and the composition of the mixture. Despite these insights, the structural and dynamic properties of water confined in highly charged ionic environments remain far from fully understood.^13–17

In this work, we aim to understand the effect of HYD-TFA ions on the structure of supercooled liquid water. To this aim, we performed Molecular Dynamics (MD) simulations at a temperature low enough to ensure that, were the water in bulk, it would adopt an LDL state. We first investigated the structural features of a solution in which ions and water are mixed. Then, we used spherical water pockets of different sizes immersed in an ionic liquid matrix of HYD-TFA to study the influence of the ions on the structural properties of the confined water. In computational works, to determine the structural features of liquid water, many order parameters have been proposed, such as the local structure index (LSI),¹⁸d₅,¹⁹ζ,²⁰q₄,²¹Ψ,²²S and δS²³ and V₄.^24,25 Recently, the node total communicability (NTC) has also been introduced.^26–28

Here, to characterize the local structure of the water molecules in the mixed solution and in pockets, we used the order parameter ζ to estimate the perturbative effect of the ionic liquid matrix on an initially ordered water configuration (i.e., an LDL-like structure). As discussed by Tanaka et al.,²⁹ζ incorporates hydrogen-bond formation and proton order-features for the characterization of local environments in liquid water. Given the nature of the water-ion interactions, such features are essential to properly describe the present binary solution. Consequently, we have chosen to use the ζ descriptor instead of other parameters that have been applied to binary systems before (e.g. LSI and q).³⁰ To the best of our knowledge, this is the first time that ζ has been used for this purpose. We also tested another order parameter, d₅, which does not explicitly include hydrogen bonds, to characterize water structure in the fully mixed solution. As shown in Fig. S1 of the supporting information (SI), the d₅ descriptor does not effectively discriminate the different structural arrangements in the solution.

Methods

Molecular dynamics simulations

We performed MD simulations of bulk water and of different HYD-TFA solutions. For bulk water, two simulations of 710 water molecules were performed: one at 180 K and one at 300 K, both at 1 bar. For the HYD-TFA solutions, different starting configurations were used. In the fully-mixed simulation, a proper number of HYD-TFA molecules were randomly inserted in a box of 710 water molecules in order to have a solution at a molar fraction χ_W = 0.844. In a second set of simulations, a spherical region of LDL-like water molecules with a given radius was generated, with a defined number of ions positioned at its surface to create an unmixed solution, while preserving the same water molar fraction as in the fully mixed system. The spherical water pockets were obtained by extracting spherical regions from an equilibrium configuration of bulk water at 180 K and 1 bar. Using this procedure, three spherical water pockets with radii of 1.25, 1.75, and 2.35 nm were generated, containing respectively 307, 745, and 1194 water molecules.

To equilibrate the water-ion interface in each system, a 20 ns-long simulation at 300 K was carried out while maintaining the positions of the water molecules within the spherical region fixed. This was followed by a 100 ps temperature annealing step, gradually lowering the temperature back to 180 K, prior to initiating the 100 ns production runs. For all MD simulations, the TIP4P/2005³¹ water model was used, as it was shown to exhibit a metastable liquid–liquid critical point under deeply supercooled conditions.^5,32–35 OPLS-AA force field parameters for HYD⁺ and TFA⁻ ions were used and generated with the LigParGen³⁶ tool with the exception of the atomic partial charges that were computed using a series of optimized HYD-TFA ion pairs as described in ref. 37 MD simulations were performed with the GROMACS package (version 5.1.2)³⁸ in the NPT ensemble using a cuboidal simulation box, the velocity rescaling temperature coupling,³⁹ and the Parrinello–Rahman barostat⁴⁰ with 2 ps relaxation times. Periodic boundary conditions were used, and the long range electrostatic interactions were treated with the particle mesh Ewald⁴¹ method with a real space cutoff of 0.9 nm. The Lennard-Jones potential was truncated at 0.9 nm. The LINCS⁴² algorithm was used to constrain bond lengths along with a 2 fs time step.

Infrared spectra calculation

To compute the infrared spectra, we used a hybrid quantum/classical methodology, the Perturbed Matrix Method (PMM),⁴³ to be used in conjunction with MD simulations. The MD-PMM has already been used to calculate the OH-stretching mode of water. In previous experiments,⁹ dilute H/D isotopic mixtures (H:D fraction of ≈ 0.03) were used to prevent coupling between the OH-stretching modes. Therefore, we calculate here the IR absorption band arising from the decoupled OH-stretching mode in HOD. A brief outline of the MD-PMM approach for the calculation of infrared spectra is provided in the SI of,⁴³ as well as details on the quantum mechanical calculations needed to apply the MD-PMM.

Order parameters calculation

All distances and hydrogen bonds in the following sections were computed using MDAnalysis.⁴⁴ Two atoms are considered hydrogen bonded if the distance between donor (D) and acceptor (A) is within 0.35 nm, and the angle H–D–A is less than 30°.²⁰

d ₅ calculation. Given a reference molecule i, d₅(i) is defined as the absolute distance d_ji, where j is the fifth nearest oxygen atom to the molecule i. We redefine the set of possible neighboring atoms to include, beyond oxygen atoms, the nitrogen atoms of hydrazinium and the oxygen atoms of trifluoroacetate, and compute d₅ based on this extended set.

ζ calculation. Given a reference molecule i, ζ(i) is defined as the difference between the distance d_j′i of the first neighbor j′ not hydrogen bonded to i and the distance d_j′′i of the last neighbor j′′ hydrogen bonded to i as reported in eqn (1).

Expanding on this concept, we consider as donor and acceptor atoms, along with the oxygen of water, the nitrogen atoms of hydrazinium and the oxygen atoms of trifluoroacetate.

Results and discussion

Before studying the structural properties of water (TIP4P/2005 model) in confinement, we first characterize reference systems such as: bulk water in the LDL-like state (180 K and 1 bar for the given water model), bulk water in the HDL-like state (300 K and 1 bar for the given water model), and a fully mixed solution of water and HYD-TFA (15.6% mol) at 180 K and 1 bar. To this aim, we use ζ and d₅ as order parameters. For the calculation of ζ in the presence of the ions, the oxygen atoms of TFA and the nitrogen atoms of HYD are included as potential HB donors and acceptors. However, the values of ζ are computed specifically for water molecules only. See the Methods section for more details. The distributions of ζ calculated for bulk water at 300 K and 180 K and along the MD simulation of the fully mixed solution are shown in Fig. 1. The LDL-like state (bulk water at 180 K) exhibits a sharp distribution centered at 0.945 nm, while the HDL-like state (bulk water at 300 K), characterized by the presence, on average, of an interstitial water molecule, features lower ζ values. These results are consistent with those shown in²⁰ for the same water model. In the fully mixed solution, the distribution appears bimodal, with values skewed toward the HDL-like structure, as shown in Fig. 1. On the other hand, the d₅ descriptor reported in Fig. S1 of the SI is not particularly informative, as it does not account for hydrogen bonds and provides no information about the four closest species to the reference molecule. Therefore, all of the following analyses are performed using the ζ order parameter.

Fig. 1 Distribution of ζ for LDL-like bulk water (180 K and 1 bar, blue), HDL-like bulk water (300 K and 1 bar, black) and water fully mixed with hydrazinium-trifluoroacetate at 15.6% mol at 180 K and 1 bar (orange). Top panel shows an example snapshot of the fully mixed water solution.

To understand the contribution of different structural organizations of water to the ζ distribution in the fully mixed solution, we characterize the local environment of each water molecule based on its hydrogen-bonded (h-bonded) and non-bonded partners. For each molecule, we assign a label in the format XXX–YYY, where XXX indicates the type of hydrogen-bonded partners (‘water’ if only water molecules are bonded, ‘HYD’ if at least one HYD molecule is present, or ‘TFA’ if at least one TFA molecule is present); YYY denotes the closest non-bonded neighbor, categorized as either ‘water’, ‘HYD’, or ‘TFA’. The results are reported in Fig. 2 in which only contributions greater than 10% of the total distribution are shown; the minor contributions can be found in the SI (Table S1 and Fig. S2).

Fig. 2 Contributions to the total distribution of ζ values for water molecules in the fully mixed solution (orange line) at 180 K, arising from distinct structural motifs. Only contributions ≥10% are shown. The sum of these four distributions is displayed alongside the total distribution (dashed and solid orange lines, respectively). The most frequently sampled structures corresponding to each of the four structural motifs are shown at the top (dashed lines represent hydrogen bonds). We label each water molecule as XXX–YYY, where XXX denotes the type of hydrogen-bonded partners (‘water’ if only water, ‘HYD’ if at least one HYD, ‘TFA’ if at least one TFA), and YYY the closest non-bonded neighbor (water, HYD, or TFA).

Concerning the water–water coordination motif in the mixed solution, the water molecules exhibit, on average, slightly lower ζ values compared to bulk water (see distribution a in Fig. 2 and the corresponding representative structure reported in panel (a). This suggests that the presence of ions perturbs the HB network even when they are not directly h-bonded in the first hydration shell. The presence of a TFA ion directly h-bonded to the water molecule (panel b in Fig. 2) shifts the ζ values to lower values, although the distribution still remains at roughly the same position as the water–water coordination motif (see distributions a and b in Fig. 2). The most pronounced disruption of the water coordination structure arises from the presence of non-bonded HYD ions within the first hydration shell (see representative structures in panels c and d of Fig. 2). In these configurations, a central water molecule loses one of its hydrogen bonds as a HYD ion intrudes into the first shell. Basically, the hydrazinium ion replaces the fourth hydrogen-bonded water molecule without establishing a hydrogen bond itself. As a consequence, the coordination parameter ζ decreases significantly, and the distribution maximum shifts to much lower ζ values (around zero). On the contrary, when HYD is the first non-bonded partner but does not intrude into the first hydration shell, i.e. when water forms four hydrogen bonds and hydrazinium is the closest non-hydrogen-bonded partner, the distribution of ζ corresponds to that obtained for the water–water case (Fig. S3 in the SI). Overall, comparison of curves a, b, c, and d in Fig. 2 shows that the difference between the ζ distributions in bulk water at 180 K and the fully mixed solution at the same temperature can be largely attributed to a decrease in the number of hydrogen bonds in the first hydration shell. A minor contribution to this difference arises from alterations in the nature of both hydrogen-bonded and non-bonded interactions, which modify the local coordination environment even when the number of hydrogen bonds remains unchanged.

Additional analyses of the hydrogen-bond network in the fully mixed solution are reported in Section S1 of the SI, including (i) classification based on the number of hydrogen bonds per molecule (Table S2 and Fig. S3, S4); (ii) contributions of hydrogen-bond patterns where the same ionic species acts as both a hydrogen-bond and non-bonded partner (Tables S3, S4 and Fig. S5); (iii) detailed statistics of hydrogen-bonding and partner distributions comparing 180 K bulk water with the fully mixed solution (Tables S5 and S6).

Overall, the HB network in our simulated fully-mixed solution differs significantly from that in pure water, and the presence of ions interacting with water appears incompatible with the maintenance of a LDL-like water structure. Therefore, elaborating on our previous hypothesis of ion-poor/LDL-like and ion-rich/HDL-like water regions, we investigate a set of solutions in which the ions/water concentration is the same as in the fully mixed solution but water is confined within a pocket surrounded by ions. The ζ distributions for the three water pockets, compared to that of the fully mixed solution, reveal distinct interaction patterns of water under confinement (Fig. 3). The low-ζ peak, associated with the presence of hydrazinium ions as non-bonded partners at the water-ion interface, is significantly reduced relative to the fully mixed solution. In contrast, the high-ζ peak, which primarily arises from water–water interactions within the confined water spheres, is more pronounced. The relative intensities of the two peaks vary slightly among the three pockets, reflecting differences in their surface-to-volume ratios. Additionally, the high-ζ peak in the distributions of the confined pockets is shifted to lower values compared to that of bulk LDL-like water. This shift is rationalized by analyzing the spatial dependence of ζ as a function of the distance from the water-ion interface.

Fig. 3 ζ distributions of water molecules confined in pockets with radii of 1.25 (black line), 1.75 (red line), and 2.35 nm (green line), shown alongside those for bulk water (blue line) and water in the fully mixed state (orange line), all systems are at 180 K. An example snapshot of a water pocket is shown as an inset in the figure.

To this end, we computed the ζ values of water molecules at decreasing distances from the water-ion interface for each of the three systems. Specifically, starting from the interface of each spherical pocket, we iteratively reduced the core radius by 0.1 nm and calculated ζ for the water molecules within each successive spherical shell (An example of a shell is shown in Fig. S6 of the SI). The corresponding distributions are reported in Fig. S7–S9 of the SI. In Fig. 4 we report the frequency at ζ = 0.0945 (i.e., the ζ value of the peak for bulk water at 180 K and 1 bar) obtained from the distributions at increasing distances from the ions. The figure shows that, as the distance from the ions increases, the structure of the water molecules converges toward that of bulk water at 180 K. This trend indicates that (i) in the deep core of the pocket the LDL-like structure is recovered; (ii) the ionic liquid has a strong perturbative effect on the structure of water, even in regions where the ions are not directly present (see Fig. S10 in the SI). The three trends in Fig. 4 indeed consistently show that the effect of ions on the structure of water extends up to 0.70–0.75 nm from the surface. Therefore, water molecules are affected by the presence of ions located beyond the second hydration shell of water (see Fig. S11 in the SI for the oxygen–oxygen and oxygen-hydrogen radial distribution functions). In addition, the radial distribution functions confirms the different organization of the local structure of water in the fully mixed system compared to bulk water. In fact, the shift in the second O_w–O_w peak and the loss of the features beyond the third O_w–H_w peak are consistent with a denser first hydration shell and a loss of long-range order, respectively. An additional support to the results obtained with ζ is given by the calculation of the infrared (IR) band arising from the decoupled OH-stretching mode in HOD of each water pocket using the PMM method. This approach has already proven effective in our previous work, where PMM-based spectra showed good agreement with experimental IR data for the same aqueous environments.¹¹ The results (Fig. S12 in the SI) are fully consistent with those obtained from the ζ analysis (see Fig. 3), showing a decrease in intensity at the peak frequency of the LDL-like bulk water signal (3300 cm⁻¹). This reduction reflects the diminished population of highly coordinated, tetrahedral water molecules under confinement.

Fig. 4 Intensity of the ζ distribution at ζ = 0.0945 nm, corresponding to the LDL peak in bulk water at 180 K, as a function of the radial coordinate measured inward from the ion/water interface, with the interface defined at zero. Data are shown for spherical shells of 0.1 nm thickness (each data point represents the midpoint of a shell). Reference intensities at ζ = 0.0945 nm for bulk water (dashed blue line) and fully mixed water (dashed dark red line) are also included.

All together, these results explain why, in the ζ distributions shown in Fig. 3, the high-ζ peak is shifted to lower values compared to that of LDL-like bulk water.

Conclusions

Our findings demonstrate that the presence of ions significantly alters the structural properties of supercooled water, as reflected in the distributions of the ζ parameter. In the fully mixed state with hydrazinium trifluoroacetate, the hydrogen-bond network is strongly disrupted, particularly due to the replacement of hydrogen-bonded water molecules with hydrazinium ions. In contrast, within confined water pockets, a less pronounced effect is observed, indicating that nanosegregation helps preserve local tetrahedral order. The effect of the ions is scale-dependent, with smaller pockets exhibiting stronger perturbations due to their higher surface-to-volume ratio and greater exposure to the ionic interface. The structural influence of the ions can be rationalized in terms of both direct and indirect effects. At the interface, the structure is directly perturbed by the presence of the ions, disrupting local hydrogen bonding. Beyond the interface, an indirect influence extends into the water phase, reaching approximately 0.7–0.75 nm from the surface. Beyond this range, particularly in the core of larger pockets, the LDL-like structure is largely preserved. Our findings are consistent with a recent study by Sciortino and co-workers,⁴⁵ which demonstrated that pure water droplets containing more than 1000 molecules, corresponding to radii of approximately 2.35 nm, as in our study, are sufficiently large to maintain cores exhibiting the structural characteristics of LDL-like water.

The calculated IR spectra further support these findings, revealing a clear correlation between ζ and the spectral balance between LDL-like and HDL-like water signatures. In particular, the attenuation of the 3300 cm⁻¹ absorption band associated with LDL-like water reflects the structural disruption captured by ζ.

Taken together, these results, based on our spherical water pocket model, illustrate how ionic solutes reshape the local structure of supercooled water and highlight the stabilizing role of sufficiently large nanosegregated domains in preserving hydrogen-bond order. This work offers new insight into the molecular-scale effects of hydrazinium trifluoroacetate in aqueous environments.

Author contributions

Nico Di Fonte: data curation, formal analysis, methodology, software, writing – original draft, writing – review & editing. Gianluca Dell’Orletta: formal analysis, methodology, writing – original draft. Isabella Daidone: conceptualization, methodology, writing – original draft, writing – review & editing. Laura Zanetti-Polzi: conceptualization, methodology, writing – original draft, writing – review & editing. All authors contributed equally to all aspects of the work.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI. Supplementary information: Initial structure files, and MD setup files. See DOI: https://doi.org/10.1039/d5cp02330k

Acknowledgements

The authors acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support (Project Code: IsCc3). This work was supported by the European Union – NextGenerationEU under the Italian Ministry of University and Research (MUR), PRIN2022-P2022MC742PNRR, CUP E53D23018440001.

References

1. P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos and J. Russo, et al. , Chem. Rev., 2016, 116, 7463–7500.
2. C. Huang, K. T. Wikfeldt, T. Tokushima, D. Nordlund, Y. Harada, U. Bergmann, M. Niebuhr, T. Weiss, Y. Horikawa and M. Leetmaa, et al. , Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15214–15218.
3. P. H. Poole, F. Sciortino, U. Essmann and H. E. Stanley, Nature, 1992, 360, 324–328.
4. J. C. Palmer, P. H. Poole, F. Sciortino and P. G. Debenedetti, Chem. Rev., 2018, 118, 9129–9151.
5. P. G. Debenedetti, F. Sciortino and G. H. Zerze, Science, 2020, 369, 289–292.
6. T. E. GartnerIII, L. Zhang, P. M. Piaggi, R. Car, A. Z. Panagiotopoulos and P. G. Debenedetti, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 26040–26046.
7. J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos and P. G. Debenedetti, Nature, 2014, 510, 385–388.
8. Z. Zhao and C. A. Angell, Angew. Chem., 2016, 55, 2474–2477.
9. S. Woutersen, B. Ensing, M. Hilbers, Z. Zhao and C. A. Angell, Science, 2018, 359, 1127–1131.
10. F. Caporaletti, D. Bonn and S. Woutersen, J. Phys. Chem. Lett., 2021, 12, 5951–5956.
11. L. Zanetti-Polzi, A. Amadei and I. Daidone, J. Chem. Phys., 2021, 155, 104502.
12. J. Bachler, I. Daidone, L. Zanetti-Polzi and T. Loerting, Phys. Chem. Chem. Phys., 2024, 26, 9741–9753.
13. H. Abe, T. Takekiyo, M. Shigemi, Y. Yoshimura, S. Tsuge, T. Hanasaki, K. Ohishi, S. Takata and J.-I. Suzuki, J. Phys. Chem. Lett., 2014, 5, 1175–1180.
14. H. Abe, T. Yamada and K. Shibata, J. Mol. Liq., 2018, 264, 54–57.
15. H. Abe, T. Takekiyo, Y. Yoshimura and A. Shimizu, J. Mol. Liq., 2019, 290, 111216.
16. F. Pabst, J. Kraus, S. Kloth, E. Steinrücken, M. Kruteva, A. Radulescu, M. Vogel and T. Blochowicz, J. Mol. Liq., 2021, 155, 174501.
17. A. Semmeq, K. Anand, A. Carof, A. Bastida and F. Ingrosso, Molecules, 2025, 30, 1678.
18. E. Shiratani and M. Sasai, J. Chem. Phys., 1998, 108, 3264–3276.
19. M. J. Cuthbertson and P. H. Poole, Phys. Rev. Lett., 2011, 106, 115706.
20. J. Russo and H. Tanaka, Nat. Commun., 2014, 5, 3556.
21. J. R. Errington and P. G. Debenedetti, Nature, 2001, 409, 318–321.
22. R. Foffi and F. Sciortino, J. Phys. Chem. B, 2022, 127, 378–386.
23. F. Martelli, H.-Y. Ko, E. C. Ouz and R. Car, Phys. Rev. B, 2018, 97, 064105.
24. N. A. Loubet, A. R. Verde and G. A. Appignanesi, Eur. Phys. J. E, 2024, 47, 61.
25. N. A. Loubet, A. R. Verde and G. A. Appignanesi, J. Chem. Phys., 2024, 160, 144502.
26. C. Faccio, M. Benzi, L. Zanetti-Polzi and I. Daidone, J. Mol. Liq., 2022, 355, 118922.
27. C. Faccio, N. Di Fonte, I. Daidone and L. Zanetti-Polzi, J. Mol. Liq., 2023, 392, 123425.
28. N. Di Fonte, C. Faccio, L. Zanetti-Polzi and I. Daidone, J. Chem. Phys., 2024, 161, 034505.
29. R. Shi and H. Tanaka, J. Chem. Phys., 2018, 148, 124503.
30. E. Duboué-Dijon and D. Laage, J. Phys. Chem. B, 2015, 119, 8406–8418.
31. J. L. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
32. J. L. Abascal and C. Vega, J. Chem. Phys., 2010, 133, 234502.
33. T. Sumi and H. Sekino, RSC Adv., 2013, 3, 12743–12750.
34. R. S. Singh, J. W. Biddle, P. G. Debenedetti and M. A. Anisimov, J. Chem. Phys., 2016, 144, 144504.
35. J. W. Biddle, R. S. Singh, E. M. Sparano, F. Ricci, M. A. González, C. Valeriani, J. L. Abascal, P. G. Debenedetti, M. A. Anisimov and F. Caupin, J. Chem. Phys., 2017, 146, 034502.
36. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
37. B. Doherty, X. Zhong, S. Gathiaka, B. Li and O. Acevedo, J. Chem. Theory Comput., 2017, 13, 6131–6145.
38. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1, 19–25.
39. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
40. M. Parrinello and A. Rahman, Phys. Rev. Lett., 1980, 45, 1196.
41. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
42. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
43. L. Zanetti-Polzi, S. Del Galdo, I. Daidone, M. D'Abramo, V. Barone, M. Aschi and A. Amadei, Phys. Chem. Chem. Phys., 2018, 20, 24369–24378.
44. R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Domanski, D. L. Dotson, S. Buchoux and I. M. Kenneyet al., MDAnalysis: a Python package for the rapid analysis of molecular dynamics simulations, Los Alamos Natl. Lab., 2019.
45. S. M. Malek, F. Sciortino, P. H. Poole and I. Saika-Voivod, Phys. Rev. Lett., 2025, 134, 138001.

---