JOURNAL OF CHEMICAL PHYSICSVOLUME 118, NUMBER 9

1 MARCH 2003

COMMUNICATIONS

Ab initio studies of quasi-one-dimensional pentagon and hexagon ice
nanotubes
J. Bai
Department of Chemistry and Physics, University of Nebraska, Lincoln, Nebraska 68588

C.-R. Su
Department of Physics and Center for Atomic and Molecular Nanosciences, Tsinghua University, Beijing,
100084, China

R. D. Parra
Department of Chemistry, DePaul University, Chicago, Illinois 60604

X. C. Zenga)
Department of Chemistry, University of Nebraska, Lincoln, Nebraska 68588

H. Tanaka
Department of Chemistry, Okayama University, Okayama 700-8530, JapanK. Koga
Department of Chemistry, Cornell University, Ithaca, New York 14853

J.-M. Li
Department of Physics and Center for Atomic and Molecular Nanosciences, Tsinghua University, Beijing,
100084, China

共Received 13 November 2002; accepted 3 January 2003兲
Ab initio plane-wave total-energy calcuation is carried out to study the relative stability of thequasi-one-dimensional 共Q1D兲 pentagon and hexagon ice nanotubes. Electronic structure
calculations indicate the two Q1D ice nanotubes have nearly the same band structures and energy
bandgap as those of proton-ordered bulk ice Ih . Ab initio molecular-orbital and density-functional
theory calculations, as well as three classical potential models of water, are also employed toinvestigate the relative stability of the pentagon and hexagon water clusters (H2 O) 30 , (H2 O) 60 , and(H2 O) 120 . Clusters of this kind can serve to bridge the gap between the small polygonal water rings
and the infinitely long Q1D polygon ice nanotubes. It is found that the polygon water prisms withthe size (H2 O) 120 begin to show the relative energetic behavior of the infinitely long polygon ice
nanotubes. © 2003 American Institute of Physics. 关DOI: 10.1063/1.1555091兴

Ab initio computer simulations1 can explore regions ofphase diagram not easily accessible by laboratory experiments. For example, the first direct observation of phasetransformation of the densest ice phase 共among thirteen crystalline phases of ice兲, ice X, was established from ab initio
molecular dynamics simulation,2 whereas the experimental
evidence of existence of the ice X is thus far only inferredfrom spectroscopy.3 Recently, classical molecularmechanics4,5 and molecular dynamics simulations6,7 of water
encapsulated in carbon nanotubes have been reported, whichsuggests possible existence of four quasi-one-dimensional
共Q1D兲 polymorphs of ice nanotube. Among these Q1D polymorphs, pentagon and hexagon ice nanotubes are the twomost stable ones.4 In the previous studies, the TIP4P potential of water8 has been used. Since the potential is derived
from fitting to experimental data for bulk water, extension of
the potential to highly confined water may be only qualita-

tive or semiquantitative. Thus, the predicted polymorphs of
Q1D ice remains to be confirmed, at least by more quantitative means, although direct experimental observation ofthese phases is the ultimate confirmation.9 In this communication, we report results of ab initio pseudopotential totalenergy calculation to further affirm the existence of polymorphs of Q1D ice nanotube. Our calculations show thatboth hexagon and pentagon ice nanotubes are metastable
solid phase in vacuum at 0 K and they are nearly isoenergetic. Moreover, calculation of the electron density of statesindicated that the two Q1D polymorphs have nearly the sameenergy bandgap as the proton-ordered bulk ice Ih . Finally,
relative stability of the finite-size pentagon and hexagon water clusters at 0 K is also examined using all-electron quantum chemistry methods as well as three additional classicalpotential models of water.
For Q1D ice nanotubes, we used the Cambridge Serial
Total Energy Package 共CASTEP兲10 to calculate the total energy per molecule and forces within the framework ofdensity-functional theory 共DFT兲. The exchange-correlationa兲Electronic mail: xzeng1@unl.edu

