Room temperature electrofreezing of water

Good Thursday friends Worldwide from Dr. TJ Gunn in Houston Texas!………

Room temperature electrofreezing of water yields a missing dense ice phase in the phase diagram

Nature Communications 10, Article number: 1925 (2019) 


Water can freeze into diverse ice polymorphs depending on the external conditions such as temperature (T) and pressure (P). Herein, molecular dynamics simulations show evidence of a high-density orthorhombic phase, termed ice χ, forming spontaneously from liquid water at room temperature under high-pressure and high external electric field. Using free-energy computations based on the Einstein molecule approach, we show that ice χ is an additional phase introduced to the state-of-the-art TP phase diagram. The χ phase is the most stable structure in the high-pressure/low-temperature region, located between ice II and ice VI, and next to ice V exhibiting two triple points at 6.06 kbar/131.23 K and 9.45 kbar/144.24 K, respectively. A possible explanation for the missing ice phase in the TPphase diagram is that ice χ is a rare polarized ferroelectric phase, whose nucleation/growth occurs only under very high electric fields.


Ice exhibits an exceptionally rich TP phase diagram due to the extraordinary adaptability of water’s hydrogen-bonding networks to different environmental temperatures (T) and pressures (P). To date, at least 17 crystalline ice phases (ice IhIc, ice II to ice XVII) have been produced in the laboratory1,2,3. A number of “computer ice” phases have also been predicted from molecular dynamics (MD) simulations and density functional theory (DFT) computations, including very-low-density porous ices (density ρ < 0.85 g cm−3) such as s-III4, s-IV5, ice ITT6, and sL;7 low-density ices (0.85 g cm−3 ≤ ρ ≤ 1.0 g cm−3) such as silica-like ice polymorphs8,9,10, ice 011, ice i, and ice i’;12 high-density ices (1.0 g cm−3 < ρ < 1.4 g cm−3) such as ice B13; and super-high-density ices (ρ > 2.0 g cm−3), which entail partial ionization14,15,16,17,18,19,20.

Among the 17 bulk ice phases observed in the laboratory, ice XI is believed to be ferroelectric21,22 and also has been suggested to exist on Uranus and Neptune23,24, although the recent theoretical calculations indicated that it would be antiferroelectric ice in nature25. The ice VIII is antiferroelectric, but it is likely to be ferroelectric in an applied electric field26. The ice XV, the hydrogen-ordered form of ice VI phase, is antiferroelectric (P1¯)(P1¯) according to experimental observation1, whereas it is predicted to be a ferroelectric Cc hydrogen-ordered structure based on local density functional approach27,28. However, Del Ben et al.29 used high-level ab initio computation and predicted that the antiferroelectric phase is indeed the ground state, suggesting that more accurate density-function approaches should be considered (see below). The ice Ic is a metastable ice crystal with hydrogen disordered, but it is predicted to be ferroelectric when all hydrogens are ordered (also named as ice XIc, space group I41md) based on computer simulations30,31. Indeed, as the water dipole moments can add up to produce a net moment or cancel each other, either ferroelectric ice polymorphs can exist for special crystalline structures. The fabrication of bulk ferroelectric ice, however, is still a challenging task, as without the assistance of dopants as catalysts, the phase-transformation time for a single-phase ferroelectric ice is estimated to be on the order of 104years23. Hence, pure bulk ferroelectric ice is rare in nature.

