Two-dimensional partitioned square ice confined in graphene/graphite nanocapillaries

The initial configuration of the MD simulation is shown in Figure 1(a), and is similar to that used in Ref. Algara-Siller et al. (2015). A long capillary is formed by two parallel graphene-based sheets consisting of different numbers of graphene monolayers. We use the most common graphene stacking ABA in our multilayer simulation systems. (We have also tested the AAA stacking for comparison to find out how the stacking affects the ice patterns. The result shows that square ice is formed independently of whether the graphene stacking is AAA or ABA. See Figure S1 in Supplementary Material for details.) Two water reservoirs with each containing 2000 water molecules are connected by the capillary. The number of graphene monolayers on each side of the water molecules is varied between $1\leq n\leq 6$. The length and the width (in $z$ and $x$ direction as shown in Figure 1(a)) of the graphene capillary are fixed at 68 Å and 56 Å, respectively. The height $h$ of the graphene capillary (i.e., the distance between the two graphene sheets in $y$ direction) is chosen to be 6.5 Å or 9.0 Å to accommodate one or two layers of water molecules, respectively. The carbon atoms in the graphene sheets are fixed during simulations, as previous studies have shown that square ice is formed independently of whether the confinement is provided by rigid or flexible graphene sheets.Algara-Siller et al. (2015); Ruiz Pestana, Felberg, and Head-Gordon (2018)

MD simulations are performed in the isothermal, isobaric ensemble (NPT), in which the temperature (298 K) and the pressure are controlled by the Nosé-Hoover thermostat and barostat, respectively. Periodic boundary conditions are applied in the three directions of the simulation box. In the equilibration run, the lateral pressure $P$ is kept at 1 bar ($10^{-4}$ GPa) for the first 5 ns, during which water molecules fill in the graphene capillary. After that, $P$ is increased to 1.1 GPa (critical lateral pressure under which the square icing patterns can form, as shown in Figure S2 in Supplementary Material) and maintain for more than 25 ns, considering the extra lateral forces acting on water molecules in experiments.Algara-Siller et al. (2015) The temperature is fixed at 298 K during the simulations. A time step of 1.0 fs is used for the velocity-Verlet integrator. Water molecules are modeled using the extended simple point charge (SPC/E) model, including the long-range Coulomb potential and the short-range Lennard-Jones potential between the interaction sites.Berendsen, Grigera, and Straatsma (1987) The SPC/E water model has been proved to be accurate in the simulations of monolayer and bilayer square ice structure by comparing with experimental results.Algara-Siller et al. (2015) The water-carbon interaction is modeled by a 12-6 Lennard-Jones potential between the carbon and the oxygen atoms. The distance and the energy parameters are $\sigma=0.328$ nm and $\epsilon=0.114$ Kcal/mol,Che and Theodorakis (2017); Gordillo and Martí (2000); Joly (2011); Werder et al. (2008) respectively, and the distance between graphene monolayers is 3.4 Å in graphite sheets. The long-range interactions are computed using the particle-particle particle-mesh (PPPM) algorithm with an accuracy of 10^-4. The simulations are carried out by using LAMMPS,Plimpton (1995) and the snapshots are rendered in VMD.Humphrey, Dalke, and Schulten (1996)

To quantitatively determine and distinguish the partitioned square icing pattern from the homogeneous square icing pattern and the liquid water, square icing parameters, $MCV_{1}$ and $MCV_{2}$, are proposed in this study:

where $L$ and $M$ are two free positive integers, $l$ is an integer that runs from 1 to $L$, $m$ is an integer that runs from 1 to $M$, and $A_{lm}$ is the number of water molecules whose coordinate and orientation satisfy:

where $z$ is the coordinate of the water molecule in $z$ direction, $Z_{\text{min}}$ and $Z_{\text{max}}$ are the coordinate minimum and maximum of the graphene capillary in $z$ direction, respectively, and $\theta$ represents the direction of the hydrogen-oxygen bonds in the water molecule, which runs from $\theta=0$ to $\theta=\pi/2$. The overbar in Eqs. (1) and (2) indicates the average, i.e.,

These square icing parameters analyze the uniformity of the distribution of water molecules in spatial position and direction of the hydrogen-oxygen bonds, so as to be able to identify the different configurations of water molecules. From the data of many MD simulations (all of our simulation results including, but are not limited to, the simulation cases presented in this article), a map of the square icing parameters is generated, as shown in Figure 1(b), in which $L=13$ and $M=18$ are chosen based on trials. The partitioned square ice (partitioned SI) can be easily distinguished from the homogeneous square ice (homogeneous SI) and liquid water when we take $MCV_{1}=0.4$ and $MCV_{2}=0.6$ as thresholds.

We carried out a series of MD simulations of water confined between two parallel graphene-based sheets, which consisted of different numbers of graphene monolayers. When there is only one graphene monolayer in the sheet, the confined water molecules form homogeneous square icing patterns in tune with the results revealed in Ref.Algara-Siller et al. (2015). From the side and the top views of the simulation results, shown in Figure 2(a) and (b), we can see the square icing patterns with hydrogen bonds in two perpendicular directions, and the water molecules in the top and the bottom layers form exactly the same structure. When we add another graphene monolayer, i.e., changing the number of graphene monolayers to two, the water molecules between them exhibit partitioned square icing patterns, as shown in Figure 2(c) and (d). We can observe two distinct regions in which the water molecules both formed obvious square icing patterns but with different directions of hydrogen bonds. Both of the structures of ice grains in the homogeneous and partitioned square patterns (except for grain boundaries and point defects) satisfy the ice rule, i.e., every water molecule is hydrogen-bonded to its four nearest intralayer neighbors, with every water molecule being a double donor and a double acceptor of hydrogen bonds, as shown in Figure 2. Most of the hydrogen bonds have switched to the in-plane configuration for the formation of the square ice structure. There are also some hydrogen atoms that lie between the two layers, which can form hydrogen bonds with oxygen atoms in adjacent layers and connect the two ice layers. The orientations of water molecules at the grain boundaries between the two regions are random, indicating that the partitioned square ice is a more unstable crystalline state with higher energy than the homogeneous square ice.

We also performed another series of MD simulations to consider the possibility of conversion from partitioned structure to homogenous square patterns by changing the pressure and the temperature. The lateral pressure $P$ was increased with a step of 0.1 GPa in the pressurization process, and the temperature $T$ was increased with a step of the 5 K in the heating process, respectively. The variation of the potential energy of the confined water during the pressurization and heating process is shown in Figure 3(a) and (b), respectively. (We also performed simulations with different pressurization rates and heating rates, which show similar results to that in Figure 3(a) and (b). See Figure S3 in Supplementary Material for details.) As the lateral pressure increases, water can transform from the liquid phase to the square ice structure when the pressure exceeds 1.0 GPa, and the square ice will decompose into amorphous phase if the pressure is approximately above 4.0 GPa, as shown in Figure 3(a). As the temperature increases, the partitioned SI structure will transform into the homogeneous SI structure at around 320 K, and the partitioned SI structure will melt into liquid phase eventually when the temperature exceeds 420 K, as shown in Figure 3(b). In some cases, the transformation from the partitioned SI structure to the homogeneous SI structure occurs in the pressurization process, as shown in Figure S2 and Movie S1 in Supplementary Material.

Hysteresis simulations are also performed through the pressurization/depressurization process by increasing the lateral pressure until 2.0 GPa and then decreasing the lateral pressure back to 1 bar ($10^{-4}$ GPa), as well as through the heating/cooling process by increasing the temperature until 550 K and then decreasing the temperature back to 300 K. The pressurization/depressurization process is similar to the cooling/heating process of the confined water, and both exhibit a large hysteresis loop, as shown in Figure 3(c) and (d), implying that the formation of the square ice is a first-order phase transition. It is worth noting that the homogeneous SI structure does not transform into the partitioned SI structure during the cooling process.

As discussed in the previous section, the number of graphene monolayers has a strong effect on the formation of partitioned square ice. To investigate the effect of the number of layers, we perform simulations with more graphene monolayers ($n=3-6$), and the icing parametric diagram is shown in Figure 4. Each case is repeated ten times by changing the initial velocity distributions of the water molecules while keeping the other settings identical (the snapshots of the ten repeated simulations for $n=2$ are shown in Figure S3). By comparing the distributions of the icing parameters for $n=1$ and $n=2$, we find that it is energetically unfavorable for the water molecules confined between two parallel graphene sheets ($n=1$) to form partitioned square icing patterns. In contrast, for the graphite ($n=2$), the partitioned square icing patterns become common. The key factor is the different effects of graphite and graphene monolayer on water molecules, caused by the short-range nature of the interaction between water molecules and the carbon substrate. The diagram of the icing parameters in Figure 4 also shows that $n=2$ does not always lead to the formation of the partitioned square icing patterns, and there is randomness in the formation of the partitioned structure.

From the distributions of icing parameters for $n$ = 2 – 6 in Figure 4, we can see that the partitioned square icing phenomenon is insensitive to the number of graphene monolayers as long as $n\geq 2$. This result can be explained by the short-range nature of the interaction between water molecules and the carbon substrate. In the force field that we use,Werder et al. (2008) the water-carbon interaction is modeled by a 12-6 Lennard-Jones potential between the carbon and the oxygen atoms. The distance and the energy parameters are $\sigma=0.328$ nm and $\epsilon=0.114$ kcal/mol, respectively, and the distance between the graphene monolayers is 3.4 Å in the graphite sheets. These force field parameters have been validated in many experiments and used in many MD simulations.Che and Theodorakis (2017); Gordillo and Martí (2000); Joly (2011); Werder et al. (2008) According to these force fields, the second layer of carbon atoms exerts a strong force on the water molecules between the graphite sheets. However, from the third layer of carbon atoms, the force exerted on the water molecules will quickly vanish. Therefore, the short-range nature of the water-carbon interactions determines the strong difference between $n=1$ and $n\geq 2$, as well as the insensitivity to $n$ when $n\geq 2$.

The variation of the potential energy of the confined water during the icing process for a typical case is shown in Figure 5(a). After the first 5 ns, the lateral pressure $P$ is increased to 1.1 GPa and crystalline structures of water molecules begin to form. Even though the water molecules gain some extra energy from the high lateral pressure, more energy is released during the formation of hydrogen bonds, and the potential energy of water molecules drops rapidly until the system gradually stabilizes. The process of square ice formation lasts about $5-10$ ns in our simulation, including an initialization stage and a growth stage as shown in Figure 5(a). In the initialization stage, water molecules change from a disordered liquid state to the locally ordered structure, while in the growth stage the locally ordered structure gradually expands until filling the whole space in the capillary. After that, the icing patterns remain stable eventually.

The variation in the potential energy shows different features for the formation of different ice structures, as shown in Figure 5(b). All cases are taken into account in Figure 5(b), and the cases with the same number of graphene monolayers and homogeneous/partitioned square icing patterns are grouped, and their potential energy is averaged per nanosecond. In the first 5 ns, the lateral pressure $P$ is kept at 1 bar ($10^{-4}$ GPa) for the equilibration. Water molecules confined between two parallel graphene sheets ($n=1$) have significantly higher potential energy than those confined between graphite sheets ($n\geq 2$). According to the force field, the second layer of carbon atoms exerts a strong attractive force on the water molecules between the graphite sheets, which is different from the repulsion caused by the first layer of carbon atoms. As a result, water molecules gain higher potential energy from the carbon sheets when $n=1$ at a fixed temperature. After the pressurization at 5 ns, water molecules gain extra energy from the high lateral pressure and, in the meantime, they release some energy to form hydrogen bonds, due to the formation of square ice. When the system gradually stabilizes after 20 ns, the water molecules for $n=1$ and the water molecules with homogeneous square icing patterns for $n\geq 2$ have almost the same difference in the potential energy compared to that before the pressurization. This is because they gain the same amount of energy from the same lateral pressure and consume the same amount of energy to form the same homogeneous square icing patterns in the pressurization process. Different from the situation for $n=1$ in which water molecules have enough potential energy to support the formation of all hydrogen bonds, water molecules for $n\geq 2$ have significantly lower potential energy, making it energetically unfavorable to form unidirectional, homogeneous icing patterns. Therefore, for $n\geq 2$, the partitioned square icing structures with disordered local regions are more energetically favorable to form.

For $n=1$ (i.e., the water molecules are confined between two parallel graphene sheets), there are two typical scenarios of square ice formation, and their distributions of water molecules during the formation of square ice are shown in Figure 6 and Movie S4. In Figure 6(a), square ice forms without partitioned patterns during the initialization stage. Therefore, in the subsequent growth stage, the square icing pattern is always homogeneous and eventually becomes stabilized. In contrast, in Figure 6(b), there are two regions with different orientations of square icing patterns during the initialization stage. However, in the growth stage, one of them gradually becomes smaller and finally disappears, resulting in the homogeneous square icing pattern in the end. These typical processes reveal that for $n=1$, partitioned square icing patterns will finally disappear during the growth stage because of the poor stability even if they are formed during the initialization stage. In the simulations, the energy gain of water molecules comes from graphene sheets and lateral pressure. The critical value of lateral pressure for ice formation can be estimated by the adhesion according to the previous studyAlgara-Siller et al. (2015) ($P_{W}=E_{W}/d\approx 0.9$ GPa, where $E_{W}\approx 30$ meV/Ų is the difference between graphene-graphene and graphene-water adhesion energies per unit area, and $d\approx 5.6$ Å is the difference in the distance between plates due to the trapped water molecules), and confirmed by calculations of the per-atom potential energy (about 1.0 – 1.1 GPa as shown in Figure S2 in Supplementary Material). Therefore, after the pressurization to 1.1 GPa for $n=1$, water molecules have enough potential energy to support the formation of all hydrogen bonds, which can drive water molecules at the grain boundaries between the two regions to change from a disordered liquid state to the ordered crystalline state, and the homogeneous square icing patterns will eventually form.

