Electrofreezing of liquid water at ambient conditions

Water is routinely exposed to external electric fields. Whether, for example, at physiological conditions, in contact with biological systems, or at the interface of polar surfaces in countless technological settings, water responds to fields on the order of a few V Å−1 in a manner that is under intense investigation. Dating back to the 19th century, the possibility of solidifying water upon applying electric fields – a process known as electrofreezing – is an alluring promise that has canalized major efforts since, with uncertain outcomes. Here, we perform long (up to 500 ps per field strength) ab initio molecular dynamics simulations of water at ambient conditions under external electric fields. We show that fields of 0.10 − 0.15 V Å−1 induce electrofreezing to a ferroelectric amorphous phase which we term f-GW (ferroelectric glassy water). The transition occurs after ~ 150 ps for a field of 0.15 V Å−1 and after ~ 200 ps for a field of 0.10 V Å−1 and is signaled by a structural and dynamic arrest and the suppression of the fluctuations of the hydrogen bond network. Our work reports evidence of electrofreezing of bulk liquid water at ambient conditions and therefore impacts several fields, from fundamental chemical physics to biology and catalysis.


Introduction
With at least 20 known crystalline forms and counting, the baroque phase diagram of water is the most complex of any pure substance [1] and is continuously under construction.Two amorphous ices, a low-density amorphous (LDA) and a high-density amorphous (HDA) ice [2], encompass a large set of sub-classes [3]; a third, medium density amorphous ice has recently been proposed [4], while a plastic amorphous ice has been suggested to exist at high pressures [5].Water is also routinely exposed to external electric fields (EFs).The range of strengths 0.1 − 1 V/ Å is particularly relevant, as it represents the range continuously produced by molecular dipoles fluctuations [6] in aqueous solutions [7][8][9] and to which water is exposed in countless technological/industrial settings [10].Recent developments have shown that the reaction rates of common organic reactions can be increased by one to six orders of magnitude upon applying external EFs [11][12][13][14][15][16], hence paving the way to the adoption of EFs as efficient catalyzers.Comparable EFs, generated by charge separation, endow microdroplets with strong and surprising catalytic power [17][18][19][20].
Historically, the possibility of manipulating water kinetics via EFs was first proposed by Dufour in 1862 [21].As experimental techniques matured over the years, such an opportunity became more tangible: the role of EFs on the heterogeneous nucleation of ice in cirrus clouds was addressed in the 1960s [22], and several other investigations followed, starting a vivid scientific debate [23][24][25][26][27][28].Recently, Ehre et al. [29] have shown that the kinetics of electrofreezing of supercooled water on pyroelectric materials is highly heterogeneous, favoring the crystallization on positively charged surfaces.Early -and pioneering -computational investigations based on classical molecular dynamics simulations also joined forces.According to these studies, liquid water undergoes electrofreezing to crystalline ice when exposed to external static EFs in the order of ∼ 0.5 V/ Å [30,31], and the effects of oscillating EFs have also been investigated [32].On the other hand, ab initio molecular dynamics (AIMD), which account for chemical reactions, and experiments have more recently shown that ∼ 0.3 V/ Å represents a threshold above which water molecules undergo dissociation into oxonium (H 3 O) + and hydroxide (OH) − ions [33][34][35][36].Seemingly, below this threshold, thermal energy and the associated large intrinsic field fluctuations taking place at the molecular scale impede the ordering of the hydrogen bond network (HBN), a necessary step for crystallization to occur.This task can instead be achieved, according to classical simulations, upon tuning the working pressure to ∼ 5 kPa and imposing external EFs of ∼ 0.2 V/ Å [37].
The application of EFs to liquid water induces a fast response of water molecules which align their dipole parallel to the field direction.On the other hand, the intrinsic cooperativity of HBs acts as a competing force, slowing down the relaxation and, in turn, driving the sample out of equilibrium.Therefore, in order to follow the response of water, it is necessary to probe the system at (non-)overlapping time windows rather than averaging over the entire simulation, and this is the paradigm we have decided to adopt (we report, in the SI, a comparison between quantities of interest computed at disjoint time windows and averaged over entire trajectories).In this study, we perform long (∼ 250 ps) AIMD simulations and show that EFs in the order of 0.10 ≤EFs≤ 0.15 V/ Å induce a structural transition to a new ferroelectric glassy state that we will call f-GW (ferroelectric glassy water).This transition occurs after ∼ 150 ps and is signaled by the freezing of the translational degrees of freedom, the suppression of the fluctuations of the HBN, and the drop in the potential energy.Our work represents the first evidence of electrofreezing of liquid water occurring at ambient conditions.

