Case Study: Water and Ice

Timothy A. Isgro, Marcos Sotomayor, and Eduardo Cruz-Chu

Illustration by Hirohito Ogasawara, from
http://www-ssrl.slac.stanford.edu/images/water4ssrl-large.jpg.
1 The Universal Solvent
Water is essential for sustaining life on Earth. Almost 75% of the Earths
surface is covered by it. It composes roughly 70% of the human body by
mass [1]. It is the medium associated with nearly all microscopic life pro-
cesses. Much of the reason that water can sustain life is due to its unique
properties.
Among the most essential and extreme properties of water is its capabil-
ity to absorb large amounts of heat. The heat capacity of water, which is the
1
highest for compounds of its type in the liquid state, measures the amount of
heat which needs to be added to water to change its temperature by a given
amount. Water, thus, is able to eectively maintain its temperature even
when disturbed by great amounts of heat. This property serves to main-
tain ocean temperatures, as well as the atmospheric temperature around the
oceans. For example, when the sun rises in the morning and a large amount
of heat strikes the surface of the Earth, the vast majority of it is absorbed
by the ocean. The ocean water, however, does not exhibit a drastic increase
in temperature that could make it inhospitable to life. In constrast, when
the sun sets in the evening and that heat is taken away, the oceans do not
become too cold to harbor life. Water also has high latent heats of vapor-
ization and melting, which measure the amount of heat needed to change a
certain amount of liquid water to vapor and ice to liquid, respectively. These
two properties also aid in the maintenance of ocean temperature.
An anomalous property of water, which is especially important for life in
cold climates, is the occurrence of its maximum density at 4

C, just above
the 0

C freezing point. Therefore, as water begins to drop in temperature

when its surface is cooled, the water at the surface will reach this highest
density. The surface water will then begin to sink, displacing the warmer
water below itself. This results in the freezing of water from the surface
downward. The eect takes place in ponds and lakes that are located in cold
climates. Their surface freezes, and the water underneath remains liquid,
protecting the organisms that live there.
Water vapor also plays a key role in determining the amount of energy
from the sun which is retained on Earth. The spectrum of water vapor is such
that the molecule is transparent to most radiation in the sunlight range, but it
absorbs a considerable amount of infrared radiation at wavelengths typically
emitted by the Earth upon release of the energy provided by the sun. Thus,
water vapor is responsible for Earths greenhouse eect and provides higher
temperatures than would otherwise be expected on Earth.
Many other properties of water are no less surprising in their relationship
to the development of life. Water maintains an identity as a universal sol-
vent, that is, it may serve to transport a multitude of chemical constituents
throughout the bodies of living organisms. Water has the ability to ionize
substances in solution which is important in establishing the proper electro-
static concentrations utilized by so many biological processes. Water also
possesses an abnormally high surface tension, a property seen in the shape
of water droplets and very important in living systems. The unique proper-
2
ties of water give it incredible versatility allowing it to foster and maintain
a stable relationship with life on Earth. In this case study we will explore
some of the properties of water mentioned above. The following les will be
used:
Figure 1: Files for the case study may be found on the Desktop of your laptop in the
folder case-study-water/ or in BioCoRE as the le case-study-water.tar.gz. If it should
become necessary for you to download the les from BioCoRE, do so and untar them with
the command tar -xzvf case-study-water.tar.gz. The folder case-study-water/ containing
all les will be created in your working directory
3
2 The Liquid and Solid Phases of Water
Water exists above 100

C as a gas, above 0

C as liquid, and below 0

C
as ice. Here we describe the latter two phases that are biologically the most
relevant.
2.1 Liquid Water
The uid phase of any substance is usually the hardest to study at the
microscopic level by theoretical means. Common approximations of non-
interacting volumeless particles employed for ideal gases are not applica-
ble for liquids, and neither are the elegant theories applied to periodic and
symmetric systems found in solid phases. Experimental and theoretical de-
scriptions of liquid phases are usually limited to macroscopic properties well
described by continuous models. However, it is well known that hydrogen
bonding between water molecules is important in establishing the structure
of water, whether it be in solid or liquid form. The attraction which exists
between molecules is indeed the reason for waters high boiling point, which
is roughly 100

C higher than it would be if water was non-polar. If that were

the case, all water on Earth would exist in gaseous form.
Although the microscopic structure of liquid water is unknown, the pair
distribution function of liquid water (see Figure 14) obtained from experi-
ments and simulations suggests that at short scales, water features an or-
der that favors hydrogen bonding. However, molecules in liquid phase have
enough kinetic energy to continuously form and break the bonds. Moreover,
charged and polar solutes present in water under physiological conditions
drastically inuence the arrangement and bonding of water molecules sur-
rounding them.
Experimental techniques to determine the structure of liquid water are
lacking since each water molecule undergoes rapid rearrangement on the or-
der of femtoseconds. The need for a better understanding of water at the
microscopic level has forced the development of computational methods that
describe the structure and dynamics of individual water molecules, and many
studies have been carried out using these techniques. Many predict locally
ordered hydrogen bonding whereby rings of molecules continually form and
break. The rings may be made up of three to seven members, rather than
the six-membered rings which exist in ice [1]. The rapidly changing nature of
liquid ice structure, and of all liquids, for that matter, makes it very dicult
4
to study and characterize.
2.2 Ice
Figure 2: Five water molecules
showing the arrangement of hy-
drogen bonding in crystalline ice.
View the le 5-molecules-ice.pdb
using VMD.
Hydrogen bonding is also responsible for the
crystal structure of ice. Upon freezing at at-
mospheric pressure, water reduces the inter-
action energy between molecules by forming
a regular network in which each molecule is
hydrogen-bonded to four others. The oxygen
atom bonds to two hydrogen atoms from other
molecules and each hydrogen atom bonds to an
oxygen atom from another molecule (Figure
2). With four hydrogen bonds per molecule,
the ice crystal lattice conformation reduces the
overall energy of the water molecules. The lat-
tice also happens to be less dense than liquid
water structure. In fact, the crystal structure
of ice leaves enough room between molecules
to t entire other water molecules. The dif-
ference in density between liquid and crystal
water is visible in Fig. 3. A movie of the crystal lattice of ice being formed
is available in the le movies/ice-crystal.mpg.
The structures of ice shown in Fig. 3 represent hexagonal ice, the type of
ice that you are most familiar with, such as ice cubes in a glass or snow on
the ground. However, it is not the only type of ice known. In fact, at least 14
phases of crystal ice are known to form depending on ambient pressure and
temperature. Hexagonal ice, or ice 1h, is the form which exists at 1 atm and
0