0021-9606/2003/118(9)/3913/4/$20.00

3913© 2003 American Institute of Physics

Downloaded 26 Mar 2003 to 128.253.229.143. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

3914

J. Chem. Phys., Vol. 118, No. 9, 1 March 2003

Bai et al.

FIG. 2. Electronic density of states 共DOS兲 of Q1D hexagon ice nanotube
共solid line兲, pentagon ice nanotube 共dotted line兲, and proton-ordered bulk ice
Ih 共dash-dotted line兲.

FIG. 1. A lateral cut of pentagon and hexagon ice nanotubes with ABABstacking 共Ref. 4兲 results in a finite-size 共a兲 pentagon and 共b兲 hexagon water
cluster, respectively. Big and small spheres represent oxygens and hydrogens, respectively, and the long thinner bonds denote the hydrogen bonding.

functional is treated in the Perdew–Wang 共PW91兲11
generalized-gradient approximation 共GGA兲. The wave functions are expanded by using a plane wave basis set with akinetic energy cutoff of 450 eV. The ion-valence electroninteractions are represented by ultrasoft pseudopotential.12
The Brillouin zone was sampled with (1⫻1⫻5) k points ofa Monkhorst–Pack grid.13 Because the hexagon and pentagon ice nanotubes are Q1D, the supercell geometry is taken
to be a tetragonal cell with the dimension L⫻L⫻L z where
the z direction is chosen to be the axial direction of the ice
nanotubes. We used a 10-molecule supercell for pentagon
and a 12-molecule supercell for hexagon ice nanotube. The
two water layers in the supercell were chosen to be the
ABAB-stacking4 共Fig. 1兲. Thus, L z is just 2a z where a z is themean lattice constant of the Q1D ice in the z direction. In thecalculations, L was chosen to be 20 Å. The CASTEP code
allows full geometry optimization for Q1D periodic systems,
i.e., L z is allowed to vary to achieve the zero-pressure condition in the axial direction. The energy criterion for geometry optimization is 5⫻10⫺6 eV/molecule. The final optimized value of L z is 5.672 Å, corresponding to theminimized energy of ⫺471.5501 eV/molecule for hexagon
ice nanotube. Those for the pentagon ice nanotube are 5.557
Å and ⫺471.4386 eV/molecule, respectively. Hence, theGGA calculation indicates that the hexagon ice nanotube isslightly more stable than the pentagon ice nanotube; the energy difference is 0.1115 eV/molecule or 10.76 kJ/mol.
We also independently examined the relative stability of
the hexagon and pentagon ice nanotubes using the Vienna
ab initio simulation package 共VASP兲.14 In this case, theexchange-correlation functional was treated in the local-

density approximation 共LDA兲 with the Ceperley–Alderpotential15 based on quantum Monte Carlo simulations. Kinetic energy cutoff is taken to be 395.7 eV. The Brillouinzone was sampled with (1⫻1⫻8) k points of a Monkhorst–
Pack grid. Again, the supercell contains 10 molecules forpentagon ice nanotube and 12 molecule for hexagon ice
nanotube. We used 20 Å⫻20 Å⫻L z supercell size. The energy criterion for geometry optimization was set at
10⫺4 eV/supercell. The optimized value of L z is 5.21 Å forboth hexagon and pentagon ice nanotubes, both are slightly
smaller than those obtained from the GGA calculation. It is
known that LDA generally overbinds water molecules,16
which leads to a smaller a z . Moreover, in contrast to the
GGA calculation, the minimized energy for pentagon ice
nanotube is actually lower than that for hexagon ice nanotube, but the energy difference is merely 0.0012
eV/molecule. That the prediction of the relative stability between the two Q1D polmorphs are opposite on the basis of
GGA and LDA suggests that the two ice polymorphs may benearly isoenergetic. To summarize, two independent ab initio
calculations show that the Q1D pentagon and hexagon ice
nanotubes are metastable phases of ice at 0 K in vacuum
since 3D bulk ice is expected to be the more stable phase at
this condition. However, the Q1D ice nanotubes can becomea stable phase relative to other confined Q1D phases in ahydrophobic pore with a diameter about 1 nm.6 The van der
Waals attraction between the the confined water and the wall,
although weak in comparison with the hydrogen bonding
interaction, can further stabilize the Q1D ice nanotubes relative to other 3D bulk phases so that the confined Q1D ice canbe in equilibrium with a bulk phase 共e.g., liquid water兲.7
Figure 2 displays the calculated electron density of states
共DOS兲 for both hexagon and pentagon ice nanotubes. For
comparison, the DOS of proton-ordered bulk ice Ih is also
shown in Fig. 2. The DOS calculation was based on PW91
GGA within the DFT. The Fermi energy was set at zero.
Interestingly, we find that the Q1D ice nanotubes have thenearly same value of the energy bandgap, E g ⬇5 eV, as thatof the bulk ice Ih . Note that the optical absorption
experiments17,18 have shown that the absorption edge which

