Molecular Simulation, Vol. 29 (10–11), October–November 2003, pp. 619–626

Formation of Quasi Two-dimensional Bilayer Ice inHydrophobic Slits: A Possible Candidate for Ice XIII?
J. BAIa, X.C. ZENGa,*, K. KOGAb and H. TANAKAca
Department of Chemistry, University of Nebraska-Lincoln, Lincoln, NE 68588, USA; bDepartment of Chemistry, Cornell University, Ithaca,
NY 14853, USA; cDepartment of Chemistry, Okayama University, Okayama 700-8530, Japan(Received October 2002; In final form February 2003)

We performed molecular dynamics simulations of water
confined to hydrophobic slit nanopores. Two five-site
potential models of water were employed: the ST2 modeland the TIP5P model. The simulations confirm ourprevious simulation results on basis of the four-site
TIP4P model of water, that is, upon cooling the confined
liquid water may undergo a first-order phase transition toeither a bilayer crystalline ice phase or to a bilayer
amorphous ice. The selection of the bilayer crystalline
phase occurs at the constant normal pressure condition
whereas the selection of the bilayer amorphous ice phaseis likely to occur at the fixed pore-width condition. This
computer-simulation generated ice form, if confirmed by
crystallographic or spectroscopic experiments, may be a
candidate for “ice XIII”, provided that purely quasi twodimensional ice forms (metastable in vacuum) can beviewed as a stand-alone solid phase of ice.
Keywords: Ice XIII; Bilayer ice; ST2; TIP5P

INTRODUCTION
Three-dimensional bulk ice displays a rich phase
diagram [1]. Under pressure, bulk ice can exhibit
many solid-phase transformations, which include avariety of metastable crystalline phases and amorphous phases [2]. To date, 13 forms of bulk crystalline
ice have been found either in nature or inlaboratories. Ice I alone has two polymorphs, one is
the hexagonal ice Ih (i.e. the normal ice) and anotheris the metastable cubic ice Ic. The latest metastableform of ice, ice XII, was uncovered in 1998 [3].
It is also known that new phases can emerge inmicroscopically confined environments [4 – 6] such

*Corresponding author.
ISSN 0892-7022 print/ISSN 1029-0435 online q 2003 Taylor & Francis Ltd
DOI: 10.1080/0892702031000103158

as slit pores. The special character of confinedsystems stems from their spatial inhomogeneity.
Forced to pack into spaces severely constricted by
two confining walls, molecules in the vicinity of thefluid – wall interface tend to arrange themselves inlayers parallel with the interface. The resulting
oscillations in local density are reflected in propertiesof the confined phase that can differ drastically fromthose of the bulk counterpart. Some of the new phasebehavior, such as wetting and layering transitions, isdue solely to the presence of the interface; others,
such as capillary freezing, can be regarded “simply”
as bulk behaviors modified by the presence of theconfining walls.
Recently, we reported computer simulation evidences of the possible existence of a new quasi twodimensional (Q2D) form of crystalline ice [7 –10] anda Q2D amorphous ice [11] for water confined to ahydrophobic nanoslit. Both types of transition—liquid water to Q2D ice and amorphous ice—arestrongly first-order. In our previous MD simulations,the four-site TIP4P water model was employed. In
these MD simulations, we used two ways to control
the confinement conditions. In one case [7 – 10], weallowed the two confining walls to be movable inorder to achieve a constant normal pressure
condition. In another case, we fixed the distance inbetween the two walls while allowing the lateral areaof the simulation cell to change in order to achieve a
constant lateral pressure condition. In the former
case, the bilayer ice crystal was observed uponcooling whereas in the latter case the amorphousbilayer ice was observed. Presumably, by fixing
the pore width gives rise to certain “frustration” in

620

J. BAI et al.

the normal direction (perpendicular to the walls) toprevent the formation of perfect crystal. Furthermore, the change in the lateral area during thesimulation seems to allow the system to be moreeasily “trapped” in the metastable amorphous state.
Could the formation of these Q2D ice andamorphous ice be an artifact of the TIP4P [12]model of water? To address this question, we have
performed a new set of MD simulations with otherpopular models of water, the five-site ST2 model [13]and TIP5P model [14]. Note that the four-site TIP4P
model has a planar structure whereas both ST2 and
TIP5P models are tetrahedral-like in geometry. The
ST2 model is amongst the earliest molecular models
of water, developed by Stillinger and Raman in the
1970s, and is known to be strongly tetrahedral. It is a
modified version of the BNS model [13]. The TIP5P isthe latest model in the TIP series developed byJorgensen and coworkers. One improvement overthe previous TIP models is that the TIP5P canquantitatively reproduce the temperature at which
the maximum density occurs.

