Formation of hot ice caused by carbon nanobrushes

Although water is a simple molecule, there are a large number of ice polymorphs: 18 ice phases have been identified experimentally.^1,2 There also exist many clathrate hydrate and filled ice phases.³ Discovery of new ice and clathrate hydrates is a challenging problem.^4–10 A route to find a new stable crystalline phase is computation of the Gibbs energies of possible structures.^11–14 The existence of ice XVI was predicted by this method prior to the experimental synthesis which was achieved by vacuum pumping of CS-II neon hydrate.⁶ 

A silica zeolite structure is expected to have the corresponding ice structure because the coordination number 4 is common to ice and zeolite.^12,15–20 In previous studies, we evaluated the stability of such zeolite-based ice structures to find candidates of novel low-density ice and found that several of them are quite stable under negative pressure.^14,17 ITT ice, which is shown in Fig. S1, is one of the stable zeolite-based ices. There are channels in this ice. Zeolite ITT has the largest channel in the zeolite structures collected in the Database of Zeolite Structure web site which provides structural information on all of the zeolite framework types that have been approved by the Structure Commission of the International Zeolite Association.²¹ The radius of the channels is suitable to accommodate a carbon nanotube (CNT) with a chirality of (6,6). This suggests that the ice ITT structure, or similar ones, may form with the aid of the CNTs. The resultant structure can be considered as a hydrate of the CNT, and it is in contrast to tubular ices inside CNTs.²² 

In this study, we perform all-atom molecular dynamics (MD) simulations of water in the presence of CNTs. The CNTs are fixed in the simulation cell with a triangular arrangement so that they match the ITT structure. This configuration is a model of a nanobrush made of CNTs. It is shown that a novel ice crystal, which is slightly different from ITT, indeed forms spontaneously when there is a space that can accommodate two layers of water molecules between nanotubes. The melting temperature of this crystal is much higher than that of ice Ih. A theoretical framework is proposed to estimate thermodynamic properties of the new crystal in which the chemical potential of the crystal is calculated using the quasiharmonic approximation.

The GROMACS 4.6 package is used for MD simulations.^23,24 The temperature is kept by the Nosé-Hoover method.^25,26 The particle mesh Ewald (PME) method is employed with a real space cutoff distance of 0.9 nm.^27,28

We employ the TIP4P/2005 model for water. This model well reproduces the phase diagram of ice phases unless the pressure is extremely high.^29,30 A recent MD study suggested that the TIP4P/2005 model also reproduces the structure of water under 2D confinement.³¹ The intermolecular Lennard-Jones (LJ) and intramolecular (bond, angle, and dihedral) parameters of CNT are taken from the OPLS-AA model (Table S1).^32,33 The standard geometric mean rule is used for unlike pair interactions.

Figure 1 shows a snapshot of MD simulations. The system consists of three types of regions: (1) open-ended CNTs with length L and water between and inside the CNTs, (2) bulk liquid water with a width of 3 nm, and (3) a vapor region with a width of 6 nm. The number of CNTs is either four or sixteen. To fix each CNT at a desired position in the simulation cell, four carbon atoms located at an edge of the CNT and other four carbon atoms at the other edge are selected and their positions are constrained using a harmonic potential with a force constant of 20 MJ mol⁻¹ nm⁻². Table S2 summarizes the system size and the number of water molecules in the cell. MD simulations are performed without barostats to keep the gap between the CNTs, d, constant. The vapor region is required to maintain the pressure at ∼0 bar in the constant volume simulations.

Rough hydrophobic surfaces, including carbon nanobrushes, can be water-repellent.³⁴ However, it is also known that water molecules can flow through narrow hydrophobic spaces such as the inside of CNTs.³⁵ We examine whether water molecules penetrate into the space between the CNTs spontaneously in preliminary MD simulations in which water molecules are absent in region 1 in the initial configuration. It is found that water molecules enter region 1 unless the gap is too narrow (d < 0.7 nm) when the temperature is 280 K and the tube length is L = 6 nm. Therefore, region 1 is prefilled with water molecules in the production MD simulations for d ≥ 0.7. Such an initial configuration is prepared by (i) generation of a random liquid configuration in the simulation cell, (ii) elongation of the cell to make a vacuum region, (iii) relaxation of liquid water using an isothermal MD simulation for 1 ns at 280 K, and (iv) insertion of CNTs in the liquid region with removing water molecules to avoid overlapping.

