Molecular dynamics analysis of amorphous ice transformations: elastic properties and nucleation

Unlike conventional ﬁrst-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 shear and compression rates. The system size effects are considered as well. The comparison with the experimental data enriches our understanding of the transitions observed. Our modelling gives new information on the nucleation mechanisms during transformations between low-density and high-density amorphous ices and, thus, helps answering the question about ﬁrst-order nature of this transition.


Introduction
Water is central to many physical, biological and industrial applications and exhibits physical properties that are qualitatively different from those of most other liquids. Despite the huge amount of theoretical and experimental research aimed at studying the many specific properties of water, it continues to be one of the most mysterious pure substances in the Universe. Supercooled water, known for its polymorphism, attracts a lot of interest. Today the phase diagram of water boasts no less than 19 modifications 1 of crystalline ice and 3 amorphous forms (not counting their derivatives) 2 , the most popular of which are the low density amorphous ice (LDA) and the high density amorphous ice (HDA). Amorphous ice practically does not occur on Earth in nature, but LDA is the most widespread form of H 2 O in the Universe, as it is a part of interstellar dust and comet nuclei. LDA can be obtained from the gas phase 3 or from liquid water 4 . Actually, the discovery by Mishima et al. 5 of solid-state amorphization (SSA) Ih → HDA in 1984 launched the entire epoch of the studies of the nature of these non-equilibrium transformations. It is noteworthy that amorphous forms of water were obtained from stable phases of all three states of aggregation 6 . Under isothermal compression of LDA or of hexagonal ice Ih, a transition occurs, which leads to the formation of HDA, and LDA can then be restored from HDA under isothermal decompression.
An important question about the nature of amorphous ices is the hypothesis of the two-liquid model of water and the existence of a second critical point. One 11 of several [11][12][13] theoretical scenarios to explain water anomalies (such as a sharp increase in its isothermal compressibility, isobaric heat capacity, coefficient of the thermal expansion upon cooling below the equilibrium freezing point) postulates the existence of a first-order-like transition between two forms of liquid water that correspond to LDA and HDA below the glass transition line (Fig. 1) and ends at the liquid-liquid critical point (LLCP). It is assumed that LLCP is located in the area called the "no man's land" corresponding to the temperature range of ∼145-200 K, where only crystalline forms of water are observed in experiments due to very rapid crystallization. With all the abundance of experimental observations [14][15][16][17][18][19] that indirectly indicate the possibility of the existence of the liquid-liquid phase transition in supercooled water, not a single unambiguous experiment has proved this thesis yet. Experimental problems under these conditions are associated with the inevitable rapid crystallization of a metastable liquid state, although the recent study by Kim et al. 20 provides a convincing evidence.
Molecular modelling plays an important role in studying the behavior of water under such conditions 7, 21-23 . It allows heating at sufficiently high rates to avoid crystallization and at the same time to study the structure in detail. Moreover, some models of water demonstrate the liquid-liquid transition, others do not 24,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].  7 . 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 experimental work 8,9 are indicated by black dotted lines; for the TIP4P/Ice model, the lines are in red dotted lines 10 . The white arrows show the paths considered in our work. Various approaches are used to study structure of various amorphous ices: an analysis of partial radial distribution functions (RDFs) 30,31 , the structural factor 32 , the classification using neural networks 23 . In this work, we rely on the partial RDF O-O and on the distribution of angles between hydrogen-bonded oxygen atoms that reflects the order in the tetrahedral structure. It was shown 32 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. Our interest in studying the behavior of elastic properties in transformations between amorphous forms of ice is caused by the experimental works, where elastic moduli were measured by the ultrasonic method 8 . It was 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.

Simulation details
In this work, we consider classical molecular dynamics modelling of amorphous ice using the TIP4P/Ice water model 10 . 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 33 with the GPU acceleration for TIP4P models 34,35 . 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 36 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 NV T (the Nose-Hoover thermostat) and NPT ensembles (the Nose-Hoover barostat). The visualization and structure analysis are carried out in the Ovito software package 37 .

Elastic moduli
Suzuki et. al 38 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 2/11 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 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 38,39 , 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 40 .
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 NV T 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  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. 2 (a)), exponential convergence is observed with decreasing compression rate. Moreover, for lower velocities, the value is closer to the experimental one 8 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 .

3/11
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. 41 . They took a similar approach using several water models, TIP4P/Ice including.