C and is thus the type of frozen water found on Earth. As implied by its
name, it is constructed of a hexagonal lattice of water molecules, a feature
which is reected by the 6-fold symmetry of snowakes. The hexagonal
lattice is commonly identied by four axes labeled a
1
, a
2
, a
3
, and c (see
Fig. 3). Note that all three a-axes are equivalent. The six lattice faces
perpendicular to the a-axes are known as the primary prism planes and
the two faces perpendicular to the c-axis are known as basal planes. (The
lattice may also be identied by three axes called a, b, and c, whereby a =
a
1
, b = a
2
, and c = c.)
5
Figure 3: Views of the crystal structure of ice (a-c) and a view of liquid water (d). a)
Along the c-axis of the lattice (showing the basal plane) with all 3 equivalent a-axes in
the plane of the page and the c-axis coming out of the page, b) With an a-axis running
left to right (showing the primary prism plane), c) Along an a-axis, d) Liquid water at
300 K. The dierence in the density between ice and liquid water is apparent. The density
of liquid water at 4

C is 1.00 g/cm
3
, while that of ice at 0

C is 0.917 g/cm
3
[2]. View the
les section02/ICES-hex.pdb (a-c) and section02/water-300K.pdb (d) using VMD.
6
Exercise 1: Hydrogen bonding. In ice 1h, each molecule is hydrogen bonded to
four others. However, in liquid water, that number varies from molecule to molecule.
(1) Do you think the number of hydrogen bonds in liquid water should be more or
less than four? Why? Load the le water-298K-1atm.pdb, a snapshot of equilibrated
water at 298 K and 1 atm, into VMD. Count the average number of hydrogen bonds
for each molecule by running the script count-Hbonds.tcl using the source command
in the TkConsole. Make sure before hand that you are in the right directory; you may
have to rst change the directory with the command cd case-study-water/section02
in the TkConsole. Does your result match what you expected?
(2) Load the VMD saved state water-298K-1atm.vmd. Use VMD to count the
number of hydrogen bonds that the three water molecules shown in color make with
other molecules: (a) green, (b) orange, (c) yellow. Use 3.0 Angstroms and 30 degrees
as the cuto for hydrogen bonding. Give the number of hydrogen bonds for each of
the three molecules along with the resids of the molecules to which it is hydrogen
bonded.
3 Antifreeze Proteins
The ability of atomic and molecular substances to change the behavior of
water by lowering its freezing point or elevating its boiling point is evident to
most people, especially on cold winter days when salt is scattered on streets
to prevent ice from forming. Nature, too, has a means of causing the freezing
point depression of water. It uses proteins and the crystalline regularity of
ice to carry out the task.
Thermal hysteresis proteins, or antifreeze proteins (APs) as they are
sometimes referred to, enable nature the same type of control over the freez-
ing of water as that aorded to street salters. However, the manner in which
water interacts with the AP diers from salt. Whereas freezing point de-
pression via adding sodium chloride to water is a colligative property (the
eect is proportional to the number of particles added), the freezing point
depression which antifreeze proteins induce is noncolligative. The reason is
that APs act by binding to ice surfaces already formed and inhibiting the
further formation of ice, rather than by altering the properties of cold liquid
water.
To understand specically how APs work, one must relate their structure
to the crystal structure of ice. In accordance with usual growth of crystalline
7
AP Type Typical Size Secondary Amino Acid
(kDa) Structure Features
Type I 3.3-4.5 -helical ThrX
2
AsxX
7
Type II 6-14 globular disulde bridges
Type III 6-14 globular none
Type IV 12 4 amphipathic Glx-rich
-helices
Table 1: Properties of the four types of antifreeze proteins. It is interesting to note that
the types exhibit notable variety in structure and composition, yet perform essentially the
same task.
solids, normal ice growth takes place along the axis with the highest atomic
density, the a-axes in this case (see Fig. 3b). APs take advantage of this fact
by binding to faces, or planes, of ice which intersect the a-axes, restricting
ice growth.
Figure 4: Ice crystal formation
in a glass capillary in the pres-
ence of winter ounder type I an-
tifreeze protein. The shape re-
ects the binding of winter oun-
der AP to the pyramidal {2021}
plane of ice. Picture from [3].
APs are found in many cold-surviving ani-
mals, such as sh and insects. Fish APs have
been the most widely studied both experimen-
tally and computationally. All APs are cate-
gorized into one of four types (Table 1). What
makes them perhaps most interesting is that
while all APs perform the same basic task,
each type is quite dierent in structure. We
will use VMD to investigate a type I AP, and
then look at a possible ice recognition motif
for a type III AP.
Type I APs are typically small proteins,
roughly 3.3-4.5 kDa in size and composed of -
helical secondary structure. Most type I APs
consist of the amino acid repeat ThrX
2
AsxX
7
which binds to ice. Based on ice etching experiments, it has been determined
that HLPC6, a winter ounder type I AP, typically binds to the {2021} plane
(equivalently {201} in a, b, c-axis notation) of ice 1h along the 0112 direc-
tion in which the ice oxygen atoms have a repeat distance of approximately
16.7

A [4]. Figure 4 shows how macroscopic ice formation is aected by the
presence of winter ounder AP. The eect of ounder AP on the formation
8
Figure 5: Winter ounder antifreeze protein HPLC6 bound to the {2021} plane of ice.
a) A {2021} plane cut into the ice crystal. The a-axis points out of the page as in Figure
3c. b) HPLC6 bound to the surface. Note the regular spacing of Thr residues shown in
vdW representation, which match the ice lattice spacing. The oxygen atoms to which Thr
residues bind are shown as black spheres. c) A front view of the binding. The a-axis
runs left to right parallel to the page. View the le section03/winter-ounder-lattice.pdb
and the saved state section03/winter-ounder-lattice.vmd using VMD.
of ice crystals reects the atomic properties of the proteins binding to the ice
surface.
Indeed, the proteins sequence and structure seems to t the {2021} motif
quite nicely. HPLC6 is a single -helix with threonine residues lining one side
of the helix and having a distance of approximately 16.7