Downloaded 26 Mar 2003 to 128.253.229.143. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 118, No. 9, 1 March 2003

corresponds to the energy bandgap, E g , is 7.8 eV for protondisordered ice Ih . Our calculated DOS for the bulk ice Ih isin very good agreement with a previous theoretical calculation of electron energy spectrum for the proton-disordered
ice Ih . 19 In that work, the tight-binding approach was used
for which a hopping matrix element was adjusted to fit theexperimental bandgap 7.8 eV. For bulk ice Ih both calculations show that the DOS has a few singularity-like peaks dueto the high degrees of degeneration at these values ofenergy.19 In addition to the bandgap, we find the two Q1D
ice nanotubes show nearly the same electron energy spectrum as that of ice Ih except some small difference in finepeaks. We therefore tentatively conclude that the overall
DOS features appear to be not very sensitive to the local
hydrogen-bonding structure so long as the entire molecularsystem has long-range positional order.
The fine electronic structural differences between the
pentagon and hexagon ice nanotubes and the bulk ice Ih canbe shown by using the Mulliken population analysis. For
bulk ice Ih , it is found that the Mulliken charge is 0.49e for
H and ⫺0.98e for O. In the case of pentagon ice nanotube,
however, O still has a charge ⫺0.98e as that of ice Ih , but H
exhibits two different charges, 0.51e for in-plane H or 0.47e
for H involved in hydrogen bonding in the axial direction.
The calculated electron density distribution also showed thatthe hydrogen bonds in the plane of pentagon exhibit aslightly denser electron density than that for the hydrogenbonds in the axial direction. Finally, in the case hexagon ice
nanotube, a variety of Mulliken charges on O and H sites
was found. In fact, H has eight different charges ranging
from 0.47e to 0.53e while O has two charges, ⫺1.00e or⫺1.01e. This result indicates that the O–O–O angles in the
hexagon ice nanotube deviate appreciably from the tetrahedral angle 109.47°, and thus gives rise to a quite different
electronic distribution compared to that of pentagon ice
nanotube.
In addition to the pseudopotential total-energy calculation for the Q1D ices, we also performed large-scalequantum-chemistry calculations to examine energetics offinite-size pentagon and hexagon water clusters consisting of30, 60, and 120 molecules. Clusters of this kind20 can serveto bridge the gap between the small polygonal waterrings21,22 and the infinitely long Q1D polygon ice nanotubes.
Both molecular-orbital and density-functional methods
共implemented in GAUSSIAN 98 software package兲23 were employed. The water clusters were generated from a lateral cutof the Q1D ice nanotubes. Thus, these clusters are either
pentagonal or hexagonal prisms 共see Fig. 1兲. The geometriesof the pentagon and hexagon water clusters (H2 O) 30 ,(H2 O) 60 , and (H2 O) 120 were fully optimized using theB3LYP/6-31⫹G共d兲, B3LYP/6-31G, and HF/6-31G levels oftheory, respectively. The geometries of the optimized clusters
were then used for single-point energy calculation at the HF/
6-31G共d,p兲 level for all clusters. At this level, the number ofbasis amounted to about 3000 for (H2 O) 120 . Furthermore,
the correlation corrections were evaluated at the MP2/631G共d兲 level, but only for the smaller (H2 O) 30 and (H2 O) 60
clusters.
Single-point energy calculations at the HF/6-31G共d,p兲

