E, and TIP5P ... - ACS Publications

Subscriber access provided by UNIV OF DURHAM

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Phase Diagrams of TIP4P/2005, SPC/E, and TIP5P Water at High Pressure Takuma Yagasaki, Masakazu Matsumoto, and Hideki Tanaka J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b04441 • Publication Date (Web): 17 Jul 2018 Downloaded from http://pubs.acs.org on July 19, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Phase Diagrams of TIP4P/2005, SPC/E, and TIP5P Water at High Pressure Takuma Yagasaki, Masakazu Matsumoto*, and Hideki Tanaka Research Institute for Interdisciplinary Science, Okayama University, Okayama, 7008530, Japan

Abstract We investigate high-pressure ice phases using molecular dynamics simulations. Spontaneous nucleation of a new crystalline solid, named ice T2, is observed in a simulation of TIP4P/2005 water at 260 K and 3.3 GPa. The phase diagram of ices VI, VII, T2 and recently reported two other hypothetical ices, ice R and ice T, is determined using the direct coexistence method and the Clausius-Clapeyron equation for TIP4P/2005, SPC/E and TIP5P water. It is found that there exists at least one pressure region in which a hypothetical ice phase is the most stable at ambient temperature with those models. Although the hypothetical ices may be metastable in reality, these ices could be of great importance toward a comprehensive understanding of the phase behaviors of water including many metastable ice polymorphs settled in the hidden area of T-P space. The unit cell of ice T2 is tetragonal with a space group of I41/acd and it contains 152 water molecules. This is probably the most entangled structure among crystals which have been found in nucleation simulations without bias.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction There exist eighteen crystalline phases of pure water.1-5 There also exist many crystalline phases of aqueous solutions: clathrate hydrates, semiclathrate hydrates, and filled ices.6-12 The variety of crystalline structures arises from the flexibility and versatility of the hydrogen bond network of water. The chair and boat forms of sixmembered rings are the most stable fragments of the hydrogen bond network, and ice at ambient pressure, Ih, consists only of these two types of rings. Other types of ring structures are found in high-pressure ices and clathrate hydrates. For example, CS-I and CS-II clathrate hydrates consist of the flat five- and six-membered rings, and there are highly distorted four-membered rings in ice VI. The number of ice phases is still increasing. Kosyakov and Shestakov, Conde et al., and Jacobson et al. suggested that empty CS-II hydrate is thermodynamically stable under negative pressure on the basis of theoretical calculations.13-15 In 2014, Falenty et al. succeeded to prepare the empty CS-II hydrate experimentally by vacuum pumping of Ne hydrate.4 This crystal was named ice XVI. In 2016, Rosso et al. reported a new ice phase, ice XVII, obtained by emptying C0 filled ice.5 Ices XVI and XVII are less dense than ice Ih and they are metastable at ambient pressure. Several other low-density ice phases have been proposed in recent simulation studies.16-18 There might be high-density metastable ices that have not been found in experimental studies. It was reported that ice VII, which is the densest molecular crystal of water, becomes a plastic crystal at very high temperature and pressure (T ~ 560 K and P ~ 8 GPa for TIP5P water) in classical molecular dynamics (MD) simulations.19 The oxygen atoms are located in a BCC arrangement while each water molecule rotates freely in the plastic phase.20-21 A recent ab initio MD study also supports the potential existence of BCC plastic ice.22 A different plastic phase, FCC plastic ice, forms at higher temperature and pressure.23-24 These plastic ice phases have not been observed experimentally, suggesting that they may be metastable ices for real water. It should be noted that neutron diffraction experiments showed that the orientation of water molecules in LiCl doped ice VII is disordered similarly to that of the BCC plastic ice phase.25


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The TIP4P/2005 model is a rigid non-polarizable four-site model for water.26 This model was developed to reproduce several properties of ices Ih, II, III, and V as well as liquid water. The phase diagram of TIP4P/2005 water is similar to that of real water for P < 1 GPa.23 Recently, we found rhombohedral and tetragonal high-density ices that can be stable phases for the TIP4P/2005 model under very high pressure.27-28 The two phases, named ice R and ice T, are partially plastic phases in which certain water molecules can rotate. The structures of these ices consist of densely packed parts that resemble the ice VII structure and distorted hydrogen bond network structures that resemble the ice VI structure. We evaluated the chemical potentials of ices R, T, VI, and VII at 300 K.28 The obtained VI/R, R/T, and T/VII coexistence pressures are 3, 7.5, and 17.5 GPa, respectively. Ice R and ice T have not been found experimentally. In this study, we perform MD simulations of TIP4P/2005 water and observe spontaneous nucleation of a new tetragonal crystalline phase at T = 260 K and P = 3.3 GPa. The new phase is named ice T2. The phase diagram of ices VI, VII, R, T, T2, and liquid water is determined using the direct coexistence method29 and the ClausiusClapeyron equation30 not only for the four-site TIP4P/2005 model but also for three-site SPC/E water31 and five-site TIP5P water.32 It is shown that all three models cannot reproduce the phase diagram of real water at high pressure.