Water molecules are classified into immobile solidlike and mobile liquidlike molecules. A trajectory is divided into time bins with a width of 0.1 ns, and $δi2=rit+Δt−rit2$ is calculated for each bin, where r_i is the coordinate vector of the oxygen atom of the ith water molecule.^30,36 We set $Δt$ = 0.01 ns and define that water molecule with $δi2<0.01 nm2$ as immobile molecules. A water molecule is considered to be hydrogen bonded with a different water molecule when the O–O distance is less than 0.35 nm and the angle between the O–O vector and one of the O–H bonds is less than 30°.³⁷ We define that two immobile water molecules belong to the same cluster when they form a hydrogen bond.

Figure 2(a) shows the time evolution of the potential energy for the large system (16 CNTs). The length of the CNTs is set to L = 6 nm, and the gap between the CNTs is d = 0.92 nm. The potential energy decreases sharply after an induction time of 17 ns at T = 290 K. Snapshots of this trajectory are shown in Figs. 3(a)–3(d) (Multimedia view). Clusters of immobile water molecules form in region 1, but they dissociate quickly. The number of molecules in the largest cluster does not exceed 200 during the induction period [Fig. 3(e)]. A large cluster consisting of ∼300 molecules forms at t ∼ 17 ns. This cluster grows as time evolves, and eventually, region 1 is fully filled with immobile water molecules at t ∼ 25 ns. This behavior indicates that this is a first order transition from the liquid phase via nucleation.

The network structure of the crystal which emerged at T = 290 K is identified as the same as that of a hypothetical structure coded dtc,^14,38 which has not been observed as either a zeolite or a hydrate crystal in reality. The same structure forms more rapidly at T = 288 and 286 K [Fig. 2(a)]. We perform three additional MD simulations at T = 286, 288, 290, and 292 K starting from different initial liquid configurations and confirm that the dtc crystal forms spontaneously in all the simulations at T = 286, 288, and 290 K (Fig. S2).

The dtc structure, which is quite similar to the ITT structure, is shown in Fig. 3(f) and Fig. S1. There are 4-, 5-, and 6-membered rings in the structure. The channel with 18 water molecules on the circumference is made of regular tiling of hexagons and is isomorphic to (9,0) CNT. All water molecules in this structure are four-connected and obey the Bernal-Fowler-Pauling ice rules.^39,40 The obtained solid can be considered as a hydrate crystal of the CNTs and a kind of composite crystal with an incommensurate structure. We call it dtc CNT hydrate. In the ideal crystal of infinitely long (6,6) CNTs and the dtc water structure, the ratio of the number of water molecules to the number of carbon atoms is 1/1.8 (excluding the water molecules inside CNTs).

The empty CS-II hydrate crystal, also known as ice XVI, can be synthesized experimentally from CS-II Ne hydrate by emptying the guest molecules.⁶ The melting point of empty CS-II hydrate is lower than that of ice I under ambient pressure. However, ice XVI becomes more stable than ice I under negative pressure because of its large molar volume.¹¹ Likewise, the dtc ice crystal without CNTs (empty dtc hydrate) occupies a region in the phase diagram of ice under negative pressure as the most stable phase.¹⁴ 