A between them.
It is thought that the protein uses this perfect structural complementarity
to ice to bind to the surface of the crystal and inhibit its growth [5]. Figure
5 shows the protein bound to a {2021} plane of the ice lattice. Note the
regularity of the threonine residues which is a perfect complement to the
lattice periodicity. The transparent surface of the protein, which is colored
by residue type, also displays an interesting feature of the protein. One
side of the protein is lined very regularly with hydrophilic residues (white
is hydrophobic), as mentioned, which may lead one to believe that binding
occurs via these residues. However, mutation experiments and computer
9
simulations have led to a second hypothesis which states that the hydrophobic
residues along the helical surface are the ones that are actually important
in the protein recognizing the ice surface [3, 6]. Only after such recognition
may the hydrophilic residues bind to the crystal. The drastic contradiction
between the two models poses a very fundamental question which remains
unanswered: Is ice more hydrophobic than water?. Computational and
experimental studies being conducted today may provide the answer.
Figure 6: Ocean eel pout type
III antifreeze protein showing sur-
face complemetarity to the pri-
mary prism face of ice. Residues
N14, T18, and Q44 are shown
in vdW representation. View the
le section03/eel-ice.pdb and the
saved state section03/eel-ice.vmd us-
ing VMD.
Type III APs dier from their type I
counterparts in several ways. They are typ-
ically larger, 6-14 kDa in size, globular in
structure, do not contain a dominant amino
acid residue like type I APs (Ala), and
lack any apparent repeat sequence of amino
acids. As a result, it is very dicult to pre-
dict how a given type III AP will bind to
ice. It is known from ice etching experi-
ments, however, that the most prominent
face to which type III APs bind is the pri-
mary prism plane {100} of ice [5], which
is shown in Fig. 3b. The antifreeze pro-
tein from ocean eel pout has been classi-
ed as a type III AP, and conserved hy-
drophilic residues N14, T18, and Q44 have
been identied as essential to antifreeze ac-
tivity. These residues make up a face of
the protein which is unusually at. As in
the case of type I APs, surface complemen-
tarity with ice is believed to be important,
too. By examining the crystal structure of
ocean eel pout AP, in relation to the pri-
mary prism crystal plane of ice, we are able
to see the surface complementarity estab-
lished by the three residues (see Fig. 6).
It may also be interesting to note that the
fourth residue (not highlighted) which adds
to the surface complementarity is A16. A
simple look at the molecules surface when paired with ice may be enough
to identify the ice binding mechanism. Indeed, energy minimization calcu-
10
lations of ocean eel pout have established that low energy congurations of
the protein-ice surface exist when lattice matching is established, just as
in type I APs [5].
While seeing a possible ice binding motif is interesting in its own right,
researchers are also interested in exactly how the protein recognizes and
attaches itself to ice, that is, the binding mechanism. In fact, ocean eel
pout AP induces steps to form in the ice surface, leading to the hypothesis
that residue N14 initiates binding of the protein to the crevice between basal
{001} and primary prism planes {100} of ice [7]. Molecular dynamics is an
ideal method to test such a hypothesis, but in order to do so one must know
how to represent, or model, the water which our protein will be solvated
in and binding to. Before giving a description of these kind of microscopic
models, we will explore in the next section some properties of liquid water.
Exercise 2: Type I AP binding. The originally proposed [8] binding mechanism for
type I AP HPLC6 had the protein binding to the primary prism plane of ice. This may
have been the most obvious choice, as ice grows perpendicular to this plane, but ice
etching experiments [4] proved it to be wrong. Try to replicate the thought process
of early AP scientists by guring out the best possible docking of HPLC6 on the
primary prism face of ice. The le 1WFA-prism-sheet.pdb contains a sheet of primary
prism face ice along with HPLC6. Use VMD translate and rotate commands in the
TkConsole, such as $HLPC selection moveby {1 2 5} and $HLPC selection move
[trans axis x 40 deg], to move the protein into a hypothetical docked position on the
ice, which takes advantage of the proteins threonine repeat distance. Note that the
t will not be as good as for {2021} binding. Emphasize the ice oxygen atoms (as
vdW spheres, for instance) with which the protein make contact. Create two VMD
snapsots from dierent angles.
4 The Wettability of Water
The importance of water in sustaining life processes is very apparent. On
the nanoscopic level, living organisms have evolved proteins which are pre-
cisely tuned to take advantage of a water environment. For cold-dwelling
organisms, the interaction between water and its antifreeze proteins prevents
freezing in extreme conditions. On the macroscopic level, some life forms
have evolved features to take advantage of another interesting property of
water: its high surface tension.
11
The wetting behavior of a surface, its ability to interact with a liquid,
can be characterized by the surface tension of the liquid as well as the features
of the solid surface. Qualitatively, good wetting behavior is observed if a
liquid spreads widely over a certain surface, and non-wetting behavior is
observed if the liquid tends to avoid the surface, forming small droplets. We
will begin with an investigation of the physics behind surface wetting and
later see how one insect has evolved features to benet from waters surface
tension in an exquisite and surprising way.
Figure 7: A picture of a water strider standing on the surface of water.
Note the extensive surface coverage of the back legs. Picture from http://www-
math.mit.edu/%7Edhu/Striderweb/adultside.jpg
4.1 The physics of surface tension
In order to understand the wetting behavior of liquids, we need to understand
the microscopic molecular organization of a homogeneous liquid. In a liquid
state, molecules are free to move in relation to one another but at the same
time are in close contact. Due to their close proximity, intermolecular forces
between molecules, which are typically irrelevant in gases, play a major role
in establishing the properties of the liquid.
A closer view into a liquid structure reveals that not all molecules within
the liquid are under the same inuence. In general, we can separate the
molecules into two groups: those moving freely in the bulk while others are
12
transiently captured at the surface. This is because each molecule at the
surface has interacting neighbors only at one side, and the resulting force is
in the direction of the bulk. The energy of a surface molecule is therefore
higher than that of a bulk molecule, and energy must be expended to move
a molecule from the interior to the surface. This energy distribution on the
surface is responsible for the phenomenon known as surface tension.
Surface tension is measured in units of force per length (typically dynes/cm),
and it quanties the amount of force required to break a liquid lm of given
length. An equally meaningful quantity is the surface free energy, mea-
sured in units of energy per area (typically ergs/cm
2
). Even though both
names qualitatively describe the same phenomenon, the term surface ten-
sion is older and consequently used more [9] . Water has a surface tension
of 72.8 dynes/cm at 20