MODELS AND METHODS
The five sites of ST2 and TIP5P model consist of
one oxygen (O), two hydrogen (H) and two virtualsites (M). For the ST2 model, the bond length of
OZM is 0.8 Å and OZH is 1.0 Å. Both the /MOM
and /HOH angles are tetrahedral-like (109.478).
The oxygen site has Lennard – Jones (LJ)parameters, e ¼ 0:07575 kcal=mol and s ¼ 3:10 A:
Two hydrogen sites carry a charge, q ¼ 0:2357e:
A negative charge 2 q is assigned to each virtualsite (M). The ST2 potential function involves ashort-range switching function that smoothly cuts
off the potential at the OZO distance (roo) of two themolecules. At distances less than roo ¼ 3:1287 A;
switching function is active but it completely cuts
 In contrast, nooff the interaction at roo ¼ 2:0106 A:
short-range switching function is involved in the
TIP5P potential function. For TIP5P, the bondlength of OZH is 0.9572 Å and OZM is 0.7 Å. The
angle /MOM is the same as ST2 but /HOH is104.528. The H site has a charge, q ¼ 0:241e: The M
site has a negative charge 2 q. The LJ parametersof O site are e ¼ 0:16 kcal=mol and s ¼ 3:12 A:
The wall is a smooth planar wall, and the wall –oxygen site interaction is taken to be the nonpolar
9 –3 LJ interaction [7,15] with the LJ parameters Thus, thee ¼ 1:1827 kcal=mol and s ¼ 2:4782 A:
wall is hydrophobic.

As in our previous MD simulations, we used two
kinds of confined systems [7,11]. In both systems,
water was confined in between two plane-parallelsmooth solid walls. In the first system (system I) a
constant normal-pressure condition was imposed [7]while the lateral area of the system was fixed. To thisend, the distance in between the two planar walls
was allowed to change. In the second system (system
II), the constant lateral-pressure condition was
imposed [11] while the distance in between the twowalls was fixed. The lateral area of the simulation cell
was allowed to vary in order to achieve the constantlateral-pressure condition. For both systems, theperiodic boundary conditions were applied in lateral
directions (x and y) in parallel with the two walls.
In system I, the simulation cell contained 192 water2
molecules.† The planar wall is a 34:16 £ 28:797 A

rectangular prism (for TIP5P) or 33:55 £ 29:06 A2
prism (for ST2). The normal pressure was fixed at
1 GPa. The MD simulation started at a temperature of
330 K. We lowered the system temperature in steps(typically 10 K in temperature change) from 330 K
until a freezing transition was observed. We thenraised the temperature sequentially (also in steps) up
to the point where a melting transition was observed.
Smaller temperature steps (1 –4 K in temperaturechange) were used near the transition point.
At each temperature, the MD run time was in
between 2 and 15 ns. Note that the MD timestep is
0.4 fs. For TIP5P model, longer MD times up to 40 ns
were used near the transition point, whereas for ST2
model it was found that 4 ns MD time was sufficient
to observe the 2D ice formation process near thetransition point, due to the strongly tetrahedral
character of ST2 model.
In system II, the distance between the twowalls was fixed. The pore space in the normal
direction (z direction) just accommodates about
two layers of water molecules. The simulation
cell contained 200 water molecules. The lateralpressure was fixed at 0.1 MPa. As for system I,
the temperature of system II was lowered insteps starting from 330 towards 250 K, and thenraised sequentially until a melting transition wasobserved.

RESULTS AND DISCUSSION
Figure 1a shows the temperature dependence of thequenched potential energy (QPE) for both ST2 and
TIP5P water in system I. The QPE was obtained byusing the steepest descent method to map a few tens†
Previously we have studied effects of system size [8,11] and periodic-boundary conditions [9,10] for the TIP4P model. The largest
system size studied consisted of 960 water molecules [11]. It was found that the final hexagonal morphology of 2D bilayer ice is notdependent on the system size.

ICE IN HYDROPHOBIC SLITS621

FIGURE 1 Temperature dependence of the QPE in (a) system I and (b) system II.

of instantaneous configurations of the confinedwater onto their inherent structure [16]. Note that
the QPE includes a vibrational energy 3RT where R is
the gas constant and T is the temperature of thesystem, but excludes the interaction potential energy
between water and the walls. In system I, as thetemperature was lowered the QPE exhibits a sharp
drop at 293 K for ST2 and at 280 K for TIP5P. In theheating process, a sharp jump in the QPE occurredat 320 K for ST2 and at 327 K for TIP5P.
This hysteresis-loop behavior demonstrates

the appearance of a strong first-order phasetransition in the confined system. Figure 1b shows
the temperature dependence of the QPE for system
II. Again, a similar hysteresis-loop behavior wasobserved. Figure 2 displays the quenched configuration of the TIP5P water prior to and after the phasetransition in the system I. As shown in Fig. 2b,c, thequenched configuration in system I is a bilayer ofhexagonal ice crystal, whereas in system II it is apolygonal amorphous bilayer ice. Both types ofphase transition were observed previously on

622

J. BAI et al.

FIGURE 2 Quenched structures of the confined water in the state of (a) liquid, (b) bilayer crystal and (c) amorphous bilayer ice. Larger
spheres denote oxygen and smaller spheres denote hydrogen.

the basis of the TIP4P model. However, the bilayer
hexagonal ice structure obtained from the two fivesite models differs somewhat from that from thefour-site model. Indeed, the hexagons in both the ST2
and TIP5P bilayer crystal appear to be much lessdistorted than those in the TIP4P bilayer crystal. In
Figs. 3 – 5 we show more detailed structural
information for the TIP5P and ST2 bilayer ice andamorphous bilayer ice.
Figure 3a,b displays a typical geometry and theassociated three OZOZO angles (a, b and g) madeby a center oxygen site with its three nearest oxygensites in the same layer of the Q2D ice crystal. Figure 3c
displays the OZOZO angle distribution for

the TIP5P and ST2 bilayer crystals. It can be seenthat the OZOZO angles in the TIP5P bilayer crystal
exhibit five sharp peaks at 115.98, 117.98, 119.58, 121.38
and 124.98. The highest peak is located at 121.38.
However, the OZOZO angles of the ST2 bilayer
crystal shows only one broad peak which is locatedat 1208. In contrast, for the TIP4P bilayer crystal three
distinctive OZOZO angles were found, which are
1088, 118.58 and 135.58. Based on the angle
distribution for TIP4P, several bilayer ice rules weresuggested. In particular, we observed that hexagons
in the TIP4P bilayer crystal exhibit two orientations,which were denoted as orientation A or B. It was
found that all hexagons in the same row adopt either

ICE IN HYDROPHOBIC SLITS623

FIGURE 3 A typical geometry and the associated three OZOZO angles (a, b and g) made by a center oxygen site with its three nearest
oxygen sites in the same layer of the Q2D ice crystal. (a) For five-site models—TIP5P : a ø b – g; ST2 : a – b – g; (b) For the TIP4P model[7]—a ¼ 1088; b ¼ 118:58 and g ¼ 135:58: (c) The OZOZO angle distribution of TIP5P and ST2 bilayer crystals.

the orientation A or B. In the TIP4P bilayer crystal themost preferred stacking for the hexagon rows
appears to be ABAB. However, for both TIP5P andST2 bilayer crystals the OZOZO angle distribution ismuch narrower than that of TIP4P. In fact, the angle
distribution peaks for TIP5P are located within 1158
and 1258 while for ST2s the peak of angle distribution
is located within 1108 and 1308.
Another local geometric information is the
distribution of the nearest-neighbor OZO distance.
For TIP5P and ST2 bilayer crystal, the distribution isshown in Fig. 4. Comparing the averaged OZO
distance in the bilayer crystal to the OZO distancein water dimer (2.68 Å for TIP5P [14] and 2.74 Å forST2 [13]) it can be found that the averaged OZO
distance in the bilayer crystal is greater. In contrast,
the averaged OZO distance in the TIP4P bilayer
crystal is nearly the same as that in the TIP4P dimer(2.74 Å). As shown in Fig. 4, the OZO distance
distribution for the two five-site models are muchbroader than that for the four-site TIP4P model [9,11].
These quantitative differences in local geometricconfiguration stems mainly from the molecular
structure differences between the four-site planar
model and the five-site tetrahedral model.
The two-dimensional radial distribution function
was also calculated to provide global structural
information of the bilayer solids. Figure 5a,b show
the OZO radial distribution of TIP5P and ST2 bilayer
crystal and amorphous bilayer ice, respectively. First,
it can be seen that for TIP5P crystal the first peak isslightly shifted from the origin, which does not occur
for the ST2 bilayer crystal. Since the first peakdescribes the relative location of the oxygen in the
opposing layer of the crystal, it appears that the twolayers of the TIP5P crystal do not exactly overlap