Results
In Fig. 1 we report the infrared (IR) spectra for bulk water without field (violet line) and for increasingly higher applied fields (0.05 V/ Å, blue; 0.10 V/ Å, orange; 0.15 V/ Å, red) computed over the time window [200 − 250] ps of the respective trajectories.In the absence of applied fields, the position of the OH stretching band -located at 3220 cm −1 -and that of low-frequency libration mode -at 560 cm −1 -are in good agreement with the experimental data [38].Upon exposure of the water sample to external EFs, we observe a contraction of the frequency range ascribed to the vibrational Stark effect [39], as also reported in Ref. [40] and -on limited frequency domains -in Ref. [41].This contraction indicates that the field imposes novel selection rules on the molecular vibrations.The largest frequency shift is associated with the OH stretching band; the corresponding red-shift is in the order of ∼ 75 cm −1 each 0.05 V/ Å, up to an EF of 0.10 V/ Å.However, a milder further red-shift in the order of 35 cm −1 occurs at a field of 0.15 V/ Å.The red-shift of the OH stretching is generally associated with stronger hydrogen bonds (HBs) [42] and the development of more "ice-like" environments [43].The reduced magnitude of the relative red-shift upon increasing the field from 0.10 V/ Å to 0.15 V/ Å, as quantified by the difference in frequencies reported in the inset of Fig. 1, suggests that the effect of the applied field becomes less intense.Moving towards lower frequencies, the weak libration+bending combination mode band at 2200 cm −1 is commonly associated with the strength of the hydrogen bond network (HBN) [44].The presence of the external EFs induces an enhancement and a concurrent slight blue-shift of this band, further suggesting that the EF causes a strengthening of the HBN.A similar effect has been reported on the IR spectra of water undergoing supercooling [45] where the strengthening of the HBN is, instead, induced by the reduction of thermal energy.A stronger effect on the vibrational spectrum occurs at lower frequencies, the signature of librational modes.The application of an external EF induces a significant blue-shift and, at 0.10 V/ Å and 0.15 V/ Å, the development of a clear new band peaked at ∼ 1000 cm −1 .This band has been ascribed to the breaking of the isotropy of molecular rotations and the preferential alignment with the field direction [40], as shown in Fig. S1 of the SI.
The picture emerging from the inspection of the IR spectra, therefore, indicates that the EF affects the topology of the HBN in several ways: the red-shift of the OH stretching band occurring upon increasing the applied EF is indicative of a strengthening of the HBs, while the blue-shift of the libration+bending combination mode band and that of the librations at lower frequencies suggests some degree of ordering of the HBN.At the same time, the appearance of a new peak in the librational band indicates an alignment of the molecular dipoles along the field direction.
The strengthening of the HBs (or their stiffening, as shown in Fig. S7) and the alignment of the molecular dipoles along the field direction mirror an enhancement of spatial correlations that also persists over time.In order to test this hypothesis, we report, in Fig. 2, the G OO (r, t), the Van Hove correlation function computed between oxygen atoms only and in the same time window [201 − 250] ps on which we have computed the IR spectra shown in Fig. 1.In Fig. 2 a) we report G OO (r, t) in the absence of external fields.We can observe that weak spatial correlation in the region ∼ 2.8 Å and ∼ 4.5 Å, corresponding to the first and second shells of neighbours, rapidly wear off in timescales of ∼ 5 − 10 ps.The application of a field of 0.05 V/ Å, reported in panel b), induces an extension of spatial correlations over slightly longer timescales.Radically stronger responses are induced by more intense fields: a field of 0.10 V/ Å (panel c)) and a field of 0.15 V/ Å (panel d)) clearly strengthen spatial correlations between ∼ 2 − 3 Å and 4 − 5 Å and extend them to timescales above ∼ 35 ps.
By projecting the partial Van Hove correlation functions on the reduced domain constituted by the spatial distances only (i.e., by removing the temporal dependence), we obtain the oxygen-oxygen radial distribution functions g OO (r).In Fig. 3 a) we report the g OO (r) computed in the time window [201 − 250] ps.Without any applied field (violet), the g OO (r) is that of bulk liquid water with a first peak located at ∼ 2.8 Å and a second peak at ∼ 4.5 Å. Adding a small EF of 0.05 V/ Å (blue) we observe an increase in the intensity of both the first and the second peaks with a reduction of the population between the first and second peaks.Upon doubling the intensity of the field and reaching 0.10 V/ Å (orange) we observe an enhanced increase in the intensity of both the first and second peaks and a further depletion of water molecules populating the interstitial region.An additional increase in the field intensity to 0.15 V/ Å (red) does not show appreciable changes in the g OO (r) with respect to the previous case, suggesting that no further major structural changes occur in the sample.Fig. S4 of the SI reports the g OO (r) computed in consecutive time windows of 50 ps starting from the beginning of our simulations, hence providing a glimpse of the dynamical structural transformations.In agreement with the profiles of the Van Hove functions, it is possible to observe that the g OO (r) for 0.05 V/ Å converges to the same profile after 50 ps (Fig. S3-a), while convergence is achieved only after 150 ps for 0.10 V/ Å (S3-b) and 200 ps for 0.15 V/ Å (S3-c).We notice, at this point, that the g OO (r) at fields of 0.10 V/ Å and 0.15 V/ Å at convergence, i.e., after 200 ps of simulation, strikingly resemble the g OO (r) of supercooled water or that of low-density amorphous (LDA) ice [3].This comparison is, instead, less accurate if one does not take into account the out-of-equilibrium nature that drives the process, and computes the radial distribution functions over the entire trajectories, as reported in Fig. S10-a as well as in several previous works.In order to rule out the effect of the simulation box, we have performed longer simulations (up to ∼ 500 ps) for systems with 256 H 2 O molecules at densities of 0.92 g/cm 3 and 0.95 g/cm 3 .Our results, reported in Fig. S8 of the SI, show that the development of a glassy-like g OO (r) is independent of the system size and density.Considering the high computational cost of performing AIMD simulations, we can not produce an equilibrated supercooled sample or an LDA via realistic quenching rates to compare the relative radial distribution functions.Therefore, in order to understand whether our g OO (r)s belong to a glassy sample or to a supercooled sample, we look at dynamical properties, namely the diffusivity measured via the mean squared displacement (MSD).Our results are reported in Fig. 3 b).We stress here that, like for the g OO (r), the MSD are computed on the time window [201 − 250] ps.It is possible to appreciate how the slope of the MSD drastically drops as soon as we introduce an EF.In the presence of a weak field of 0.05 V/ Å (blue) the sample is still liquid, although the mobility is strongly reduced compared to the case without field (violet).Upon increasing the field to 0.10 V/ Å (orange) and to 0.15 V/ Å (red) the MSD profiles indicate that water's translational degrees of freedom are confined to molecular vibration and to the rattling within the cage of the local neighbourhood.Computing the MSD over wider time windows implies accounting for the contribution of water molecules still in the liquid phase, hence artificially increasing the slope of the MSD, as shown in Fig. S10-b.We posit that this might be one of the reasons why the f-GW phase has been overlooked in previous studies.
In Fig. 4 we report the profile of the potential energy computed performing single point calculations on 1000 configurations randomly chosen within the time window [201−250] ps.Panel a) reports the profile as a function of the chosen molecular configurations.Without any applied field (violet) the potential energy fluctuates around the dashed violet line.Upon introducing a field of 0.05 V/ Å (blue) we observe a decrease in potential energy for almost all configurations, with an average value (dashed blue line) sitting below the case of water without field.Stronger drops in potential energy occur in the presence of EFs of 0.10 V/ Å (orange) and of 0.15 V/ Å (red).In panel b) we report the average potential energy -relative to the zero-field case in kcal/molas a function of the field strength.The drop in potential energy is clearly visible and shows how EFs of 0.10 V/ Å and 0.15 V/ Å drag the system into lower potential energy basins [46,47].It is worth noticing that the reduction in potential energy occurs along with a reduction of ∼ 14% of the entropy, as reported in Ref. [48] for the same system and numerical setups.The amorphization of liquid water involves a sensitive change in the fluctuations and topology of the HBN [49], which can be quantitatively inspected via the ring statistics.Therefore, in order to confirm that the structural and dynamical changes induced by the EFs indeed prompt a rearrangement in the HBN, we compute P (n), the normalized probability of having a ring of length n ∈ [3,10] in time windows of 50 ps.In hexagonal/cubic ice at 0 K and without defects, the P (n) is centered at n = 6, indicating that only hexagons are present.In Fig. 5 we report P (n) for strengths of 0.05 V/ Å, 0.10 V/ Å, and 0.15 V/ Å.Each case is reported against the P (n) determined in the absence of the EF (cyan circles).In the case of 0.05 V/ Å during the first 50 ps (black circles, panel a)), we can observe that the topology of the HBN overlaps almost perfectly with that of liquid water.Upon increasing the simulation time, the topology of the HBN responds to the presence of the field by increasing the number of hexagonal and heptagonal rings while reducing the number of longer rings.Overall, the response of the HBN to the presence of a weak field resembles the transformation of the HBN topology upon cooling [49,50].Upon doubling the field intensity to 0.10 V/ Å (panel b)), the topology of the HBN drastically changes even within the first 50 ps of simulation.In particular, we observe an increase in hexagonal and heptagonal rings with a corresponding decrease in longer rings.At consecutive simulation time windows, we observe a further sharpening of the P (n) with a considerable increase of hexagonal rings and a depletion of octagonal and longer rings.The topology of the HBN within the last 50 ps of our simulation is remarkably similar to that of LDA (obtained from classical simulations [49]).A similar behaviour occurs when we apply a field of 0.15 V/ Å (panel c)): the HBN reacts to the presence of the field increasing the population of hexagonal and heptagonal rings while decreasing the population of longer rings.Upon increasing the simulation time, the topology of the HBN further increases the population of hexagonal rings while decreasing longer rings, including heptagonal rings.
The gradual rearrangement of the topology of the HBN described above occurs on slower timescales compared to the alignment of water's dipole moment (see Fig. S3) and clearly shows that, although single water molecules react very quickly to the presence of EFs, the overall network of bonds reorganizes itself into new steady configurations on longer times, as also partially reported in Ref. [51].Such time-dependence, key in our investigation, can be seen in the gradual build-up of four-coordinated water molecules shown in Fig. S9.This gradual build-up in time leads to an increase in four-coordinated molecules up to 15% from the early stages of the simulation.Such an increase in the percentage of four-coordinated environments also induces a gradual enhancement of the local order.We report, in Fig. S6, P (I), the probability distribution of the local structure index I estimated on consecutive windows of 50 ps.It is possible to appreciate the development of bimodality in the later stages of our simulations for fields of 0.10 V/ Å (middle panel) and 0.15 V/ Å (lower panel).The lower panel of Fig. S6 reports a comparison between P (I) computed in the time window [201 − 250] ps for 0.15 V/ Å and for LDA at T = 200 K obtained from classical molecular dynamics simulations.Despite the differences in simulation techniques, the local structure of liquid water under EF strongly resembles that of LDA.The information collected so far indicates that our samples gradually readjust to the presence of external EFs.The slow evolution in time involves (i) the gradual development of four-folded configurations interacting via stronger HBs, (ii) the congruent development of more ordered local environments, (iii) the slow reduction of translational and rotational degrees of freedom, (iv) a drop in the potential energy, and (v) the gradual rearrangement of the HBN topology towards configurations richer in hexagonal rings.Eventually, after exposing the samples of liquid water to a field of 0.15 V/ Å for ∼ 150 ps, we observe a complete freezing of translational degrees of freedom, hence suggesting that our sample might be glass.Although the definition of glassy water is precise (molecular relaxation time exceeding 100 s or the shear viscosity reaching to 1013 poise), our simulations are too short to access these quantities.On the other hand, it has been recently shown that the transition to glass upon quenching liquid water is clearly signaled by the damping in the fluctuations of the HBN topology [49], which we here evaluate and report in Fig. 6 for the three cases in presence of the EF and against the fluctuations computed in liquid water without EF (cyan circles).For all cases, we determine σ(n) in time windows of 50 ps.In the case of 0.05 V/ Å, we can observe that, with respect to the case in the absence of the field, the fluctuations are strongly damped but for hexagonal and pentagonal rings, which fluctuate in a comparable measure.Upon increasing the simulation time, the fluctuations of the HBN are reduced for all cases but for the hexagonal rings, which become increasingly enhanced with the simulation time.Considering that the sample is liquid (although with strongly reduced diffusion), we posit that the diffusion occurs via changes in the HBN mostly involving hexagonal rings.At 0.10 V/ Å and 0.15 V/ Å, we observe a drastic suppression of the fluctuations of the HBN, to values well below those of the liquid.Such marked reduction of the fluctuations is responsible for the suppression of long-range density fluctuations occurring in correspondence with the transition to glassy water [49,52], a characteristic that differentiates liquid water from glassy states [53].Therefore, our findings are strongly indicative of a transition to a glass.

Discussion
In this work, we have performed long ab initio simulations of bulk water at ambient conditions in the presence of applied external electric fields (EFs) in the range 0.05 ≤EFs≤ 0.15 V/ Å.We have inspected the out-of-equilibrium process at disjoint time windows and recorded the results on each window.In the presence of an EF of 0.05 V/ Å, the dipoles align along the direction parallel to the EF while the diffusivity becomes sluggish.Overall, the inspected quantities computed within the last ∼ 150 ps of simulation are stable in time, indicating that the system is genuinely a liquid.Upon increasing the EF to 0.10 V/ Å and to 0.15 V/ Å, we observe a transition to a new ferroelectric glass that we call f-GW (ferroelectric glassy water).The amorphization occurs after ∼ 150 − 200 ps and is signaled by the freezing of the translational degrees of freedom and a drop in the potential energy, indicating that the sample has reached a metastable basin on the potential energy landscape.The evolution in time of the radial distribution functions and of other structural descriptors report an enhancement of the first and the second shells of neighbours along with a drastic depletion of the entries populating the space between them, as expected in the low-density glassy water.Similarly, the hydrogen bond network (HBN) undergoes a progressive structural reorganization favoring hexagonal motifs and a corresponding suppression of its fluctuations, as expected in the low-density glassy water state [49].
Our work represents the first evidence of electrofreezing of liquid water at ambient conditions, a task that has been attempted since 1862 [21].The new f-GW phase can be unveiled only by accessing and isolating late portions of long AIMD simulations, and is a new tile in the complex phase diagram of water.Therefore, it enriches our understanding of the physics of this complex material.Nonetheless, the conditions explored in this work are ubiquitous in industrial and natural settings, fields that can potentially benefit from this work.For example, water is routinely exposed to natural EFs comparable to the ones explored in this work when at the interface with enzymes, proteins, and biological membranes, defining the biological functionality and stabilizing such complex structures.
We infer that an experimental validation of our finding and the realization of the f-GW phase might be relatively straightforward by exploiting modern experimental settings.Many laboratories are nowadays capable of quantifying the field strengths generated in the proximity of emitter tips [10,13,54] -such those established by STM and AFM apparatus -, which fall in the same range required to transition to the new f-GW.Nonetheless, we posit that lower fields may induce electrofreezing to f-GW on longer time scales, accessible to accurate interaction potentials such as, e.g., MB-Pol [55] or Neural Network potentials [56].

