Physica A 292 (2001) 87–101

www.elsevier.com/locate/physa

Computer simulation of water–ice transition inhydrophobic nanopores
Jan Slov%aka , Hideki Tanakab; ∗ , Kenichiro Kogac ,Xiao C. Zengd

a E. Hata Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Acad. Sci.,

165 02 Prague 6, Czech Republicb Department of Chemistry, Faculty of Science, Okayama University, 3-1-1 Tsushima-naka,

Okayama 700-8530 Japanc Department of Chemistry, Fukuoka University of Education, Fukuoka 811-4192, Japand Department of Chemistry and Center for Materials & Analysis, University of Nebraska, Lincoln,

Nebraska 68588, USA

Received 23 October 2000

Abstract
A series of molecular dynamics simulations is performed in order to examine in more detailthe results of a previous simulation which shows that a thin 2lm of water, when con2ned in ahydrophobic nanopore, freezes into a bilayer ice crystal composed of two layers of hexagonalrings. Three simulations are carried out and each starts with a di4erent initial con2guration but
has the same number of molecules and the area density. Using a previously introduced solid-likecluster de2nition, we monitor the dynamic process of crystallization. We 2nd that only in one
case the con2ned water completely freezes into perfect bilayer ice whereas in other two cases,
an imperfect crystalline structure consisting of hexagons of slightly di4erent shapes is observed
and this imperfection apparently hinders the growth of perfect bilayer crystal. After adjusting
the area density to match spatial arrangements of molecules, the latter two systems are able tocrystallize completely. As a result, we obtain three forms of bilayer crystal di4ering in the area
density and hexagonal rings alignment. Further analyses of these bilayer crystals provide moreinsightful explanation on the in:uence of the boundary condition and the simulation-cell size onc 2001 Elsevier Science B.V. All rights
the diversity of possible crystallographic structures. reserved.
PACS: 64.70.Ja; 61.20.Ja
Keywords: Water; Nanopore; Nucleation; Bilayer ice∗ Corresponding author. Fax: +81-86-251-8497.

E-mail address: htanakaa@cc.okayama-u.ac.jp (H. Tanaka).

c 2001 Elsevier Science B.V. All rights reserved.
0378-4371/01/$ - see front matter 
PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 5 7 9 - 3

88

J. Slov9ak et al. / Physica A 292 (2001) 87–101

1. Introduction
Nanoscale systems such as nanoclusters, nanowires, and nanoplates are currently
attracting considerable attention because of their great potential to make important
contribution in nanotechnology. The design and preparation of nanostructured materials including, for example, quantum dots, single electron transistors [1], thin layered microstructures, reagent 2lms for biosensors, and devices for optoelectronics [2],require knowledge and control of nanoarchitectures from the very early stages ofself-organization. This requirement includes the control of nucleation growth, morphology and dissolution of crystals.
The large surface area-to-volume ratio of nanoscale materials also provides an opportunity for fundamental studies of the e4ects of nanoscale size and surface interactions on the phase transitions. In this paper, we focus on the phase transition of athin 2lm of water con2ned in hydrophobic slit nanopores by molecular dynamics (MD)computer-simulation method. Tools of computer simulation have been proved to be very
useful to obtain deeper insight into the complex behavior of con2ned :uids. In particular, computer simulations for water enable us to understand on the molecular level
the co-operativity of hydrogen bonding, solvation, and hydrogen bond network rearrangement dynamics. Moreover, computer simulations also allow us to study metastablestates of supercooled water (like low- and high-density amorphous [3]) under conditions that are either expensive or diMcult to achieve in laboratory experiments. Some
new phases and new features of metastable water have been found from computer
simulations.
Recently, we performed an MD simulation of liquid–solid phase transition of acon2ned TIP4P water [4,5] in a narrow slit composed of two parallel hydrophobic
walls [6,7]. The width of the slit is about one nanometer, just enough to accommodate
two layers of water molecules at certain density. When lowering the temperature of thesystem at a 2xed normal pressure, the water undergoes a phase transition from a liquid
to a solid. It is found [6] that the resulting solid phase is a crystalline ice, namely,
a bilayer ice crystal. The structure of the bilayer ice resembles none of the known
structures of the ice polymorphs. Nevertheless, the ice crystal still retains basic featuresof bulk ices. For example, every molecule is hydrogen bonded to four-nearest-neighbormolecules. Each layer of ice is composed of slightly distorted hexagonal rings (butstill very :at compared to those in ice Ih or Ic), and the two layers are completely inregistry.
The original simulations were performed at a normal pressure of 50, 150 MPa
and 1 GPa, respectively, and in all three cases, the resulting solid phase has exactly
the same bilayer crystalline structure. Computationally, it is more diMcult andtime consuming to detect the phase transition at lower normal pressures because thetransition temperature would be much lowered with decreasing the pressure. It isalso found that once the bilayer ice is formed, it will remain stable even after removing the two con2ning walls, at least within the timescale of computer
simulation.