The onset temperature of crystallization, 290 K, should be lower than the solid/liquid coexistence temperature because the nucleation time at the coexistence temperature is longer than the simulation time due to a very low thermodynamic driving force and the presence of the free energy barrier that is much higher than the thermal energy. We perform MD simulations starting from a dtc structure to determine the melting temperature that is defined as the average between the highest temperature at which the dtc structure persists for the simulation time of 50 ns and the lowest temperature at which the crystal melts completely. This melting temperature is the same as the coexistence temperature because the phase transition occurs without nucleation due to the presence of the interface of the solid in region 1 and liquid water in region 2. The potential energy shown in Fig. 2(b) indicates that the coexistence temperature is ∼303 K. This is much higher than the melting point of ice Ih for the employed force field model, ∼250 K.²⁹ It is known that the dissociation temperatures of semiclathrate hydrates and other water-rich hydrate crystals can be higher than the ice point under ambient pressure. Tetrabutylammonium fluoride (TBAF) hydrate is probably the most stable semiclathrate hydrate. However, the difference in the melting temperature between TBAF hydrate and ice I, ∼30 K,⁴¹ is much smaller than the corresponding value for the dtc crystal with CNTs, ∼50 K. Therefore, dtc CNT hydrate can be called “hot ice” (the term hot ice often means sodium acetate hydrate that resembles ice I but melts at a much higher temperature).

Elevation of the melting point caused by confinement has been observed for water in carbon nanotube and between two graphene sheets where the van der Waals interaction with the carbon sheet stabilizes the water structure.^42–46 A substantial difference from these systems is that dtc is neither 1D nor 2D crystal but a 3D crystal. Another large difference is the small number of unstable four-membered rings:¹⁴ the ratio of the number of four-membered rings to that of water molecules for dtc, 1D ice in nanotube, and 2D ice between graphene sheets is 30%, 100%, and 75%, respectively. This is a reason for the extremely high melting temperature of dtc CNT hydrate.

As shown in Fig. S3, the melting temperature of the 4-CNT system (L_x × L_y × L_z = 3.468 × 3.003 × 18.000 nm) is the same as that of the larger 16-CNT system (6.936 × 6.007 × 18.000 nm). Hereafter, we present the results obtained from MD simulations of the smaller system.

Figure 4(a) shows the melting temperature of the solid plotted against the length of the CNTs, L. The gap between the CNTs is d = 0.92 nm. The melting temperature increases with increasing L and almost levels off at L ∼ 4 nm. As shown in Fig. S4, there exist fluctuations of the solid/liquid interface. The length scale of the fluctuations reaches 1.5 nm at 300 K. Thus, the dtc crystal melts more easily with shorter L. The melting temperature is also dependent on the gap d. Figure 4(b) shows that the melting temperature is maximized at d = 0.94 nm.

The CNTs are fixed in the simulation cell using a harmonic potential. Does dtc CNT hydrate spontaneously form without such constraints? In other words, is it possible to synthesize the crystal by simply cooling the aqueous solutions of the CNT? To assess this possibility, we calculate the free energy change using the gap d as the reaction coordinate. The increase in d involves the increase in the cell dimensions along the x and y directions (the CNTs are parallel to the z axis). The resultant free energy change can be expressed as

where p_ii is the diagonal element of the pressure tensor. Both p_xx and p_yy are nonzero even when d is very large because of the presence of the water/vapor interfaces.⁴⁷ To eliminate this contribution, we calculate the difference between w with and without the CNTs,

Figure 4(c) shows the free energy profile Δw plotted against d at T = 280 K and 320 K (the diagonal elements of the pressure tensor at 280 K are shown in Fig. S6). The tube length is L = 6 nm. There are two minima at d = 0.34 and 0.94 nm at both temperatures. The minimum of Δw at d = 0.94 nm corresponds to the maximum of the melting temperature shown in Fig. 4(b). The CNTs are in contact with each other at d = 0.34 nm. The free energy at d = 0.94 nm is much higher than that at d = 0.34 nm, indicating that aggregation of the CNTs occurs in the absence of the constraints. In order to synthesize the dtc structure experimentally, methods to prevent the aggregation, such as fixing short CNTs on a base or placing some spacers between CNTs to keep a regular separation, are required.