Numerical simulations
We performed ab initio molecular dynamics (AIMD) simulations using the software package CP2K [57], based on the Born-Oppenheimer approach.The external electric fields (EFs) are static, homogeneous and directional (i.e., along the z-axis).The implementation of external EFs in Density Functional Theory (DFT) codes can be achieved via the modern theory of polarization and Berry's phases [58][59][60].In particular, owing to the seminal work carried out by Umari and Pasquarello [61], nowadays AIMD simulations under the effect of static EFs with periodic boundary conditions are routinely performed.The reader who is interested in the implementation of EFs in atomistic simulations can refer to the following literature: Refs.[58,59,[61][62][63][64][65][66].The main simulation here presented consists of a liquid water sample containing 128 H 2 O molecules arranged in a cubic cell with side parameter a = 15.82Å, so as to reproduce a density of 0.97 g•cm −3 .Furthermore, additional simulations were executed on bigger cubic cells composed of 256 water molecules and having edges of 20.05 Å and 20.26 Å.In such a case, lower densities of 0.95 g•cm −3 and 0.92 g•cm −3 were simulated, respectively.To minimize undesirable surface effects, the structures were replicated in space by employing periodic boundary conditions.We applied static and homogeneous EFs of intensities equal to 0.05 V/ Å, 0.10 V/ Å, and 0.15 V/ Å from a zero-field condition in parallel simulation runs.The maximum field strength of 0.15 V/ Å was chosen to prevent water splitting known to occur at larger field intensities [33][34][35][36]67].In the zero-field case we performed dynamics of 50 ps whereas, for each other value of the field intensity, we ran dynamics of at least 250 ps.Besides, as for the simulations of the lower-density states only a single field intensity of 0.15 V/ Å was simulated -in addition to the fieldless cases -for time-scales of ∼ 500 ps (ρ = 0.95 g•cm −3 ) and ∼ 450 ps (ρ = 0.92 g•cm −3 ).This way, we accumulated a global simulation time approaching 2 ns, whilst a time-step of 0.5 fs has been chosen.
Wavefunctions of the atomic species have been expanded in the TZVP basis set with Goedecker-Teter-Hutter pseudopotentials using the GPW method [68].A planewave cutoff of 400 Ry has been imposed.Exchange and correlation (XC) effects were treated with the gradient-corrected Becke-Lee-Yang-Parr (BLYP) [69,70] density functional.Moreover, in order to take into account dispersion interactions, we employed the dispersion-corrected version of BLYP (i.e., BLYP+D3(BJ)) [71,72].The adoption of the BLYP+D3 functional has been dictated by the widespread evidence that such a functional, when dispersion corrections are taken into account, offers one of the best adherence with the experimental results among the standard GGA functionals [73].It is well-known, indeed, that neglecting dispersion corrections leads to a severely over-structured liquid (see, e.g., Ref. [74] and references therein).Moreover, a nominal temperature slightly higher than the standard one has been simulated in the main simulations to better reproduce the liquid structure (i.e., T = 350 K).Furthermore, the additional simulations at lower density regimes were executed at a lower (supercooling) temperature of T = 250 K (see the SI for the respective results).
Albeit the BLYP+D3 functional represents e reasonably good choice, computationally more expensive hybrid functionals, such as revPBE0, when simulated along with the quantum treatment of the nuclei performs excellently well for water, as demonstrated by Marsalek and Markland [75].However, since sufficiently large simulation boxes are necessary to track structural transitions, the inclusion of the nuclear quantum effects is beyond the scope of the present work.Moreover, IR absorption line shapes of liquid water (and ice) are overall reproduced remarkably well by standard AIMD simulations, which include by their nature the explicit quantum adiabatic response of the electrons [76].In addition, the adherence of the IR and of the Raman spectra evaluated by some of us [40] under zero-field conditions with recent experimental results [38,77] justifies a posteriori the classical treatment of the nuclei.As a consequence, the dynamics of ions was simulated classically within a constant number, volume, and temperature (NVT) ensemble, using the Verlet algorithm whereas the canonical sampling has been executed by employing a canonical-sampling-throughvelocity-rescaling thermostat [78] set with a time constant equal to 10 fs.IR spectra have been determined by means of the software TRAVIS [79] (see the SI for further information).