Computational Details The GROMACS 4.6 package is employed for MD simulations.33-34 The temperature and pressure are maintained using a Nosé-Hoover thermostat and a Berendsen barostat.35-36 The particle mesh Ewald method is used for long-range Coulomb interactions with a real space cutoff length of 0.9 nm.37-38 We consider a cubic simulation cell containing randomly located 27000 TIP4P/2005 water molecules. An equilibration simulation is performed for 1 ns at T = 300 K and P = 1 bar. Nucleation simulations are started from the final configuration of the equilibration run. The simulation time is 100 ns. We also perform single-phase MD simulations with smaller simulation cells. The number of water molecules is 2160, 3456, 2625, 3240, 2432, and 2500 for ice VI, VII, R, T, T2, and liquid water, respectively. The crystal configurations are generated using



the GenIce tool.39 The length of the simulations is 3 ns and the last 1 ns is used to calculate the oxygen-oxygen radial distribution function and the orientational autocorrelation function. The crystal/liquid coexistence temperature is obtained at several pressures using the direct coexistence method.29, 40-41 We prepare rectangular cells of the crystals using the GenIce tool.39 The number of water molecules is 4320, 9216, 4200, 5184, and 4864 for ice VI, VII, ice R, ice T, and ice T2, respectively. An isochoric and isothermal MD simulation is performed for 100 ps at T = 2000 K with constraints that fix the coordinates of a half of molecules in the cell to make a crystal/liquid coexistence state. Then, isobaric and isothermal simulations are performed for several temperatures at each pressure. We also calculate the melting temperature of ice II using the direct coexistence method for SPC/E water because it was reported that this water model overestimates the thermodynamic stability of ice II in the high pressure region focused in this study, P > 1 GPa.42 The number of molecules is 4320 for the coexistence of ice II and liquid water.

Results and Discussion Figure 1a presents the time evolution of the potential energy for simulations of 27000 TIP4P/2005 water molecules at T = 260 K and P = 3.3 GPa. We find a decrease in the potential energy for five out of eight trajectories. The density increases at the same time as shown in Figure 1b. We divide a trajectory into time bins with a width of 1 ns and calculate the following quantity for each bin = 〈 + ∆ − 〉 1 where ri is the coordinate vector of the oxygen atom of i-th water molecule.43-44 We set ∆t = 0.1 ns. Water molecules with δ2 < 0.004 nm2 are assumed to be immobile molecules. Figure 2 shows snapshots of a simulation at elapsed times of 8, 16, 24, 32, 48, and 64 ns. The hydrogen bonds between the immobile molecules are represented by black lines. We observe formation and annihilation of clusters of immobile molecules. For example, all clusters found at t = 8 ns disappear until t = 16 ns. A cluster forms at t = 16 ns (marked with a circle) and grows gradually with increasing time. The size of the cluster begins to increase rapidly at t ~ 40 ns. Eventually, the simulation cell is filled


Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with immobile water molecules. This behavior is consistent with the time evolution of the potential energy shown by the red curve in Figure 1a.

Figure 1. Time evolution of the (a) potential energy and (b) density of eight nucleation simulations starting from eight different liquid configurations. The temperature and pressure are 260 K and 3.3 GPa, respectively. The TIP4P/2005 model is employed.



Page 6 of 33

Figure 2. Snapshots of a nucleation simulation at 260 K and 3.3 GPa. The potential and the density of this trajectory are shown by red curves in Figure 1. The hydrogen bonds between the immobile molecules are represented by black lines.

A large part of the final structure of the trajectory shown in Figure 2 seems to be a single crystal of unknown structure (as shown below, the radial distribution function of this solid is different from that of any other high-pressure ice). In order to identify the basic repeating unit (i.e., unit cell) of this structure, we apply a pattern matching technique which evaluates the similarity between local structures of selected areas in the simulation cell.39 This method is useful to detect crystalline structures without a priori knowledge of the structure. The procedure is given in Ref. 28. We obtain a bunch of translational vectors connecting two molecules having similar local molecular configurations for a spherical region of radius r = 0.8 nm from each molecule. The primary vectors of the unit cell can be determined from this bunch of vectors. The unit cell is tetragonal and has body-centered symmetry, and its dimensions are 2 × 2 × 0.75 nanometers. This result is then utilized to obtain the average molecular density distribution in the unit cell and to determine the space group of the crystal. The positions of water molecules in a unit cell are listed in Table 1. We name this crystal ice

T2.

Table 1. Oxygen atom fractional coordinate for six types of water in the ice T2 structure. The space group is I41/acd with a = 1.9912 nm and c = 0.7575 nm.

x

y

z

W1

0.000

0.000

0.250

W2

0.109

0.109

0.750

W3

0.074

0.674

0.903

W4

0.061

0.235

0.706

W5

0.026

0.111

0.043

W6

0.185

0.204

0.034


Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The structure of ice T2 is shown in Figure 3. There are six types of water molecules. The space group of ice T2 is I41/acd and the unit cell contains 152 molecules. Such crystalline structure has not been reported.45 The unit cell of ice T2 is larger than that of any ice that can be prepared experimentally (the unit cell of ice XVI consisting of 136 molecules is the largest one).4 The Ice T2 structure is probably the most entangled among crystals which have nucleated spontaneously in MD simulations.46-48

Figure 3. Unit cell of ice T2. Water molecules of types W1 to W6 (cf. Table 1) are depicted as red, blue, cyan, green, yellow, and gray spheres. Hydrogen bonds are shown with white cylinders.

The determined unit cell structure is used as a template of pattern matching for detection of the crystalline structure emerged in the MD trajectory.27-28,

39

A local

structure is considered to be matched to the template when the root mean square displacement of atomic coordinates between them is less than 0.07 nm. The blue dots in Figure 4 are water molecules of ice T2 detected by the pattern matching. The comparison between Figures 2 and 4 demonstrates that the immobile molecules indeed form the ice T2 structure. The cluster at t = 16 ns is not detected by the pattern matching because it is smaller than the unit cell.



Figure 4. Results of the pattern matching using the T2 unit cell. The trajectory is the same as that shown in Figure 2. The detected T2 unit cells and water molecules in them are shown by boxes and dots. The box color is different for different cell orientations. Concentric patterns of T2 unit cells are seen clearly.

Nucleation of ice T2 does not occur at T = 260 K and P = 3.2 GPa within the simulation time of 100 ns. Multiple nuclei form without any induction time at P = 3.6 GPa, indicating that the free energy barrier for crystallization becomes small enough to be quickly overcome by thermal fluctuations at this pressure (Figure S1). We also perform MD simulations at a higher temperature, T = 290 K, and find that ice R forms at P ~ 5 GPa. The radial distribution functions (RDF) of high-pressure ices are shown in Figure 5. The RDF of ice T2 is very similar to those of ices R and T, but the first minimum of T2 is farther than those of ices R and T, and a characteristic dent at 0.8 nm found for ices R and T is missing in the RDF of T2.


Page 8 of 33

Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Oxygen-oxygen radial distribution functions of high-pressure ices at 300 K and 6 GPa.

The orientational autocorrelation functions are shown in Figure 6a. The rotational relaxation time of ice T2 is much longer than those of partially plastic ice R and ice T, suggesting that ice T2 does not possess plasticity. Figure 6b shows the orientational autocorrelation functions of the six types of water molecules in ice T2. The relaxation time of W6 molecules (gray spheres in Figure 3) is shorter than those of the other types of water molecules because each of W6 water molecules can form either three hydrogen bonds or one bifurcated hydrogen bond together with the three hydrogen bonds.

Figure 6. Orientational autocorrelation functions of (a) ices VI, VII, R, T, T2 and liquid water and (b) six types of water molecules in ice T2 at T = 300 K and P = 6 GPa.



We determine the melting temperature of ice T2 at several pressures using the direct coexistence method.29 Figure 7 shows the time evolution of the potential energy of the T2/liquid coexistence system at P = 3.2 GPa. The potential energy decreases at T = 324 K because of crystal growth while the potential energy increases at T = 328 K because of dissociation of the crystal. We assume that the melting temperature of ice T2 is 326 K at this pressure. The melting temperature is also calculated for ice VI, R, and T. The filled symbols in Figure 8 show the obtained melting temperatures.

Figure 7. Time evolution of the potential energy of the T2/liquid coexistence system at P = 3.2 GPa for TIP4P/2005 water.

Figure 8. Crystal/liquid (solid) and crystal/crystal (dotted) coexistence curves for TIP4P/2005 water calculated using the Clausius-Clapeyron equation. The filled symbols show the melting temperatures determined from the direct coexistence simulations.


Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A crystal/liquid coexistence curve is obtained from a coexistence point using the Clausius-Clapeyron equation.23, 30 The molar enthalpy and volume of phase i, hi and vi, and those of phase j, hj and vj, are obtained from two independent single-phase MD simulations performed at a coexistence point, Tco and Pco. An increment in the temperature, ∆T, is chosen. We set ∆T = 2 K or –2 K in this study. Then, the coexistence pressure at the new temperature is calculated from

∆ =

ℎ − ℎ ∆ . 2 −

When the slope of coexistence curve, dPco/dTco, is very large, we set ∆P = 0.5 GPa and obtain ∆T using Eq. (2). The T2/liquid coexistence curve calculated from a coexistent point at Tco = 326 K and Pco = 3.2 GPa is shown by the solid cyan curve in Figure 8. The coexistence curve passes through two other coexistence points determined by the direct coexistence method. The other crystal/liquid coexistence curves are also consistent with the results of the direct coexistence simulations. As shown in Figure 1, nucleation of ice T2 is observed in five out of eight trajectories within the simulation time of 100 ns at T = 260 K and P = 3.3 GPa. This temperature is approximately 70 K lower than the melting temperature of ice T2. The potential energy and the density change abruptly after the nucleation and no precursors (i.e. step-down or decay of the potential energy, emergence of crystallites of metastable phases, etc.) are observed prior to the nucleation. The nucleation rate of ice T2 at P = 3.3 GPa is higher than that of ice Ih under ambient pressure.49 This can be explained from the fact that the melting temperature is much higher for ice T2 (~330 K) than for ice Ih (~250 K for TIP4P/2005 water) and thus dynamics of water molecules is much faster in the nucleation of ice T2 than in the nucleation of ice Ih when they are compared at the same supercooling. However, it remains unsolved why such a complex structure can nucleate from liquid spontaneously and quickly without the help of formation of precursor states. It is possible that the surface free energy between the crystal and liquid, which is a dominant factor for the nucleation rate, is lower for ice T2 than for ice Ih because of the low symmetry of T2. Evaluation of the crystal/liquid surface free energy is time consuming and beyond the scope of this study. Further studies are needed to determine the main causes of the rapid nucleation.



The intersection of the VI/liquid and T2/liquid coexistence curves in Figure 8, T = 309 K and P = 2.7 GPa, is the VI/T2/liquid triple point. The VI/T2 coexistence curve is calculated from this triple point. Similarly, the T2/R and R/T coexistence curves are calculated from the T2/R/liquid triple point at T = 327 K and P = 3.2 GPa and the R/T/liquid triple point at T = 679 K and P = 12.6 GPa. It is impossible to realize the VII/liquid coexistence for TIP4P/2005 water because either transformation from VII to BCC plastic ice or formation of ice T at the VII/liquid interface occurs very quickly in MD simulations.19,

23, 28

However, it is possible to estimate the VII/T coexistence

pressure at a selected temperature from the enthalpies of the two crystals with an assumption that the entropy of ice VII is the same as that of ice VI28 (The difference in entropy between two non-plastic proton-disordered ice phases is very small. This is evident from the slope of coexistence curve, dPco/dTco = ∆Sco/∆Vco, in the phase diagram of real water shown in Figure 9d). The VII/T coexistence curve is obtained from the coexistence pressure at 300 K estimated under this assumption. Figure 9a shows the phase diagram of TIP4P/2005 water constructed from the crystal/crystal and crystal/liquid coexistence curves. The stable regions of ice T2, R, and T exist between those of ice VI and VII at ambient temperature. Note that ices T2, R, and T have been discovered because they form spontaneously in MD simulations. There is always the possibility that an ice phase is more stable than the phases considered to construct the phase diagram but it has not been found because of slow nucleation (the nucleation of the most stable phase is not always the fastest, which is known as Ostwald’s step rule).


Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9. Phase diagram of water for the (a) TIP4P/2005, (b) SPC/E, and (c) TIP5P models. The experimental phase diagram is given in panel (d) for comparison. We do not consider ice VIII, which is the hydrogen-ordered counterpart of ice VII, to construct the phase diagrams because it was reported that the transition temperature between ice VII and ice VIII is quite low (< 100 K) in classical molecular simulations.24



We also determine the phase diagrams of SPC/E and TIP5P water. As shown in Figure 9b, ice T2 and ice R are stable high-pressure crystalline phases for SPC/E water. There exists a region of ice II although the transition pressure between ice II and ice VI is ~ 0.6 GPa in experiments. The stable region of ice VII exists in the very high pressure regime, P > 10 GPa, for SPC/E and TIP4P/2005 water while ice VII forms P ~ 2 GPa in experiments. This indicates that the thermodynamic stability of ice VII is significantly underestimated for these two water models. Figure 9c presents the phase diagram of TIP5P water. Four point-charge sites are located in a tetrahedral arrangement in the TIP5P model. This charge distribution stabilizes the ice VII structure which is interpenetrating diamond lattices. However, there still exists the region of hypothetical ice T2 between the regions of ice VI and ice VII in the phase diagram. Vega and colleagues calculated the phase diagram of high pressure ice phases using free energy calculations based on the Einstein crystal model and the Clausius-Clapeyron equation.24,

42

Figures S3a and S3b present the phase diagram of TIP4P/2005 water

calculated in the present study and that by Vega et al. Only the liquid/VI coexistence curve is shared in the two phase diagrams because Vega et al. did not consider ices R, T, and T2. The location of the liquid/VI coexistence curve is almost the same for the two studies. Figures S3c and S3d present the phase diagram of SPC/E water in this study and that of Vega et al. The liquid/II and II/VI coexistence curves in Figure S3c are quite similar to those in Figure S3d. It is often observed that a metastable solid phase quickly forms and the most stable solid phase cannot nucleate because of the slow dynamics of molecules in the metastable solid. For example, the transition from diamond to graphite never occurs at ambient condition. Therefore, it is possible that ice R, T, or T2 is indeed more stable than ice VII at a certain thermodynamic state for real water but the emergence of ice R, T, or T2 is prohibited by quickly formed ice VII. However, such kinetic phenomenon can occur only when the thermal energy is quite low compared with the binding energy between molecules. The temperature of ~300 K is not so low for hydrogen bonding of water. Therefore, probably, ices R, T, and T2 are metastable phases in reality. TIP4P/2005, SPC/E, and TIP5P are non-polarizable rigid water models. One may guess that these approximations cause the difference between the calculated and experimental phase diagrams. Experimental studies have shown that the length of the 14 ACS Paragon Plus Environment

Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


intramolecular O-H bond increases with increasing pressure.50-51 However, this elongation is evident only at very high pressure, P > 20 GPa. Therefore, it is unlikely that the emergence of ice R, ice T, and ice T2 is due to the rigidity of the models. The three water models tend to underestimate the stability of ice VII. Ice VII is a highly symmetric crystal and all water molecules are equivalent in it.1 It is expected that the introduction of polarizability, which can stabilize molecules at anisotropic condition such as the surface,52-53 would not improve the underestimation of the stability of ice VII. A possible explanation for the poor description of the phase diagram is that the 12-6 Lennard-Jones potential is employed in the three models although the repulsive r–12 term has no theoretical justification. The r–12 term has been used in classical potential models in order to reduce computational cost. The repulsive interaction should be softer than the r–12 term. It is likely that the difference between the true repulsive interaction and the r–12 term can be ignored at ambient pressure but becomes noticeable for high pressure ices.54 If so, it is worth trying to parameterize a water potential model for high pressure using a softer repulsive interaction such as the exp-6 Buckingham function.

Conclusions High-pressure ices have been examined using MD simulations. It is found that TIP4P/2005 water crystallizes spontaneously at 260 K and 3.3 GPa. The crystal structure, named ice T2, is extraordinarily complex as a crystal of a pure substance: the unit cell is tetragonal with a space group of I41/acd and it contains 152 water molecules. A part of the hydrogen bond network is quite distorted in ice T2. This crystal has a variety of coordination geometries likewise other two hypothetical high-pressure ice phases, ice R and ice T. To the best of our knowledge, the ice T2 structure is the most entangled among crystals that have been observed in nucleation MD simulations without bias.46-48 It is challenging to design a particle so that its fluid phase crystallizes into a complex structure.55 Formation of complex crystal structures occurs quite rarely in a single component liquid of a simple molecule.56 Most complex crystals of real materials



consist of multiple components.57 Engel et al. reported the formation of a quasicrystal in a single component atomic system using a well-designed artificial interparticle interaction.58 Although the TIP4P/2005 water model was designed to reproduce properties of the real material and there is no intentional bias to make complex crystals, the water model yields the unexpected complex structure. Water molecules find their correct position and orientation for the crystal growth of ice T2 in a certain region of TP space although the repeating unit of this ice is very long. Coulomb interactions in a hydrogen bonding system are not so long-ranged. How do they convey the structural information to the next unit cell? As a purely mathematical problem, it is interesting to examine how much complex structure a simple particle can form spontaneously, and the relation between the complexity of a crystal and its nucleation rate.58-60 The phase diagram of ices at high pressures determined by the direct coexistence method and the Clausius-Clapeyron equation shows that the most stable phase at 300 K alters from liquid water to ices VI, T2, R, T, and VII with increasing pressure for the TIP4P/2005 water model. Ice R and ice T2 appear in the phase diagram of SPC/E water. The stable region of ice T2 also emerges in the phase diagram of TIP5P water. The enhanced stability of the experimentally undiscovered high-pressure ices may originate from the repulsive r–12 term of the Lennard-Jones interaction. Use of a more realistic potential function for the repulsive interaction might improve the phase behavior of water under high pressure in simulations. Recent experimental studies have shown that the stability of ices can be controlled by inclusion of solutes.25, 61-62 For example, the stable region of ice II disappears when a certain amount of NH4F is doped.61 The exotic and complex ice R, ice T, and ice T2 as well as unexpected even more complex ones might be produced experimentally with the help of dopants. Metastable phases can affect properties of stable phases. For example, the liquidliquid phase transition hypothesis suggests that there exists a critical point of two metastable liquid phases in the deeply supercooled region and large fluctuations of liquid water at higher temperatures is attributed to the presence of this critical point.63-67 It is important to find a new metastable phase because it may shed new light on the anomalous behavior of water. Computer simulations will play an important role for the discovery of metastable phases because they allow us to examine water under extreme 16 ACS Paragon Plus Environment

Page 16 of 33

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


or unreal conditions (e.g., very high pressure, negative pressure, and biased states due to unreal interactions) so that metastable phases can be detectable.

Supporting Information A movie of the nucleation process at T = 260 K and P = 3.3 GPa, results of nucleation simulations, a comparison between the phase diagrams of the present study and Vega et al., and a CIF file of ice T2.

Author Information Corresponding Author Phone: +81-(0)86-251-7846 E-mail: [email protected]

Acknowledgments The present work was supported by a grant of MORINO FOUNDATION FOR MOLECULAR SCIENCE and MEXT as "Priority Issue on Post-Kcomputer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use) using computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp180204). Calculations were also performed on the computers at Research Center for Computational Science, Okazaki, Japan.



References 1.

Petrenko, V. F.; Whitworth, R. W., Physics of Ice; Oxford University Press: Oxford, 1999.

2.

Salzmann, C. G.; Radaelli, P. G.; Hallbrucker, A.; Mayer, E.; Finney, J. L., The Preparation and Structures of Hydrogen Ordered Phases of Ice. Science 2006, 311, 1758-1761.

3.

Salzmann, C. G.; Radaelli, P. G.; Mayer, E.; Finney, J. L., Ice XV: A New Thermodynamically Stable Phase of Ice. Phys. Rev. Lett. 2009, 103, 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, 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.

Sloan, E. D.; Koh, C. A., Clathrate Hydrates of Natural Gases; CPC Press: Boca Raton, 2008.

7.

Vos, W.; Finger, L.; Hemley, R.; Mao, H.-K., Novel H2-H2O Clathrates at High Pressures. Phys. Rev. Lett. 1993, 71, 3150-3153.

8.

Fowler, D. L.; Loebenstein, W. V.; Pall, D. B.; Kraus, C. A., Some Unusual Hydrates of Quaternary Ammonium Salts. J. Am. Chem. Soc. 1940, 62, 11401142.

9.

McMullan, R.; Jeffrey, G. A., Hydrates of the Tetra N‐Butyl and Tetra I‐ Amyl Quaternary Ammonium Salts. J. Chem. Phys. 1959, 31, 1231-1234.

10.

Feil, D.; Jeffrey, G. A., The Polyhedral Clathrate Hydrates, Part 2. Structure of the Hydrate of Tetra Iso-Amyl Ammonium Fluoride. J. Chem. Phys. 1961, 35, 1863-1873.

11.

Bai, J.; Angell, C. A.; Zeng, X. C., Guest-Free Monolayer Clathrate and Its Coexistence with Two-Dimensional High-Density Ice. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 5718-5722.

12.

Bai, J.; Zeng, X. C., Polymorphism and Polyamorphism in Bilayer Water Confined to Slit Nanopore under High Pressure. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 21240-21245.


Page 18 of 33

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13.

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, 49-51.

14.

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, 10298-10307.

15.

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, 034510.

16.

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. Sci. Adv. 2016, 2, e1501010.

17.

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 GuestFree Form. Chem. Phys. Lett. 2017, 671, 186-191.

18.

Matsui, T.; Hirata, M.; Yagasaki, T.; Matsumoto, M.; Tanaka, H., Communication: Hypothetical Ultralow-Density Ice Polymorphs. J. Chem. Phys. 2017, 147, 091101.

19.

Takii, Y.; Koga, K.; Tanaka, H., A Plastic Phase of Water from Computer Simulation. J. Chem. Phys. 2008, 128, 204501.

20.

Himoto, K.; Matsumoto, M.; Tanaka, H., Lattice- and Network-Structure in Plastic Ice. Phys. Chem. Chem. Phys. 2011, 13, 19876-19881.

21.

Himoto, K.; Matsumoto, M.; Tanaka, H., Yet Another Criticality of Water. Phys. Chem. Chem. Phys. 2014, 16, 5081-5087.

22.

Hernandez, J. A.; Caracas, R., Proton Dynamics and the Phase Diagram of Dense Water Ice. J. Chem. Phys. 2018, 148, 214501.

23.

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, 543-555.

24.

Aragones, J. L.; Vega, C., Plastic Crystal Phases of Simple Water Models. J. Chem. Phys. 2009, 130, 244504.



25.

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 H-Bonds. Sci. Rep. 2016, 6, 32040.

26.

Abascal, J. L. F.; Vega, C., A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505.

27.

Mochizuki, K.; Himoto, K.; Matsumoto, M., Diversity of Transition Pathways in the Course of Crystallization into Ice VII. Phys. Chem. Chem. Phys. 2014, 16, 16419-16425.

28.

Hirata, M.; Yagasaki, T.; Matsumoto, M.; Tanaka, H., Phase Diagram of TIP4P/2005 Water at High Pressure. Langmuir 2017, 33, 11561-11569.

29.

Ladd, A. J. C.; Woodcock, L. V., Triple-Point Coexistence Properties of the Lennard-Jones System. Chem. Phys. Lett. 1977, 51, 155-159.

30.

Kofke, D. A., Direct Evaluation of Phase Coexistence by Molecular Simulation Via Integration Along the Saturation Line. J. Chem. Phys. 1993, 98, 4149-4162.

31.

Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P., The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269-6271.

32.

Mahoney, M. W.; Jorgensen, W. L., A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910-8922.

33.

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, 435-447.

34.

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, 1701-1718.

35.

Nosé, S.; Klein, M. L., Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055-1076.

36.

Allen, M. P.; Tildesley, D. J., Computer Simulation of Liquids. Clarendon press: Oxford, 1987.

37.

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, 10089-10092.


Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


38.

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G., A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593.

39.

Matsumoto, M.; Yagasaki, T.; Tanaka, H., Genice: Hydrogen-Disordered Ice Generator. J. Comput. Chem. 2018, 39, 61-64.

40.

García Fernández, R.; Abascal, J. L. F.; Vega, C., The Melting Point of Ice Ih for Common Water Models Calculated from Direct Coexistence of the SolidLiquid Interface. J. Chem. Phys. 2006, 124, 144506.

41.

Yagasaki, T.; Matsumoto, M.; Tanaka, H., Mechanism of Slow Crystal Growth of Tetrahydrofuran Clathrate Hydrate. J. Phys. Chem. C 2016, 120, 3305-3313.

42.

Sanz, E.; Vega, C.; Abascal, J. L. F.; MacDowell, L. G., Phase Diagram of Water from Computer Simulation. Phys. Rev. Lett. 2004, 92, 255701.

43.

Vatamanu, J.; Kusalik, P. G., Molecular Dynamics Methodology to Investigate Steady-State Heterogeneous Crystal Growth. J. Chem. Phys. 2007, 126, 124703.

44.

Yagasaki, T.; Matsumoto, M.; Tanaka, H., Adsorption of Kinetic Hydrate Inhibitors on Growing Surfaces: A Molecular Dynamics Study. J. Phys. Chem. B 2018, 122, 3396-3406.

45.

Gražulis, S.; Daškevič, A.; Merkys, A.; Chateigner, D.; Lutterotti, L.; Quirós, M.; Serebryanaya, N. R.; Moeck, P.; Downs, R. T.; Le Bail, A., Crystallography Open Database (Cod): An Open-Access Collection of Crystal Structures and Platform for World-Wide Collaboration. Nucleic Acids Res. 2012, 40, D420D427.

46.

Moon, C.; Taylor, P. C.; Rodger, P. M., Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003, 125, 4706-4707.

47.

Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T., Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326, 1095-1098.

48.

Yi, P.; Locker, C. R.; Rutledge, G. C., Molecular Dynamics Simulation of Homogeneous Crystal Nucleation in Polyethylene. Macromolecules 2013, 46, 4723-4733.

49.

Espinosa, J. R.; Sanz, E.; Valeriani, C.; Vega, C., Homogeneous Ice Nucleation Evaluated for Several Water Models. J. Chem. Phys. 2014, 141, 18C529.



50.

Besson, J. M.; Pruzan, P.; Klotz, S.; Hamel, G.; Silvi, B.; Nelmes, R. J.; Loveday, J. S.; Wilson, R. M.; Hull, S., Variation of Interatomic Distances in Ice VIII to 10 Gpa. Phys. Rev. B 1994, 49, 12540-12550.

51.

Song, M.; Yamawaki, H.; Fujihisa, H.; Sakashita, M.; Aoki, K., Infrared Observation of the Phase Transitions of Ice at Low Temperatures and Pressures up to 50 Gpa and the Metastability of Low-Temperature Ice VII. Phys. Rev. B 2003, 68, 024108.

52.

Jungwirth, P.; Tobias, D. J., Specific Ion Effects at the Air/Water Interface. Chem. Rev. 2006, 106, 1259-1281.

53.

Yagasaki, T.; Saito, S.; Ohmine, I., Effects of Nonadditive Interactions on Ion Solvation at the Water/Vapor Interface: A Molecular Dynamics Study. J. Phys. Chem. A 2010, 114, 12573-12584.

54.

Vega, C.; Abascal, J. L., Simulating Water with Rigid Non-Polarizable Models: A General Perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663-19688.

55.

Engel, M.; Trebin, H.-R., Self-Assembly of Monatomic Complex Crystals and Quasicrystals with a Double-Well Interaction Potential. Phys. Rev. Lett. 2007, 98, 225505.

56.

Johnston, J. C.; Kastelowitz, N.; Molinero, V., Liquid to Quasicrystal Transition in Bilayer Water. J. Chem. Phys. 2010, 133, 154516.