For $n\geq 2$ (i.e., the water molecules are confined between two parallel graphite sheets), there are also several typical scenarios of square ice formation, and their distributions of water molecules during the formation of square ice are shown in Figure 7 and Movies S5 & S6. In Figure 7(a), square ice forms without partitioned patterns during the initialization stage (same as the graphene case shown in Figure 6(a)), leading to a homogeneous square icing pattern that becomes stable eventually. In contrast, in Figure 7(b),here are some regions with square icing patterns in different orientations during the initialization stage. However, these small regions gradually enlarge and stabilize in the growth stage instead of becoming smaller and disappearing, finally leading to the partitioned square icing patterns. Even if more than two regions with different icing orientations are formed during the initialization stage in some particular cases, as shown in Figure 7(c), a certain local region gradually enlarges and stabilizes while the other eventually disappears, forming the two-region partitioned icing patterns in the end. As analyzed in the previous section, water molecules have significantly lower potential energy when $n\geq 2$, which is not enough for all the water molecules at the grain boundaries to form the hydrogen bonds and change their disordered liquid state to the ordered crystalline state. Simulation results reveal that the randomness of the partitioned structure for $n\geq 2$ originates from the initialization stage. Because of the short-range nature of the interaction between the water molecules and the carbon substrate, the effect of the carbon substrate determines the stability of the partitioned structure during the growth stage. The difference in stability eventually determines the probability of forming partitioned square icing patterns.

The mechanisms of square ice formation generalized from the typical scenarios are further confirmed and quantified by analyzing the time evolution of the number of ice molecules, as shown in Figure 8. After the stabilization, the total number of water molecules between the two parallel graphene-based sheets is approximately 500, most of which form the ordered structure, and the rest are grain boundaries or point defects. In Figure 8(a), the number of ice molecules for case A rises suddenly after pressurization (at 5 ns) and then remains almost constant, corresponding to that square icing pattern is always homogeneous and eventually becomes stabilized (Figure 6(a)). The number of ice molecules for case B is significantly lower than case A after pressurization (at 5 ns), which means there are multiple regions with different orientations of square icing patterns and water molecules at the grain boundaries are disordered during the initialization stage. Then it rapidly increases until about 10 ns, corresponding to that, regions with different orientations become a homogeneous square icing pattern in the growth stage (Figure 6(b)). Similarly, in Figure 8(b), the numbers of ice molecules for cases A, B, and C correspond to the typical scenarios shown in Figure 7(a), (b), and (c), respectively. The numbers of ice molecules for cases B and C are always lower than that of case A because the orientations of water molecules at the grain boundaries between the two regions are random, indicating that the partitioned structure is stable during the growth stage for $n\geq 2$. The results show that carbon substrate of graphene and graphite determined the stability of the partitioned structure during the growth stage.

Besides the force by the carbon sheets on the water molecules, the confinement dimension is another important factor for the formation of the square icing patterns, as shown in Figure S6 in the Supplementary Material. In all of the above simulations, the height of the graphene capillary $h$ (the separation between the two sheets in $y$ direction) is fixed at 9.0 Å, which is appropriate for the formation of the bilayer square ice. To consider the effect of the confinement dimension, we reduced $h$ to 6.5 Å, which is appropriate for the formation of the flat monolayer square ice. We also varied the number of graphene monolayers $n$ to find out the possibility of partitioned square icing patterns. However, we did not observe the partitioned square icing patterns for the single layer of water molecules. To further verify the numerical simulation, we repeated the simulation for another 9 cases by using different initial velocity distributions, and still did not obtained the partitioned square icing patterns.

A square icing parameter map is produced for these simulation results, as shown in Figure 9. From the results of single-water-layer cases, we can see that it is energetically unfavorable for the single-layer water molecules confined between two parallel graphite sheets to form partitioned square icing patterns, indicating that the interaction between different layers of water molecules is a dominant factor in partitioned square icing. Further analysis of the results for $h=9.0$ Å reveals that both the top and the bottom layers of water molecules have the same behavior during the crystallization process and form identical square icing patterns after stabilization. This feature indicates that there is strong interaction between the layers of the water molecules that maintains these aligned structures. Like the hydrogen bonds between water molecules in the same layer, the interaction between the top and the bottom layers will consume part of the energy. When there is only one layer of water molecules between graphite sheets, the energy consumed for the interaction between layers is not needed. As a consequence, the water molecules have more potential energy to support the formation of all hydrogen bonds in the graphene capillary, even if they are confined between two parallel graphite sheets ($n\geq 2$).