Network topology
In order to probe the topology of the hydrogen bond network (HBN), we employed ring statistics, a theoretical tool that has proven to be instrumental in investigating the network topology in numerically simulated network-forming materials.The ring statistics is only one of many graph-based techniques to investigate network topologies and, in the case of water, it has helped in understanding the connections between water anomalies and thermodynamic response functions [50,80] as well as the properties of glassy water [52].We construct rings by starting from a tagged water molecule and recursively traversing the HBN until the starting point is reached or the path exceeds the maximal ring size considered (10 water molecules in our case).The definition of hydrogen bond follows Ref. [81].We do not distinguish between the donor-acceptor character of the starting water molecule.Supplementary information.A Supporting Information (SI) file with additional analyses and results accompanies the current work.
water: The synergistic action of several factors.The Journal of chemical physics 150 (9), 094506 (2019) [81] Luzar, A., Chandler, D.: Hydrogen-bond kinetics in liquid water.Nature 379(6560), 55-57 (1996) 1 Additional results Infrared spectra shown in Fig. 1 of the main text have been determined by means of the software TRAVIS [1,2] from the centers of the Maximally Localised Wannier Functions (MLWFs) [3,4] calculated on the fly during the ab initio molecular dynamics (AIMD) simulations.Molecular dipoles from MLWFs centers can be determined as: where e is the electron charge, r i is the position vector of the MLWF center i, Z j is the atomic number of the nuclei j whilst R j is the position vector of this latter.This way, the IR spectra at the investigated field intensities were computed as the Fourier transform of the molecular dipole autocorrelation function along the last 50 ps of the respective simulation trajectories.
In order to track molecular reorientations under the field action, we compute the distributions of the angle θ formed between the instantaneous water molecular dipole vectors and the field direction (i.e., z-axis), Fig. S1.Interestingly, whilst the field is capable of reorienting a large fraction of water dipoles already at 0.05 V/ Å, the electrostatic potential gradient producing this field strength does not induce a net suppression of the translational degrees of freedom of the molecules, as shown in Fig. 3 of the main text.The enhancement of the water dipoles at increasingly high fields is also visible from the dipole distributions reported in Fig. S2-a, showing a progressive shift towards larger magnitudes and a slight narrowing of the distributions.Fig. S2-b reports the profile of the dipole moment with the field strength.It is possible to recognize a linear regime holding up to a strength of 0.10 V/ Å.Thus, the transition from the liquid to the f-GW phase is also marked by the breakdown of the linear response regime to external electric fields (EFs).Additionally to this analysis, it is worth monitoring the temporal dependence of the P(θ) distributions at disjoint time windows, a procedure allowing for disclosing the dynamical response of the sample.As displayed in Fig. S3, the field-induced reorientation of the molecular dipoles takes place on fast timescales and achieves saturation within the first 50 ps of dynamics at all field intensities, with the exception of the weakest field (Fig. S3-a), where nonetheless the convergence of the dipolar response is reached in less than 100 ps.
Fig. S4 reports the oxygen-oxygen radial distribution functions computed at consecutive, disjoint time windows of 50 ps.At 0.05 V/ Å, the g OO (r) converges to a steady profile after 50 ps (Fig. S4-a), while convergence is achieved only after 150 ps for 0.10 V/ Å (Fig. S4-b) and 200 ps for 0.15 V/ Å (Fig. S4-c).The g OO (r) computed within the last 50 − 100 ps for 0.10 V/ Å and 0.15 V/ Å resemble the g OO (r) of a low-density amorphous (LDA) (see main text).To shed some light on the dynamical reorganization of the water structure induced by the external field, we have evaluated the oxygen-oxygen radial distribution function differences between the last 50 ps and the first 50 ps time frames of each simulation, as reported in Fig. S4-d.Whereas at 0.05 V/ Å structural differences between the initial and the final time windows appear to be small -as also visible in Fig. S4-a -, a field of intensity equal to 0.10 V/ Å induces much larger global reorganizations towards more structured molecular correlations in the system (Fig. S4-d, yellow curve).On the other hand, the evidence that these differences are smaller in the sample exposed to a 0.15 V/ Å field (Fig. S4-d, red curve) has to be ascribed to a faster initial reorganization taking place since the first 50 ps of dynamics, whereas longer timescales (∼ 200 ps) are somehow needed for bringing to completion the structural transition in the simulated sample, as shown in Fig. S4-c.In Fig. S5 we report P (q), the distribution of the tetrahedral order parameter q defined as q = 1 − 3 8 where ψ jk is the angle formed between the oxygen atoms of the water molecule under consideration and its nearest neighbour oxygen atoms j and k.The tetrahedral order parameter q was originally proposed by Chau and Hardwick [5] and subsequently rescaled by Errington and Debenedetti [6] so that the average value of q varies from 0 for an ideal gas to 1 for a regular tetrahedron.From Fig. S5 it is possible to observe that the samples in the absence of a field and in the presence of a field of 0.05 V/ Å show a very similar tetrahedral character.A major change occurs at stronger fields signaling the transition to the more ordered f-GW phase.Similar conclusions can be drawn upon inspecting the local structure index (LSI) [7,8], an insightful order parameter that can be employed to characterize the LDL and HDL molecular environments and defined as the inhomogeneity on the distribution of radial distances where ∆ j+1,j = r j+1 − r j is the distance between particles within a cutoff distance of 3.7 Å from a reference molecule and ⟨∆⟩ is the average overall neighbours of a molecule within the given cutoff.The LSI, therefore, provides a convenient quantitative measure of the fluctuations in the distance distribution surrounding a given water molecule within a sphere defined by a radius of 3.7 Å.In doing so, the index I measures the extent to which a given water molecule is surrounded by well-defined first and second coordination shells.In Fig. S6, we report the LSI computed for the three EFs here inspected at time windows of 50 ps.It is possible to observe the development of hints of a bimodal distribution in the cases of 0.10 V/ Å and 0.15 V/ Å in correspondence with the transition to f-GW.This can be also appreciated from the lower panel of Fig. S6, reporting the LSI computed in the last time window [201−250] ps for 0.15V/ Å and for the LDA simulated via classical molecular dynamics at T = 200 K.The latter has been obtained upon quenching liquid water from T = 300 K to T = 200 K at a quenching rate of 1 K/ns, as reported in Refs.[9][10][11][12] Somewhat related to the local and global degree of order of the H-bond network, is its kinetics.In particular, we performed a structural analysis of the H-bond network and identified a H-bond through the following geometric conditions (that must be simultaneously fulfilled): two water molecules are considered as H-bonded if R (OO) ≤ 3.5 Å and ∠O-H• • •O≤ 30 • , where R OO is the instantaneous distance between the oxygen atoms.From this, we calculated the time autocorrelation function of H-bonds as: where the indices i and j run on all pairs of first-neighbour molecules which at t 0 were H-bonded, t 0 being the time at which the measurement process begins; s ij = 1 if the criterion for the presence of a H-bond is fulfilled, s ij = 0 otherwise.The results were averaged over hundreds of initial configurations.Fig. S7 shows the continuous (Fig. S7-a) and intermittent (Fig. S7-b) autocorrelation functions c(t) of the H-bonds for different field intensities.Within the intermittent definition of c(t), a given Hbond is allowed to cleave within timescales ≤ 5 fs to account for bond fluctuations.Thus, within this latter fast timescale, we always assign to s ij a value equal to 1 when considering the intermittent autocorrelation function (see eq. ( 4)).The application of a 0.05 V/ Å field induces only a relatively moderate -with respect to the zerofield case -slow down of the dynamics of the H-bond network recorded by means of the continuous c(t) function (Fig. S7-a).Instead, significantly more drastic effects are recorded upon applying fields of 0.10 and 0.15 V/ Å.Interestingly, the changes produced on the H-bond network kinetics by these field regimes qualitatively resemble those induced by a sizable (∼ 40 K) decrease of the temperature [13].This is also visible from the intermittent H-bond autocorrelation function displayed in Fig. S7b.Although the H-bond characteristic time recorded at zero field (violet curve) is extended by the application of a field strength of 0.05 V/ Å (blue curve), a visible decay of the intermolecular correlations within the timescales of our simulations is recorded at the latter regime.Instead, a field of 0.10 V/ Å (orange curve) and a field of 0.15 V/ Å (red curve) clearly strengthen the H-bond persistence over sizably longer timescales.These results are fully consistent with the picture emerging from the partial Van Hove correlation functions shown in the main text (Fig. 2).
In Fig. S8 we report the oxygen-oxygen radial distribution function computed with a larger simulation box of 256 water molecules, at densities of 0.92 g•cm −3 and 0.95 g•cm −3 and at a temperature of 250 K (panels (a) and (b), respectively) for a sample without EF and a sample with a field of 0.15 V/ Å.These simulations reach ∼ 500 ps.It is possible to observe the development of an f-GW-like g OO (r) at both densities, indicating that the transition to f-GW reported in our work is not an artifact of small simulation boxes and that takes place for different densities.In Fig. S9 we report d = 4, the percentage of four-folded water molecules at consecutive time windows.The blue stripe corresponds to the case of liquid water in the absence of EFs.In the presence of 0.05 V/ Å (red squares), d = 4 increases from ∼ 50% to ∼ 53% within the first 50 ps of the simulation, and keeps gradually increasing reaching a maximum of ∼ 56% in the last two time windows.Upon increasing the field to 0.10 V/ Å (green diamonds) we can observe that d = 4 computed within the first 50 ps is roughly the same as the case for 0.05 V/ Å computed on the same time window.On the other hand, the d = 4 linearly increases by ∼ 6% in the second and in the third time window.The d = 4 reaches then a plateau in correspondence with the last two time windows.Upon increasing the field strength to 0.15 V/ Å, we observe a sudden increase in the d = 4 to ∼ 56% within the first time window.Further increases occur at the later stages of the simulation as for the cases previously inspected with lower field strengths.It is worth noticing that the percentage of four-coordinated water molecules for 0.10 V/ Å and 0.15 V/ Å is almost indistinguishable towards the end of the simulation.
To further stress the relevance of the sampling at disjoint time windows, we report in Fig. S10 a series of structural and dynamical observables determined at the highest field intensity here explored (i.e., 0.15 V/ Å) and calculated over the whole 250-ps-long trajectory and on the last 50 ps of dynamics.The g OO (r) at 0.15 V/ Å determined over the whole trajectory exhibits smaller (higher) peaks (dips) with respect to the same quantity calculated over the last 50 ps of dynamics of the same trajectory, as displayed in Fig. S10-a.Interestingly, the importance of sampling at consecutive time frames is even more visible from dynamical rather than structural properties.In fact, to adequately evaluate the field effects on the translational degrees of freedom, the sampling at consecutive time windows here adopted appears to be necessary for the timescales affordable by ab initio simulations.As shown in Fig. S10-b, indeed, the mean squared displacement of the oxygen atoms determined over the whole trajectory at 0.15 V/ Å witnesses the mixing of different translational regimes, a circumstance leading to an underestimation of the EF-induced damping effect.This is not only true for translational but also for rotational degrees of freedom, which are intimately related to the dynamics of the H-bond network.By direct comparison of the continuous (Fig. S10-c) and intermittent (Fig. S10-d) H-bond autocorrelation functions determined over different timescales (see legends), diverse H-bond characteristic times emerge.All these findings prove that relevant information on the effects produced by the application of external fields on liquid water can be unveiled only by accessing and isolating late portions of long ab initio simulations.Exclusively by adopting this strategy it is possible to catch the electrofreezing effect induced by the field on the roto-translational degrees of freedom of water and, presumably, of other H-bonded systems.