Let us consider the equilibrium between dtc CNT hydrate and bulk water at temperature T and (external) pressure p₀ by the quasiharmonic approximation. The method for calculation of the Helmholtz free energy is closely related to what have been applied to clathrate hydrates.⁴⁸ It provides an insight into the stable temperature ranges of dtc CNT hydrate being in equilibrium with bulk liquid water and also a rough estimation of the axial cell dimensions, l, and the normal pressures, p_n.

In order to modify the method applicable to the equilibrium, we first set up the thermodynamic restrictions specific to the present system. Water molecules are confined in a rectangular cell in which the CNTs of length L are placed. Unlike the MD simulations, the length of CNTs is the same as the cell dimension, l. Periodic boundary conditions are applied in all the directions. The number of the CNTs is fixed to an integer m that does not violate the crystallographic condition for the dtc structure. The water molecules inside the cell are assumed to be in equilibrium with bulk water. This equilibrium is guaranteed by equivalence between the chemical potential of water in the cell and that of pure liquid water or ice. Instead of variation of the number of molecules in the cell, the cell dimension is subject to change. The Helmholtz free energy is a function of temperature, T, the number of water molecules in the cell, N_w, the axial cell dimension, l, and the area of the plane perpendicular to the nanobrush per nanotube, a (the area is determined by the gap d). The cell volume is given by lma. The normal pressure exerting on the planes perpendicular to the CNTs, p_n, is calculated as

Similarly, the lateral pressure, p_l (in the radial direction), exerting on the planes parallel to the CNTs is given as

The thermodynamic potential can be expressed by

For the new thermodynamic potential, we have

and

where μ is the chemical potential of water in the dtc structure. This provides a way to calculate μ as

The equilibrium condition is that the chemical potential of dtc CNT hydrate, μ, is equal to the chemical potential of bulk water, μ = μ₀(T, p₀). That is, we solve the equation

combined with Eq. (3) to find an appropriate $l$ value. In any stable system for a given $l$⁠, the normal pressure calculated by Eq. (9) must be positive and the following inequality should be satisfied:

Whenever the above inequality holds, dtc CNT hydrate could be in equilibrium with bulk water. There must be, however, a temperature above which a positive p_n cannot be obtained and the dtc crystal becomes metastable. It is the upper limit of the melting temperature of dtc CNT hydrate, T₀. Practically, the temperature, T₀, is calculated in reference to Eqs. (9) and (10) as

where l₀ is the optimum l to minimize A(T₀, l, a, N_w) under given (T₀, a, N_w), i.e., p_n = 0. Below T₀, dtc CNT hydrate can be in equilibrium with bulk water by raising the normal pressure, p_n. There must be a temperature at which the p_n value agrees with that of liquid water inside the cell. It is the melting point of dtc CNT hydrate, which must be lower than T₀. Exploration of the melting point via the calculation of the chemical potential of liquid water in the cell is beyond the scope of the present study, but it is obtained from the present MD simulations.

It is noted that even the normal pressure, p_n, can be different from the bulk one, p₀. The discrepancy arises from the interfacial tension between water molecules and CNTs, which depends on the morphology of water and the effective surface area. Alternately, the system can be regarded as a solution of the CNT where water inside is separated from bulk water by a semipermeable membrane that inhibits transfer of the solute. It gives rise to a kind of osmotic pressure and causes a difference between internal and external pressure. In our theoretical method, such grandcanonical view is switched to the isothermal-isobaric one to handle crystalline states as described above.

The Helmholtz free energy at T and V for dtc CNT hydrate of N_w water molecules is evaluated in terms of the energy at quenched structures, U_q (l, a, N_w), the vibrational free energy, F_v (T, l, a, N_w), and the residual entropy, S_r (N_w), as

Here, we assume that dtc CNT hydrate is fully hydrogen-disordered. According to Pauling’s approximation,³⁹ S_r (N_w) is given as

where k_B is the Boltzmann constant. If the anharmonic vibrational free energy can be neglected, the vibrational free energy is approximated as