Ice nanotubes

3915

FIG. 3. Energy difference 共per molecule兲 ⌬E⫽E hex ⫺E pen vs the number of
H2 O molecules in hexagon and pentagon water clusters.

level show that the pentagon water clusters are energeticallymore stable than the hexagon water clusters. Figure 3 showsthat the energy difference per molecule versus the size of theclusters. For (H2 O) 30 the hexagon cluster is 0.244 kJ/mol
less stable than the pentagon one. For the (H2 O) 60 and(H2 O) 120 the energy difference increases to 0.370 kJ/mol and
0.346 kJ/mol, respectively. The latter result indicates that therelative energy difference begins to decrease for the water
clusters with the size (H2 O) 120 . Zero-point vibration energy
difference was also calculated for the two (H2 O) 30 clusters,
which is 0.0251 kJ/mol. Although the hexagon water cluster
has a slightly lower zero-point vibration energy, the difference is about one order of magnitude smaller than the singlepoint energy difference between pentagon and hexagon water
clusters. Thus, at the HF/6-31G共d,p兲 level the pentagon water
clusters are more stable than the hexagon water clusters. At
the MP2/6-31G共d兲 level, the energy calculations also indicated that the pentagon water clusters are more stable thanthe hexagon water clusters. For (H2 O) 30 the hexagon cluster
is higher in energy by 0.370 kJ/mol, whereas for (H2 O) 60 thehexagon cluster is higher by 0.908 kJ/mol. The correlation
correction at the MP2 level increases the relative energy difference between the pentagon and hexagon water clusters,
but the qualitative results obtained at the Hartree–Fock levelare not changed. For infinitely long Q1D ice nanotubes, on
the other hand, our ab initio pseudopotential total-energy calculations suggest that the pentagon and hexagon ice nanotubes are likely isoenergetic. We therefore conclude that one
can start to see the relative energetic behavior of ‘‘bulk’’Q1D ice for polygon water prisms at the size (H2 O) 120 .
Finally, we examined the relative stability of three hexagon and pentagon water clusters 关 (H2 O) 60 , (H2 O) 120 and(H2 O) 180] by using three additional potential models of water: TIP5P,24 SPC/E,25 and Dang–Chang polarizable
models.26 The optimized structures for every water model
were obtained with the steepest-descent method. Zero-point
vibration energy4,27 was evaluated using normal-mode analysis, but only for TIP5P and SPC/E models. The vibrationenergy was then added to the potential energy of the corresponding optimized water cluster structure. Figure 3 also displays the energy difference between hexagon and pentagon

Downloaded 26 Mar 2003 to 128.253.229.143. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

3916

water clusters for all three models of water, respectively. For
both nonpolarizable models 共TIP5P and SPC/E兲, the hexagon
clusters were predicted to be more stable than the pentagon
clusters. However, with the Dang–Chang polarizable model,
the pentagon clusters were predicted to be slightly morestable, in agreement to the prediction based on the allelectron quantum-chemistry calculation. We thus conclude
that the inclusion of polarizibility in the water model appears
to be important in predicting the relative stability of large
water clusters.
One of the authors 共X.C.Z.兲 is grateful for valuable discussions with Professor Ken Jordan and Professor Guang-Tu
Gao. This work is supported by the U.S. NSF, by the NSFC
and Chinese Ministry of Education and Ministry of Scienceand Technology, and by the JSPS and Japan Ministry of Education. Part of the computation was performed on Research
Computing Facility at the University of Nebraska-Lincoln.
1