J. BAI et al.

624

FIGURE 4

The OZO bond length distribution of TIP5P and ST2 bilayer crystals.

with each other. The OZH· · ·O hydrogen bonding in
between the two layers has typically an angle
about 108 with respect to the normal (z) direction.
On the other hand, the two layers of ST2 crystal
nearly overlap with each other, a feature similar to
the TIP4P bilayer crystal. Second, it can be seen thatthe location of the major peaks in the radial
distribution functions is about the same for boththe TIP5P and ST2 bilayer crystal, indicating theoverall structure of the two bilayer crystals are about
the same. However, the major peaks of TIP5P crystal
are further split into a number of sharp peaks
whereas those of ST2 crystal are much broader withless fine features. Moreover, the radial distributionfunction of the TIP5P and ST2 amorphous bilayer ice
also resembles each other except the second major
peak for the TIP5P amorphous ice is further split into
three sharp peaks.
To summarize, our new sets of MD simulations
with two five-site water models confirms ourprevious simulation on the basis of the four-site
models [7,11], that is, in system I the bilayer ice
crystal forms upon cooling whereas in system II the
amorphous bilayer ice forms. A question arises: why
do the two systems result in different polymorphous
structure of bilayer solids? The answer to thisquestion is still inconclusive. We have foundpreviously that if the confined water can undergo afreezing transition into system I, the bilayer crystal islikely to form with some defects [9,10]. It seems that
there exists some special numbers of water molecules when used in the simulation cell to form
perfect bilayer crystals. When the number of defects
is small, they can be removed by repeating several

heating and cooling processes to the defected crystal.
But this does not work for system II. Can it bepossible to form a perfect crystal in system II undercertain conditions? To explore this possibility, we
performed an MD simulation with 192 ST2 molecules in system II at the same simulation conditions
and with 4 ns MD time at each temperature, except
we abruptly changed the temperature from 295 to270 K (25 K change in temperature). After theamorphous bilayer ice was observed, the amorphous
ice was heated from 270 K directly to either 280, 290,
295 or 300 K. When it was below 300 K, we only
observed that some polygons in the amorphous icechanged into some other polygons. When it was at
300 K, however, the polygonal structure of theamorphous ice began to rearrange and eventually ahexagonal network was formed with a few defects.
At this point, we increased the temperature further
from 300 to 310 K and performed the MD simulation
for about 0.8 ns, and then lowered the temperature
directly to 280 K before the structure reached a stableconfiguration at 310 K. In this case, a perfect bilayer
crystal was formed at 280 K after about 2 ns MD timein system II. The average potential energy of thecrystal is about 0.387 kcal/mol lower than that of theamorphous ice. Therefore, once the crystal wasformed, we did not detect any evidence oftransformation from the bilayer crystal to amorphous bilayer ice during further heating and coolingprocesses. These results show that the perfecthexagonal network is very stable despite the factthat the formation of bilayer crystal is a rare event insystem II. Similar behavior was found for TIP5P in
system II.

ICE IN HYDROPHOBIC SLITS

625

FIGURE 5 Two-dimensional radial distribution of oxygen of (a) TIP5P and (b) ST2 bilayer crystal and amorphous bilayer ice.

CONCLUSION
We have performed MD simulations of water
confined to hydrophobic nanoslits to investigate
the first-order phase transition phenomena. Two
five-site water models, the TIP5P and ST2, were
employed. Our new set of simulation confirms ourprevious simulations with the TIP4P model. It isfound that, upon cooling, the confined water transforms into a bilayer ice crystal in system I and into an
amorphous bilayer ice in system II. From the OZO