where h is the Planck constant, ν_j stands for the jth normal mode frequency, and $…$ indicates the average over the generated hydrogen-disordered configurations. The number of water molecules is fixed to 640, and we generate 10 hydrogen-disordered structures using the GenIce tool.⁴⁹ 

Calculation of the chemical potential of bulk water is required to consider the equilibrium. The bulk region is ice Ih when the temperature is lower than the melting point, T_m0. The chemical potential of bulk ice is calculated similarly from the Helmholtz free energy, A(T, V, N_w), by changing the independent variable from the volume, V, to the bulk pressure, p₀.

The chemical potential of bulk liquid water is obtained as follows. The chemical potential of liquid water, μ_l (T, p₀), at temperature T (=T_m0 + ∆T) can be estimated roughly as

from the chemical potential of ice, μ_i(T_m0, p₀), the entropy difference between ice and liquid water, ∆s_m, and the heat capacity difference at constant pressure, ∆C_p. The entropy difference is obtained from the heat of fusion, ∆H_m, as ∆s_m = ∆H_m/T_m0 and the heat capacity difference is simply calculated from the enthalpies of both ice Ih and liquid water by performing MD simulations at several temperatures close to T_m0. Ice Ih of TIP4P/2005 water is known to melt at 252.1 K.²⁹ The anharmonic vibrational free energies may not be negligible. However, we can expect that the magnitudes of the anharmonic free energies in ice Ih and the dtc crystal (without CNTs) are hardly different from each other. A more extensive discussion on the legitimacy of ignoring the anharmonic free energy is given elsewhere.⁴⁸ It is also noted that in the harsh confinement as imposed by the CNTs of this type, the anharmonic free energy is rather small compared with the bulk phase and the chemical potential evaluated by the above method may lead to a little lower value.

Figure 5 plots the chemical potentials at p_n = 0, equivalently the minimum values of the Helmholtz free energies per mole of water for dtc CNT hydrate, A(T, l₀, a, N_w), at various nearest neighbor distances, d (rather than a), as a function of temperature. The chemical potential of bulk water at p₀ = 0.1 MPa is also drawn along with that of dtc ice without CNTs. The chemical potential of dtc ice without CNTs is higher than that of liquid water by 3–5 kJ mol⁻¹ depending on the temperature. The existence of CNTs affects the chemical potential of water in two opposite ways. The CNTs shift the vibrational frequencies of water molecules to higher side, which increases the chemical potential.⁵⁰ They act as an external force field and decrease the chemical potential. The latter contribution is much larger, and the resultant chemical potential at p_n = 0 with CNTs (solid lines) is much lower than that without CNTs (magenta dotted line).

In Fig. 5, the temperature of the intersection of the chemical potential curve of bulk water (light blue dotted line) and that of dtc CNT hydrate (solid lines), T₀, increases with d. Once it reaches the optimal size, i.e., d = 0.94 nm, the temperature gradually decreases as the distance increases. The optimum gap distance from this calculation agrees well with that from the MD simulations. Each intersection corresponds to the upper limit of the melting point for a given d. Those are higher by about 35 K than the corresponding melting points obtained from the MD simulations. This indicates that the normal pressure, p_n, at the coexisting state of liquid water and dtc CNT hydrate are fairly high.

The axial interlayer distance in the dtc crystal, which is proportional to l, is depicted in Fig. 6 against T for several gap values. The corresponding values obtained from the MD simulations at 260 K are also marked. As the gap decreases, the dtc crystal expands axially. This occurs due to squeezing of water molecules by the lateral compression. The trend obtained from the theoretical calculation agrees, although only semiquantitatively, with the MD simulations.

The temperature dependence of the normal pressure is illustrated in Fig. 7. As is expected from the above argument, the normal pressure increases with decreasing temperature. The slope becomes gentler with the gap size. The pressure is higher for a gap having higher T₀, which is also anticipated from Eq. (8). It is noted that the axial length in the statistical mechanical calculation is approximately l = 3 nm, but the surface effect normal to the nanobrush is removed owing to the periodic boundary condition applied to the nanobrush-dtc ice complex without bulk water in the statistical mechanical calculations. The discrepancy in the distance between the MD simulations and the statistical mechanical calculations seems to arise in part from the finite axial length in the MD simulation and the surface areas normal to the nanobrush, which reduce the normal pressure.