C and 1 atm.
To minimize the free energy of the system, liquids tend to form shapes
that reduce their surface. As a result, drops and bubbles are spherical, since
the sphere has the lowest surface/volume ratio. Even though gravity breaks
that ideal geometry, causing liquids to t into the container shape, the nal
form always tends to minimize the exposed surface [10] .
By 1800, Laplace realized that the spherical shape also leads to a pressure
dierence between the two sides of the surface and stated a simple equation
that correlates pressure dierence, surface tension, and the radius of the
sphere:
P
in
= P
ex
+ 2
T
r
(1)
where P
in
is the internal pressure, P
ex
is the external pressure, T is the
surface tension, and r is the radius of the sphere [9]. The equation states
that the pressure inside of a curved surface is always higher than the pressure
outside. The pressure dierence is cancelled when the curvature ratio tends
to be innite, i.e., the surface is at.
4.2 Wetting a surface
When a liquid is in contact with an inert solid phase, the liquid wets the
surface. Liquid molecules at the solid-liquid interface are now in a dier-
ent environment than the ones that are either in the bulk or at the exposed
surface. Those molecules feel two kinds of forces: cohesive forces acting be-
tween like molecules and adhesive forces acting between dierent molecules.
13
Figure 8: The wetting behavior of a water droplet depends on the surface characteristics.
Over a hydrophilic surface, the water droplet spreads. Over a hydrophobic surface, it
tries to avoid contact and forms a bubble. This gure was created from the trajectories
associated with Fig. 10.
The balance between cohesive and adhesive forces determines the wetting
properties of the surface (Fig. 8). When the cohesive forces of water can
be counterbalanced by the adhesive forces of the substrate, a liquid droplet
tends to spreads over the surface. When the cohesive forces of water are
stronger than the adhesive forces of the substrate, the droplet tries to avoid
the surface, keeping its spherical shape and reducing the surface tension.
Think about a tiny water droplet resting on an inert solid surface, such
as your desk. Within the droplet, water molecules are held together by
cohesive forces. At the interface between the water droplet and the desk,
adhesive forces emerge. However, cohesive forces are stronger and surface
tension will still hold the droplet in a roughly spherical shape (unless your
desk is quite exotic).
4.3 Water Contact Angle
We now introduce a method to evaluate the balance between the solid surface-
liquid water interactions by using one of the most sensitive, beautiful, and
cheapest instruments in the world: a small water drop.
Consider a water drop on a smooth and clean surface. The droplet will
either spread completely over the surface, partially, or not at all [11]. This
14
Figure 9: Water droplet resting on an inert surface. The contact angle is dened by the
angle between the solid surface and the tangent to the surface liquid at the solid-liquid
interface. View the movie section04/WCA.mpg.
15
behavior may be quantied by the water contact angle (WCA), the angle
between the tangent to the surface of the liquid and the tangent to the
surface of the solid at the solid-liquid interface (Fig. 9). If the contact angle
between a liquid and a solid is close to zero, the surface is considered wet.
If the angle is close to or greater than 90

, the liquid does not wet the

surface. Wetting is analogous to hydrophobicity in that it reects the degree
of interaction of certain molecules with water. For a WCA greater than 150

,
the surface is considered superhydrophobic.
Macroscopically, the water contact angle of a droplet in equilibrium with
the surface is described by Youngs equation:
cos =

SV

SL

LV
(2)
where is the contact angle, and the symbols denote the surface tensions
of each interface: solid-vapor, solid-liquid, liquid-vapor. Youngs equation
was stated in 1805 and describes the equilibrium relationship between the
interfaces [12].
Microscopically, the water contact angle depends on the solids surface
topology as well as the specic solid-water interactions. We will explore
one of the most popular surfaces, silicates, also known as glasses. Exposed
functional groups establish the accessible regions that interact with other
molecules and thus determine the silica wetting properties. At the silica
surface, we have two functional groups: siloxanes ( -O-Si-O- ) and silanols
( -Si-OH ). The silanol group, with its hydroxyl group, is the dominant
factor for hydrophilicity, attracting water to the surface. As more silanols
are exposed, the surface becomes more hydrophilic. Conversely, a complete
dehydroxylated silica surface, without silanols, is less hydrophilic, and would
exhibit a WCA of around 44

[13].
Molecular dynamics allows one to explore the behavior of a silica surface
in contact with water (Fig. 10). One can use the experimental values for the
silica WCA to calibrate silica force elds, the set of equations and parameters
that dene the molecular interactions. Once the force eld is calibrated, one
can explore many scientic questions involving water-silica interactions, like
water permeation (Fig. 11).
16
Figure 10: Force eld (FF) parameters determine the wetting behavior in molecular
dynamics simulations. For the same system two dierent silica force elds were evaluated.
The top gure corresponds to the starting set-up: 699 water molecules in a cube of 18

A
and a dehydroxylated silica slab of 595920

A. The silica for the two dierent force elds
are pictured blue [14] and green [15]. Snapshots correspond to dierent time during the
simulations: 20, 50 and 100 ps. View the simulations by loading the le section04/Pore.psf
and either section04/HydroPHOBIC.Pore.dcd or section04/HydroPHILIC.Pore.dcd into it
using VMD.
4.4 Walking on water
The mechanical motion of the water strider is a charming example of how
nature exploits the surface tension of water to an advantage. Water striders
are small insects that live on the surface of ponds and quiet rivers. They
are typically a few centimeters in length, although some varieties like Gi-
gantometra gigas have legs up to 20 cm long. As their name suggests, these
insects can walk on water! The surface of water is quite an unsuitable place
for many insects, so water striders take advantage of this to populate, feed,
and breed themselves.
Three pairs of long hydrophobic legs are the secret of the striders mo-
17
Figure 11: Permeation of water through a silica nanopore for two dierent force elds
(c.f. Fig. 10) [14, 15]. The top snapshot shows the initial set-up of the system: a silica
pore of 20

A diameter and 37

A height and two water boxes, each box with 6062 water
molecules. Remaining snapshots correspond to dierent times during the simulations,
namely, at 50, 100 and 190 ps. View the simulations by loading the le section04/Surf.psf
and either section04/HydroPHOBIC.Surf.dcd or section04/HydroPHILIC.Surf.dcd into it
using VMD.
bility. Precisely engineered to support their light weight, these legs have
thousands of little hairs which are covered in wax. A careful combination
of topology and structural features in the legs allow these insects to stand
and walk without piercing the water. Recent studies [16] proposed that
the particular arrangement of hair at the legs surface is the key factor for
understanding their exceptional hydrophobicity. Accurate measurements of
the WCA reveals a superhydrophobic surface, with an angle of 167.6