Can a pure bulk ferroelectric ice be produced in the laboratory by other means? The answer to this question is still highly sought today. One possible way could be through the application of an ultrahigh electric field for the electrofreezing of water. Electrofreezing is known to play an important role in many natural processes, ranging from tropospheric dynamics to frost damage in cells32,33,34,35,36,37. In addition, several known ice structures have already been determined from computer simulations of electrofreezing under ultra-high electric fields10,38,39,40,41. Svishchev and Kusalik39,41 showed from their MD simulations that a polar crystal with the structure of cubic ice Ic can be achieved via electrofreezing of supercooled liquid water. In a later work, they demonstrated a MD simulation of the formation of a quartz-like ice polymorph through electrofreezing10. This quartz-like ice structure was originally proposed by Bernal and Fowler42 as a type of dense ice polymorph. Stutmann40 investigated the effects of an ultra-high electric field (tens of V nm−1) on bulk water. When the electric field magnitude was 30 V nm−1, the MD simulation showed that liquid water transforms into a crystal-like structure40, which can be either defective polar cubic ice or amorphous ice. Recently, Hu et al.38 provided simulation details on the behaviour of glassy water in external electric fields, including the formation of a body-centred-cubic (bcc) ice phase at 77 K. This bcc ice phase is polarized ferroelectric ice VIII, as determined by its lattice constant of 3.19 ± 0.17 Å and oxygen–oxygen (O–O) radial distribution function (RDF)43,44. However, Saitta et al.45 showed from ab initio MD simulation that the threshold strength of electric field that makes water molecules dissociate is 3.5 V nm−1. As such, the electric fields considered in aforementioned simulations were considerably higher than the threshold.

In this study, we report the formation of a previously unreported ice structure, termed ice χ, which can be observed to form spontaneously in the MD simulation of liquid water at room temperature and under an electric field below the threshold strength. The field-direction-dependent result indicates that ice χ is a rare ferroelectric phase. DFT calculations also indicate that the ice χ is dynamically stable even in zero field. Most importantly, our free-energy computation shows that ice χ is not merely a new crystalline structure but a “missing ice” phase in the contemporary TP phase diagram of ices and ice χ belongs to the family of high-density ices (1.0 g cm−3 < ρ < 1.4 g cm−3). In the newly obtained phase diagram of ice with the TIP4P/2005 water model, ice χ is located in the high-pressure region between ice II and ice VI at low temperatures, and to the left of ice V at relatively high temperatures. The electric-field-induced crystallization of liquid water may serve as an alternative approach to attain new phase structures of water, particularly the ferroelectric ices.


Dissociation of water under an intense electric field

The external electric field strengths (E) used in our classical MD simulations (see below) are in the range of 0–3.5 V nm−1. In a previous study, Saitta et al.45 showed from ab initio MD simulation that the threshold strength of electric field that can lead to dissociation of water molecule is 3.5 V nm−1. It is noteworthy that Saitta et al.45 used the Perdew–Burke–Ernzerhof (PBE) functional and the Berry theory approach to the description of an external electric field (implemented in Quantum Expresso package). Here, instead of the PBE functional, we employed the state-of-the-art dispersion-corrected vdW-DF2 exchange-correlation functional (also called the rPW86-vdW2 functional) for the ab initio MD simulation (see Supplementary Methodsfor simulation details). We note that the vdW-DF2 functional has been shown to be particularly accurate for computing relative energies and transition pressures for known phases of ice (see Computational Method section). Four independent ab initio MD simulations were performed, two with the electric field of 5.0 V nm−1, while another two with the electric field of 10.0 V nm−1 (Supplementary Movies 14 and Supplementary Fig. 1). Our ab initio MD simulations show that although dissociation of water molecule can be clearly seen, within 1 ps simulation time, in bulk liquid water at 270 K and 10.0 V nm−1, dissociation of water was not observed within 5 ps simulation time for the liquid water at 270 K and 5.0 V nm−1; thus, such an event would be very unlikely to occur when the water was under the electric field <3.5 V nm−1.

Spontaneous formation of ice χ under high electric field

