trunk/fennellDissertation/iceChapter.tex

\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}

Water has proven to be a challenging substance to depict in
simulations, and a variety of models have been developed to describe
its behavior under varying simulation
conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04}
These models have been used to investigate important physical
phenomena like phase transitions and the hydrophobic
effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models
available, it is only natural to compare the models under interesting
thermodynamic conditions in an attempt to clarify the limitations of
each of the models.\cite{Jorgensen83,Jorgensen98,Baez94,Mahoney01} Two
important property to quantify are the Gibbs and Helmholtz free
energies, particularly for the solid forms of water, as these predict
the thermodynamic stability of the various phases. Water has a
particularly rich phase diagram and takes on a number of different and
stable crystalline structures as the temperature and pressure are
varied. This complexity makes it a challenging task to investigate the
entire free energy landscape.\cite{Sanz04} Ideally, research is
focused on the phases having the lowest free energy at a given state
point, because these phases will dictate the relevant transition
temperatures and pressures for the model.

The high-pressure phases of water (ice II-ice X as well as ice XII)
have been studied extensively both experimentally and
computatuionally. In this chapter, standard reference state methods
were applied in the {\it low} pressure regime to evaluate the free
energies for a few known crystalline water polymorphs that might be
stable at these pressures. This work is unique in the fact that one of
the crystal lattices was arrived at through crystallization of a
computationally efficient water model under constant pressure and
temperature conditions. Crystallization events are interesting in and
of themselves;\cite{Matsumoto02,Yamada02} however, the crystal
structure obtained in this case is different from any previously
observed ice polymorphs in experiment or simulation.\cite{Fennell04}
We have named this structure Ice-{\it i} to indicate its origin in
computational simulation. The unit cell of Ice-$i$ and an axially
elongated variant named Ice-$i^\prime$ both consist of eight water
molecules that stack in rows of interlocking water tetramers as
illustrated in figure \ref{fig:unitCell}A,B. These tetramers form a
crystal structure similar in appearance to a recent two-dimensional
surface tessellation simulated on silica.\cite{Yang04} As expected in
an ice crystal constructed of water tetramers, the hydrogen bonds are
not as linear as those observed in ice I$_\textrm{h}$; however, the
interlocking of these subunits appears to provide significant
stabilization to the overall crystal. The arrangement of these
tetramers results in open octagonal cavities that are typically
greater than 6.3\AA\ in diameter (see figure
\ref{fig:protOrder}). This open structure leads to crystals that are
typically 0.07 g/cm$^3$ less dense than ice I$_\textrm{h}$.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/unitCell.pdf}
\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
elongated variant of Ice-{\it i}.  For Ice-{\it i}, the $a$ to $c$
relation is given by $a = 2.1214c$, while for Ice-$i^\prime$, $a =
1.7850c$.} 
\label{fig:iceiCell}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/orderedIcei.pdf}
\caption{Image of a proton ordered crystal of Ice-{\it i} looking 
down the (001) crystal face. The rows of water tetramers surrounded by
octagonal pores leads to a crystal structure that is significantly
less dense than ice I$_\textrm{h}$.}
\label{fig:protOrder}
\end{figure}