distance and OZOZO angle distributions of thebilayer crystals we find that hexagons are much lessdistorted and thereby do not exhibit any special
orientations like those in the TIP4P bilayer ice.
We also find that ST2’s hexagonal network is more
flexible compared to the TIP4P’s and TIP5P’s. Finally,
although the formation of a bilayer crystal is notpreferred in system II, it is still conceivable togenerate a bilayer crystal via multiple heating andcooling processes (somewhat like a simulated
annealing). Once the bilayer crystal is formed, no626

J. BAI et al.

reverse process, namely, from the bilayer crystal toamorphous bilayer ice was observed, indicating thatthe bilayer crystal phase is indeed more stable at low
temperatures even in system II. In conclusion, thecomputer-simulation generated Q2D bilayer ice may
be considered as a possible candidate of the “Ice
XIII” or whatever the next Roman numerals would
be available at the time, provided that the 2D ice
phases in vacuum can be recognized as a stand-alone(metastable) solid ice. Alternatively, one could
introduce Roman numerals specifically for 2D ice
phases. In any case, ultimate confirmation andnaming of the computer-generated 2D ice structure
requires crystallographic or at least spectroscopic
experiments. For now, we name the 2D bilayer ice the“Nebraska ice” because according to the Otoe Indian
the word “Nebraska” means “flat water” [17].
Acknowledgements
X.C. Zeng is supported by US NSF and Office ofNaval Research. K. Koga and H. Tanaka aresupported by the Japan Society for the Promotionof Science (JSPS) and the Japan Ministry of
Education.

References[1] Petrenko, V.F. and Whitworth, R.W. (1999) Physics of Ice(Oxford, New York).
[2] Mishima, O. and Stanley, H.E. (1998) “The relationshipbetween liquid, supercooled and glassy water”, Nature 396,
329.

[3] Lobban, C., Finney, J.L. and Kuhs, W.F. (1998) “Thestructure of a new phase of ice”, Nature 391, 268.
[4] Gubbins, K.E., Sliwinska-Bartkowiak, M. and Suh, S.-H.
(1996) “Molecular simulation of phase transitions in pores”,
Mol. Simulation 17, 333.
[5] Gelb, L.D., Gubbins, K.E., Radhakrishnan, R. and SliwinskaBartkowiak, M. (1999) “Phase separation in confinedsystems”, Rep. Prog. Phys. 62, 1573.
[6] Schoen, M. (1993) Computer Simulation of Condensed Phases inComplex Geometries (Springer, Heidelberg).
[7] Koga, K., Tanaka, H. and Zeng, X.C. (1997) “Freezing ofconfined water: a bilayer ice phase in hydrophobic
nanopores”, Phys. Rev. Lett. 79, 5264.
[8] Slovák, J., Tanaka, H., Koga, K. and Zeng, X.C. (1999)“Confined water in hydrophobic nanopores: dynamics offreezing into bilayer ice”, Phys. Rev. E 60, 5840.
[9] Slovák, J., Tanaka, H., Koga, K. and Zeng, X.C. (2001)“Computer simulation of water–ice transition in hydrophobic nanopores”, Physica A 292, 87.
[10] Slovák, J., Tanaka, H., Koga, K. and Zeng, X.C. (2003)“Computer simulation of bilayer ice: structures and thermodynamics”, Physica A 319, 163.
[11] Koga, K., Tanaka, H. and Zeng, X.C. (2000) “First-ordertransition in confined water between high density
liquid and low density amorphous phases”, Nature 112,8910.
[12] Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impay, R.W.
and Klein, M.L. (1983) “Comparison of simple
potential functions for simulating liquid water”, J. Chem.
Phys. 79, 926.
[13] Stillinger, F.H. and Rahman, A. (1974) “Improved simulationof liquid water by molecular dynamics”, J. Chem. Phys. 60,
1545.
[14] Mahoney, M.W. and Jorgense, W.L. (2000) “A five-site model
for liquid water and the reproduction of the density anomaly
by rigid, nonpolarizable potential functions”, J. Chem. Phys.
112, 8910.
[15] Steele, W.A. (1974) Interaction of Gases with Solid Surfaces(Pergamon Press, Oxford).
[16] Stillinger, F.H. and Webber, T.A. (1984) “Packing structures
and transitions in liquids and solids”, Science 225, 983.
[17] Bradley, D. (2000) “Water, water. . .”, Alchemist, www.chemweb.com/alchem/articles/985883672133.html.