!
Standing on the water surface is a basic mechanics problem. In order to
remain aoat, the weight of the strider must be supported by a force in the
opposite direction. That compensating force can come from two sources: the
buoyancy force like in a swimming boat (via Archimedes principle) and/or
surface tension. In the case of the water strider, the weight is essentially
counterbalanced by surface tension, and buoyancy is irrelevant, since the
volume of water displaced by the strider is very small. As the strider stands
on the surface, their hydrophobic legs push the water downwards, generating
a concave curvature surface around each leg. Thus, the vertical component
of the surface tension force is in the upward direction and is proportional to
18
the surface area of the water-leg interface.
Water striders propel themselves using their middle legs, achieving speeds
of 1.5 m/s [17], while the back legs are used to steer and brake. The middle
legs can apply force ve times larger than the striders weight. Conveniently,
a single superhydrophobic leg can stand 15 times the striders weight before
piercing the surface [16].
Superhydrophobic surfaces are not exclusive properties of water striders.
Lotus and rice leaves, to keep them from getting soaked, also have superhy-
drophobic surfaces [18], formed by capturing air to avoid contact with water.
Nature has tuned the properties of these surfaces with a meticulous com-
bination of non-wetting materials and nanostructures to the benet of the
organism.
5 Modeling Water at the Molecular Level
We have learned already about the relevance of water in many aspects of our
everyday life. In particular, its role in biological processes is fundamental,
but sometimes not completely understood at the atomic level. How can one
explore the microscopic behavior of water and thereby predict or explain its
role in biological systems? Experiments, of course, allow the measurement of
many properties of water like density, diusion coecient, heat capacity, and
melting and boiling temperatures. However, theoretical and computational
modeling, especially taking into account current advances in computer sim-
ulations, provides insight into systems in which water plays a fundamental
role at the microscopic level.
The computer is becoming a microscope that can see scales experiments
cannot reach, but, of course, the accuracy of the computational description
depends on the underlying theoretical model of water. A model needs to re-
produce the known behavior of water and, at the same time, predict unknown
properties. While continuum models eectively reproduce macroscopic prop-
erties of water, the discrete nature of water molecules and their interaction
and inuence on the dynamics of molecular structures needs to be taken into
account every time a biological system, i.e., a protein, derives function from
interaction with individual water molecules. This is indeed the case very
often.
The water molecule seems, at rst glance, simple and easy to model.
It has two hydrogen atoms and one oxygen atom with overall charges of
19
+1 e for each hydrogen and -2 e for the oxygen atom. However, a great
deal of information needs to be accounted for in a model. Where are these
charges located? Are these charges distributed in a spherical, uniform, and
symmetrical arrangement? How far are the hydrogens from the oxygen?
What is the angle formed by the hydrogens and the oxygen? Is this a exible
or rigid molecule? How does it interact with the surrounding molecules?
All these questions and many others need to be addressed for a faithful
description of water and its role in living cells.
5.1 Explicit Water Models
With some simplifying assumptions, a basic theoretical water model can be
constructed [19]. Imagine that the charges can be considered as point charges
located at the center of each atom. Furthermore, assume the distance be-
tween atoms H and O and the HOH angle are xed. Actually, if we use the
location of the nuclei to place the charges, the dipole moment will be exces-
sively high. By shifting the position of the negative charge from the oxygen
atom towards the hydrogen atoms along the HOH angle bisector, we can de-
crease the dipole moment and improve the model (see Fig. 12a). In addition,
we need to take into account the fact that molecules actually do not overlap.
We can incorporate this property into our model by dening a Lennard-Jones
like spherical repulsive potential centered in the molecule. If two molecules
get too close to each other, forces originated from this repulsive potential
will push them apart. The rigid theoretical model we just have constructed
was actually proposed in 1933 by Bernal and Fowler, and can be considered
the forerunner of subsequent three-point-charge models such as single point
charge (SPC) and transferable intermolecular potentials (TIPS), both widely
used in simulations of biomolecular systems [20]. In 1949, Rowlinson pro-
posed another theoretical model of water using a very similar approach to
that used by Bernal and Fowler, but the negative charge was split above and
below the molecular plane at the oxygen center (see Fig. 12b) in order to re-
produce the quadrupole moment of the water molecule [19]. Accordingly, the
Rowlinson model can be considered as the forerunner of subsequent multi-
point charge models used in computer simulations such as TIP4P, TIP5P,
and others [19].
Once the basic characteristics of the theoretical model are set (number of
charges, geometrical arrangement, type of interaction potentials, etc.), the
parameterization process that follows can become extremely complicated.
20
Figure 12: Water molecule models. a) Three site model of a water molecule in which
the hydrogen atoms (white) and the oxygen atom (red) are linked through rigid or exible
bonds. The positive charges are located at the nuclei of the hydrogen atoms. The negative
charged is located along the HOH angle bisector (blue sphere). Current models used in
computer simulations, SPC and TIP3P, are based on this kind of arrangement. b) Four
site model of a water molecule. Hydrogen atoms (white) and the oxygen atom (red) are
linked through rigid or exible bonds. In this case, the negative charge is split in two (blue
spheres) and located above and below the oxygen nuclei. Use the available VMD state
section05/models.vmd to further explore the water molecules presented in this gure.
Small changes in the parameters can have relevant implications for the macro-
scopic properties of water, such as density, diusion coecient, and specic
heat. Moreover, current theoretical models used in simulations of water in-
clude separate spherical dispersion and repulsion terms for the hydrogens
and oxygen, as well as exible harmonic bonds between the water molecule
atoms. These are called intra-molecular interactions. In some cases, even the
non-additive nature of water molecule interactions is included as molecular
polarizability. All these new properties present a challenge for both the force
eld parameterization and the simulation itself, which needs to deal with
more complicated and costly computations. Many of these models are still
under development and certainly fail to reproduce all properties of water as
determined by experiments [20].
Two of the most popular models of water used in simulations of biomolec-
ular systems are SPC and TIP3P. Both use three point charges but their
respective parameters are slightly dierent. For instance, the equilibrium
HOH angle is 109.47

in the SPC model, and 104.52

in the TIP3P model.

The SPC parameters have been chosen to reproduce energy and pressure of
liquid water at ambient conditions, as well as the second peak in the radial
21
distribution function of oxygen atoms (see Section 6.2). The TIP3P parame-
ters have been chosen to yield reasonable structural and energetic results for
both gas-phase dimers and pure liquids. The values of the dipole moment,
dielectric constant, and self diusion coecients for some popular models
used in computer simulations are provided in Table 5.1.
Model Dipole Dielectric Self Diusion
moment, D constant 10
5
cm
2
/s
SPC 2.27 65.0 3.85
TIP3P 2.35 82.0 5.19
TIP4P 2.18 53.0 3.29
Exp. 2.95 78.4 2.30
Table 2: Computed values for the dipole moment, the dielectric constant,
and self diusion coecient for SPC, TIP3P, and TIP4P models of water (see
http://www.lsbu.ac.uk/water/models.html). Experimental values are given as a reference.
Exercise 3: Modeling Water. Using VMD and the saved state sec-
tion05/models.vmd answer the following questions:
(1) What is the distance between the oxygen and hydrogen atoms?
(2) What is the angle formed by hydrogens and oxygen?
(3) By looking into the le section05/par all27 prot lipid.inp get the parameters uti-
lized for the TIP3P water model. Be sure to include bond, angle, and Lennard-Jones
terms. Hint: look for the words TIP3P or TIPS3P.
5.2 Implicit Water Models
Usually, if not always, biomolecular systems need to be described in wet
environments, where the role of the solvent may be essential. However, sim-
ulation of the solvent and the biomolecular system at the microscopic level
can be prohibitive. We already described the parameterization diculties of
a molecular model of water. Furthermore, the cost of the computation of the
interactions of all atoms is largely increased by the presence of an explicit
solvent. Large reductions in computing time can be obtained if description of
bulk solvent eects can be incorporated into the empirical potential function
describing the biomolecular system through an implicit solvent model. This
22
kind of implicit model needs to describe two important eects of a polar sol-
vent: the solvation of charges embedded in a polar medium, and the damping
of the electrostatic interaction between charges of the solute due to polariza-
tion by the medium. Most of the implicit models are based on the solution of
the Poisson-Boltzmann equation or on the generalized Born approach. Other
models include integral equations, stochastic methods, or screened Coulomb
potentials (SCP) that allow calculation of electrostatic properties in proteins.
Solvation energies and conduction properties of wide biological channels can
be reasonably estimated using implicit solvent models over time scales that
current simulations using explicit water models simply cannot reach. Al-
though these implicit models permit the exploration of biomolecular systems
over longer time scales, essential details may not be captured. For instance,
hydrodinamical eects and the interaction of water with proteins forming
and breaking water-protein hydrogen bonds during unfolding and folding is
essential, and, thus, a description of the events at the atomic level is required.
6 Simulating Water
As one might be able to infer from Section 5, reproducing the macroscopic
properties of water with either explicit or implicit models is not a simple
task. In this section, we will test the TIP3P model of water with molecular
dynamics simulations and analyze the results to give a feel for how these
models are used in practice.
6.1 Simulation: Melting Ice
One of the most peculiar and challenging points of computational water
modeling is the inability to reproduce freezing in simulated water molecules
when the temperature is decreased below the experimental freezing point.
Freezing of water in a simulation has been reported only one time in the
literature [21], however, multiple simulations on the order of microseconds
were necessary to see a single freezing event. The work used the TIP4P
model for water-water interaction.
In this section, we will not attempt to freeze liquid water, rather, we will
try the opposite. Beginning with a lattice of ice, we will simulate the melting
process.
We begin the simulation with the periodic ice lattice described in Section
23
2.2 and minimize the lattice for 1000 steps. The minimization simply allows
the lattice to adjust slightly to accomodate inaccuracies in our force eld,
which was actually not parameterized for ice. After minimization, we begin
a process of heating, starting at an initial temperature of 100 K and incre-
menting by 1 K every 2 ps. This heating protocol is somewhat unrealistic,
since it corresponds to a temerature change of 500 billion K in one second.
Ideally, the process should be much slower to allow equilibration after every
temperature jump, but one is often constrained to short time simulations in
molecular dynamics due to lack of computer power. Moreover, the heating
protocol will still reveal interesting melting behavior.
We also run the simulation with periodic boundary conditions. Thus, we
are eectively simulating an innite ice lattice, so that water molecules on
the edge will not cause premature melting. The initial boundary conditions
must be chosen to take the exact dimensions of the ice lattice into account.
The results of the simulation reveal an interesting melting mechanism (see
les section06/ICES.psf and section06/melting-sim/melt-100-grad2-01.dcd,
or the movie melting.mpg). The ice lattice is preserved for the rst 100 ps
of simulation when the temperature reaches 150 K. At this point, a slight
disruption in the lattice occurs and expands along a plane perpendicular to
lattice vector b until it stabilizes around 150 ps, when the temperature has
reached 200 K. The simulation continues until the sheared melting region
begins to expand when the simulation at 290 ps reaches 245 K. Melting
occurs beginning with the shear plane and expanding rapidly over the next
20 ps (during which the temperature is raised by another 10 K) until the
entire system is disordered and melted.
Two interesting issues arise from the results of the simulation. First, we
must ask if the shearing which occurs is a natural phenomenon or if it is
related to an inaccuracy of our forceeld. There is no denitive answer with
the data available, however, one could run a second simulation whereby the
water molecules have a dierent value for their angular spring constant. This
might change the manner in which the shearing occurs or may prevent that
it occurs altogether.
The second issue which should be noticed is that the melting temper-
ature is lower than the experimental value. This reects the diculty of
reproducing macroscopic properties in the modeling of water. Water force
eld parameters are generally looking to reproduce quantities for liquid wa-
ter at ambient room temperature since many of the chemical and biological
processes to be studied occur at that condition. One might also argue, in
24
Figure 13: Ice melting simulation. Water molecule congurations are shown at three
times and temperatures during the melting simulation: a) 20 ps, 110 K; b) 200 ps, 200
K; c) 400 ps, 300 K. Note the near perfect ice lattice at the beginning of the simulation
which begins to melt due to lattice shearing. View the melting simulation with les
section06/ICES.psf and section06/melting-sim/melt-100-grad2-01.dcd using VMD.
fact, that obtaining a melting temperature which is 20-30 K dierent than
experiment is quite good!
Exercise 4: Melting ice. The melting simulation is available in the direc-
tory section06/melting-sim/ as NAMD input and output les. Load the psf sec-
tion06/ICES.psf and into it the dcd section06/melting-sim/melt-100-grad2-01.dcd in
VMD and view the melting process.
Changing the angle spring constant of water may have an eect on the melting process
shown. What changes in the process do you expect for a larger spring constant or for
a smaller one?
6.2 Analysis: Pair Distribution Function*
In the previous simulation, we were able to observe the transition of ice to
water. We knew denitively that the solid phase existed at the beginning
of the simulation and that the liquid phase existed at the end. But what
happened during the melting process when it may have been unclear if a true
liquid phase existed? Moreover, what about species dierent from water such
as polymer gels, whose very nature make it dicult to distinguish a solid or
liquid phase? Is there some way of quantifying the ordering to tell us which
phase exists?
One way to observe the structural arrangement of molecules is through
computing the pair distribution function. The pair distribution function
25
g(r
i
, r
j
), or g(r) with r = |r
i
r
j
| for a homogeneous and isotropic material,
is dened as the probability of nding a pair of atoms at a distance r,
relative to the probability expected for a completely random distribution at
the same density [22]. In the canonical ensemble the actual formula is
g(r
1
, r
2
) =
N(N 1)

2
Z
NV T
_
dr
3
dr
4
. . . exp(V(r
1
, r
2
, . . . r
N
)) (3)
Here N is the total number of atoms, is the density of the system, Z
NV T
is the canonical partition function, = 1/kT, and V is the potential en-
ergy governing the motion of the molecules. The pair distribution function
can be obtained through a Fourier transformation of the structure factor ob-
tained by x-ray and neutron diraction patterns. Then, the experimentally
determined values of the pair distribution function can be compared to the
values computed using molecular models. An equivalent denition when us-
ing systems of identical atoms and useful for computer simulations is given
by
g(r) =
V
N
2
_

j=i
(r r
ij
)
_
(4)
where . . . denotes an ensemble average.
In order to gain a more intuitive knowledge of the pair distribution func-
tion, we analyze g(r) computed for the oxygen atoms of liquid water as
presented in Fig. 14. The value of g(r) is 0 for r 2.5

A, meaning that
the probability of nding two oxygen atoms at a distance r 2.5

A is zero.
0 2 4 6 8 10
r ()
0
1
2
3
4
g
(
r
)
Figure 14: Pair distribution
function for oxygens of liquid wa-
ter (TIP3P model) computed us-
ing NAMD and VMD.
Does this make sense? Yes! Water molecules
do not overlap, and hydrogen atoms share
some room with oxygen atoms. Therefore, we
should not nd oxygen atoms too close to each
other. The value of the pair distribution func-
tion suddenly increases for r increasing beyond
2.5

A and reaches a maximum of 3.8 at
r 3

A. This maximum value of g(r) means
that if we pick a random oxygen atom in the
system, and count the number of oxygen atoms
surrounding it at a distance r 3

A, then we
will nd a number that is approximately three
times larger than what we would have found
in a completely random distribution of oxygen
26
atoms at the same density. Is this reasonable?
Yes! Hydrogen bonds between adjacent wa-
ter molecules favor structured conformations
of the molecules at short distances, even in
liquid water. The value of the pair distribution function then decreases and
uctuates around one, reecting the fact that at longer scales the structure
of water molecules is lost, and the distribution of oxygen atoms is fairly ran-
dom, as expected for liquid water. The pair distribution for oxygen atoms
of water using the TIP3P water model described above compares reasonably
well with experimental values. The theoretical model simulated seems to
reproduce at least one of many properties of real water!
Exercise 5: Order in Water. Many structural properties can be tested using the
pair distribution function. Using the provided scripts (section06/calcpdf.tcl and sec-
tion06/pdf.tcl) along with VMD and the provided trajectories try to answer the fol-
lowing questions:
(1) How does the pair distribution function look like for liquid water when computed
using hydrogen atoms only? Hint: Modify the script calcpdf.tcl in order to select
hydrogens (name H1 H2) instead of oxygens (name OH2). The script is already set
to use a trajectory of liquid water (melt-300-01.dcd). In order to use it, open the
Tk Console in VMD and type source calc.pdf (make sure that you are in the right
directory by using the commands pwd and cd) and then plot the resulting list.
How does the plot compares with Fig. 14?
(2) How does the pair distribution look like when computed on a crystal of water?
Hint: Modify the script calcpdf.tcl so as to use the trajectory melt-10-01.dcd. Note
that the computation of the pair distribution function is very demanding and it may
take over 10 minutes to go over all the frames of the provided trajectories. While you
wait, check the phases of water provided in the trajectories by opening another session
of VMD, loading the le ICES.psf, and into it the trajectories melt-300-01.dcd and
melt-10-01.dcd.
6.3 Theory: Temperature in Molecular Dynamics*
Determining the temperature of a substance is a simple task in practice,
as long as one has a thermometer present. Complications may arise, how-
ever, when the amount of substance you want to measure is very small. For
instance, determining the temperature of a microscopic object would be a
considerable challenge. An ever greater challenge arises at the atomic level,
27
when the very concept of temperature becomes skewed. Temperature is a
macroscopic property. It is a single number that measures the average ki-
netic energy of a large system as a whole. How then are we to determine an
accurate temperature for an atomic system, one in which the kinetic energy
may vary wildly from atom to atom? Better yet, does temperature even
have a meaning for very small systems like a single molecule of water? These
questions need to be considered as one attempts an atomic representation of
physical phenomena.
Molecular dynamics make use of the equipartition theorem of statistical
mechanics [23] to determine the temperature of a simulation at any time. For
a system of N atoms, the equipartition theorem equates the thermal energy
of the system to the average kinetic energy of its atoms and accounts for any
constraints on the atoms:
(3N N
c
)
2
k
B
T = E
kin
(5)
from which follows
T =
2
(3N N
c
)k
B
E
kin
. (6)
Here, N
c
is the total number of internal constraints, such as xed bonds,
which may be the case for hydrogen atoms in a simulation, or xed bond
angles, which may be the case for some water models, as discussed earlier.
Note that 3N N
c
is the total number of degrees of freedom of the system.
Furthermore, note N
c
also includes constraints that may be placed on the
system to conserve momentum. In a system in which the momentum is xed
but no other constraints are placed on the atoms, N
c
= 3, since the entire
center of mass velocity is constrained to oset the atomic momenta.
With relationship (6), we can use the kinetic energy of the system as a
thermometer. The kinetic energy of a system of atoms at any time t is
dened as
E
kin
(t) =
1
2
N

i=1
m
i
|v
i
(t)|
2
. (7)
In this equation the sum is over all atoms of the simulated sample. One can
then express the instantaneous temperature of the system as
T(t) =
1
(3N N
c
)k
B
N

i=1
m
i
|v
i
(t)|
2
. (8)
28
Figure 15: The left panel shows the instantaneous temperature computed for a periodic
box of 11763 water molecules over a 50 ps molecular dynamics simulation. The right panel
shows the corresponding distribution of temperatures sampled during the simulation (blue
histogram) and the theoretically derived distribution (black).
Furthermore, one can use the average kinetic energy of the system
E
kin
=
_
1
2
N

i=1
m
i
|v
i
|
2
_
(9)
to dene the average temperature of the system as
T =
2
(3N N
c
)k
B
_
1
2
N

i=1
m
i
|v
i
|
2
_
. (10)
In the above equations (9) and (10), the brackets denote a temporal average.
Figure 15 shows the time dependence of T(t) for a simulated water system.
One can recognize that T(t) uctuates around its average value T. The
right panel of Fig. 15 presents the sampled distribution of T(t) values.
One can recognize from Fig. 15 that the temperature exhibits a Gaussian
distribution around T with a width
T
p(T) =
1
_
2
2
T
exp
_

(T T)
2
2
2
T
_
. (11)
29
Here, p(T) is the probability of determining a system temperature T from the
sample employing Eq. 8. The reason why the temperature readings produce
a rather wide distribution is actually the nite size of the simulated system,
as we will demonstrate below. In fact, the larger the system, the narrower
the distribution of temperatures. This decrease in uctuation with size may
be understood if one thinks about a coin ip. On average, the probability
of obtaining heads or tails is 50%. This is not the case, however, after you
examine a single ip. The uctuation from the average value is 50%! Upon
increasing the number of ips, however, the deviation decreases as 1/
_
N
flips
.
In order to derive the distribution (11) from rst principles, we begin with
the fact that the velocity of each atom is distributed according to the well
known Maxwell distribution. For each Cartesian component of the velocity,
e.g., for the x-component of the velocity of the jth atom, v
jx
, the Maxwell
distribution is
p(v
jx
) =
_
m
j
2k
B
T
exp
_