Preparatory stage
Preparation of LDA We obtain LDA by fast isobaric cooling of liquid water at a rate of 10 K/ns. 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. It was also shown 42 that the structure of the obtained amorphous ice changes weakly with a variation in the cooling rate in the range of 0.1-1 K/ns, the density deviations are within statistical error. As a result, we get a system at T = 77 K and P = 0.01 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 30 , but the first peak is overestimated, which was already shown earlier 10 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 g/cm 3 to 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.01 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.0-3.5 Å 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, SSA of hexagonal ice occurs already at 1.1 GPa, and here the compression rate plays a certain role. 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 8 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 .

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 (a)) 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 8, 9 . a) b) In fact, in MD we have implemented two methods for obtaining HDA: from hexagonal ice and from LDA (point 4 in Fig. 4 (a)). Comparing the radial distribution functions of the resulting structures ( Fig. 4 (b)), we can conclude that the HDA structure is independent of the method of its formation for these two cases considered. By analogy with a 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. 5. 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

Heating-induced transformations
Before we start discussing isobaric heating of amorphous ices, it is necessary to understand what rate of heating makes it possible to obtain an adequate description at a not very high computational costs. To do this, we carry out several processes of HDA heating at constant pressure (0.05 GPa) at different rates in the range Q heat = 0.25-50 K/ns, and Fig. 6 (a) shows how the density behaves with increasing temperature. At a heating rate of less than 7 K/ns, the curves almost converge, and for further analysis we will consider heating at a rate Q heat = 7 K/ns. Thus, the temperature dependencies of density of LDA and HDA under isobaric heating at the same rate is shown in Fig. 6 (a). 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. S3 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 and the tetrahedral structure change can be seen in more detail in Supplementary Fig. S3 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 phases, 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 8 , which was attributed by authors precisely to the HDA-LDA transition. In our modelling, however, the specific volume increases gradually (Fig. 6 (b)). 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 8 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 temperatures above 225 K to practically zero values is consistent with the fact that both systems transform into the liquid phase.

Nucleation during LDA ↔ HDA transitions
We made an attempt to see nucleation, using the fact that in HDA phase for an oxygen atom, the probability of detecting a neighboring oxygen atom at a distance of 3.0-3.5 Å increases. We classify molecules with a reduced probability of finding a neighbor at a distance of 3.0-3.5 Å as an LDA-like state and molecules with an increased probability as HDA. First, in the process of the LDA → HDA transition under isothermal compression (T = 77 K), the number of molecules characteristic of HDA increases with increasing pressure; in addition, growing HDA clusters are formed, which is shown in Fig. 8 (a). The size of the critical nucleus can be estimated as 5 molecules. That is, the nucleation takes place.
Next, let us consider the inverse transformation HDA → LDA during isothermal decompression to negative pressures. By analogy with the previous case, the formation of growing LDA clusters is observed as well ( Fig. 8 (b)).
But during isobaric heating (P = 0.05 GPa) of HDA, tiny LDA clusters of 2-3 molecules are distributed homogeneously in the simulation box and their number gradually increases without the formation and growth of critical nuclei (see Supplementary  Fig. S4 online)). Therefore, in this case there is no sign of nucleation during HDA → LDA transition. 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 nucleation and growth events. The absence of nucleation 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 8 : 1. There is a satisfactory agreement with experimental results 30 for the RDFs of LDA and HDA. The extra height of the first peaks of RDFs in the TIP4P/Ice molecular dynamics tells that the O-O and H-O bonds are too rigid. This fact, however, is a well-known feature of classical water models that do not take into account nuclear quantum effects at low temperatures 43, 44 . 2. During the isothermal SSA of Ih into HDA, the amorphization pressure depends on the compression rate, but no significant size effect was observed. The MD results are varying (very high) compression rates show a reasonable convergence to the experimental results.
3. The non-equilibrium modelling of pressure-induced transformations between LDA and HDA reflects the main features observed experimentally: the jump in density at the LDA → HDA transition and the hysteresis. The reverse transition of the HDA → LDA during decompression is observed at negative pressures. The RDFs of HDA structures obtained from Ih and from LDA are indistinguishable. 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 8 . The results of the corresponding MD calculations of bulk and shear moduli of HDA are in a good overall agreement with the experimental data 8 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 45,46 ). 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 41 ). It has been shown 47, 48 that including nuclear quantum effects in MD modelling improves the prediction of the Ih bulk modulus.
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 nucleation of a new phase as soon as a high degree of metastability is achieved. There are clear pictures of nucleation and growth events in the system of 2880 molecules for both transitions. This fact is an important argument in favor of the first order nature of the LDA ↔ HDA transition.
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 and then to liquid water. Contrary to the experimental data, in MD modelling there is no rapid increase of specific volume. The absence of nucleation 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 equilibruim line. Moreover, the system approaches the region near the second critical point 7 . 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. 8 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 that 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 nucleation and growth mechanism at high levels of metastability that confirms the first order nature of the LDA ↔ HDA transition.

Data availability
The datasets obtained and analysed during the current study available from the corresponding author on reasonable request.