We have investigated the phase behavior of water in the presence of CNTs. The CNTs are fixed in the simulation cell so that they can stabilize the ice lattice structure that has the ITT zeolite topology. It is found that a crystalline structure forms spontaneously in MD simulations although it is different from ITT. The topology of the crystalline structure is the same as that of a hypothetical zeolite structure, dtc. This crystal, we call dtc CNT hydrate, forms via nucleation and growth, and a hysteresis is found between the freezing and melting processes, suggesting that this is a first-order phase transition. The melting temperature of dtc CNT hydrate can be tuned by the tube length and the gap between the tubes. The free energy of the CNTs with water molecules forming the dtc structure is much higher than that of the CNTs that are in contact with each other. This indicates that the aqueous solutions of CNTs do not form dtc CNT hydrate by cooling because of aggregation of the CNTs. However, there is a possibility that dtc CNT hydrate is experimentally realized by the use of nanobrush made of short CNTs.

Thermodynamic properties of dtc CNT hydrate are examined using the quasiharmonic approximation considering the equilibrium between dtc CNT hydrate and bulk liquid water or ice. This provides an insight into the stable temperature ranges of dtc CNT hydrate and a rough estimation of the axial cell dimensions, l, and the normal pressures, p_n. The chemical potential of dtc ice without CNTs is higher than that of liquid water by 3–5 kJ mol⁻¹ depending on the temperature. However, the interactions between water molecules and CNTs reduce the chemical potential of water so that the equilibrium with bulk water can be observed. We can estimate the upper bound of the melting temperature of dtc CNT hydrate and the most favorable gap size, which agree reasonably with what we obtain from the MD simulations. It is also found that the interlayer distance decreases with decreasing temperature, which is due mainly to the abrupt increase in the normal pressure rather than simple thermal expansion.

In a recent paper, Kumar et al. suggested ways to find a new guest-free zeolitic ice: (1) incorporation of attractive guests, (2) application of pressure, and (3) relief of the strain in the water network via incorporation of other hydrogen-bonded molecules or ions.⁵¹ Our approach can be classified in the first category. A solid surface can stabilize water and promote formation of zeolitic ices. Several recent studies showed that the rate of ice nucleation can be controlled by the geometry of the solid surface.^52–60 Our result suggests that the structure and the melting point of ice can also be controlled by the solid surface. Phase transitions of water have been used for separation processes and thermal energy storage.^61,62 The nanobrush may be useful for this purpose because the melting point of ice can be controlled so that it is close to the room temperature.

We have focused only on the CNT of chirality (6,6) which can be included in the 18-ring channel. In a previous paper, we proposed a new type of ultralow-density ice named aeroice, which is expected to be one of the most stable ice phases under certain thermodynamic conditions.¹⁴ The aeroice structure is made by extending the prismatic part of zeolitic ice.¹⁷ It is possible to design a dtc-like ice structure having larger cylindrical channels by extending the bilayer wall part of the dtc structure in a way similar to the generation of the aeroice structures. There also exist zeolitelike network frameworks having larger channels such as bzd, dtd, and gcd.³⁸ These structures, or unexpected ones, may form in MD simulations in the presence of CNTs with larger radii.

The supplementary material contains two tables summarizing the potential parameters used in MD simulations and the number of molecules in typical simulation cells and six images showing the ITT and dtc structures, potential energies of four independent MD simulations, the potential energy of the small 4-CNT system, fluctuation of the dtc/water interface, the tetrahedrality parameter, and the pressure anisotropy.

The present work was supported by a grant of the MORINO FOUNDATION FOR MOLECULAR SCIENCE, JSPS KAKENHI, Grant No. 16K05658, 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 No. hp180204). MD simulations were also performed on the computers at the Research Center for Computational Science, Okazaki, Japan.