m
j
2k
B
T
v
2
jx
_
(12)
and the average kinetic energy
jx
=
1
2
m
j
v
2
jx
is accordingly

jx
=
_
m
j
2k
B
T
_
+

dv
jx
1
2
m
j
v
2
jx
exp
_

m
j
2k
B
T
v
2
jx
_
. (13)
Introducing the integration variable y =
_
m
j
/2k
B
T v
jx
leads to

jx
=
k
B
T

_
+

dy y
2
exp
_
y
2

. (14)
The integral arising here is equal to
1
2

. Hence, we nd

jx
=
1
2
k
B
T (15)
which is the expected result. The average of the energy squared is

2
jx
=
_
m
j
2k
B
T
_
+

dv
jx
_
1
2
m
j
v
2
jx
_
2
exp
_

m
j
2k
B
T
v
2
jx
_
. (16)
Introducing the same integration variable y, we nd

2
jx
=
(k
B
T)
2

_
+

dy y
4
exp
_
y
2

. (17)
30
The integral arising here is equal to
3
4

. Hence, we nd

2
jx
=
3
4
( k
B
T )
2
. (18)
The results described for
jx
and
2
jx
permit us to determine the mean
square deviation of the kinetic energy uctuations. First, one derives the
well known result
(
jx
(t)
jx
)
2
= (
2
jx
(t) 2
jx
(t)
jx
+
jx

2
) =
2
jx

jx

2
.
(19)
Then, using Eqs. (15) and (18), we can conclude

2
jx

jx

2
=
1
2
(k
B
T)
2
. (20)
From this follows the mean square uctuation of the kinetic energy of the
entire simulated system with (3N N
c
) degrees of freedom

2
E
= E
2
kin
E
kin

2
=
(3N N
c
)
2
(k
B
T)
2
. (21)
One can dene similarly the temperature uctuation
T
, and using the de-
nition of T as given in Eqs. (9) and (10), one concludes

2
T
= T
2
T
2
=
4
(3N N
c
)
2
k
2
B

2
E
(22)
or

2
T
=
2T
2
3N N
c
. (23)
From our derivation we conclude that over the course of a simulation, the
temperature will change even when the system is in equilibrium. Deviations
in atomic interactions will cause changes in atomic velocities, so an equili-
brated system whose temperature is 300 K at one time may later have a
temperature of 297 K at another time and a temperature of 305 K still later.
The above result states that the temperature distribution shown in Fig. 15
should be characterized through a width that depends on temperature and
particle number. As one can see in Fig. 15, the Gaussian distribution with
a width predicted by Eq. (23) matches the simulated distribution. One con-
cludes therefore that the observed temperature uctuation in the simulation
is a nite size eect, with relative uctuations on the order of 1/

3N N
c
.
31
References
[1] D. Voet, J. G. Voet, and C. W. Pratt. Fundamentals of Biochemistry.
John Wiley & Sons, Inc., New York, 1999.
[2] D. R. Lide, editor. CRC Handbook of Chemistry and Physics. CRC
Press, Boca Raton, 74th edition, 1994.
[3] P. L. Davies, J. Baardsnes, M. J. Kuiper, and V. K. Walker. Struc-
ture and function of antifreeze proteins. Phil. Trans. R. Soc. Lond. B,
357:927935, 2002.
[4] C. A. Knight, C. C. Cheng, and A. L. DeVries. Adsorption of alpha-
helical antifreeze peptides on specic ice crystal surface planes. Biophys-
ical Journal, 59:409418, 1991.
[5] J. D. Madura, K. Baran, and A. Wierzbicki. Molecular recognition
and binding of thermal hysteresis proteins to ice. Journal of Molecular
Recognition, 13:101113, 2000.
[6] S. M. McDonald, A. White, P. Clancy, and J. W. Brady. Binding of an
antifreeze polypeptide to an ice/water interface via computer simulation.
AIChe, 41(4):959973, 1995.
[7] Z. Jia, C. I. DeLuca, H. Chao, and P. L. Davies. Structural basis for
the binding of a globular antifreeze protein to ice. Nature, 384:285288,
1996.
[8] A. L. DeVries. Antifreeze proteins and glycopeptides in cold-water shes.
A. Rev. Physiol., 45:245260, 1983.
[9] A. W. Adamson. Physical Chemistry of Surfaces. John Wiley & Sons,
3rd edition, 1976.
[10] P. Atkins and J. de Paula. Atkins Physical Chemistry. Oxford Univer-
sity Press, 2002.
[11] F. Heslot, N. Fraysse, and M. Cazabat. Molecular layering in the spread-
ing of wetting liquid drops. Nature, 338:640642, 1989.
[12] P. Roura and J. Fort. Local thermodynamic derivaiton of youngs equa-
tion. Journal of Colloid and Interface Science, 272:420429, 2004.
32
[13] R. Lamb and N. Furlong. Controlled wettability of quartz surfaces.
Journal of the Chemical Society-Faraday Transactions, 78:6173, 1982.
[14] A. Brodka and T. W. Zerda. Properties of liquid acetone in silica pores:
Molecular dynamics simulation. Journal of Chemical Physics, 104:6319
6326, 1996.
[15] J. Hill and J. Sauer. Molecular mechanics potential for silica and zeolite
catalysts based on ab initio calculations. 1. dense and microporous silica.
Nature, 338:640642, 1989.
[16] X. Gao and L. Jiang. Water-repellent legs of water striders. Nature,
432:36, 2004.
[17] D. L. Hu, B. Chan, and J. W. M. Bush. The hydrodynamics of water
strider locomotion. Nature, 424:663666, 2003.
[18] L. Feng, S. Li, Y. Li, H. Li, L. Zhang, J. Zhai, Y. Song, B. Liu, L. Jiang,
and D. Zhu. Super-hydrophobic surfaces: From natural to articial.
Advanced Materials, 14:18571860, 2002.
[19] J. L. Finney. The water molecule and its interactions: the interaction
between theory, modelling, and experiment. Journal of Molecular Liq-
uids, 90:303312, 2001.
[20] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and
M. L. Klein. Comparison of simple potential functions for simulating
liquid water. JCP, 79:926935, 1983.
[21] M. Matsumoto, S. Saito, and I. Ohmine. Molecular dynamics simula-
tion of the ice nucleation and growth process leading to water freezing.
Nature, 416:409413, 2002.
[22] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford
University Press, New York, 1987.
[23] Kerson Huang. Statistical Mechanics. Wiley Eastern Limited, New
Delhi, 1988.
33