Results from our previous study indicated that Ice-{\it i} is the
minimum energy crystal structure for the single point water models
investigated (for discussions on these single point dipole models, see
the previous work and related
articles\cite{Fennell04,Liu96,Bratko85}). Our earlier results only
considered energetic stabilization and neglected entropic
contributions to the overall free energy. To address this issue, we
have calculated the absolute free energy of this crystal using
thermodynamic integration and compared to the free energies of ice
I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common low-density ice
polymorphs) and ice B (a higher density, but very stable crystal
structure observed by B\`{a}ez and Clancy in free energy studies of
SPC/E).\cite{Baez95b} This work includes results for the water model
from which Ice-{\it i} was crystallized (SSD/E) in addition to several
common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
field parametrized single point dipole water model (SSD/RF). The
axially elongated variant, Ice-$i^\prime$, was used in calculations
involving SPC/E, TIP4P, and TIP5P. The square tetramer in Ice-$i$
distorts in Ice-$i^\prime$ to form a rhombus with alternating 85 and
95 degree angles. Under SPC/E, TIP4P, and TIP5P, this geometry is
better at forming favorable hydrogen bonds. The degree of rhomboid
distortion depends on the water model used but is significant enough
to split the peak in the radial distribution function which corresponds
to diagonal sites in the tetramers.

\section{Methods}

Canonical ensemble ({\it NVT}) molecular dynamics calculations were
performed using the OOPSE molecular mechanics package.\cite{Meineke05}
The densities chosen for the simulations were taken from
isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and
200K. Each model (and each crystal structure) was allowed to relax for
300ps in the {\it NPT} ensemble before averaging the density to obtain
the volumes for the {\it NVT} simulations.All molecules were treated
as rigid bodies, with orientational motion propagated using the
symplectic DLM integration method described in section
\ref{sec:IntroIntegration}.

We used thermodynamic integration to calculate the Helmholtz free
energies ({\it A}) of the listed water models at various state
points. Thermodynamic integration is an established technique that has
been used extensively in the calculation of free energies for
condensed phases of
materials.\cite{Frenkel84,Hermans88,Meijer90,Baez95,Vlot99} This
method uses a sequence of simulations during which the system of
interest is converted into a reference system for which the free
energy is known analytically ($A_0$). This transformation path is then
integrated in order to determine the free energy difference between
the two states:
\begin{equation}
\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda 
)}{\partial\lambda}\right\rangle_\lambda d\lambda,
\end{equation}
where $V$ is the interaction potential and $\lambda$ is the
transformation parameter that scales the overall
potential. Simulations are distributed unevenly along this path in
order to sufficiently sample the regions of greatest change in the
potential. Typical integrations in this study consisted of $\sim$25
simulations ranging from 300ps (for the unaltered system) to 75ps
(near the reference state) in length.

For the thermodynamic integration of molecular crystals, the Einstein
crystal was chosen as the reference state. In an Einstein crystal, the
molecules are harmonically restrained at their ideal lattice locations
and orientations. The partition function for a molecular crystal
restrained in this fashion can be evaluated analytically, and the
Helmholtz Free Energy ({\it A}) is given by
\begin{eqnarray}
A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
[\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
)^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
(-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
)^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
\label{eq:ecFreeEnergy}
\end{eqnarray}
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
\ref{eq:ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
$K_\mathrm{\omega}$ are the spring constants restraining translational
motion and deflection of and rotation around the principle axis of the
molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
minimum potential energy of the ideal crystal. In the case of
molecular liquids, the ideal vapor is chosen as the target reference
state.

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/rotSpring.pdf}
\caption{Possible orientational motions for a restrained molecule. 
$\theta$ angles correspond to displacement from the body-frame {\it
z}-axis, while $\omega$ angles correspond to rotation about the
body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
constants for the harmonic springs restraining motion in the $\theta$
and $\omega$ directions.}
\label{waterSpring}
\end{figure}

Charge, dipole, and Lennard-Jones interactions were modified by a
cubic switching between 100\% and 85\% of the cutoff value (9\AA). By
applying this function, these interactions are smoothly truncated,
thereby avoiding the poor energy conservation which results from
harsher truncation schemes. The effect of a long-range correction was
also investigated on select model systems in a variety of manners. For
the SSD/RF model, a reaction field with a fixed dielectric constant of
80 was applied in all simulations.\cite{Onsager36} For a series of the
least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
simulations were performed with longer cutoffs of 12 and 15\AA\ to
compare with the 9\AA\ cutoff results. Finally, results from the use
of an Ewald summation were estimated for TIP3P and SPC/E by performing
calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
mechanics software package.\cite{Tinker} The calculated energy
difference in the presence and absence of PME was applied to the
previous results in order to predict changes to the free energy
landscape.

\section{Results and discussion}

The free energy of proton ordered Ice-{\it i} was calculated and
compared with the free energies of proton ordered variants of the
experimentally observed low-density ice polymorphs, I$_\textrm{h}$ and
I$_\textrm{c}$, as well as the higher density ice B, observed by
B\`{a}ez and Clancy and thought to be the minimum free energy
structure for the SPC/E model at ambient conditions.\cite{Baez95b} Ice
XI, the experimentally-observed proton-ordered variant of ice
I$_\textrm{h}$, was investigated initially, but was found to be not as
stable as proton disordered or antiferroelectric variants of ice
I$_\textrm{h}$. The proton ordered variant of ice I$_\textrm{h}$ used
here is a simple antiferroelectric version that has an 8 molecule unit
cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules
for ice B, 1024 or 1280 molecules for ice I$_\textrm{h}$, 1000
molecules for ice I$_\textrm{c}$, or 1024 molecules for Ice-{\it
i}. The larger crystal sizes were necessary for simulations involving
larger cutoff values.

\begin{table}
\centering
\caption{HELMHOLTZ FREE ENERGIES FOR SEVERAL ICE POLYMORPHS WITH A 
VARIETY OF COMMON WATER MODELS AT 200 KELVIN AND 1 ATMOSPHERE}
\begin{tabular}{ l  c  c  c  c }
\toprule
\toprule
Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or $i^\prime$) \\
& kcal/mol & kcal/mol & kcal/mol & kcal/mol \\
\midrule
TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\
TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\
TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\
SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\
SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\
SSD/RF & -11.51(4) & NA & -12.08(5) & -12.29(4)\\
\bottomrule
\end{tabular}
\label{tab:freeEnergy}
\end{table}

Table \ref{tab:freeEnergy} shows the results of the free energy
calculations with a cutoff radius of 9\AA. It should be noted that the
ice I$_\textrm{c}$ crystal polymorph is not stable at 200K and 1 atm
with the SSD/RF water model, hense omitted results for that cell. The
free energy values displayed in this table, it is clear that Ice-{\it
i} (or Ice-$i^\prime$ for TIP4P, TIP5P, and SPC/E) is the most stable
state for all of the common water models studied.

With the free energy at these state points, the Gibbs-Helmholtz
equation was used to project to other state points and to build phase
diagrams.  Figures \ref{fig:tp3phasedia} and \ref{fig:ssdrfphasedia}
are example diagrams built from the free energy results. All other
models have similar structure, although the crossing points between
the phases exist at slightly different temperatures and pressures. It
is interesting to note that ice I does not exist in either cubic or
hexagonal form in any of the phase diagrams for any of the models. For
purposes of this study, ice B is representative of the dense ice
polymorphs. A recent study by Sanz {\it et al.} goes into detail on
the phase diagrams for SPC/E and TIP4P in the high pressure
regime.\cite{Sanz04}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/tp3PhaseDia.pdf}
\caption{Phase diagram for the TIP3P water model in the low pressure 
regime. The displayed $T_m$ and $T_b$ values are good predictions of
the experimental values; however, the solid phases shown are not the
experimentally observed forms. Both cubic and hexagonal ice $I$ are
higher in energy and don't appear in the phase diagram.}
\label{fig:tp3phasedia}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdrfPhaseDia.pdf}
\caption{Phase diagram for the SSD/RF water model in the low pressure 
regime. Calculations producing these results were done under an
applied reaction field. It is interesting to note that this
computationally efficient model (over 3 times more efficient than
TIP3P) exhibits phase behavior similar to the less computationally
conservative charge based models.}
\label{fig:ssdrfphasedia}
\end{figure}

\begin{table}
\centering
\caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) 
temperatures at 1 atm for several common water models compared with
experiment. The $T_m$ and $T_s$ values from simulation correspond to a
transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
liquid or gas state.}
\begin{tabular}{ l  c  c  c  c  c  c  c }
\toprule
\toprule
Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
\midrule
$T_m$ (K)  & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\
$T_b$ (K)  & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\
$T_s$ (K)  & - & - & - & - & 355(3) & - & -\\
\bottomrule
\end{tabular}
\label{tab:meltandboil}
\end{table}

Table \ref{tab:meltandboil} lists the melting and boiling temperatures
calculated from this work. Surprisingly, most of these models have
melting points that compare quite favorably with experiment. The
unfortunate aspect of this result is that this phase change occurs
between Ice-{\it i} and the liquid state rather than ice
I$_\textrm{h}$ and the liquid state. These results are actually not
contrary to previous studies in the literature. Earlier free energy
studies of ice I using TIP4P predict a $T_m$ ranging from 214 to 238K
(differences being attributed to choice of interaction truncation and
different ordered and disordered molecular
arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
Ice-{\it i} were omitted, a $T_m$ value around 210K would be predicted
from this work. However, the $T_m$ from Ice-{\it i} is calculated at
265K, significantly higher in temperature than the previous
studies. Also of interest in these results is that SSD/E does not
exhibit a melting point at 1 atm, but it shows a sublimation point at
355K. This is due to the significant stability of Ice-{\it i} over all
other polymorphs for this particular model under these
conditions. While troubling, this behavior turned out to be
advantageous in that it facilitated the spontaneous crystallization of
Ice-{\it i}. These observations provide a warning that simulations of
SSD/E as a ``liquid'' near 300K are actually metastable and run the
risk of spontaneous crystallization. However, this risk changes when
applying a longer cutoff.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
\caption{Free energy as a function of cutoff radius for (A) SSD/E, (B) 
TIP3P, and (C) SSD/RF. Data points omitted include SSD/E:
I$_\textrm{c}$ 12\AA\, TIP3P: I$_\textrm{c}$ 12\AA\ and B 12\AA\,
and SSD/RF: I$_\textrm{c}$ 9\AA . These crystals are unstable at 200 K
and rapidly convert into liquids. The connecting lines are qualitative
visual aid.}
\label{fig:incCutoff}
\end{figure}

Increasing the cutoff radius in simulations of the more
computationally efficient water models was done in order to evaluate
the trend in free energy values when moving to systems that do not
involve potential truncation. As seen in figure \ref{fig:incCutoff},
the free energy of all the ice polymorphs show a substantial
dependence on cutoff radius. In general, there is a narrowing of the
free energy differences while moving to greater cutoff
radius. Interestingly, by increasing the cutoff radius, the free
energy gap was narrowed enough in the SSD/E model that the liquid
state is preferred under standard simulation conditions (298K and 1
atm). Thus, it is recommended that simulations using this model choose
interaction truncation radii greater than 9\AA\ . This narrowing
trend is much more subtle in the case of SSD/RF, indicating that the
free energies calculated with a reaction field present provide a more
accurate picture of the free energy landscape in the absence of
potential truncation.

To further study the changes resulting to the inclusion of a
long-range interaction correction, the effect of an Ewald summation
was estimated by applying the potential energy difference do to its
inclusion in systems in the presence and absence of the
correction. This was accomplished by calculation of the potential
energy of identical crystals with and without PME using TINKER. The
free energies for the investigated polymorphs using the TIP3P and
SPC/E water models are shown in Table \ref{tab:pmeShift}. TIP4P and
TIP5P are not fully supported in TINKER, so the results for these
models could not be estimated. The same trend pointed out through
increase of cutoff radius is observed in these PME results. Ice-{\it
i} is the preferred polymorph at ambient conditions for both the TIP3P
and SPC/E water models; however, there is a narrowing of the free
energy differences between the various solid forms. In the case of
SPC/E this narrowing is significant enough that it becomes less clear
that Ice-{\it i} is the most stable polymorph, and is possibly
metastable with respect to ice B and possibly ice
I$_\textrm{c}$. However, these results do not significantly alter the
finding that the Ice-{\it i} polymorph is a stable crystal structure
that should be considered when studying the phase behavior of water
models.

\begin{table}
\centering
\caption{The free energy of the studied ice polymorphs after applying 
the energy difference attributed to the inclusion of the PME
long-range interaction correction. Units are kcal/mol.}
\begin{tabular}{ l  c  c  c  c }
\toprule
\toprule
Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or Ice-$i^\prime$) \\
\midrule
TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\
SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\
\bottomrule
\end{tabular}
\label{tab:pmeShift}
\end{table}

\section{Conclusions}

The free energy for proton ordered variants of hexagonal and cubic ice
$I$, ice B, and recently discovered Ice-{\it i} were calculated under
standard conditions for several common water models via thermodynamic
integration. All the water models studied show Ice-{\it i} to be the
minimum free energy crystal structure in the with a 9\AA\ switching
function cutoff. Calculated melting and boiling points show
surprisingly good agreement with the experimental values; however, the
solid phase at 1 atm is Ice-{\it i}, not ice I$_\textrm{h}$. The
effect of interaction truncation was investigated through variation of
the cutoff radius, use of a reaction field parameterized model, and
estimation of the results in the presence of the Ewald
summation. Interaction truncation has a significant effect on the
computed free energy values, and may significantly alter the free
energy landscape for the more complex multipoint water models. Despite
these effects, these results show Ice-{\it i} to be an important ice
polymorph that should be considered in simulation studies.

Due to this relative stability of Ice-{\it i} in all manner of
investigated simulation examples, the question arises as to possible
experimental observation of this polymorph.  The rather extensive past
and current experimental investigation of water in the low pressure
regime makes us hesitant to ascribe any relevance of this work outside
of the simulation community.  It is for this reason that we chose a
name for this polymorph which involves an imaginary quantity.  That
said, there are certain experimental conditions that would provide the
most ideal situation for possible observation. These include the
negative pressure or stretched solid regime, small clusters in vacuum
deposition environments, and in clathrate structures involving small
non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
our predictions for both the pair distribution function ($g_{OO}(r)$)
and the structure factor ($S(\vec{q})$ for ice I$_\textrm{c}$ and for
ice-{\it i} at a temperature of 77K.  In a quick comparison of the
predicted S(q) for Ice-{\it i} and experimental studies of amorphous
solid water, it is possible that some of the ``spurious'' peaks that
could not be assigned in HDA could correspond to peaks labeled in this
S(q).\cite{Bizid87} It should be noted that there is typically poor
agreement on crystal densities between simulation and experiment, so
such peak comparisons should be made with caution.  We will leave it
to our experimental colleagues to determine whether this ice polymorph
is named appropriately or if it should be promoted to Ice-0.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
\caption{Radial distribution functions of Ice-{\it i} and ice 
I$_\textrm{c}$ calculated from from simulations of the SSD/RF water
model at 77 K.}
\label{fig:gofr}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{./figures/sofq.pdf}
\caption{Predicted structure factors for Ice-{\it i} and ice 
I$_\textrm{c}$ at 77 K.  The raw structure factors have been
convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
width) to compensate for the trunction effects in our finite size
simulations. The labeled peaks compared favorably with ``spurious''
peaks observed in experimental studies of amorphous solid
water.\cite{Bizid87}}
\label{fig:sofq}
\end{figure}

Revision:	2978
Committed:	Sun Aug 27 19:34:53 2006 UTC (18 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	23199 byte(s)
Log Message:	fixing up some chapters