Fig. 1
Fig. 1 Infrared (IR) absorption spectra of liquid water determined at zero field (violet line) and under different field intensities as detailed in the legend.Arrows are guides for the eye qualitatively following the field-induced modifications of the bands.In the inset, we report the vibrational Stark effect of the OH stretching band.Data are computed on the time window [201 − 250] ps.

Fig. 2
Fig. 2 Partial Van Hove correlation functions between the oxygen atoms (i.e., G OO (r, t)) as a function of the time and of the intermolecular distance in the absence of the field (a) and in presence of static electric fields with intensities equal to 0.05 (b), 0.10 (c), and 0.15 V/ Å (d).Data are computed on the time window [201 − 250] ps.

Fig. 3
Fig. 3 (a) Oxygen-oxygen radial distribution functions at different electric field strengths (see legend).Dashed arrows qualitatively depict field-induced modulation of the hydration shells.(b) Mean squared displacement (MSD) of the oxygen atoms at various field intensities (see legend).In the inset, a logarithmic plot of the self-diffusion coefficient of the oxygen atoms as a function of the field strength.Data are computed on the time window [201 − 250] ps.

Fig. 4
Fig. 4 Potential energy computed via single point calculations on 1000 configurations randomly chosen within the time window [201 − 250] ps of each simulation.(a) Profile of the potential energy in a.u.for water without field (violet), water in the presence of 0.05 V/ Å (blue), water in the presence of 0.10 V/ Å (orange), and water in the presence of 0.15 V/ Å (red).Dashed lines represent the average value.(b) Average potential energy relative to the zero-field case in kcal/mol as a function of the field strength.

Fig. 5
Fig. 5 Probability distribution P (n) of having a ring of length n ∈ [3, 10] computed at different time windows during our simulations.The upper panel refers to the applied field E = 0.05 V/ Å, the middle panel to the applied field E = 0.10 V/ Å, the lower panel to the applied field E = 0.15 V/ Å.The cyan circles refer to the zero-field case.The black squares refer to the first 50 ps, red diamonds to the time window 51 − 100 ps, blue upper triangles to the time window 101 − 150 ps, the left magenta triangles to the time window 151 − 200 ps, and the green lower triangles to the window 201 − 250 ps.The dashed arrows emphasize the change in P (n) at consecutive time windows.

Fig. 6
Fig. 6 Fluctuations σ(n) computed on the ring statistics at different time windows during our simulations.The upper panel refers to the applied field of E = 0.05 V/ Å, the middle panel to the applied field of V = 0.10 V/ Å, and the lower panel to the applied field of E = 0.15 V/ Å.The cyan circles refer to the case of no field.The black squares refer to the first 50 ps, red diamonds to the time window 51 − 100 ps, blue upper triangles to the time window 101 − 150 ps, the left magenta triangles to the time window 151 − 200 ps, and the green lower triangles to the window 201 − 250 ps.The dashed arrows emphasize the P (n) change at consecutive time windows.

Fig. S1
Fig. S1 Distributions of the θ angle formed between the instantaneous molecular dipole vectors during the last 50 ps of the respective trajectories and the EF direction for diverse field intensities applied along the z-axis.

Fig
Fig. S2 (a) Distributions of the magnitude of the water dipoles extracted from the last 50 ps of the respective simulations at different field intensities determined from the MLWFs centers.(b) Average water dipole and associated standard deviation extracted from the distributions in (a).It is noteworthy pointing out the interruption of the linear response regime for field strengths producing water electrofreezing.

Fig. S3
Fig. S3 Distributions of the θ angle at 0.05 V/ Å (a), 0.10 V/ Å (b), and 0.15 V/ Å (c) measured at consecutive time windows.The black solid curves refer to the first 50 ps, dashed red lines to the time window 51 − 100 ps, dotted blue curves to the time window 101 − 150 ps, the dashed-dotted magenta lines to the time window 151 − 200 ps, and the dashed-dotted-dotted green curves to the window 201 − 250 ps.

Fig. S4
Fig. S4 Oxygen-oxygen radial distribution functions determined for disjoint time frames of 50 ps each determined from the ab initio molecular dynamics simulations conducted at 0.05 V/ Å (a), 0.10 V/ Å (b), and 0.15 V/ Å (c) field strengths.The black solid curves refer to the first 50 ps, dashed red lines to the time window 51 − 100 ps, dotted blue curves to the time window 101 − 150 ps, the dashed-dotted magenta lines to the time window 151 − 200 ps, and the dashed-dotted-dotted green curves to the window 201 − 250 ps.In panel (d), the point-by-point difference of the oxygen-oxygen radial distribution functions determined in the final and the initial time frames (i.e., during the last 50 and first 50 ps, respectively) for diverse field intensities is shown.

Fig. S6
Fig. S6 Upper panel: P (I) computed for the case at 0.05 V/ Å, 0.10 V/ Å, and 0.15 V/ Å at disjoint time windows of 50 ps.Lower panel: comparison between the LSI computed for 0.15V/ Å and the LSI computed for LDA at T = 200 K from classical molecular dynamics.

Fig. S7
Fig. S7 Log-log continuous (a) and linear-scale intermittent (b) H-bond autocorrelation functions calculated during the last 50 ps of the respective simulations.

Fig. S8
Fig. S8 Oxygen-oxygen radial distribution functions in the absence of the field (violet solid lines) and in the presence of a field strength equal to 0.15 V/ Å from ab initio molecular dynamics simulations conducted on boxes containing 256 H 2 O molecules in the supercooled regime (T = 250 K), and for densities of 0.92 g•cm −3 (a) and 0.95 g•cm −3 (b).

Fig. S9
Fig.S9Percentage of four-coordinated water molecules computed on consecutive time windows and compared against the case of bulk water in the absence of EFs (blue stripe).Red circles refer to the case of 0.05 V/ Å, green diamonds to 0.10 V/ Å, blue triangles to 0.15 V/ Å.

Fig
Fig. S10 (a) Oxygen-oxygen radial distribution functions, (b) oxygen mean squared displacement, (c) log-log continuous and (d) intermittent autocorrelation functions of the H-bonds determined at an EF intensity of 0.15 V/ Å calculated over the whole trajectory (black lines) and during the last 50 ps of dynamics (red lines).