J. Slov9ak et al. / Physica A 292 (2001) 87–101

89

In a previous study [7], e4ects of the simulation cell size and those of the initialconditions were not investigated systematically. The aim of this work is twofold: toexamine e4ects of the initial con2guration on the time evolution of the crystal phase andto elucidate the extent to which the size of simulation cell and the periodic boundary
conditions a4ect the process of crystallization and the shape of 2nal crystal.

2. Simulation and structure analysis
2.1. Simulation
We perform three MD simulations of water con2ned within two plane parallel structureless walls, at 2xed temperature (T =253 K) and 2xed normal pressure (Pz =1 GPa).
The water molecules interact with each other via TIP4P potential and with the two walls
via the 9-3 Lennard–Jones potential. The long-range water–water potential is smoothly
Q The simulation cell is a rectangular prism with lateral dimensions
truncated at 8:665 A.
Q and Ly = 65:684 A;
Q thus, keeping the area density A = N=A = N=Lx Ly
Lx = 65:01 A
2xed. The con2ned system consists of a 2xed number of water molecules N = 896.
The distance between the walls, Lz , is about 1 nm. Periodic boundary conditions are
imposed in the x and y directions. To generate distinct initial con2gurations, we use
con2gurations stored during a simulation at T = 257 K at which the system is wellequilibrated and never shows a tendency towards freezing (at least not on computer
simulation time-scale). We select three distinct con2gurations in the course of simulation, which are separated by a time lag of about 2 ns. Next, we assign new velocities to the molecules according to Boltzmann distribution and let the system develop
into equilibrium state at the same temperature. Then, the temperature is lowered to253 K and the system is let to reach a new equilibrium state again. In this way, we
generate three con2gurations that will be used as the starting con2guration for MD
simulations.
It has been shown [8,9] that the induction time and the size of critical nucleus
depend strongly on the number of particles of the system. For a three-dimensional
Lennard–Jones system, Swope and Andersen concluded that at least 15 000 atoms
are required to describe correctly the kinetics of crystal growth from a melt. For thequasi-two-dimensional con2ned water, it is unclear how large the system size should
be in order to remove artifacts due to the 2nite system size and the periodic boundary
condition. Other than e4ects of the system size, the rectangular shape of the simulation
cell would enforce a particular arrangement of molecules and thus would also play acertain role in crystallization.
2.2. Structure analysis
Every instantaneous con2guration generated in the MD simulation entails a vibrational displacement from a local minimum of the potential energy hypersurface. The

90

J. Slov9ak et al. / Physica A 292 (2001) 87–101

Fig. 1. A typical example of spatial arrangement of water molecules in the clusters based on our cluster
de2nition. The solid lines depict the hydrogen bonds in each layer, the dashed ones are HB connecting
molecules in upper and lower layers.

