Article Cite This: Langmuir 2017, 33, 11561-11569
pubs.acs.org/Langmuir
Phase Diagram of TIP4P/2005 Water at High Pressure Masanori Hirata,† Takuma Yagasaki,‡ Masakazu Matsumoto,‡ and Hideki Tanaka*,‡ †
Graduate School of Natural Science and Technology and ‡Research Institute for Interdisciplinary Science, Okayama University, Okayama 700-8530, Japan S Supporting Information *
ABSTRACT: We report a new ice phase that forms spontaneously at the interface between ice VII and liquid water in molecular dynamics simulations of TIP4P/2005 water. The new phase is structurally quite similar to an ice phase originally found to be a precursor in the course of the homogeneous nucleation of ice VII in liquid water. A part of the water molecules in these ice phases can rotate easily because the number of hydrogen bonds in them is less than four, and thus they can be regarded as partial plastic phases. A rough estimate suggests that these phases are thermodynamically more stable than either ice VI or ice VII for 3 GPa < P < 18 GPa at T = 300 K. Although the partial plastic phases would be metastable states at any point in the phase diagram of real water, they might be realized experimentally with the aid of dopants and/or solid surfaces.
■
INTRODUCTION There exist many ice polymorphs. Currently, 17 ice phases have been reported experimentally.1−5 They are distinguished by the network topology of the hydrogen bonds. Low-pressure ice I has a porous structure consisting only of six-membered rings of hydrogen-bonded water molecules with a tetrahedral angle of 109.5°. Other ring structures appear in high-density ice polymorphs: ice II contains flat six-membered rings in addition to the chair-type six-membered rings, ice III includes five- and seven-membered rings, and ice VI consists of two interpenetrating lattices that contain highly distorted fourmembered rings. Ice VII, which is the highest-density molecular crystal, consists of two interpenetrating lattices of cubic ice. Classical molecular dynamics (MD) and Monte Carlo simulations can be used to find a new ice phase. A theoretical study by Kosyakov et al.6 as well as simulation studies of Jacobson et al.7 and Conde et al.8 suggested that guest-free structure II clathrate hydrate is a stable ice phase at negative pressures. Although this ice is metastable under ambient pressure, it was realized experimentally by Falenty et al.4 several years later. They prepared this ice, named ice XVI, by 5 days of vacuum pumping on neon hydrate at 142 K. The thermodynamic properties of ice XVI were investigate by using the quasiharmonic method.9 Huang et al. recently suggested that there exist more stable ice phases than ice XVI under deeply negative pressures using the combination of DFT calculations and classical molecular simulations.10,11 Classical molecular simulations have also suggested the existence of several new ice structures at high pressures. Takii et al. discovered that there is a plastic phase in the phase diagrams of TIP5P, TIP4P, and SPC/E water.12 In the plastic ice phase, each water molecule rotates freely while its center-of-mass © 2017 American Chemical Society
position is bound to the lattice site of ice VII. Aragones et al. reported a different plastic ice phase using the TIP4P, SPC/E, and TIP4P/2005 water models.13,14 The oxygen atoms are located in a BCC arrangement in the plastic phase found by Takii et al., and they are located in an FCC arrangement in the plastic phase of Aragones et al. Thus, these plastic phases are called BCC and FCC plastic phases, respectively. Himoto et al. found the emergence of critical phenomena including the tricriticality for the phase transition between the BCC plastic phase and ice VII.15,16 Mochizuki et al. observed a metastable ice phase having the same unit cell as methane A17 in the course of the crystallization of ice VII.18 The space group of the unit cell is R3, so we call it ice R in this study. In the unit cell of ice R, one molecule is coordinated by 12 water molecules at almost equivalent distances and the water molecule orients randomly while the others do not: this ice is a partial plastic phase. So far, the existence of the plastic ice phases at high pressure has not been affirmed experimentally. This might suggest whether the plastic ice phases are stable or metastable states in real systems. However, metastable states can be realized by, for example, dopants and temperature control. Indeed, neutron diffraction studies have shown that the hydrogen bond network of ice VII is strongly disturbed by lithium halide, and the orientation of water molecules in doped ice VII is similar to that in the BCC plastic ice phase. Moreover, Special Issue: Tribute to Keith Gubbins, Pioneer in the Theory of Liquids Received: May 26, 2017 Revised: August 2, 2017 Published: August 10, 2017 11561
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir
Figure 1. Snapshots of the MD simulation of the liquid/ice VII coexistence system at T = 300 K and P = 6 GPa in the sequence of the elapsed time. Local structures are identified by the pattern-matching technique. Orange dots indicate the water molecules whose local configuration matches the 2 × 2 × 2 extended unit cell of ice VII including 16 molecules. Other water molecules are not drawn for clarity. Boxes indicate the local configuration that matches the unit cell of ice T. Unit cells of different orientations are shown by boxes with different colors. For blue boxes, the c axis of the unit cell is parallel to the z axis (vertical direction in this picture) of the simulation cell. Ice VI-like structures also appear sporadically in the liquid region, but they do not reach the critical size. See also the video in the Supporting Information. (a) t = 0 ps. The initial configuration. (b) t = 340 ps. Ice T nucleates on ice VII and quickly covers the interface between the liquid and ice VII. (c) t = 1140 ps. The liquid region is replaced by distorted ice T. (d) t = 4740 ps. Ice T erodes ice VII. (e) t = 6240 ps. The simulation cell is filled with a single crystal of ice T.
distinct from the partial plastic ice phase found by Mochizuki et al., ice R.18 The thermodynamic stabilities of ices T, R, VI, and VII are roughly evaluated with several assumptions.
NMR studies suggested the presence of a plastic state of water molecules in mesoporous silica.19,20 It should be noted that not only stable phases but also metastable phases are worth scrutiny because a metastable phase may have different properties from the stable phase under a given condition, and there are industrially useful metastable phases such as the martensite phase of carbon steels, the diamond phase of carbon, and amorphous phases of silica. There are many classical force fields for water. The TIP4P/ 2005 model is one of the most reliable ones.21,22 This model was originally developed on the basis of the TIP4P model to reproduce various properties of liquid water under ambient conditions, the densities of ice II and ice V, and the thermal stability of ice III with Ewald summation. It was demonstrated that the TIP4P/2005 model well reproduces the phase diagram of ice polymorphs in a very wide range of pressure and temperature.13 This model has been widely used to examine liquid water,23−28 electrolyte solutions,29,30 ice polymorphs,8,10 and clathrate hydrates.31−34 There might be other stable ice phases of TIP4P/2005 water that are candidates for metastable phases of real water. The order−disorder phase transition of ice VII to the BCC plastic phase occurs easily in MD simulations because the oxygen atoms need not move. It is, however, generally difficult to observe the structural phase transition from a solid phase to a different solid phase without rare event sampling methods such as transition path sampling because molecules are tightly bound in solid phases and the homogeneous nucleation rate is quite low.35,36 However, nucleation from a metastable solid state to a more stable solid one is expected to be strongly enhanced under the coexistence of a metastable liquid phase because such a process occurs promptly at the interface with the liquid phase. With this facile nucleation due to the intervention of the liquid phase, such a transition is observed without any sophisticated treatment. In this study, we perform MD simulations of the coexistence of liquid water and various ice phases using the TIP4P/2005 model and find that an unidentified crystalline phase emerges at the ice VII/liquid interface at P ≈ 6 GPa and T ≈ 300 K. The obtained crystal, we call ice T, is similar to but
■
COMPUTATIONAL DETAILS We compare four types of ice structures: ice VII, ice VI, ice R, and ice T. All of the ice structures are hydrogen-disordered and do not have net polarization. We prepare the target ice crystals in a rectangular cell with a ratio of x/y/z = 1:1:2. The system consists of 131 072 TIP4P/2005 water molecules for ice VII, 81 920 for ice VI, 72 576 for ice R, and 144 000 for ice T. To make the ice/liquid two-phase coexisting state, an MD simulation at constant volume and temperature (NVT) is performed for 100 ps at T = 1000 K with fixed atomic positions of half of the molecules in the simulation cell. The very high temperature of 1000 K is chosen to ensure complete melting of the unfixed part, and the volume is unchanged to avoid vaporization. Then, the system is annealed to the target temperature by a short (typically 300 ps) NVT MD. The ice/ liquid interfaces move slightly during the annealing because of the growth or melting of the ice part. The resultant configuration, in which almost half of the system is liquid water while the rest remains the solid phase, is used as the initial configuration of the production runs of MD simulations at constant temperature and pressure (NPT). It is possible to prepare a solid and a liquid separately and combine them into a single cell.37−39 However, this protocol is less efficient because the solid/liquid interface prepared by this method is usually quite unstable and an additional simulation is required to relax. The protocol used in this study is better and has been frequently employed to prepare liquid−solid coexistence systems in recent simulation studies.40,41 The temperature and pressure are kept constant by the Nosé−Hoover method42,43 and the Berendsen method,44 respectively. Long-range Coulomb interactions are treated with the particle mesh Ewald (PME) method. 45 The GROMACS 4.646,47 package is employed for all of the MD simulations. 11562
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir
of the crystal. Then, we examine the symmetry operations for the unit cell of the crystal to determine the space group. It is found that the unit cell is tetragonal (space group I41) and contains 72 water molecules. We refer to this tetragonal structure as ice T. Table 1 lists the positions of the oxygen
We also calculate the properties of the four types of ice structures in the single-phase systems. The system size is the same as that of the corresponding ice/liquid coexisting system.
■
RESULTS AND DISCUSSION Emergence of the New Ice Structure. We perform a MD simulation of the liquid water/ice VII two-phase coexistence at T = 300 K and P = 6 GPa. In the phase diagram of TIP4P/2005 presented by Aragones et al.,14 ice VII is the most stable phase under this condition. Thus, it is expected that ice VII grows as time evolves. However, we find that a different solid phase emerges at the liquid/solid interface as shown in Figure 1. The new solid phase grows rapidly and eventually occupies the whole system. The solid structure is observed in a certain range around the above thermodynamic condition. A video of the trajectory is given in the Supporting Information. We identify the structure of the new solid phase. Experimentally, the crystal structure of an unknown phase is determined on the basis of the diffraction pattern. In contrast, real-space configurations are obtained from MD simulations. It is possible to convert the real-space configuration to the diffraction pattern. However, this process leads to an error due to the incomplete crystallization. Therefore, we employ a realspace pattern-matching technique to determine the unit cell. We evaluate the configurational similarity of the neighborhoods of a reference oxygen atom a and a given oxygen atom b as follows. First, we define the R neighborhood of an atom a, BR(a), as a set of atom positions whose distance from a is shorter than R. A sphere of radius R = 8 Å typically contains ≈100 oxygen atoms. We regard this point cloud, BR(a), as a template. Then, the template is slid to the position of a targeted oxygen atom, b, in the simulation cell. For each point i in the slid template, there is a corresponding atom i′ that is closest to point i. (Figure 2) The configurational similarity is defined as
Table 1. Oxygen Atom Fractional Coordinates of the Ice T Structure at 6 GPaa
a
Ns
∑ (Δr ia − Δr bi′)2
(1)
i=1
y
z
0.000 0.000 0.334 0.334 0.166 0.166 0.251 0.173 0.751 0.827
0.000 0.000 0.928 0.032 0.928 0.032 0.142 0.173 0.358 0.327
0.000 0.500 0.026 0.526 0.721 0.221 0.872 0.499 0.872 0.249
The space group is I41 with a = 13.28 Å and c = 7.38 Å.
Figure 3. Schematics of the ice T structure. (a) A unit cell of ice T is shown in an orthogonal projection. The unit cell consists of four square cylinders. The molecules on the vertical walls of the cylinders are drawn with blue circles. Each cylinder contains a double helix (pink and yellow circles). The positions of the blue circles are very close to the corresponding lattice points of ice VII. Black lines are hydrogen bonds. The helix is drawn in only one of the four square cylinders. The hydrogen bonds between the walls and the helices are omitted. Some molecules are removed, and molecular positions are somewhat modified to clarify the structure and symmetry. (b) Projection onto the ab plane. A unit cell of ice T (colored) is overlaid on the BCC lattice of ice VII (open circles). Each circle indicates the position of an oxygen atom. Helices wind in opposite directions alternately. (c) Similarity between the unit cells of ice R and ice T. The c-axis lengths are close to each other. The unit cell of ice R is rotated by 90/4° along the c axis in order to match the molecular positions in the unit cells. The unit cell of ice R is almost cubic, and its space group is R3. The dotted line indicates the 3-fold symmetry axis, ⟨111⟩.
the mean square displacements of the point in the slid template from that of the corresponding atom 1 Ns
x
O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
atoms in the unit cell of ice T. Figure 3a illustrates the structure of ice T. The unit cell consists of four square cylinders, and
Figure 2. Schematic for the method to determine the primitive translational vectors. The R neighborhood of atom a (filled circles) is slid to target atom b, and the displacements (red lines) of the atoms from the corresponding points in the template are evaluated.
Da(b) =
atom
Δrai
where Ns is the number of points in BR(a), is the position of point i relative to point a in the template, and Δrbi′ is the position of atom i′ relative to target atom b. Quantity Da(b) becomes a small value when atoms a and b are located at a translationally equivalent position in the unit cell of a crystal. This procedure is repeated for many combinations of a and b for a single snapshot to assess the primitive translational vectors 11563
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir each cylinder contains a double helix of water molecules. There is no direct hydrogen bond between the two strands. Contrary to ice VII and ice VI, ice T does not have the interpenetrating double hydrogen bond network. Each water molecule in ice VII makes an undistorted tetrahedral configuration with four neighboring hydrogen-bonded molecules. In contrast, in ice T, the hydrogen bond network is strongly distorted, especially for the water molecules of the helices. Once the unit cell of ice T is determined, we can make use of it to detect the crystal structure that spontaneously forms and grows in the system by pattern matching. The computational cost of the pattern matching is huge because the orientations of the nucleated crystals are not known in advance and the system size is very large. We employ the three-dimensional geometric hashing technique to speed up the matching.48 This method is also used for the detection of other types of ice. Figure 1 demonstrates that ice T is well distinguished from ice VII and liquid water by the pattern-matching technique. There is no crystalline motif in the liquid region at t = 0 ps (Figure 1a) because the liquid part is equilibrated at the very high temperature of 1000 K in the preparation process of the initial configuration. Ice T nucleates at the liquid/solid interface at t = 340 ps and grows in the liquid phase as shown in Figure 1b,c. The whole liquid region becomes ice T at t = 1.1 ns. Then, the conversion from ice VII to ice T takes place. As shown in Figure 3b, the c axis of ice T is parallel to the z axis of the simulation cell, and the angle between the a axis and the x axis is 45°. The lattice constants of ice T are 3(21/2)L, 3(21/2)L, and 2L, where L is the lattice constant of ice VII: ice T is a 18 × 18 × 2 reconstruction of ice VII. Because of the structural similarity between the two ice phases, ice VII transforms easily to ice T. The water molecules in the unit cell are classified into four categories and shown by four different colors in Figure 3. The positions of the water molecules on the cylinder walls, shown by the dark-blue (O1 and O2 in Table 1) and light-blue (O3, O4, O5, and O6) circles, are very close to those in ice VII. The water molecules of a helix shown by the yellow circles (O8 and O10) are also only slightly displaced from the original location in ice VII. Water molecules of the other helix shown by the pink circles (O7 and O9) sequentially move to the vicinal space. The transition from ice VII to ice T can occur rapidly with the zipper action along the helices. Ice T is even more similar to ice R, which appears as an intermediate state in the simulations of the homogeneous nucleation of ice VII in the liquid phase using the TIP4P/2005 model.18 Figure 3c demonstrates that the positions of water molecules in ice T are almost the same as those in ice R when the two structures are properly overlapped. In this figure, the structure of ice R is rotated by 90/4° around the c axis. Both the unit cells of ices T and R contain the BCC packing part which is a characteristic structural motif of ice VII. The unit cell of ice V, which consists of 28 water molecules, is larger than that of any other experimentally reported ice polymorph at positive pressures. The unit cell of ice T is about 3 times larger than that of ice V. Although the unit cell of ice XVI consisting of 136 water molecules is even larger, this ice phase cannot form in liquid water because it is thermodynamically stable only at negative pressures.4,5 We perform single-phase simulations of liquid water and ices VI, VII, T, and R at 6 GPa and 300 K. Figure 4 presents the oxygen−oxygen radial distribution functions (RDF). The RDF of ice VI exhibits two peaks at 2.67 and 3.25 Å in the shortrange region. The former arises from the oxygen atoms in the
Figure 4. Oxygen−oxygen RDFs at P = 6 GPa, T = 300 K. The green, orange, blue, red, and black curves are the RDFs of ice VI, ice VII, ice T, ice R, and liquid water, respectively.
same sublattice, and the latter is due to the oxygen atoms in different sublattices; the former and latter correspond to the hydrogen-bonded pairs and unbound but close-packed pairs, respectively. The first peak of ice T is also split into two components. The minor peak at 2.77 Å comes from the strongly hydrogen-bonded pairs. The major peak at 2.97 Å corresponds to the unbound but close-packed pairs or the stretched and weakly hydrogen-bonded pairs. The strong hydrogen bonds are shown by the black solid lines in Figure 3a. The water molecules of the helix form the weak hydrogen bonds with the water molecules on the cylinder wall. The RDF of ice R is quite similar to that of ice T. The RDF of liquid water also closely resembles that of ice T, although the component at 2.97 Å is less prominent. In the RDFs of ice R, ice T, and liquid water, the first peak accompanies a subpeak or a shoulder, but their separation is less clear than that for ice VI. In ice VII, the oxygen atoms are located in the BCC arrangement. The peak at 2.83 Å and the shoulder at 3.25 Å correspond to the center−corner and corner−corner distances in the cubic unit cell of the BCC lattice. Ice VII is the densest structure among them, and indeed the unbound neighboring molecules get even closer; the sharp peak at 2.83 Å without any shoulder indicates that the distances of the hydrogen-bonded neighbors (four out of the eight center−corner distances) and those of the unbound neighbors (the other four center−corner distances) are equivalent. We calculate the rotational correlation function of water molecules defined as R (t ) =
1 Nw
Nw
∑ ⟨ui(t )·ui(0)⟩ i
(2)
where Nw is the number of water molecules and ui(t) is the unit vector parallel to the molecular dipole vector of the ith water molecule at time t. Figure 5a compares the rotational correlation functions of ice VI, VII, R, and T at P = 6 GPa and T = 300 K. The rotational relaxation is faster in ice VII than in ice VI because a rather facile exchange of the hydrogen bond partner is expected due to the equivalent hydrogen bond distance for all of the pairs and the fact that the O−O distance for hydrogen-bonding pairs is longer for ice VII. The rotational relaxation is even faster for ice T due to the distorted structure that gives rise to the peak at 2.97 Å in the RDF. Figure 5b presents R(t) for the four types of water molecules in ice T classified according to Figure 3. Function R(t) decays very quickly for the water molecules of the double helix (pink and yellow) because their hydrogen bonds are distorted. This result 11564
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir
potential energy of the system is given by half of the averaged interaction energy between a water molecule and all other water molecules.) We accumulate the attractive and repulsive water−water interaction energies separately. The attractive interaction energy mainly arises from the water molecules forming a hydrogen bond with the central water molecule. Table 2b shows that the attractive energy is higher for partial plastic ice T and ice R than for nonplastic ice VI and ice VII. The repulsive interaction energy of ice VII is much higher than that of the other three ice phases as a result of the presence of the four non-hydrogen-bonding neighboring molecules, and this results in the highest potential energy shown in Table 2a. The repulsive interaction energy is lower for less-dense ices because the distance between the central molecule and water molecules that do not form hydrogen bonds with the central molecule can be longer in less-dense structures. Phase Diagram. The dominant phase at a given temperature and pressure is determined by the Gibbs energy, G = Ek + Ep + PV − TS, where Ek and Ep are the kinetic and potential energies, respectively. Table 2a suggests that ice VI is thermodynamically stable in the low-pressure region because of the low potential energy, whereas ice VII is stable in the high-pressure region because of the high density that results in a low PV value. Ice R and ice T are expected to be stable in the intermediate region of pressure. There are several ways to find the most stable phases under given thermodynamic conditions and draw the phase diagram. One is to directly calculate the Gibbs energies of various phases by using, for example, the thermodynamic integration technique.49 The Gibbs energy of ice VI and ice VII can be calculated with an assumption that the residual entropy is given by the Pauling approximation.50 The Gibbs energies of the BCC and FCC plastic ice phases also can be calculated because there is no residual entropy.12,14 The Gibbs energies of ice R and ice T contain the residual entropies that are different from the Pauling value as a result of the partial plasticity. It is unclear how to evaluate the residual entropy of the partial plastic ice phases. Therefore, we take a different approach in investigating the phase diagram. We perform MD simulations of the ice/ liquid coexistence for various ice phases to evaluate the stability of the ice phase relative to the liquid. Sometimes a different ice phase forms in the two-phase coexisting system. This means that the newly formed phase is more stable than either the preexisting ice or liquid phases. The coexistence simulations provide a lot of information that can be used to estimate the location of phase boundaries on the P−T plane. Table 3a summarizes the most dominant phase at t = 10 ns in the MD simulations starting from the VII/liquid two-phase coexistence at various temperatures and pressures. In most cases, one phase prevails in the whole system within several nanoseconds. It is found that the crystalline phase dissociates in the low-pressure region. Ice VII grows in the low-temperature, high-pressure region. The solid phase also grows in the hightemperature and -pressure region; however, the solid becomes either BCC or FCC plastic ice. Because the FCC arrangement is denser than that of BCC, the former is more stable at higher pressures. Ice T spontaneously forms and grows at intermediate pressures of P ≈ 5 GPa in the low-temperature region of T ≤ 350 K. There is no coexistence line between ice VII and liquid water in the phase diagram of TIP4P/2005 water for 300 K ≤ T ≤ 500 K due to the presence of ice T and the plastic BCC phase.
Figure 5. (a) Rotational correlation functions of ice VI (green), VII (orange), T (blue), and R (red) at P = 6 GPa, T = 300 K. (b) Rotational correlation functions of the four types of water molecules in ice T. The curves are colored according to the classification in Figure 3. The blue and cyan curves are R(t) for the water molecules on the wall, and the yellow and magenta curves are for the water molecules of the double helix.
demonstrates that ice T is a partial plastic phase. Figure 5a shows that R(t) of ice R is very similar to that of ice T. This is because ice R is also a partial plastic phase.18 Table 2a lists the potential energies and the densities of the four ice structures at P = 6 GPa and T = 300 K. The two partial Table 2. Potential Energies and Densities of Ice Structures and Interaction Energies between a Water Molecule and Surrounding Molecules (a) Potential Energies and Densities of Ice VI, R, T, and VII at P = 6 GPa and T = 300 K PE/kJ mol−1
ρ/g cm−3
ice VI −50.11 1.550 ice R −46.75 1.637 ice T −46.52 1.641 ice VII −45.08 1.678 (b) Attractive, Repulsive, and Total Interaction Energies between a Water Molecule and Surrounding Molecules Located within 4 Å of the Central Molecule for Ice VI, R, T, and VII at P = 6 GPa and T = 300 K ice ice ice ice
VI R T VII
attractive/kJ mol−1
repulsive/kJ mol−1
total/kJ mol−1
−107.38 −101.32 −101.57 −108.01
27.66 34.70 35.67 44.65
−79.72 −66.62 −65.90 −63.36
plastic ice phases, ice R and ice T, are quite similar for both the potential energy and the density. The potential energies of the partial plastic ice phases are higher than that of less-dense ice VI and lower than that of denser ice VII. Table 2b shows the interaction energies between a water molecule and the surrounding molecules located within 4 Å from the central molecule. (r ≈ 4 Å is a common minimum position of the RDFs of the four ice structures.) The total short-range interaction energies of the partial plastic ice phases also lie in between those of ice VI and ice VII, indicating that the difference in the energetic stability can be explained without considering long-range interactions. (Note that the absolute values are larger in Table 2b than in Table 2a because the 11565
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir
5 GPa for TIP4P/2005 water reported by Aragones et al.14 We observe that ice R spontaneously forms at T = 300 K and P = 6 GPa. Figure 6 shows snapshots of the simulation, and the video is supplied in the Supporting Information. Because ice R and ice VI do not have common molecular configurations in their unit cells, ice R does not nucleate at the interface. Rather, nuclei of ice R randomly form in the bulk liquid region (Figure 6b). As a result, the liquid region is filled with crystal grains of different orientations (Figure 6c). Then, the region of ice R grows into that of ice VI (Figure 6d), and ice VI disappears eventually (Figure 6e). Table 3a shows that ice T is the final state in the VII/liquid coexistence system under the same thermodynamic conditions, T = 300 K and P = 6 GPa. We cannot determine which is more stable at this point for ice R and ice T from the liquid/solid coexistence simulations. At intermediate pressures, P = 4 and 5 GPa, the final state is ice T in Table 3a and ice VI in Table 3c. Again, we cannot determine the most stable phase for these thermodynamic conditions. The formation of ice R at T = 300 K and P = 6 GPa and the structural similarity between ice R and ice T imply that the thermal stability of ice R is comparable to that of other ice phases discussed in this study. We perform MD simulations of the ice R/liquid coexistence system for 10 ns. The final states are summarized in Table 3d. Comparison between parts b and d of Table 3 indicates that the melting temperature of ice R is higher than that of ice T. This means that ice R occupies a certain region in the phase diagram of TIP4P/2005 water as a stable phase, although this ice was first assumed to be a metastable state that emerges in the nucleation simulations of ice VII in liquid water.18 It would be possible to compare the thermodynamic stabilities of two solid phases by using MD simulations of the coexistence of the two solid phases. However, it is difficult to prepare the initial configuration because the structure of the solid/solid interface is unknown. This problem can be circumvented by the insertion of a liquid layer between the two solids, but it requires a larger simulation cell. In addition, there are many choices in combinations of two solids. Although simulations of all of the combinations are not necessarily required at each temperature and pressure, their computational cost is quite heavy. Unlike solid/solid coexistence conditions, liquid/solid coexistence conditions can be easily determined by two-phase coexistence simulations. Table 3c shows that the coexistence pressure of ice VI and liquid water is P ≈ 2.5 GPa at T = 300 K. Table 3d indicates that the coexistence pressure of ice R and liquid water is also P ≈ 2.5 GPa at T = 300 K. That is, this point is the triple point for ice VI, ice R, and liquid water. The Gibbs energies of the three phases are the same at the triple point. The entropy difference between ice VI and ice R at this point is given by
Table 3. Final States of the Solid/Liquid Coexistence Simulationsa (a) Starting from Ice VII/Liquid Coexistence P/GPa 10 8 7 6 5 4 3
300 K
400 K
450 K
500 K
7 B B B L L
F B B L
F B B L
400 K
450 K
500 K
7 7 7 T T T T T T L L (b) Starting from Ice T/Liquid Coexistence
P/GPa 10 8 7 6 5 4 3
300 K
350 K
T T T T L
T T L
T T L
T T L L
T L (c) Starting from Ice VI/Liquid Coexistence
P/GPa 7 6 5 4 3 2 P/GPa 10 9 8 7 6 5 4 3 2
350 K
300 K
350 K
400 K
R L R L 6 L 6 L 6 L L (d) Starting from Ice R/Liquid Coexistence 300 K
350 K
400 K
L L
450 K
500 K
R R L
R L
R R R R R R R R L
R R R R L
R R R L
a
Blank cells indicate that the system does not settle into a single phase within the simulation time or the simulation was not carried out. 7, 6, F, B, R, T, and L indicate ice VII, ice VI, FCC plastic ice, BCC plastic ice, partial plastic ice R, partial plastic ice T, and liquid water, respectively. Letters in bold font indicate the spontaneous formation of a solid phase that is different from the initial phases.
We perform MD simulations for 10 ns starting from the ice T/liquid water coexistence. The final states are shown in Table 3b. Although ice T nucleates at the VII/liquid interface, ice VII does not form at the T/liquid interface. This indicates that the Gibbs energy of ice VII is higher than that of ice T or liquid water. Although ice VI should be the most stable phase in the lowpressure and -temperature region, this phase does not form in the simulations of the VII/liquid and T/liquid coexistence systems. We therefore perform MD simulations of the VI/ liquid coexistence for 40 ns to assess the stability of ice VI. (The simulation time is longer than those for the VII/liquid and T/liquid systems because the growth rate of ice VI is relatively slow.) The final states are shown in Table 3c. Our simulations roughly reproduce the melting line of ice VI below
ΔS = SR − S VI =
HR − HVI T
(3)
The enthalpies of the two ice phases are calculated from the single-phase simulations. The obtained ΔS value is 0.89kB. The entropy of ice R is larger than that of ice VI because of the partial plasticity of ice R. We assume that the entropy of ice T is the same as that of ice R. We also assume that the entropy of ice VI is the same as that of ice VII at a given pressure and temperature. Figure 7 plots the Gibbs energies of ice VI, T, and VII relative to that of ice R calculated with these assumptions against pressure. The Gibbs energy curve of low-pressure ice VI 11566
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir
Figure 6. Nucleation and growth of ice R in the liquid/ice VI coexistence system at 6 GPa and 300 K as time elapses. Green dots indicate the water molecules whose local configuration matches the unit cell of ice VI. Other water molecules are omitted. Boxes indicate the local configuration that matches the unit cell of ice R. Unit cells of different orientations are shown by boxes with different colors. See also the video in the Supporting Information. (a) t = 0 ns. The initial configuration. (b) t = 2.5 ns. The first crystal nucleus appears in the liquid region. (c) t = 9 ns. The liquid region is filled with crystallites of ice R. (d) t = 13 ns. Ice R erodes ice VI. (e) t = 28 ns. The largest crystallite grows even larger, but it is still polycrystalline.
stable than ice T, probably because the ice R structure is less commensurate with ice VII than ice T as shown in Figure 3c. In the high-temperature region, T = 400, 450, and 500 K, the transition from ice VII to either the BCC or FCC plastic phase occurs while the formation of ice T is not observed (Table 3a). This indicates that the full plastic phases are more stable than the partial plastic ice T phase in this region. The comparison between parts a and d of Table 3 shows that ice R is the most stable at (T, P) = (400 K, 5 GPa) and (450 K, 6 GPa). Thus, the order of the most stable ice phases with increasing pressure is ice R, BCC, and FCC in the high-temperature region. Ice R is stable at low pressures because of the low potential energy due to the relatively strong hydrogen bonding, whereas the BCC and FCC phases are stable at higher pressures because of the higher densities. Because the entropy of the partial plastic ice R phase is smaller than that of the full plastic BCC phase, ice R is expected to disappear at a temperature higher than 500 K. Two partial plastic phases, ice T and ice R, have not been reported experimentally, although they are stable ice phases in MD simulations of TIP4P/2005 water under certain conditions. This indicates that the TIP4P/2005 model overestimates the thermodynamic stability of these ice phases. It is difficult to elucidate what is missing or overemphasized in the force field model among many possibilities. For example, the Lennard-Jones function (r−12) is used in the TIP4P/2005 model, but other functions such as the Buckingham potential (exponential function) may be better to describe the shortrange repulsion between oxygen atoms at very high pressures. Multibody interactions due to the polarization and chargetransfer effects are also not explicitly included in the TIP4P/ 2005 model. One may think that quantum mechanical simulations can be used to examine the high-pressure ice phases. If the ice phase is stable at low temperatures including T = 0 K, then it is possible to analyze the properties of the ice phase based on only the optimized structure. Because the effect of thermal motions is essential to the entropic stability of the high-pressure plastic ice phases, the analysis based on the optimized structure is of no use and long MD simulations are required. Current computational power is still insufficient to perform such long simulations with quantum mechanical force fields. The TIP4P/2005 model will likely continues to play a major role in studies on ice polymorphs because its computational cost is relatively low as a result of simplicity of the model and other classical water models is less reliable at high pressures.14
Figure 7. Relative Gibbs energies of ice VI (green), T (blue), and VII (orange) with respect to ice R (red) at T = 300 K.
intersects that of high-pressure ice VII at P = 6 GPa. This value is in agreement with the VI/VII coexistence curve in the phase diagram of TIP4P/2005 water reported by Aragones et al. and shown in Figure 8 in spite of the crude assumptions.14 The
Figure 8. Pressure ranges where ice T (blue) or ice R (red) is the most stable phase at 300 K are overlaid on the phase diagram of TIP4P/ 2005 water reported in ref 14 (gray solid lines and labels). Colored labels are the stable phases suggested by the liquid/solid coexistence simulations. See Table 3 for the labels.
Gibbs energy of ice R is lower than that of the other ice phases for 3 < P < 7 GPa, and ice T is the most stable phase for 7 < P < 17.5 GPa. As expected, the order of the most stable ice phases against pressure, VI, R, T, and VII, is consistent with the order of the density (and thus PV in the Gibbs energy) shown in Table 2a. The spontaneous formation of ice T on the interface of ice VII found at P = 4, 5, and 6 GPa is due to the thermodynamic driving force of GVII − GT ≈ 2.5 kJ mol−1. Ice R does not form at these pressures although ice R is more 11567
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Langmuir
■
As shown in Figure 8, ice T and ice R occupy wide regions in the phase diagram of TIP4P/2005 water. It is possible that the partial plastic phases are indeed stable phases for real water but the regions in the phase diagram are quite narrow. If so, the partial plastic phases may be detected experimentally by careful scanning of temperature and pressure near the VI/VII coexistence curve. It is more feasible that ice T and ice R are metastable states in reality. Even in this case, these ice phases might be observed experimentally by using dopants similarly to the orientationally disordered BCC structure found in LiCldoped ice VII.51 The interfacial free energy between the partial plastic phases and ice VII is expected to be low because the partial plastic phases are structurally quite similar to ice VII as shown in Figure 3. Because the entropy is higher for the partial plastic phases than for ice VII, the interfacial free energy between the liquid phase and the partial plastic phases would be lower than that between the liquid phase and ice VII. The relation between these interfacial free energies suggests that there might be a layer of the partial plastic structures on the interface of ice VII in liquid water in analogy with the quasi-liquid water layer on ice Ih in the gas phase.52
■
Article
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.7b01764. Descriptions of Supporting Information videos (PDF) Nucleation of ice T at the interface between ice VII and liquid water (6 GPa, 300 K). One second in the video corresponds to 50 ps in the simulation. (MPG) Nucleation of ice R in the coexistence of ice VI and liquid water (6 GPa, 300 K). One second in the video corresponds to 50 ps in the simulation. (MPG)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Hideki Tanaka: 0000-0003-2099-8990 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The present work was supported by JSPS KAKENHI grant number JP16K17857 and MEXT as “Priority Issue on PostKcomputer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use).
CONCLUSIONS
We have investigated high-pressure ice phases using MD simulations of the TIP4P/2005 water model. It is found that a new ice phase, named ice T, spontaneously forms at the interface between ice VII and liquid water at P ≈ 6 GPa and T ≈ 300 K. The structure of ice T is similar to that of ice R, which was originally proposed as a precursor state of the homogeneous nucleation of ice VII in liquid water. Both ice T and ice R are partial plastic phases and can be considered to be sibling phases. Ices T and R occupy certain regions in the phase diagram of TIP4P/2005 water as thermodynamically stable phases. However, these ice phase have not been observed experimentally. This implies that they are metastable phases in reality. Water molecules can form packed structures as well as network structures at high pressures. Ice VI is characterized by the highly distorted hydrogen bond network, and ice VII is the densest ice phase. The partial plastic phases emerge between ice VI and ice VII in the phase diagram because the plasticity arises from the trade-off between the network distortion and the packing. There might be other partial plastic phases in this region where the tendencies of water molecules to form network structures and to form packed structures compete. To the best of our knowledge, this is the first study showing the possibility of the existence of crystals belonging to space group I41 with 72 particles in a unit cell. The structure of ice R is the same as that of a high-pressure phase of methane, methane A.17 This fact and the similarly between ice T and ice R suggest that there might exist high-pressure phases of methane or other tetrahedral analogues in which the molecules are located in the ice T arrangement. We have evaluated the existence of new ice phases of the TIP4P/2005 water model around the triple point of ices VI and VII and liquid water. It is necessary to survey other new ice phases around other triple points. Accurate free-energy calculations to determine the stable regions of the new ice phases in the phase diagram are also required.
■
REFERENCES
(1) Petrenko, V. F.; Whitworth, R. W. Physics of Ice; Oxford University Press, 2002. (2) Salzmann, C. G.; Radaelli, P.; Hallbrucker, A.; Mayer, E. The Preparation and Structures of Hydrogen Ordered Phases of Ice. Science 2006, 311 (5768), 1758−1761. (3) Salzmann, C. G.; Radaelli, P.; Mayer, E.; Finney, J. Ice XV: A New Thermodynamically Stable Phase of Ice. Phys. Rev. Lett. 2009, 103 (10), 105701. (4) Falenty, A.; Hansen, T. C.; Kuhs, W. F. Formation and Properties of Ice XVI Obtained by Emptying a Type sII Clathrate Hydrate. Nature 2014, 516 (7530), 231−233. (5) del Rosso, L.; Celli, M.; Ulivi, L. New Porous Water Ice Metastable at Atmospheric Pressure Obtained by Emptying a Hydrogen-Filled Ice. Nat. Commun. 2016, 7, 13394. (6) Kosyakov, V. I.; Shestakov, V. A. On the Possibility of the Existence of a New Ice Phase Under Negative Pressures. Dokl. Phys. Chem. 2001, 376 (4-6), 49−51. (7) Jacobson, L. C.; Hujo, W.; Molinero, V. Thermodynamic Stability and Growth of Guest-Free Clathrate Hydrates: a Low-Density Crystal Phase of Water. J. Phys. Chem. B 2009, 113 (30), 10298−10307. (8) Conde, M. M.; Vega, C.; Tribello, G. A.; Slater, B. The Phase Diagram of Water at Negative Pressures: Virtual Ices. J. Chem. Phys. 2009, 131 (3), 034510. (9) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Anomalous Thermodynamic Properties of Ice XVI and Metastable Hydrates. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93 (5), 054118. (10) Huang, Y.; Zhu, C.; Wang, L.; Cao, X.; Su, Y.; Jiang, X.; Meng, S.; Zhao, J.; Zeng, X. C. A New Phase Diagram of Water Under Negative Pressure: the Rise of the Lowest-Density Clathrate S-III. Science Advances 2016, 2 (2), e1501010−e1501010. (11) Huang, Y.; Zhu, C.; Wang, L.; Zhao, J.; Zeng, X. C. Prediction of a New Ice Clathrate with Record Low Density: a Potential Candidate as Ice XIX in Guest-Free Form. Chem. Phys. Lett. 2017, 671, 186−191. (12) Takii, Y.; Koga, K.; Tanaka, H. A Plastic Phase of Water From Computer Simulation. J. Chem. Phys. 2008, 128 (20), 204501.
11568
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
Article
Langmuir (13) Aragones, J. L.; Conde, M. M.; Noya, E. G.; Vega, C. The Phase Diagram of Water at High Pressures as Obtained by Computer Simulations of the TIP4P/2005 Model: the Appearance of a Plastic Crystal Phase. Phys. Chem. Chem. Phys. 2009, 11 (3), 543−555. (14) Aragones, J. L.; Vega, C. Plastic Crystal Phases of Simple Water Models. J. Chem. Phys. 2009, 130 (24), 244504. (15) Himoto, K.; Matsumoto, M.; Tanaka, H. Yet Another Criticality of Water. Phys. Chem. Chem. Phys. 2014, 16, 5081−5087. (16) Matsumoto, M.; Himoto, K.; Tanaka, H. Spin-One Ising Model for Ice VII−Plastic Ice Phase Transitions. J. Phys. Chem. B 2014, 118 (47), 13387−13392. (17) Maynard-Casely, H. E.; Bull, C. L.; Guthrie, M.; Loa, I.; McMahon, M. I.; Gregoryanz, E.; Nelmes, R. J.; Loveday, J. S. The Distorted Close-Packed Crystal Structure of Methane A. J. Chem. Phys. 2010, 133 (6), 064504. (18) Mochizuki, K.; Himoto, K.; Matsumoto, M. Diversity of Transition Pathways in the Course of Crystallization Into Ice VII. Phys. Chem. Chem. Phys. 2014, 16 (31), 16419−16425. (19) Liu, E.; Dore, J. C.; Webber, J. B. W.; Khushalani, D.; Jähnert, S.; Findenegg, G. H.; Hansen, T. Neutron Diffraction and NMR Relaxation Studies of Structural Variation and Phase Transformations for Water/Ice in SBA-15 Silica: I. the Over-Filled Case. J. Phys.: Condens. Matter 2006, 18 (44), 10009−10028. (20) Webber, J. B. W.; Dore, J. C.; Strange, J. H.; Anderson, R.; Tohidi, B. Plastic Ice in Confined Geometry: the Evidence From Neutron Diffraction and NMR Relaxation. J. Phys.: Condens. Matter 2007, 19 (41), 415117. (21) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123 (23), 234505−234505−12. (22) Noya, E. G.; Menduiña, C.; Aragones, J. L.; Vega, C. Equation of State, Thermal Expansion Coefficient, and Isothermal Compressibility for Ices Ih, II, III, v, and VI, as Obtained From Computer Simulation†. J. Phys. Chem. C 2007, 111 (43), 15877−15888. (23) Vega, C.; Abascal, J. L. F.; Nezbeda, I. Vapor-Liquid Equilibria From the Triple Point Up to the Critical Point for the New Generation of TIP4P-Like Models: TIP4P/Ew, TIP4P/2005, and TIP4P/Ice. J. Chem. Phys. 2006, 125 (3), 034503. (24) Pi, H. L.; Aragones, J. L.; Vega, C.; Noya, E. G.; Abascal, J. L. F.; Gonzalez, M. A.; McBride, C. Anomalies in Water as Obtained From Computer Simulations of the TIP4P/2005 Model: Density Maxima, and Density, Isothermal Compressibility and Heat Capacity Minima. Mol. Phys. 2009, 107 (4−6), 365−374. (25) Abascal, J. L. F.; Vega, C. Widom Line and the Liquid−Liquid Critical Point for the TIP4P/2005 Water Model. J. Chem. Phys. 2010, 133 (23), 234502. (26) Overduin, S. D.; Patey, G. N. An Analysis of Fluctuations in Supercooled TIP4P/2005 Water. J. Chem. Phys. 2013, 138 (18), 184502. (27) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Spontaneous LiquidLiquid Phase Separation of Water. Phys. Rev. E 2014, 89 (2), 020301. (28) Limmer, D. T.; Chandler, D. Corresponding States for Mesostructure and Dynamics of Supercooled Water. Faraday Discuss. 2013, 167 (0), 485−498. (29) Benavides, A. L.; Aragones, J. L.; Vega, C. Consensus on the Solubility of NaCl in Water From Computer Simulations Using the Chemical Potential Route. J. Chem. Phys. 2016, 144 (12), 124504. (30) Aragones, J. L.; Rovere, M.; Vega, C.; Gallo, P. Computer Simulation Study of the Structure of LiCl Aqueous Solutions: Test of Non-Standard Mixing Rules in the Ion Interaction. J. Phys. Chem. B 2014, 118 (28), 7680−7691. (31) Yagasaki, T.; Matsumoto, M.; Andoh, Y.; Okazaki, S.; Tanaka, H. Dissociation of Methane Hydrate in Aqueous NaCl Solutions. J. Phys. Chem. B 2014, 118 (40), 11797−11804. (32) Conde, M. M.; Vega, C. Determining the Three-Phase Coexistence Line in Methane Hydrates Using Computer Simulations. J. Chem. Phys. 2010, 133 (6), 064507. (33) Liang, S.; Kusalik, P. G. Nucleation of Gas Hydrates Within Constant Energy Systems. J. Phys. Chem. B 2013, 117 (5), 1403−1410.
(34) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Adsorption Mechanism of Inhibitor and Guest Molecules on the Surface of Gas Hydrates. J. Am. Chem. Soc. 2015, 137 (37), 12079−12085. (35) Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition Path Sampling and the Calculation of Rate Constants. J. Chem. Phys. 1998, 108 (5), 1964. (36) Zahn, D. Atomistic Mechanism of NaCl Nucleation From an Aqueous Solution. Phys. Rev. Lett. 2004, 92 (4), 040801. (37) Fernández, R. G.; Abascal, J. L. F.; Vega, C. The Melting Point of Ice Ih for Common Water Models Calculated From Direct Coexistence of the Solid-Liquid Interface. J. Chem. Phys. 2006, 124, 144506. (38) Morris, J. R.; Song, X. The Melting Lines of Model Systems Calculated From Coexistence Simulations. J. Chem. Phys. 2002, 116 (21), 9352−9358. (39) Conde, M. M.; González, M. A.; Abascal, J. L. F.; Vega, C. Determining the Phase Diagram of Water From Direct Coexistence Simulations: the Phase Diagram of the TIP4P/2005 Model Revisited. J. Chem. Phys. 2013, 139 (15), 154505. (40) Smirnov, G. S.; Stegailov, V. V. Melting and Superheating of sI Methane Hydrate: Molecular Dynamics Study. J. Chem. Phys. 2012, 136, 044523. (41) Myshakin, E. M.; Jiang, H.; Warzinski, R. P.; Jordan, K. D. J. Phys. Chem. A 2009, 113 (10), 1913−1921. (42) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511. (43) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (45) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: an N· Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089−10092. (46) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (47) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (48) Wolfson, H. J.; Rigoutsos, I. Geometric Hashing: an Overview. IEEE Comput. Sci. Eng. 1997, 4 (4), 10−21. (49) Meijer, E. J.; Frenkel, D.; LeSar, R. A.; Ladd, A. J. C. Location of Melting Point at 300 K of Nitrogen by Monte Carlo Simulation. J. Chem. Phys. 1990, 92 (12), 7570−7575. (50) Pauling, L. The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement. J. Am. Chem. Soc. 1935, 57 (12), 2680−2684. (51) Klotz, S.; Komatsu, K.; Pietrucci, F.; Kagi, H.; Ludl, A. A.; Machida, S.; Hattori, T.; Sano-Furukawa, A.; Bove, L. E. Ice VII from Aqueous Salt Solutions: From a Glass to a Crystal with Broken HBonds. Sci. Rep. 2016, 6 (1), 32040. (52) Sazaki, G.; Zepeda, S.; Nakatsubo, S.; Yokomine, M.; Furukawa, Y. Quasi-Liquid Layers on Ice Crystal Surfaces Are Made Up of Two Different Phases. Proc. Natl. Acad. Sci. U. S. A. 2012, 109 (4), 1052− 1055.
11569
DOI: 10.1021/acs.langmuir.7b01764 Langmuir 2017, 33, 11561−11569