Bai et al.

J. Chem. Phys., Vol. 118, No. 9, 1 March 2003

M. C. Payne, M. P. Teter, D. C. Allen, T. A. Arias, and J. D. Joamopolous,
Rev. Mod. Phys. 64, 1045 共1992兲; D. K. Remler and P. A. Madden, Mol.
Phys. 70, 921 共1990兲.
2
C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev.
Lett. 69, 462 共1992兲.
3
A. F. Goncharov, V. V. Struzhkin, M. S. Somayazulu, R. J. Hemley, and H.
K. Mao, Science 273, 218 共1996兲; K. Aoki, H. Yamawaki, M. Sakashita,and H. Fujihisa, Phys. Rev. B 54, 15673 共1996兲; Ph. Pruzan, E. Wolanin,
M. Gauthier, J. C. Chervin, B. Canny, and D. Häusermann, J. Phys. Chem.
B 101, 6230 共1997兲.
4
K. Koga, R. D. Parra, H. Tanaka, and X. C. Zeng, J. Chem. Phys. 113,5037 共2000兲.

5K. Karapetian and K. D. Jordan, unpublished results on molecular mechanics study of TIP4P ice in carbon nanotubes at 0 K.
6
K. Koga, G. T. Gao, H. Tanaka, and X. C. Zeng, Nature 共London兲 412,
802 共2001兲.
7
K. Koga, G. T. Gao, H. Tanaka, and X. C. Zeng, Physica A 314, 462
共2002兲.
8
W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L.
Klein, J. Chem. Phys. 79, 926 共1983兲.
9
Y. Maniwa, H. Kataura, M. Abe, S. Suzuki, Y. Achiba, H. Kira, and K.
Matsuda, J. Phys. Soc. Jpn. 共in press兲.
10
The CASTEP software code is distributed and maintained by Accelyrs Inc.
11
J. P. Perdew and Y. Wang, Phys. Rev. B 46, 12947 共1992兲.
12
D. Vanderbilt, Phys. Rev. B 41, 7892 共1990兲.
13
H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 共1976兲.
14
G. Kresse and J. Hafner, Phys. Rev. B 47, 558 共1993兲; G. Kresse and J.
Furthmüller, ibid. 54, 11169 共1996兲.
15
D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 共1980兲.
16
D. R. Hamann, Phys. Rev. B 55, R10157 共1997兲.
17
M. Seki, K. Kobayashi, and J. Hakahara, J. Phys. Soc. Jpn. 50, 2643
共1981兲.
18
K. Kobayashi, J. Phys. Chem. 87, 4317 共1983兲.
19
V. F. Petrenko and I. A. Ryzhkin, Phys. Rev. Lett. 71, 2626 共1993兲.
20
C. J. Tsai and K. D. Jordan, J. Phys. Chem. 97, 5208 共1993兲.
21
K. Kim, K. D. Jordan, and T. S. Zwier, J. Am. Chem. Soc. 116, 11568
共1994兲.
22
K. Liu, J. D. Cruzan, and R. J. Saykally, Science 271, 929 共1996兲.
23
M. J. Frisch et al., GAUSSIAN 98, Revision A.9, Gaussian, Inc. Pittsburgh,
PA, 1998.
24
M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 共2000兲.
25
H. J. C. Brendsen, J. R. Griegera, and T. P. Straatsma, J. Phys. Chem. 91,
6269 共1987兲.
26
L. X. Dang and T. M. Chang, J. Chem. Phys. 106, 8149 共1997兲.
27
A. Pohorille, L. R. Pratt, R. A. LaViolette, M. A. Wilson, and M. A.
MacEroy, J. Chem. Phys. 87, 6070 共1987兲.

Downloaded 26 Mar 2003 to 128.253.229.143. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp