More than 300 kinds of porous ice structures derived from zeolite frameworks and space fullerenes are examined using classical molecular dynamics simulations. It is found that a hypothetical zeolitic ice phase is less dense and more stable than the sparse ice structures reported by Huang *et al.* [Chem. Phys. Lett. **671**, 186 (2017)]. In association with the zeolitic ice structure, even less dense structures, “aeroices,” are proposed. It is found that aeroices are the most stable solid phases of water near the absolute zero temperature under negative pressure.

## INTRODUCTION

Water is known to have many polymorphs. Seventeen ice polymorphs, including metastable ones, have been found experimentally.^{1–6} It is uncommon that a pure substance has such a lot of crystal structures. The variety of ice polymorphs arises from the fact that water molecules prefer the tetrahedrally directed network topology due to the hydrogen bonds.

Some of the ice polymorphs were discovered very recently. They are of particular interest because they are less dense than normal ice and are supposed to be stable under negative pressure.^{7} Falenty *et al.* succeeded to produce empty sII clathrate hydrate by vacuum pumping of neon hydrate in 2014, and this ice was named ice XVI.^{5} Existence of this porous ice was originally predicted theoretically,^{8} and its properties were surveyed by computer simulations.^{9,10} The same methodology was applied to hydrogen-filled ice C_{0} to produce ice XVII.^{6} Many hypothetical ices have been produced by computer simulations.^{10–12} Conde *et al.* predicted that empty sH clathrate hydrate is more stable than any other known ice phases under extreme negative pressures.^{10} Recently, Huang *et al.* proposed that two ultralow-density ice phases occupy wide areas in the phase diagram of water at negative pressures.^{11,12} They searched for possible stable ice structures by a Monte Carlo packing algorithm and evaluated thermodynamic stabilities of them by the combination of DFT calculations and classical molecular dynamics (MD) simulations with the TIP4P/2005 water force field. The network structures of the two stable low-density ice phases are the same as the zeolite frameworks of RHO and FAU.^{13} We call them ice RHO and ice FAU in this paper. The network structures of ice XVI and the empty sH clathrate hydrate are the same as the MTN and DOH zeolite frameworks, and thus we refer to them as ice MTN and ice DOH, respectively.

Although the vapor phase is thermodynamically more stable than any other phase on approaching to the limit of zero pressure, solid phases can be metastable at low temperatures in the negative pressure region. A survey of metastable ice phases would be helpful to understand the complex phase behavior of water in confined geometries.^{14–22} Consideration on possible metastable phases is also important because they might affect kinetics of phase transitions.^{23–25} Although various porous ice structures have been proposed, the comprehensive knowledge of them is still lacking. We pose a question to disambiguate the problem: what is the most stable solid phase under negative pressure?

Hypothetical ice structures derived from zeolite frameworks can be good candidates for the stable ice phases under negative pressure.^{26} Although more than 200 of zeolite frameworks have been reported, most of them have not been considered in the early studies on porous ices.^{9–11,27} Some of ice analogues might be more stable than either ice MTN, DOH, RHO, or FAU in certain regions in the phase diagram at negative pressures. There are many clathrate hydrate structures other than several well-known ones such as sII and sH, and they also can be candidates for stable negative pressure ice.^{28–31}

In this paper, we examine exhaustively ice structures of 219 zeolite and 84 clathrate hydrate frameworks using classical MD simulations. It is found that a zeolitic ice is more stable than either ice MTN, DOH, RHO, or FAU at deeply negative pressures. We also demonstrate that it is possible to design very sparse ice structures of arbitrary density from several zeolitic ices, and they can be more stable than any of zeolitic ices and empty clathrate hydrates under negative pressure at low temperatures.

## METHOD

The GenIce tool^{32} is used to generate zeolitic ice structures. First, SiO_{2} zeolite structures are taken from the Zeolite Database.^{13} The oxygen atoms are removed and the silicon atoms are replaced by oxygen atoms. Then, hydrogen atoms are placed so that the structure obeys the ice rule and the net dipole moment of the system becomes zero. We do not consider proton-ordered structures unless otherwise noted.

Clathrate hydrates are made of polyhedral cages. Typical clathrate hydrate structures are space-filling tilings of four types of polyhedra, called Kasper polyhedra, consisting only of five- and six-membered rings. Such structures are called “space-fullerenes”^{29} associated with the fact that Buckminsterfullerene structures are made of hexagons and pentagons. The positions of the cages in a clathrate hydrate structure are isomorphic to the atomic positions of the corresponding Frank-Kasper type alloy structure.^{33} This means that there are many potential structures for clathrate hydrates although only a few structures have been found experimentally.^{30} We evaluate all the 84 space-fullerene structures proposed by Dutour Sikirić,^{29} including the hypothetical ones, as candidates of low-density ice structures. There are a few structures that belong to both zeolite and hydrate groups. For example, clathrate hydrate structure sI (CS-I), which is a dual of the A15 Frank-Kasper structure, is also called zeolite MEP. Some of the clathrate hydrates, such as sVII (zeolitic ice SOD) and sH (zeolitic ice DOH), are not classified into space fullerenes because either of their hydrogen bond network contains rings other than pentagons or hexagons. Correspondence between several classifications is tabulated in Table S1 of the supplementary material. The space-fullerene ice structures are also generated by the GenIce tool.^{32}

The number of water molecules ranges from ∼1000 to ∼10 000 because both the size and shape of the unit cell are different for each ice structure. The TIP4P/2005 water force field model is employed.^{34} The temperature, *T*, is fixed at 77 K with a Nosé-Hoover thermostat in all MD simulations.^{35,36} The ice structures generated by GenIce are first equilibrated by a 0.1 ns isothermal-isochoric simulation. Then, an isothermal-isobaric MD simulation is performed at a pressure of *p* = 1 bar for 1.1 ns for further relaxation of the ice structure. The pressure is kept constant by a Berendsen barostat.^{37} The following 2 ns isothermal-isobaric MD simulation is used to calculate the potential energy, *E*, and the molar volume of the system, *V*.

## RESULTS AND DISCUSSION

There are 219 four-coordinated structures in the zeolite DB. However, seventeen of the 219 zeolitic ice structures quickly collapse even at 77 K because of the presence of unstable motifs such as three-membered rings of hydrogen bonds. Figure 1 plots the potential energies of the 202 of 219 zeolitic ices at 1 bar and 77 K against the molar volume. The volume of the zeolitic ices ranges from 20 to 35 cm^{3} mol^{−1}. There is a tendency that the potential energy is lower for denser structures. We find that the volume of ITT ice is higher than any other zeolitic ice, indicating that this ice occupies a region in the phase diagram of water under deeply negative pressure as the most stable phase because of the *pV* term in the Gibbs energy. Figure 1 also shows the relation between the potential energy and the volume for 75 space fullerene ices (nine out of 84 space fullerene ice structures quickly collapse at 77 K because they contain highly distorted five or six membered rings). The range of the volume of the space fullerene ices is quite narrow because the structural variety arises from the combination of only four cage types. Among the 84 space fullerene structures, C15, which is also known as sII in Jeffrey’s nomenclature or the MTN zeolite framework, is the most stable phase under negative pressure because of not only its lowest potential energy but also the highest volume. Thermodynamic properties of all the structures are listed in Table S2 of the supplementary material.

The structures of some zeolitic ices are shown in Fig. 1. Ice FAU, proposed originally by Huang *et al.*,^{12} consists of truncated octahedron (4^{6}6^{8}) connected by hexagonal prisms (4^{6}6^{2}). If the truncated octahedron and the prisms are regarded as vertices and edges, the topology of the super-network is identical to ice Ic. Ice RHO, which was also introduced by Huang *et al.*,^{11} is slightly denser than ice FAU. The structure of ice RHO is a simple-cubic lattice of truncated cuboctahedra (4^{12}6^{8}8^{6}) connected by octagonal prisms (4^{8}8^{2}). Ice ITT is the lowest density ice in Fig. 1. The potential energy of ITT ice is lower than those of ices FAU and RHO. Ice ITT is made of the stacking of double layers of water molecules with the octadecagonal channels in the triangular arrangement. They penetrate toward the c axis. In a double layer, cubes (4^{6}) connect 3-fold symmetric units. The 4^{6} cubic cluster is known as the most stable structure of the water 8-mer, and the 3-fold unit resembles one of the most stable structures of the water 20-mer (4^{9}5^{6}).^{38} The layers are connected by the 3-fold symmetric units with the 3-membered rings of hydrogen bonds.

We discuss the phase diagram of water at low temperatures and negative pressures. It is possible to determine roughly the most stable phase at a given pressure from Fig. 1 without further simulations. The Gibbs energy of ice phase *i*, *G*_{i}, is given by

where *U* and *S* are the internal energy and the entropy; the former is the sum of the mean kinetic and the potential energy, *E*. In general, the isothermal compressibility, *κ*, is very small for solid phases. We assume that *κ* = 0 for all the ice phases for simplicity. This enable us to replace *U*_{i}(*p*) and *V*_{i}(*p*) with *U*_{i}(*p* = 1 bar) and *V*_{i}(*p* = 1 bar). We also assume that the entropy is independent of ice structures (note that the effect of the *TS* term on *G* is expected to be small even when the entropy depends on ice structures because of the low temperature). Under these assumptions, the equilibrium pressure for ice *i* and ice *j*, at which the Gibbs energies of the two phases are identical, is expressed by

This means that the equilibrium pressure is approximated by the negative of the slope of the line segment connecting the two points corresponding to ices *i* and *j* in Fig. 1. The most stable ice is ice Ih at 0 bar. Figure 1 indicates that ice MTN (ice XVI) becomes more stable than ice Ih at −3.4 kbar. The slopes of lines connecting ice Ih and ice structures other than MTN are larger than that for the Ih–MTN line, indicating that the coexisting pressures of ice Ih and the other ices are lower than −3.4 kbar and thus they are metastable with respect to the two stable ice phases at this pressure, Ih and MTN. That is, the lower convex hull of the points in Fig. 1 indicates the most stable phases with decreasing pressure and the transition pressures. The phase transition pressure from MTN to RHO and that from MTN to FAU are *p* = −7.2 kbar and −5.0 kbar, respectively. These values are close to the values reported by Huang *et al.*, −5.2 kbar and −4.5 kbar at 77 K, in spite of our crude assumptions. Huang *et al.* suggested that the FAU structure is the most stable under deeply negative pressures considering various structures of 230 different space groups. However, Fig. 1 shows that the transitions from MTN to ITT occur at a higher pressure of −3.9 kbar. ITT is the most stable phase among the structures compared in Fig. 1 in the deeply negative pressure region, although the transition pressures might be varied from the present values when more rigorous methods, such as the quasi-harmonic technique,^{39} are employed to calculate the free energies.

There are octadecagonal channels in the framework of ITT, and they are likely to be able to accommodate carbon nanotubes of about 8 nm diameter, which can also contain a chain of water molecules inside.^{40} It might be possible to produce the ITT ice framework experimentally by freezing liquid water containing nanotubes. Some tricks may be required to prevent strong aggregation of the nanotubes. For example, confinement of the system in a narrow space, where attractive interactions between the wall and the nanotubes may help making a proper space to form the ice framework between the nanotubes. Other zeolitic ice structures may also be observed in confined spaces with or without guest molecules.

Is it possible to design stable low-density ice structures other than the zeolitic and space fullerene types? A clue is the structural feature of ices FAU, RHO, and ITT: all these structures have large pores surrounded by the super-network made of prismatic edges and polyhedral vertices. FAU, RHO, and ITT are made of hexagonal prismatic, octagonal prismatic, and cubic edges, respectively (Fig. 1). Such prismatic ice structures are mechanically stable at low temperatures.^{41} They are thermodynamically stable in carbon nanotubes,^{42,43} and in fact it is possible to produce them experimentally.^{44} The prismatic motifs are also often found as stable structures of water nanoclusters.^{38,45} By elongation of the length of the prismatic edges of FAU, RHO, or ITT, we can obtain less dense structures without changing the super-network topology. In the limit, it is possible to assume that all water molecules belong to the prisms, and therefore the potential energy of the ice structure can be approximated by that of the ice nanotubes.

Let us consider the case of ice FAU. The super-network structure of FAU, which is isomorphic to ice Ic, consists of hexagonal ice nanotubes. The hexagonal ice nanotube is one of the lowest energy structures confined in cylindrical pores.^{41} If all the prisms are extended, we obtain an aerogel-like structure^{46} which is still topologically isomorphic to ice Ic. We call it “aeroice” *n*xFAU, where *n* is the number of repetitions of hexagonal prisms at each super-edge. The aeroice structure 1xFAU coincides with the original ice FAU structure, and 0xFAU is identical to the empty clathrate hydrate sVII (i.e., zeolitic ice SOD in Table S1 of the supplementary material). We generate the aeroice structures of *n* = 2, 4, 8, 16, and 32. The hexagonal ice-nanotube is hydrogen-ordered.^{41} Figure 2 presents the relation between the volume and the potential energy of the aeroice structures at 1 bar and 77 K. The aeroice structures do not dissociate in MD simulations, indicating that they are at least mechanically stable at this condition. The molecular volume increases dramatically by elongation of the prisms, while the potential energy converges to −50.6 kJ/mol. Approximately, when *n* is doubled, the number of water molecules in the unit cell is also doubled, while the volume of the unit cell is octupled, and therefore the molecular volume is quadrupled. In the same manner as Fig. 1, the negative of the slope of the line connecting ice Ih and an aeroice indicates the equilibrium pressure between the two phases. The equilibrium pressure is higher for the aeroice structures than for ice XVI (MTN) even when *n* = 2. The least dense aeroice we evaluate, 32xFAU, exhibits the equilibrium pressure of only −0.05 kbar. The equilibrium pressure increases with increasing *n*, and it converges to zero. A similar result would be obtained when we select ice RHO or ITT as the base structure of aeroice instead of ice FAU. We can generate arbitrarily light aeroices. For example, 100xFAU is as light as ambient air.

At *T* = 0, the most stable negative pressure ice phase is neither MTN nor ITT but is the aeroice structure with *n* > 1. The sparse aeroice structures suggest that they seem to dissociate very easily due to thermal fluctuations. In other words, the sparse aeroice structure is entropically unfavorable. Aeroices of larger *n* would be less stable and occupy regions of lower temperatures in the phase diagram at a given negative pressure. The phase diagram under negative pressure would be tessellated into small regions of aeroices with different *n*, and the aeroice phase of *n* = 1, i.e., zeolitic ice FAU, is only one of them. It might be possible to draw the phase diagram of ices including aeroice structures using molecular simulations; however, it requires a huge computational resource due to the extreme variety of aeroices (various *n* values and various zeolite frameworks which can be used for the base structure of aeroices). Therefore, we do not determine coexisting lines between negative pressure ices in the present study.

## CONCLUSIONS

We have examined three types of low-density ices, space fullerene ices, zeolitic ices, and aeroices. The structures of space fullerene ices are the same as real or hypothetical empty clathrate hydrates. It is found that the empty sII hydrate is more stable than any other space fullerene ice. The zeolitic ice structures are prepared from SiO_{2} zeolite structures given in the zeolite database. Ice ITT, which is an analogue to zeolite ITT, is more stable than ice FAU which was proposed as the most stable phase under deeply negative pressures in a previous paper.^{12} The third type of low-density ices, aeroices, are obtained from a type of zeolitic ice that consists of prismatic edges and polyhedral vertices by elongation of the prismatic edges. The aeroice structures do not collapse in MD simulations at *T* = 77 K and *p* = 1 bar, indicating that they are at least mechanically stable at the liquid nitrogen temperature. Aeroices can be more stable than any zeolitic ices at certain thermodynamic conditions under negative pressure. As the prismatic edges get longer, the density of the aeroice phase approaches zero while the potential energy converges to a certain negative value. As a result, the transition pressure from ice Ih to the aeroice phase becomes a quite small negative value.

## SUPPLEMENTARY MATERIAL

See supplementary material for the properties of all the structures in this paper, the method to interpolate the properties of aeroices, and the correlations between the potential energy and structural properties.

## ACKNOWLEDGMENTS

The present work was supported by JSPS KAKENHI Grant Nos. JP16K17857 and JP16K05658 and MEXT as “Priority Issue on Post-Kcomputer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use).