Molecular dynamics analysis of elastic properties and new phase formation during amorphous ices transformations

Unlike conventional first-order phase transitions, the kinetics of amorphous-amorphous transitions has been much less studied. The ultrasonic experiments on the transformations between low-density and high-density amorphous ice induced by pressure or heating provided the pressure and temperature dependencies of elastic moduli. In this article, we make an attempt to build a microscopic picture of these experimentally studied transformations using the molecular dynamics method with the TIP4P/Ice water model. We study carefully the dependence of the results of elastic constants calculations on the deformation rates. The system size effects are considered as well. The comparison with the experimental data enriches our understanding of the transitions observed. Our modeling gives new information about the formation mechanisms of new phase clusters during the transition between low-density and high-density amorphous ices. We analyse the applicability of the term “nucleation” for these processes.

Molecular modelling plays an important role in studying the behavior of water under such conditions 17-21 . It allows heating at sufficiently high rates to avoid crystallization and at the same time to study the structure in detail. Moreover, some water models demonstrate the liquid-liquid transition [22][23][24] , others do not 25 . The TIP4P class of water models provides a well balanced instrument that allows predictive modelling of water-based systems at reasonable computational expenses and, therefore, at relatively large length and time scales (see e.g. [26][27][28][29]. Various approaches are used to study structure of the amorphous ices: an analysis of partial radial distribution functions (RDFs) 31,32 , the vibrational properties 33 , the hydrogen bond network 21 , the dynamical order parameter 34 and classification using neural networks 35 . In this work, we rely on the partial RDFs O-O and on the distribution of angles between hydrogen-bonded oxygen atoms that reflects the order in the tetrahedral structure. It was shown 33 that the tetrahedral structure is distorted in the Ih → LDA → HDA → liquid water sequence.
One of the tasks of this work is to form a microscopic understanding of the close relationship between elastic properties and the nature of a dynamic disorder present in amorphous ices. Knowledge of elastic properties combined with direct density measurements in amorphous ices is very useful for understanding their transformation mechanisms. The three-stage scenario of the heating induced transformation HDA → LDA discovered by Gromnitskaya et al. 30 clearly demonstrates complex kinetics in amorphous networks, which indicates the complex nature of glassy water. It was also shown that elastic softening precedes amorphous-amorphous and crystal-amorphous transformations. But the theoretical interpretation of these results is still absent and has not been considered in any molecular modelling studies yet.

Methods
Simulation details. In this work, we consider classical molecular dynamics modelling of amorphous ice using the TIP4P/Ice water model 36 that, for example, together with the TIP4P/2005 variant have been recently used for LLCP calculations 17 and for supercooled and amorphous water studies 19,37,38 . Positions of particles are determined from the solution of the classical Newtonian equations of motion. The simulations are carried out for a system of N = 2880, N = 23040 and N = 77760 water molecules, in a cubic box with periodic boundary conditions using the LAMMPS software package 39 with the GPU acceleration for TIP4P models 40,41 . The cutoff radii for the Coulomb and Lennard-Jones interactions are r c,Coul = 10Å and r c,LJ = 12Å , respectively. The long-range electrostatic interactions are treated using the Particle-Particle-Particle-Mesh algorithm 42 with an accuracy of 10 −5 . The SHAKE algorithm controls the bond lengths and angles for rigid molecules. The integration time step is 2 fs. The length of molecular dynamics trajectories is from 1 ns to 100 ns. The calculations are performed in NVT (the Nose-Hoover thermostat) and NPT ensembles (the Nose-Hoover barostat). Visualization and analysis are carried out with Ovito 43 .
Elastic moduli. Suzuki et. al 44 note that it happens that for the same material, different values of elastic moduli are obtained, depending on the methods used to calculated them. In addition, it may turn out that the same method is in good agreement with experimental data for substances with a structure of one type, but for another structure it already poorly reproduces what is obtained in practice. This clearly demonstrates the limits of applicability of one method or another. Meanwhile, it is important how well the model used describes the The T-P diagram of amorphous and liquid water. The separation line between a low density liquid and a high density liquid ends at the liquid-liquid critical point (LLCP) 17 . This point is located in an area called "no man's land", marked in yellow, where only crystalline forms of ice are experimentally observed. Transition lines based on the experiment 30  www.nature.com/scientificreports/ properties of a real system, in comparison with alternative ones, especially for water, for which more than a hundred different models have been proposed, but each is used for specific purposes 26 . Unlike other works on molecular dynamics modelling 44,45 , which study elastic properties with special emphasis on structural phase transitions, where elastic moduli are calculated using the fluctuation formula of molecular dynamics, here we calculate the elastic moduli directly 46 . The definition of the isothermal bulk modulus of elasticity is given by: where P is the pressure, V is the volume, ρ is the density of the system. To find the bulk modulus, after equilibration in the NVT ensemble for at least 2 ns, it is necessary to compress the simulation box by changing the linear dimensions by ±1% (while the deformation is linear) in each of the three directions and find the slope coefficient P(ρ). The shear modulus is by definition: where τ xy is the shear stress, γ xy is the shear strain.
To calculate the shear modulus G, by analogy with the bulk modulus B, we carry out a shear strain of ±1 − 2% in one direction. G is found as the slope coefficient τ (γ ) . It is important that our system is isotropic with respect to the shear modulus. We have carried out three independent deformations in different directions and made sure that the result does not depend on the direction of the shear (see Supplementary Fig. S1 online). The values of the elastic moduli for each state (P, T) are obtained by averaging over 3-5 independent straining simulations.
In addition, the strain rate can be varied, but this affects the values of the elastic moduli. Fig. 2 shows the values of B and G for calculations with different rates of relative deformation under the same conditions. For the bulk modulus (Fig. 2a), exponential convergence is observed with decreasing compression rate. Moreover, for lower velocities, the value is closer to the experimental one 30 B exp = 11 GPa , where the elastic moduli are found by the ultrasonic method at frequencies of 5 MHz. Actually, the ultrasonic method involves the measurement of adiabatic characteristics, so that the slight difference in values is probably due to the fact that we consider the isothermal modulus. Thus, there is not much benefit in using very low compression rates (less than 10 9 s −1 ). The behavior of the strain rate dependence of the shear modulus is different. No convergence is observed in this case, and there is a significant discrepancy (about 30%) with the experimental value G exp = 5 GPa . Further, shear modulus calculations are carried out at a shear deformation rate of 5 · 10 7 s −1 .
We are able to calculate the bulk modulus of amorphous ice with a deviation of less than 10% from the experimental value, while for hexagonal ice we have a deviation of about 70%, this is consistent with the work Moreira et al. 47 . They took a similar approach using several water models TIP4P/Ice including.

Preparatory stage
Preparation of LDA. In general, water is a poor glass former, and in practice a cooling rate of more than 0.01 K/ns is required to avoid crystallization. We obtain LDA by fast isobaric cooling of liquid water at a rate of 10 K/ns. It was also shown 48 that the structure of the resulting amorphous ice changes slightly when the cooling www.nature.com/scientificreports/ rate changes by 1-2 orders of magnitude, and the density deviations are within the statistical error. As a result, we get a system at T = 77 K and P = 0.1 MPa . with a density of 0.95 g/cm 3 , whereas in the experiment under such conditions the density is 0.94 g/cm 3 . As can be seen in Fig. 3, the radial distribution function is in good qualitative agreement with the experimental one 31 , but the first peak is overestimated, which was already shown earlier 36 for the TIP4P/Ice liquid state.

Preparation of HDA.
To obtain HDA, hexagonal ice at a temperature of 77 K and a pressure of 0.1 MPa is isothermally compressed (NVT ensemble, compression is performed by decreasing the size of the computational cell linearly with time) from a density of 0.95-1.31 g/cm 3 . After equilibration the system pressure is 1.4 GPa. Upon further release of the pressure 0.01 MPa, the amorphous structure is retained, the system remains metastable. An amorphous form is actually obtained that differs from LDA, its density at T = 77 K and P = 0.1 MPa is 1.15 g/cm 3 . Moreover, RDF of HDA has noticeable differences (Fig. 3). The second HDA's peak has a broader shape and a slight splitting at 4 Å. Also, there is an increased probability of finding water molecules at an O-O distance of 3.1-3.7 Å from the central water molecule. We repeated the calculation for a larger system, but no significant size effect is observed (see Supplementary Fig. S2 online).
In the experiment 5 , SSA of hexagonal ice occurs already at 1.1 GPa, but the compression rate plays a certain role. In the simulation, we consider isothermal compression with different compression rates in the range of 0.1-5 GPa/ns. In Supplementary Fig. S2 online dependencies of P(ρ) for different compression rates are given. The amorphization pressure decreases with decreasing compression rate. At the same time, our lowest compression rate q = 0.1 GPa/ns (in terms of density, this corresponds to a compression rate of 0.05 g/cm 3 /ns ), which is several orders of magnitude higher than the experimental rate 30 10 −12 GPa/ns . It is problematic to mimic the compression rate in molecular dynamics as in the experiment, nevertheless, the amorphization pressure obtained by us in the MD model is extrapolated to the corresponding experimental value well .

Results
Pressure-induced transformations. In this section of our study, we carry out isothermal ( T = 77 K ) compression/decompression of amorphous ices. For the initial system, we take the LDA obtained earlier (point 1 in Fig. 4) and compress it at a rate of 0.1 GPa/ns to a pressure of 2 GPa (point 3). Then the pressure is gradually released, until negative values (point 6), and compression begins again. Thus, in our model, transformations between amorphous forms of ice occur in a cycle that resembles "hysteresis", as in experiments 30,49 .
In fact, in MD we have implemented two methods for obtaining HDA: from hexagonal ice and from LDA (point 4 in Fig. 4). Comparing the radial distribution functions of the resulting structures (see Supplementary  Fig. S3 online), we can conclude that the HDA structure is independent of the method of its formation for these two cases considered.
By analogy with a conventional phase transition, the transition from LDA to HDA implies a change in local structure and should lead to a change of elastic characteristics. The dependence of the bulk modulus of elasticity on pressure during HDA decompression is shown in Fig. 5a. We also calculate the bulk modulus of LDA at the beginning and at the end of the compression, but this method does not allow finding intermediate values, since at 0.5-1.5 GPa the system is unstable. Nevertheless, one can be convinced that upon compression of LDA, the elastic modulus increases to values similar to HDA. During HDA decompression, the result of our model is in good agreement with the experiment 30 for the pressure dependence of bulk modulus. For the shear modulus, we have a slightly larger quantitative difference (Fig. 5b).   Fig. 4.2a. And as it is not difficult to see, at T > 210 K the densities of the systems become almost the same. Moreover, at a temperature of 240 K, the structures of heated LDA and HDA are identical, which is indicated by the similarity of the RDF and the distribution of angles between oxygen atoms (see Supplementary  Fig. S5 online), that is, we can say that the transformation of LDA and HDA into liquid water has occurred.
A more interesting question is what happens to amorphous ices before the transition to the liquid water phase? Trying answering this question, we compare the structure of the systems at 77 K and 180 K. How the RDF www.nature.com/scientificreports/ and the tetrahedral structure change can be seen in more detail in Supplementary Fig. S5 online. Up to 180 K, an increase in the peak of the angular distribution is observed, from which it can be concluded that, before the liquid phase, the HDA passes into an LDA-like state. Since the RDF is an ensemble-averaged characteristic of the system, the intermediate states between 90 and 180 K are more likely to relate to a mixture of two states, where HDA first prevails, and then LDA. It should be mentioned that under the experimental conditions at HDA heating, a rather sharp jump of 15% in the density (or specific volume) occurs at a temperature of 130-140 K 30 , which was attributed by authors precisely to the HDA-LDA transition. In our modelling, however, the specific volume increases gradually (Fig. 4.2b). This fact is can be associated with both the high heating rate and the use of the barostat, which facilitates the processes of local structural rearrangements. Besides, the behavior of the temperature dependencies of specific volume for systems with 2880 and 77760 molecules is practically the same.
As for the temperature dependencies of the elastic moduli (Fig. 7), we obtain a fairly good quantitative agreement with experiment 30 for the bulk modulus; however, we do not observe qualitatively no-monotony behavior in the model, as well as crystallization into cubic ice at 160 K with which these peculiarities are associated. The calculation of the shear modulus does not give such a good quantitative agreement. Nevertheless, both characteristics reflect softening of the amorphous network of LDA and HDA. The relatively sharp decrease in the bulk modulus of LDA at 190 K is most likely precedes further transformation into liquid water. And the decrease in the shear modulus at a temperatures above 225 K to practically zero values is consistent with the fact that both systems transform into the liquid phase.   12,21 ). But then one should expect that a process analogous to nucleation takes place during LDA ↔ HDA transitions. Using the results of our MD calculations, we make an attempt to describe the mechanism of the formation of HDA clusters during the compression of LDA and the formation of LDA clusters the reverse process. The selection rules for new phase detection plays a crucial role in such an analysis. In this work we use quite a simple algorithm based on the difference of RDFs for LDA and HDA structures. We use the fact that in HDA phase for an oxygen atom, the probability of detecting a neighbouring oxygen atom at a distance of 3.1-3.7 Å increases (Fig. 3). Let us denote the number of neighbours in the spherical layer 3.1 < r < 3.7 Å as N descr . The algorithm for detection of HDA/LDA clusters works as follows: 1. N descr is calculated for each oxygen atom. 2. If N descr 11 , then the atom corresponds to the HDA structure. If N descr 5 , then the atom corresponds to the LDA structure. Molecules with intermediate values 6 N descr 10 are excluded from the analysis. We will refer to this selection rule as the more stringent criteria (the less stringent criteria will be introduced below).
3. An HDA/LDA-atom i is included in a cluster if there is an HDA/LDA-atom j and |r i − r j | < 3.0 Å. The upper part of Fig. 4 shows the results of the application of the proposed selection algotithm. First, in the process of the LDA → HDA transition under isothermal compression ( T = 77 K ), the number of particles characteristic of HDA increases with increasing pressure; in addition, growing HDA domains are formed, which is shown in Fig. 4 (see also Supplementary Fig. S6 online). At a pressure of about 0.6 GPa a noticeable growth of domains begins. And the size of the smallest such cluster that accompanies intensive growth can be visually estimated as 5-20 molecules (this corresponds to a spherical domain of 4-9 Å). The transition pressure can be taken as 0.74 GPa, at which cluster growth is the most intense (we mark the corresponding value on the T-P diagram (Fig. 1). By analogy with the previous case, in the inverse transformation HDA → LDA during isothermal decompression formation of growing LDA clusters is observed as well. This supports the experiment of Tonauer et al. 50 , in which the HDA → LDA transformation during decompression at 140 K was considered and the "nucleation" phenomenon was revealed. The radius of a (spherical) LDA seed was estimated as 3-8 Å 50 .
The intermediate values 6 N descr 10 introduce a certain ambiguity to the selection algorithm considered. Following Fig. 3, molecules in both the LDA-like and the HDA-like local environments can have N descr in this range. In order to address the sensitivity of the algorithm to the choice of the threshold values, we consider below the case of N descr 10 for HDA and N descr 6 for LDA molecules. On the one hand, the less stringent selection rule ( N descr 10 for HDA and N descr 6 for LDA molecules) means that more particles are involved in the analysis of cluster formation. On the other hand, the more stringent criteria ( N descr 11 for HDA and N descr 5 for LDA molecules) makes it possible to assign a molecule to one or another structure with a higher probability. As can be seen from the distributions in Fig. 3, LDA has a small fraction (about 3%) of molecules with N descr = 10 , and HDA has a small fraction (about 5%) molecules with N descr = 6.
The following differences between the less stringent and the more stringent criteria can be mentioned: 1. In the less stringent case the "unassigned" molecules with 7 N descr 9 form small clusters, but in the more stringent case the molecules with 6 N descr 10 form a common background (see Supplementary Fig. S7  online). www.nature.com/scientificreports/ 2. In the less stringent case the number of molecules in the new phase clusters is higher (see Fig. 8). This effect is more pronounced for LDA clusters in HDA under decompression. For example, the snapshots on Fig. 8 show how the number of molecules of HDA clusters increases when the selection rule is changed from the more stringent variant to the less stringent variant. Clusters selected via N descr 10 are bigger and some big clusters are results of merging several small clusters selected via N descr 11 . However, the increase in volume of each cluster is less pronounced that the increase in its number of molecules. 3. In the more stringent case the new phase clusters vanish to 1-4 molecules close to the hypothetical LDA-HDA equilibrium where metastability disappears. In the less stringent case the clusters of ∼ 10 HDA-like molecules under LDA compression and ∼ 20 LDA-like molecules for HDA decompression are observed. This observation suggests that the more stringent selection rule suits better for the detection of the new phase formation (that is why we use for the kinetic spinodal shown on Fig. 1).
During isobaric heating ( P = 0.05 GPa ) of HDA, the number of LDA clusters gradually increases without the evident growth of domains (see the variation of the selected LDA-cluster in Supplementary Fig. S8 online). Therefore, in this case, we see fluctuations-based growth of LDA state, and there are no arguments to speak about a cluster growth (nucleation-like) process. In order to consider possible system size effects, we perform the same modelling for the system of 77760 molecules, but the nature of the transition remains qualitatively the same without any visible cluster growth features. The absence of the formation of clusters that trigger growth events during this isobaric heating can be explained by the nearly constant level of metastability along the isobar that goes parallel to the LDA-HDA equilibrium line.

Discussion
In this work using the MD modelling with the TIP4P/ice potential we have confirmed that the model gives a consistent description of LDA and HDA amorphous ices and have reproduced the key experimental results on LDA ↔ HDA transitions 30 : 1. There is a satisfactory agreement with experimental results 31   The lower snapshot shows the oxygen atoms with N descr 11 (the largest cluster has 19 molecules), seven HDA-clusters are highlighted by the color and the surface mesh (a number next to a cluster indicates how many molecules it contains). The top snapshot for the same time moment shows the oxygen atoms with N descr 10 (the largest cluster has 33 molecules) and the HDA-clusters that includes the clusters shown on the lower snapshot (additional atoms are highlighted in dark color). www.nature.com/scientificreports/ 4. We have analyzed the compression rate and share rate dependencies in the MD calculations of elastic moduli. We show that the bulk modulus converges quickly at decreasing compression rate. However, no convergence has been detected for the shear modulus at decreasing shear rate. Therefore, the calculations have been performed with the shear rate close to the corresponding ultrasonic frequency in the experiment 30 . The results of the corresponding MD calculations of bulk and shear moduli of HDA are in a good overall agreement with the experimental data 30 for the case of isothermal compression and for the case of isobaric heating. There is a very good quantitative agreement for the bulk modulus (the accurate prediction of the liquid water bulk modulus at room temperatures in TIP4P/2005 model has been reported recently as well 23,53 ). However, the MD results for the shear modulus are about 50% lower that the experimental data. Here it is worth mentioning that the agreement for the TIP4P/Ice result for the bulk modulus of Ih is about 70% higher that the experimental value (similarly as it was shown previously 47 ). It has been shown 54,55 that including nuclear quantum effects in MD modelling improves the prediction of the Ih bulk modulus.

5.1
In the framework of the presented model, we have developed a nearest neighbours analysis method that distinguishes clusters of HDA in LDA and vice versa. Using this method, we have shown that both the LDA → HDA and the HDA → LDA transformations at isothermal compression/decompression proceed via the formation and growth of new phase clusters as soon as a high degree of metastability is achieved. There are clear pictures of domains formation and growth events in the system of 2880 molecules for both transitions. Fig. 4 shows that at the pressure of thermodynamic LDA-HDA equilibrium there are clusters of the competing structure formed by 1-4 molecules. These clusters can be regarded as local "subcritical" fluctuations. In the region of higher metastability the size of such cluster increases. The mechanism of their growth and thermodynamic analysis requires a separate study. Here we want to emphasize that the new amorphous structure (LDA in HDA and vice versa) form as localized domains. We can call a kinetic spinodal the narrow range of pressures when the size of a larges cluster grows rapidly at compression of LDA (or decompression of HDA). These pressures obtained in our calculations for four temperatures are shown on Fig. 1 as a boundary of kinetic stability of the respective amorphous ices (and supercooled water at temperatures close to LLCP). The location of such a kinetic spinodal should depend on the compression/decompression rate. We can notice that the shape of a non-equilibrium kinetic spinodal obtained in this study are in a reasonable agreement with the LLCP parameters determined recently for TIP4P/Ice model 17 .
5.2 Both the more stringent and the less stringent variants of the N descr threshold values used for LDA/HDA detection show that the transition between LDA and HDA proceeds via formation of localised clusters of new phase. In the amorphous solids as we consider and at low temperatures any diffusion of molecules is practically not possible. We observe two mechanisms of growth of these new phase clusters: the attachment of individual molecules and the merging of initially separated clusters. Hence, on the basis of our observations we can say that the formation of new phase in the LDA ↔ HDA transformations has evident aspects of the nucleation process typical to first-order transitions.
5.3 The most problematic case for the MD modelling considered is the isobaric heating of HDA. In this, case we observe in MD modelling a gradual transition of HDA to LDA-like state and then to liquid water. Contrary to the experimental data, in MD modelling there is no rapid increase of specific volume. The absence of the formation of clusters that trigger growth events can be explained by the nearly constant level of metastability during this isobaric heating process: the system moves on the T-P diagram in parallel to the LDA-HDA equilibrium line. Moreover, the system approaches the region near the second critical point 17 . It is known that the size of critical nucleus becomes larger when the system approaches its critical point. Therefore, we can make a hypothesis that the proper explanation of the experimental data on the isobaric heating requires a theory for HDA-LDA-liquid water transformations at larger time and length scales well beyond the MD scale.

Conclusions
The TIP4P/Ice molecular dynamics modelling of amorphous ices and LDA ↔ HDA transformations gives consistent results that have been compared with the ultrasonic measurements of elastic constants in HDA, LDA and during the LDA → HDA isothermal transformation and during the HDA → LDA isobaric transformation by Gromnitskaya et al. 30 MD modelling gives the correct qualitative description of the pressure and temperature dependencies of bulk and shear moduli of HDA and LDA in these processes. The TIP4P/Ice water model provides a quantitatively accurate description for the HDA bulk modulus, but for the HDA shear modulus the absolute values are significantly lower than the experimental values.
Molecular dynamics modelling captures the essential features of the LDA ↔ HDA transition observed experimentally. Therefore, the model is able to supplement the experiment and allows studying in details the reverse process of HDA → LDA isothermal transformation that takes place at negative pressures. Using this model, we have shown that both types of isothermal LDA ↔ HDA transformations proceed via the mechanism of formation of new phase clusters at high levels of metastability.

Data availability
We provide a link to the most significant data that illustrate the observed "nucleation" events during isothermal LDA ↔ HDA transitions https:// github. com/ garku lic/ Nucle ation-during-the-LDA-HDA-trans ition. The rest of the data requires a deep dive into the analysis process. We are ready to provide other information to the interested readers on a reasonable request.