We performed numerous MD simulations, each with the same initial configuration of 560 TIP4P/2005 water molecules in the liquid state with the temperature controlled at 270 K (20 K above the bulk melting temperature Tm of 249 ± 2 K46 for the TIP4P/2005 model), while considering numerous pressures (P) in the range of 0.001–10 kbar (NPTensemble). At P = 5 kbar and E = 2.3 V nm−1, a previously unreported ice phase emerges spontaneously, as shown in Fig. 1a, b and Supplementary Movie 5. All hydrogen atoms in solid ice χ are oriented along the z-axis (see Fig. 1b), i.e., the direction of the external electric field. Thus, ice χ exhibits strong polarization and can be classified as a polar ice. During the liquid-to-solid transition, the diffusion coefficient decreases sharply from 0.46 × 10−5 to 1.38 × 10−9 cm2 s−1. Figure 1c shows the computed O–O pair correlation function, indicating that the nearest and second nearest distances of oxygen atoms in ice χ are 0.276 and 0.334 nm, respectively. The latter distance is much shorter than that in the normal ice Ih (0.495 nm)47. The well-separated peaks and valleys suggest that ice χ has long-range crystalline order. Figure 1d shows a notable decrease in the computed potential energy per water molecule (~ 2.9 kJ mol−1) when the phase transformation occurs from liquid water to the crystalline ice χ, reflecting a strong first-order transition from liquid to ice χ.

Fig. 1

As water molecules have a permanent dipole moment, in an external electric field, water molecules tend to orient along the dipole while maximizing the number of hydrogen bonds with neighbouring water molecules. Hence, in ice χ, the parallel arrangement of water dipoles is the most energetically favourable under high external electric field with the additional driving force of potential energy. The polarization energy is given by ∆W = E × p, where E is the electric field and p is the dipole moment per water molecule. The phase transition is favourable when the polarization energy of the water molecule is considered. Here, the threshold field strength in the electrofreezing simulation is 2.3 V nm−1. Thus, a rough estimate of the polarization energy per water molecule is 10.2 kJ mol−1. The total potential energy U = Uint − ∆W is approximately −63.6 kJ mol−1, where Uint is the interaction potential per water molecule. Compared with the initial liquid water, the potential energy (per molecule) difference is approximately −14.0 kJ mol−1, in which the polarization energy accounts for ~ 73%. Therefore, the high polarization energy from the strong electric field can make the potential energy difference of the system greater than the value of T × ∆S, where S is the entropy. As a result, a phase transition occurs from liquid water to ice χ. Under the electric field, the dipole orientations of water molecules tend to be along the direction of electric field, while the water molecules can still adapt to form the hydrogen-bonding network. Their interplay could lead to a different and yet a more stable solid state.

Interestingly, our MD simulation also shows that depending on the temperature, ice χ remains stable even after switching off the external field. Under the external electric field, when the field orientation is initially along the z-axis and then reversed against the z-axis, a strong hysteresis loop is observed at pressure of 5 kbar and temperature of 200 K (see Fig. 2), confirming that ice χ is not only a polar ice but also a ferroelectric ice. In addition, the permanent electric dipole moment per water molecule is 2.211 D, close to that of a TIP4P/2005 water molecule (2.305 D)48.

Fig. 2

In Fig. 3a, a semi-quantitative pressure vs. electric field strength (PE) phase diagram of the TIP4P/2005 water model at 270 K is plotted. At relatively low electric field (E < 1.0 V nm−1), the water remains in the liquid state. When the electric field is high (E > 2.5 V nm−1), a variety of ice polymorphs are observed, depending on the pressure of the system. The PE phase diagram can be divided into three regions. At P < 3.5 kbar, a polar ice arises, whose topological structure is the same as the previously reported (non-polar) ice B13 (Supplementary Movie 6 and Supplementary Fig. 2). Thus, we term this ice phase polar ice B. At 3.5 kbar ≤ P ≤ 9 kbar, liquid water transforms to solid ice χ. At P > 9 kbar, very-high-density amorphous (VHDA) ice is the more stable solid phase (Supplementary Fig. 3a). The intermolecular RDFs of this amorphous ice are plotted in Supplementary Fig. 3b. The first sharp peak location of amorphous ice is almost the same as that of the O–O RDF of ice χ, which denotes the nearest distance between oxygen atoms in the bulk ice. It is noteworthy that the RDF of ice χ exhibits much longer range correlation than that of the VHDA. As shown in Fig. 3b, the mass density of these ice polymorphs also increases as the pressure increases.