potential energy of an instantaneous con2guration (hereafter called I-structure) has
two contributions: the potential energy of the minimum-energy structure, and the thermal excitations. The minimum-energy structures are referred to as quenched structures(Q-structures [10 –12]) and can be obtained using the steepest descent method. The
Q-structures provide useful information about the structure in con2guration space andthe geometry of molecular arrangement in solid state.
We characterize the solid state and solid-like regions based on the appearance of thedomain of hexagonal cells. The hexagonal cell is considered to be an elementary unitof the bilayer crystalline structure. It is composed of a pair of coupled hexagonal ringsof hydrogen bonded molecules, i.e., a hexagonal ring in one layer with its counterpart
in another. A single hexagonal cell is the smallest unit having solid-like features andwe will regard it as the smallest embryo of the solid phase. The crystal grows via
attaching more and more hexagonal cells. We call the two hexagonal cells neighboring
if they share four molecules and connected if there is a chain of neighboring cells.
Thus, a solid cluster is de2ned as a maximally connected set of hexagonal cells. Fig. 1
shows an example of a con2guration during freezing process, which contains fourclusters of sizes 8, 5, 1. With this de2nition of solid cluster, we need not impose any
restrictions on the same shape of hexagons and their mutual arrangement.
3. Results and discussion
In a previous study, we performed an MD simulation for a system of 896 watermolecules con2ned to a hydrophobic slit at T = 253 K [7]. Here, we perform two additional MD simulations at the same temperature and with the same number of molecules
but starting from di4erent initial con2gurations of molecules. For the sake of analysis,

J. Slov9ak et al. / Physica A 292 (2001) 87–101

91

Fig. 2. The maximum cluster size as a function of time for simulation S1 (a) and simulations S2 (solidline) and S3 (dashed line) (b).

we store every thousandth con2guration with the time step being 5 × 10−16 s. These
instantaneous con2gurations are subsequently quenched and the resulting Q-structures
are used for structural analysis. Hereafter, we will refer to our three simulations as S1,S2 and S3.
During the simulations, we monitor the size distribution of solid-like clusters. Fig. 2
shows the time evolution of the maximal cluster size for the three independent

92

J. Slov9ak et al. / Physica A 292 (2001) 87–101

Fig. 3. The 2nal crystalline structures for simulation S1 (a), simulation S2 (b) and simulation S3 (c) (xyprojection).

simulations. Fig. 2(a) depicts the result of S1 [7]. The con2ned water gradually freezes
and 2nally forms a crystal composed of 224 hexagonal cells. The crystal appears at
about t = 10:5 ns and apparently is stable during the next 2.5 ns, only occasionally
breaking up few of the hydrogen bonds. A snapshot of molecular arrangement in oneof the 2nal con2gurations is shown in Fig. 3(a). The crystallization processes in S2and S3 are found to be quite di4erent from that in S1. We can see in Fig. 2(b) thatthe solid embryo in S2 grows very quickly during the 2rst 6 –7 ns but it never ends up
in a complete crystal within the simulation time scale. A closer look into the con2gurations reveals that the hexagonal cells line up in slightly di4erent directions comparedto those in S1 and have slightly di4erent shapes: the hexagons do not have exactly
the same inner angles as those in S1. This yields heterogeneity in spatial distribution

J. Slov9ak et al. / Physica A 292 (2001) 87–101

93

of molecules. In fact, there is a liquid-like hole in the solid-like system, which consumes too many molecules to generate the missing hexagonal cells in the spot. Thus,
formation of a perfect bilayer crystal cannot be completed unless some major structural rearrangements take place in the solid region. However, the hole is too small
to cause serious damage in the already formed crystalline structure such as complete
disintegration of solid region or change in the alignment of hexagons. Once the rowsof hexagonal cells are connected over the simulation cell, it is very diMcult to breakthis connection. If certain number of molecules are removed from the liquid-like hole,
the process of crystallization would complete successfully. In the case of S2, we observed a crystal composed of 218 hexagonal cells (872 molecules) which seems to be
perfectly stable after its formation. In the case of S3, a liquid-like hole appears again
but it comprises less molecules than those needed to form the missing hexagonal cells.
Thus, to complete the crystallization, additional molecules should be put into the hole.
Adding is technically more diMcult than removing molecules and we will not comeinto too much details. After adding certain number of molecules, we obtain a crystal
composed of 228 hexagonal cells (912 molecules) which again is very stable in thecourse of the simulation.
Fig. 2(b) shows the time dependence of maximal cluster size in S2 and S3. The
number of molecules is adjusted at t =16 ns for S2 (from 896 to 872) and at t =14:5 ns
for S3 (from 896 to 912). We observe quite a sharp decrease in the maximal cluster
size at t = 14:5 ns for S3. It depends on the way how the additional molecules are put
into the system. In S2 it is possible to remove some molecules from the hole and thus
the maximal cluster size remains the same whereas in S3 we have to ‘blow up’ thehole, that is, to reduce the already frozen region. In this way, we are able to change themaximal cluster size. Fig. 3(b) and (c) show the xy-projections of the 2nal crystals.
Fig. 4 is given to provide a better picture of how the molecules are arranged. The
picture corresponds to S2 but looks nearly same as those in S1 and S3.
One may ask the way we handle the system size because it is unnatural to suddenly
change the number of particles in the system. We remark that the main purpose toemploy this method is to be able to prepare various solid states and then analyze
them. We 2nd that at a certain point of simulation the system will not form a regular
crystal within the time scale of simulation. Once the freezing process is stuck, adjusting
the number of molecules and restarting the simulations o4ers a quick way, thoughunnatural, to induce formation of regular crystals. We can show that the crystals thus
obtained are indeed stable so that their existence is justi2ed despite the use of unnatural
means.
Another question to be asked is: do the molecules align in a certain way from thevery beginning of a simulation? The answer is very unlikely. Then what is the decisive
point in a simulation at which the freezing takes place? It appears that the boundary
condition play a very important role in crystallization and this somewhat dictates the
crystalline forms that 2t into the simulation cell. If the surface area of the simulation cell is not 2xed, would the results be di4erent?, e.g., using an NPT ensemble?
Would the results be di4erent if di4erent boundary conditions are chosen? To answer

94

J. Slov9ak et al. / Physica A 292 (2001) 87–101

Fig. 4. A typical example of spatial arrangement of molecules in a 2nal con2guration in simulation S2.

these questions, we have to conduct a more detailed structural analysis of the resulting
crystals.
First, let us recall the bilayer ice rules formulated previously [6]. Note that in theQ
TIP4P bilayer ice, each hexagon has a side length (O–O separation) of 2:73 ± 0:02 A◦◦◦
and three di4erent angles:  = 108 ;  = 117:5 , and  = 134:5 (all within the error◦
bar ±2 ). The 2rst ice rule is that the HOH angle is never superposed over the  and angles. The second rule is that when one OH bond of a water molecule is normal
to the hexagonal lattice plane, the other OH bond can point only along a direction
adjacent to the  angle edges. All the hexagons are congruent but with two di4erent
chiralities A and B. The hexagons with the same chirality always line up in a row. The
bilayer crystal is composed of alternating rows of A’s and B’s. Under the two bilayer
ice rules, the total number of possible arrangements of molecular orientations can becounted exactly, which is W bilayer = 2N=4 . Thus, the residual entropy of the bilayer ice
crystal is given by
S bilayer = kB In W bilayer = kB ln 2N=4 ;(1)

J. Slov9ak et al. / Physica A 292 (2001) 87–101

95

Fig. 5. Geometry of the hexagonal rings in the bilayer ice crystal.
Table 1
Arrangement and geometry of hexagonal rings in a crystal. A and B denote two types of hexagons which
generally occur in the resulting structure. l stands for longer sides and s for shorter sides of hexagons
Simulation

A

A

A

B

B

B

l

s

I
II
III

117.5
112.5
115.5

108.5
110.0
105.5

134.5
137.5
139.0

108.0
115.0
107.5

117.5
112.5
115.5

134.5
132.5
137.0

2.73
2.78
2.72

2.73
2.75
2.72

where kB is Boltzmann constant. Indeed, we 2nd these rules also hold for the bilayer
crystals generated in S2 and S3 when the crystallization completes. In all three cases,
the crystals is composed of alternating rows of type-A hexagons and type-B hexagons.
The two hexagons di4er not only in chirality but also in the inner angles and side
lengths. The direction of the rows is also di4erent. But there are always only two
types of hexagons. Fig. 5 schematically shows the geometry of the two hexagons. The
numerical values of inner angles and side lengths are given in Table 1. Another way
to state the two bilayer ice rules can be as follows. (1) The HOH angle can only besuperposed over angles A ; B ; (2) when one OH bond of water molecule is normal
to the hexagonal lattice plane, the other OH bond can point only along a direction
adjacent to the A and B angle edges. In the simulations, we actually observed some
exceptions from the bilayer ice rules. Since there are very few of them (less than 5defects in the whole crystal in S1 and S2 and no defects in S3), we consider these
exceptions non-essential. The residual entropy can be estimated with high precisionusing Eq. (1).
Note also that not every hexagonal ring has the exact shape as depicted in Fig. 5,
which is averaged over large number of con2gurations. All the angles given in◦
Q The distriTable 1 are within error ±2 and the side lengths within error ±0:02 A.
bution of inner angles of hexagons for all the three structures is shown in Fig. 6
and the comparison of O–O distance distributions is shown in Fig. 7. Figs. 8 and 9
show the binding-energy distributions and pair interaction-energy distribution. They

96

J. Slov9ak et al. / Physica A 292 (2001) 87–101

Fig. 6. Comparison of hexagon inner-angle distributions in simulation S1 (solid line), S2 (dashed line) andS3 (dotted line).

Fig. 7. Comparison of O–O distance distribution in simulation S1 (solid line), S2 (dashed line) and S3(dotted line).

J. Slov9ak et al. / Physica A 292 (2001) 87–101

97

Fig. 8. Comparison of the binding energy distributions in simulation S1 (solid line), S2 (dashed line) andS3 (dotted line).

Fig. 9. Comparison of the pair interaction energy distributions in simulation S1 (solid line), S2 (dashed line)
and S3 (dotted line).

98

J. Slov9ak et al. / Physica A 292 (2001) 87–101

are quite similar. The average potential energy per molecule is E1 = −52:79 kJ=mol for
S1, E2 = −52:59 kJ=mol for S2, and E3 = −52:92 kJ=mol for S3. Although both S1and S2 crystals are more or less regular, the S1 crystal possesses more imperfections.
Every two hexagons in a row of hexagonal cells are translationally invariant in
S2, S3 crystals. This is not the case for S1 crystal, which can be clearly seen from
Fig. 3(a).
We also found all three structures, once formed, seem to be very stable. It seems
that there is certain :exibility in O–O–O angles of hydrogen bonded triplets, as wellas in O–O distance of a pair of hydrogen bonded molecules. Potentially there areseveral possible directions in which rows of hexagons can be lined up. It would be anenormous task to classify all such structures but apparently there are many of them,each one di4ering in hexagonal shapes and even hydrogen bonded pair distances.
Our simulations indicate that for a given simulation-cell size, there are certain ‘magic’
numbers of molecules for which some special initial liquid con2gurations exist and will
freeze into regular crystalline structure without any defects. There is a certain number of
non-equivalent crystalline structures in the simulation cell of the given size. Increasing
the size of simulation cell, the total number of crystalline structures that 2t into the cellis expected to be larger. For illustration: If we increase the lateral dimensions of thesimulation cell 2 times (i.e., the area 4 times), then by taking 4 replicas of the smaller
cell we can immediately generate one crystalline structure in the larger cell. But notvice versa. In the larger cell, there are generally other crystalline structures which do
not exist in the smaller one. If we increase signi2cantly the size of simulation cell,the number of possible structures will increase but as for their hexagons, they will becloser in shape, so the di4erent structures would show up mainly due to the di4erentalignment of rows of hexagons.
Let us take a closer look at the early stages of crystallization. In the early stages ofcrystallization our analysis shows that the :exibility in hexagon angles is much higher
than that in O–O distance. Fig. 10 gives an example of an early con2guration in S2.
The alignment of hexagonal cells in one part of a cluster resembles that in one of S1
crystal and in another part resembles that in one of S2 crystal. At t = 6 − 7 ns theembryo of crystal begins to look like a portion cut from the 2nal crystal, i.e., at about
this time the alignment of hexagonal rings is established. Are their shapes establishedat 2nal stages as well? The angle analysis shows that even at the end of S2 and justbefore removing molecules from the system, there exist only two distinct angles inthe distribution, which means that there appears only one type of hexagon with inner◦◦◦angles  = 111 ,  = 111 and  = 138 . This distribution is a result of averaging,
and actually all the peaks in the distribution are very broad and not all the hexagons
have the same angles. There is quite a large portion of hexagons with all the three
angles di4erent. But the number of hexagons with  =  slightly prevails, which gives◦
the broad peak around the value of 111 . On the other hand, at the end of regularS3 simulation, the angle distribution is exactly the same as that one of S1. The angle
distributions become those depicted in Fig. 6 only after the number of molecules inthe system is adjusted.

J. Slov9ak et al. / Physica A 292 (2001) 87–101

99

Fig. 10. An example of an intermediate state during the simulation S2 (t = 4 ns).

If we let a system crystallize at an arbitrary density, generally we will obtaina crystal with some defects. These defects can move around throughout the system, change their size, and interact with each other. But with the periodic boundary condition, they will never reach a boundary and travel out of the system. If a
free boundary condition is employed, the defects could reach a boundary after certain time and take part in extending crystalline structure at the boundary after somemovement.
Is that possible to reduce the in:uence of periodic boundary condition using someother simulation method? To what extent does that method be appropriate? Further
studies are needed to answer these questions. However, we can make a few generalcomments here to elucidate the physics behind the simulations. Our system is inhomogenous so that instead of bulk pressure, in general, we need to consider stresses. We
can still de2ne the lateral pressure Px (=Py ) ≡ L−1
z (@F=@A)NTLz and the normal pressure
Pz ≡ A−1 (@F=@Lz )NTA . Here, F is the Helmholtz free energy and A the surface areaof the slit. The volume is then V = ALz . In the case of a rectangular simulation cell,
A = Lx Ly . Because of the strong system anisotropy, generally Px = Pz . Note that in realexperiments for the study of phase transition, one usually controls stresses. To mimic
the real experiments, an NPx Pz T simulation seems a better choice. In fact, the quantity
playing the role of bulk pressure is Px . This is because the free energy is extensivewith respect to A, i.e.,F(N; A; Lz ; T ) = Af(N=A; Lz ; T )

(2)

100

J. Slov9ak et al. / Physica A 292 (2001) 87–101

but not with respect to Lz . This shows that Px and A are quantities analogous to the bulk
pressure and bulk volume for bulk system. In this work we use NAPz T simulation. The
constraints can be implemented in computer simulation but they are not very natural.
More realistic constraints would partially reduce the problems associated with periodic
boundary condition and small-system size, although Px may be di4erent from the bulk
pressure in equilibrium with water between two walls.

4. Concluding remarks
A series of three MD simulations for a thin 2lm of water con2ned to a hydrophobic
slit is performed in order to examine the e4ects of the initial con2gurations on theevolution process of the crystal embryo and to study the in:uence of simulation cellsize and periodic boundary conditions on crystallization process. Upon freezing thesystem forms a bilayer ice crystal composed of hexagonal cells. In general, two typesof hexagonal cells are seen, each aligning in rows with di4erent orientations. In one
simulation, the system freezes completely while in other two simulations incomplete
structures containing a liquid-like region show up. These regions contain less than
100 molecules, i.e., less than one-tenth of the total number (896 molecules) in thesystem. It is conceivable to make the system crystallize completely by altering thenumber of molecules in the system. Once formed, the incomplete crystals do not change
signi2cantly. The direction in which the rows of hexagonal cells line up also stays 2xed
mainly due to the cell shape and the periodic boundary condition.
The in:uence of periodic boundary condition on nucleation has been reported manytimes for systems with simple intermolecular interaction (e.g. [13 –15]). It has been
pointed out that for smaller systems, the observed rapid crystallization after some induction time is not related to the formation of critical nucleus but just an artifact of
periodic boundary condition. For three-dimensional systems, the minimal system size
to avoid the artifacts has been estimated to be about 15 000 particles. A main di4erence between our system and previously studied simple systems is that we do not have
all the information about the bilayer crystalline structure in a quantitative fashion. For
Lennard–Jones solid, it is well known that one possible crystalline structure is fcc sothat one can clearly identify the defects in crystalline structure after the crystallization is completed. For bilayer ice, however, the inner angles of hexagons often di4er
marginally. A question naturally arises: Is the A type hexagon in S2 the ‘right’ oneor the A type in S3 the ‘right’ one? Unless a very long simulation is performed for alarge enough bilayer ice crystal with either free boundary condition or with periodic
boundary condition, we do not know the exact values of inner hexagon angles. A detailed study of the nucleation and a reliable estimation of critical nucleus size must
await a much larger scale simulation. The system size and shape to be employed also
depend on the size of critical nucleus since the smaller the critical nucleus the lessin:uence the cell size and shape would have.

J. Slov9ak et al. / Physica A 292 (2001) 87–101

101

Acknowledgements
HT and KK are supported by grants from JSPS and Ministry of Education and XCZ
was supported by grants from NSF, ONR, and a fellowship from JSPS.
References[1] A.P. Alivisatos, Science 271 (1996) 933.
[2] H. Rapaport, J. Phys. Chem. B 104 (2000) 1399.
[3] O. Mishima, H.E. Stanley, Nature 396 (1998) 329.
[4] W.L. Jorgensen, J. Chem. Phys. 77 (1982) 4156.
[5] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983)
926.
[6] K. Koga, X.C. Zeng, H. Tanaka, Phys. Rev. Lett. 79 (1997) 5262.
[7] J. Slov%ak, K. Koga, H. Tanaka, X.C. Zeng, Phys. Rev. E 60 (1999) 5833.
[8] J.D. Honeycutt, H.C. Andersen, J. Phys. Chem. 90 (1986) 1585.
[9] W.C. Swope, H.C. Andersen, Phys. Rev. B 41 (1990) 7042.
[10] F.H. Stillinger, T.A. Weber, J. Phys. Chem. 87 (1983) 2833.
[11] F.H. Stillinger, T.A. Weber, Science 225 (1984) 983.
[12] I. Ohmine, H. Tanaka, P.G. Wolynes, J. Chem. Phys. 89 (1988) 5852.
[13] B.J. Alder, T.E. Wainwrigt, J. Chem. Phys. 33 (1960) 1439.
[14] M.J. Mandell, J.P. McTague, J. Chem. Phys. 66 (1977) 3070.
[15] J.D. Honeycutt, H.C. Andersen, Chem. Phys. Lett. 108 (1984) 535.