57.

Chai, P.; Abramchuk, M.; Shatruk, M., Synthesis, Crystal Structure, and Magnetic Properties of Giant Unit Cell Intermetallics R117co52+∆sn112+Γ (R = Y, La, Pr, Nd, Ho). Crystals 2016, 6, 165.

58.

Engel, M.; Damasceno, P. F.; Phillips, C. L.; Glotzer, S. C., Computational SelfAssembly of a One-Component Icosahedral Quasicrystal. Nat. Mater. 2015, 14, 109-116.

59.

Dotera, T.; Oshiro, T.; Ziherl, P., Mosaic Two-Lengthscale Quasicrystals. Nature 2014, 506, 208-211.

60.

Haji-Akbari, A.; Engel, M.; Keys, A. S.; Zheng, X.; Petschek, R. G.; PalffyMuhoray, P.; Glotzer, S. C., Disordered, Quasicrystalline and Crystalline Phases of Densely Packed Tetrahedra. Nature 2009, 462, 773-777.


Page 22 of 33

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


61.

Shephard, J. J.; Slater, B.; Harvey, P.; Hart, M.; Bull, C. L.; Bramwell, S. T.; Salzmann, C. G., Doping-Induced Disappearance of Ice Ii from Water’s Phase Diagram. Nat. Phys. 2018, 14, 569-572.

62.

Shin, K.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A., Managing Hydrogen Bonding in Clathrate Hydrates by Crystal Engineering. Angew. Chem. Int. Ed. Engl. 2017, 56, 6171-6175.

63.

Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E., Phase-Behavior of Metastable Water. Nature 1992, 360, 324-328.

64.

Mishima, O.; Stanley, H. E., The Relationship between Liquid, Supercooled and Glassy Water. Nature 1998, 396, 329-335.

65.

Xu, L. M.; Kumar, P.; Buldyrev, S. V.; Chen, S. H.; Poole, P. H.; Sciortino, F.; Stanley, H. E., Relation between the Widom Line and the Dynamic Crossover in Systems with a Liquid-Liquid Phase Transition. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 16558-16562.

66.

Yagasaki, T.; Matsumoto, M.; Tanaka, H., Spontaneous Liquid-Liquid Phase Separation of Water. Phys. Rev. E 2014, 89, 020301.

67.

Palmer, J. C.; Martelli, F.; Liu, Y.; Car, R.; Panagiotopoulos, A. Z.; Debenedetti, P. G., Metastable Liquid-Liquid Transition in a Molecular Model of Water. Nature 2014, 509, 385-388.



TOC Graphic


Page 24 of 33

(a) -49 -50 -51 0

20

40 60 Time / ns

80

100

40 60 Time / ns

80

100

(b)

1.51 Density / g cm-3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Potential energy / kJ mol-1

Page 25 of 33

1.50 1.49 1.48 1.47 0

20

Figure 1 Yagasaki Matsumoto Tanaka

ACS Paragon Plus Environment


Page 26 of 33

8 ns

16 ns

24 ns

32 ns

48 ns

64 ns

Figure 2 Yagasaki Matsumoto Tanaka ACS Paragon Plus Environment

Page 27 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


b

a



Page 28 of 33

8 ns

16 ns

24 ns

32 ns

48 ns

64 ns



Page 29 of 33

16 12 g(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


VII T

8

R T2

4

VI Liquid

0 0.2

0.4

0.6 r / nm

0.8

1.0



(a) 1.0

VI VII T2

C(t)

0.8 0.6

R

0.4

T

0.2 0

1.0

Liquid

0

50

100 t / ps

150

W1 W2 W3

0.9

W4

0.8 0.7

200

(b)

C(t)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 33

W5 W6

0

50

100 t / ps

150

200


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Potential energy / kJ mol-1

Page 31 of 33

-46 328 K

-47 -48 -49

324 K

-50 0

1

2 Time / ns

3

4



5

3

T/Liquid

T2/R

4 P / GPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 33

R/Liquid

VI/T2 T2/Liquid

2 1 240

VI/Liquid 280

320 T/K

360

400



Page 33 of 33

(a) TIP4P/2005 20

VII T

P / GPa

10

R

5

T2 VI 1 200

Liquid

225

250

275

300

325

T/K

(b) SPC/E 20

VII

P / GPa

10

R

5

T2

VI

Liquid

II 1 200

225

250

275

300

325

T/K

(c) TIP5P 20

P / GPa

10

VII

5

T2 VI 1 200

225

Liquid 250

275

300

325

T/K 20

(d) Experimental

10 P / GPa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

VII

VIII

VI 1 200

225

250

275

Liquid 300

325

T/K

ACS Paragon Plus Environment Figure

9 Yagasaki Matsumoto Tanaka

E, and TIP5P ... - ACS Publications

Recommend Documents