Fig. 3

Computed structural properties of ice χ based on DFT

Next, we optimized the structure of ice χ obtained from the MD simulation using the first-principles DFT method with vdw-DF2 functional. The optimized structure is an orthorhombic crystal with space group Fdd2 (Supplementary Fig. 4). The lattice parameters of the unit cell are a = 24.34 Å, b = 12.53 Å, and c = 4.32 Å. The fractional coordinates are given in Supplementary Table 1. The top view of ice χ in the z-axis direction is displayed in Supplementary Fig. 4, where ice χ can be viewed as a column of fused oval octagons (each octagon has a centre). Based on the local surrounding environment, the 56 water molecules in the unit cell are divided into four different types (T1T2T3, and T4) with a population ratio of 1:2:2:2. All molecular dipoles are aligned along the z-axis, resulting in polar ice χ. As Supplementary Fig. 4d shows, the primitive cell of ice χ includes 14 water molecules and the corresponding lattice parameters are a = 6.63 Å, b = 12.36 Å, c = 13.68 Å, α = 28.90°, β = 64.36°, and γ = 86.74°. The phonon dispersion and density of states are computed by using the density functional perturbation theory (DFPT) method49, confirming dynamic stability of the ice χ (Supplementary Fig. 5).

To investigate stability competition of ice phases with different mass densities, we took the high-density ice VI50 as a reference and considered the low-density ice XI22, ice B13, polar ice B, and high-density ice II51 for the purpose of comparison. For the polar ice B obtained from our MD simulations, its fraction coordination based on DFT optimization (using the dispersion-corrected vdW-DF2 functional, also called the rPW86-vdW2 functional; see below) is given in Supplementary Table 2. The equilibrium volume of the unit cell, average nearest-neighbouring O–O distance, mass density, and lattice cohesive energy from our DFT calculation are summarized in Table 1, and were compared with available experimental data. For ice XI, the calculated mass density (0.927 g cm−3) is very close to the experimental value (0.93 g cm−3)22 and the calculated O–O distance (2.755 Å) between neighbouring water molecules is slightly longer than the experimental value (2.735 Å)22. The calculated lattice cohesive energy (Elatt) of ice XI differs from the experimental value52 by only 1.67 kJ mol−1. For ice II and ice VI, their mass densities of 1.178 and 1.313 g cm−3, and their average O–O distances of 2.785 and 2.815 Å, obtained from our DFT calculations, are also very close to the corresponding experimental values of 1.180 and 1.310 g cm−3, and 2.770 and 2.800 Å, respectively50,51. Overall, the vdW-DF2/DFT functional reasonably describes the intermolecular hydrogen-bonding interactions of ices, as shown in our previous work4. For ice χ, the mass density of 1.272 g cm−3 is between that of ice II (1.178 g cm−3) and ice VI (1.313 g cm−3). The average O–O distance of ice χ is 2.785 Å, which is comparable to that of ice II (2.785 Å) and ice VI (2.815 Å). As shown in Table 1, for ices XI, II, χ, and VI, the lattice cohesive energy decreases with increasing mass density. However, ice B and polar ice B have lower mass densities (1.082 g cm−3 and 1.072 g cm−3) than ice II but also lower cohesive energies.Table 1 Structural data on the ice polymorphsFull size table

Relative stability among various ice polymorphs at 0 K

To compare the relative stability of ferroelectric ice χ with neighbouring ice phases in the phase diagram, we calculated their enthalpies (based on vdw-DF2 computation) under different pressures at 0 K without including the zero-point energy correction. Figure 4 depicts the relative enthalpy (with ice VI being the reference) vs. pressure for ice XI, ice II, ice χ, ice B, and polar ice B. The point at which two curves cross marks the transition pressure between the two corresponding phases at 0 K. At 0 kbar < P < 1.32 kbar, low-density ice XI is the most stable phase with the lowest enthalpy. Next, the higher-density ice II becomes more favourable at P > 1.32 kbar. The transition pressure is close to the 2 kbar value that was previously obtained by Conde et al.53 based on the TIP4P/2005 model. At P > 6.66 kbar, ferroelectric ice χ replaces ice II as the most stable ice polymorph. For P ≥ 20.28 kbar, ice VI becomes more stable than ice χ. As shown in Fig. 4, if ice χ is “missing”, ice II would transform directly into ice VI at P = 9.79 kbar, consistent with the transition pressure of 10 kbar obtained by Conde et al.53. Although neither is the most stable phase, as indicated by the enthalpy curves, ice B and polar ice B have nearly the same stability. Again, our DFT calculations demonstrate that ice χ is one of the most stable high-density ices in the high-pressure region at zero temperature, along with the known phases of ice II and ice VI. In addition, we calculated enthalpies of these ice structures at different pressures and 0 K, using the strongly constrained and appropriately normed54 functional (Supplementary Fig. 6). The results are consistent with vdW-DF2/DFT computation, demonstrating that ice χ is a highly stable high-density ice in the high-pressure region and at zero temperature (more details see Supplementary Information).

Fig. 4

Free-energy computation of TP phase diagram

Lastly, to examine the stability of ice χ at temperatures much higher than 0 K, we performed free-energy calculations using the Einstein molecule approach method. In our previous study, we confirmed that the TIP4P/2005 water model can reasonably simulate the realistic TPphase diagram of water/ice4. Aragones et al.55,56 also showed that the TIP4P/2005 water model can describe the relative energy, critical temperature, and surface tension of liquid water and ice phases well. The TP phase diagram of water/ice is plotted in Fig. 5. Four ice polymorphs, namely, ice Ih (or hydrogen-disordered ice XI), ice II, ferroelectric ice χ, and ice VI, arise in sequence with increasing pressure at low temperature. Extrapolation of the phase boundaries to 0 K gives the corresponding transition pressures of 1.92, 5.14, and 29.73 kbar, respectively, compared with 1.32, 6.66, and 20.28 kbar predicted from the above DFT computations. The predicted transition pressure of 1.92 kbar at 0 K for ice Ih to ice II is in excellent agreement with that of the 2 kbar value previously obtained by Conde et al.53. Different from the previous TP phase diagram, the coexistent line between ice II and ice VI disappears, whereas ferroelectric ice χ occupies a region between ice II and ice VI at low temperature and part of the region of ice V at relatively high temperature. As a result, two new triple points emerge: one for ices II, V, and χ at 6.06 kbar and 131.23 K, and the other for ices χ, V, and VI at 9.45 kbar and 144.24 K. Ice B and polar ice B do not appear in the TP diagram, as they have higher Gibbs free energies. It is noteworthy that the free-energy calculations show that ice χ has the lowest free energy among ice χ, ice XIII, and ice XV in the low-temperature and high-pressure region, indicating that ice χ is more stable than ice XIII and ice XV. Overall, both MD simulations at finite temperature and DFT calculations at 0 K strongly support the existence of ferroelectric ice χ in the TP phase diagram at high pressures.

Fig. 5


We predict a new ferroelectric ice χ in the phase diagram of water. Ferroelectric ice χ has a high mass density of 1.27 g cm−3. The ferroelectric ice χ is also proven to be dynamically stable on the basis of phonon-spectrum DFT computation. In the TP phase diagram of water/ice, ferroelectric ice χ emerges in the high-pressure region, located between ice II and ice VI at low temperatures and occupying some domains of ice V at relatively high temperatures, leading to two triple points at P = 6.06 kbar and T = 131.23 K, and at P = 9.45 kbar and T = 144.24 K, respectively. The appearance of ice χ in the TP phase diagram of water/ice suggests that the ferroelectric ice χ entails high thermodynamic stability. Identification of this ferroelectric ice phase not only reveals a “missing” ice polymorph in the high-pressure region of the state-of-the-art TP phase diagram of water but also provides more precise temperature/pressure conditions for seeking the elusive ferroelectric ice. In light of the requirement of ultrahigh electric field, whether this predicted ice χχ can be produced in the laboratory via electrofreezing of liquid water remains to be an open question.


MD simulations

All MD simulations are performed in the isothermal-isobaric (NPT) ensemble, with an external electric field applied along the z-axis. The MD simulations are undertaken with the GROMACS (GROningen MAchine for Chemical Simulations) package57. For all MD simulations, the initial supercell is a cubic box containing 560 TIP4P/2005 water molecules48 and the temperature is controlled at 270 K (corresponding to room temperature for the TIP4P/2005 model, as the melting point of the bulk ice Ih based on the TIP4P/2005 model is approximately 250 K)46. To map out a semi-quantitative PE phase diagram, we set a series of pressure and electric field values, including P = 0.001, 0.01, 0.05, 0.1, 1.0, 2.0, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 9.5, and 10.0 kbar, whereas the electric field strength is varied from 0 to 3.5 V nm−1 by an increment of 0.1 V nm−1. The external electric field induced an additional force, Fi = qi·E, where Fi is the force induced by the electric field, qi is the charge of each atom, and E is the applied electric field. For all simulations, the time step is set to 2 fs. The equilibration MD run lasted for at least 10 ns and, in some cases, lasted for 200 ns. T and P are controlled by a Nosé–Hoover thermostat58 and a Parrinello–Rahman barostat59, respectively. A cutoff of 1.0 nm is adopted for the LJ interactions and the long-range electrostatic interactions are treated by the slab-adapted Ewald sum method60.

DFT calculations

The relative energies of ice χ and selected ice phases are computed using the DFT methods implemented in the VASP 5.3.5 software package61. The electron-ion interactions are described by the projector augmented wave potential62. To account for the intermolecular dispersion interactions, the exchange-correlation interaction is described by the dispersion corrected vdW-DF2 exchange-correlation functional63 (also called the rPW86-vdW2 functional). We note that the vdw-DF2 functional has been previously shown by Santra et al.61 to be particularly accurate for computing relative energies and transition pressures for known phases of ices. The electron wavefunction is expanded by a plane-wave basis up to 700 eV. The Brillouin zones are sampled by k-point grids with a uniform spacing of 2ππ × 0.04 Å−1. To confirm dynamic stability of ice χ, the phonon dispersion is computed by using the DFPT method49 as implemented in the VASP 5.4.

Monte Carlo/MD simulations of the TP phase diagram

The TP phase diagram of water and ice polymorphs (including ice Ih, ice II, ice V, ice VI, hydrogen-ordered ice XIII and ice XV, previously predicted polar ice B, and predicted ice χ from this study) was derived based on the Einstein model for crystals and the TIP4P/2005 water potential. First, to obtain reliable configurations of the ice polymorphs, isothermal-isobaric Monte Carlo simulations at temperatures from 1 to 200 K (with 10 K increments) and pressures from 1 to 24 kbar (with 1 kbar increments) are performed using a homemade code. For each candidate phase, the configurations from the Monte Carlo simulations are used to calculate the free energy on the basis of the Einstein molecule approach with the GROMACS programme57. At each TP condition, the Ewald sum method with a real-space cutoff of 8.5 Å is adopted to treat the electrostatic interactions and the pair potential is truncated at 8.5 Å. For ice Ih, ice V, and ice VI, the effect of hydrogen disorder is considered in the free-energy computation and their Pauling entropies S/NkB are taken as ln(3/2), 0.3817, and ln(3/2), respectively64. Once the free energy at a reference point is determined, the thermodynamic integration method can be used to evaluate the free energy under other thermodynamic conditions. Specifically, an initial coexistent point is located by equating the chemical potentials of two phases at a given temperature and pressure65. Next, the Gibbs–Duhem integration based on the trapezoid predictor-corrector formulas is performed to compute the phase boundaries66.