trunk/fennellDissertation/iceChapter.tex

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

As discussed in the previous chapter, 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,Jorgensen98b,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
computationally. 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. 

While performing a series of melting simulations on an early iteration
of SSD/E, we observed several recrystallization events at a constant
pressure of 1 atm. After melting from ice I$_\textrm{h}$ at 235~K, two
of five systems recrystallized near 245~K. Crystallization events are
interesting in and of themselves;\cite{Matsumoto02,Yamada02} however,
the crystal structure extracted from these systems 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:iceiCell}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 initial studies 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}). These 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 and Thermodynamic Integration}

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
200~K. Each model (and each crystal structure) was allowed to relax for
300~ps 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:IntroIntegrate}.


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,Hermens88,Meijer90,Baez95a,Vlot99}.  This
method uses a sequence of simulations over which the system of
interest is converted into a reference system for which the free
energy is known analytically ($A_0$).  The difference in potential
energy between the reference system and the system of interest
($\Delta V$) is then integrated in order to determine the free energy
difference between the two states:
\begin{equation}
 A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
\end{equation}
Here, $\lambda$ is the parameter that governs the transformation
between the reference system and the system of interest.  For
crystalline phases, an harmonically-restrained (Einstein) crystal is
chosen as the reference state, while for liquid phases, the ideal gas
is taken as the reference state. Figure \ref{fig:integrationPath}
shows an example integration path for converting a crystalline system
to the Einstein crystal reference state.
\begin{figure}
\includegraphics[width=\linewidth]{./figures/integrationPath.pdf}
\caption{An example integration path to convert an unrestrained 
crystal ($\lambda = 1$) to the Einstein crystal reference state
($\lambda = 0$). Note the increase in samples at either end of the
path to improve the smoothness of the curve. For reversible processes,
conversion of the Einstein crystal back to the system of interest will
give an identical plot, thereby integrating to the same result.}
\label{fig:integrationPath}
\end{figure}

In an Einstein crystal, the molecules are restrained at their ideal
lattice locations and orientations. Using harmonic restraints, as
applied by B\'{a}ez and Clancy, the total potential for this reference
crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
\begin{equation}
V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
\frac{K_\omega\omega^2}{2},
\end{equation}
where $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.  These spring constants are typically calculated from
the mean-square displacements of water molecules in an unrestrained
ice crystal at 200~K.  For these studies, $K_\mathrm{v} =
4.29$~kcal~mol$^{-1}$~\AA$^{-2}$, $K_\theta\ =
13.88$~kcal~mol$^{-1}$~rad$^{-2}$, and $K_\omega\ =
17.75$~kcal~mol$^{-1}$~rad$^{-2}$.  It is clear from
Fig. \ref{fig:waterSpring} that the values of $\theta$ range from $0$
to $\pi$, while $\omega$ ranges from $-\pi$ to $\pi$.  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{equation}
\begin{split}
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] \\ 
        &- 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],
\end{split}
\label{eq:ecFreeEnergy}
\end{equation}
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
potential energy of the ideal crystal.\cite{Baez95a} The choice of an
Einstein crystal reference state is somewhat arbitrary. Any ideal
system for which the partition function is known exactly could be used
as a reference point as long as the system does not undergo a phase
transition during the integration path between the real and ideal
systems.  Nada and van der Eerden have shown that the use of different
force constants in the Einstein crystal does not affect the total
free energy, and Gao {\it et al.} have shown that free energies
computed with the Debye crystal reference state differ from the
Einstein crystal by only a few tenths of a
kJ~mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can
lead to some uncertainty in the computed melting point of the solids.
\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{fig:waterSpring}
\end{figure}

In the case of molecular liquids, the ideal vapor is chosen as the
target reference state.  There are several examples of liquid state
free energy calculations of water models present in the
literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
typically differ in regard to the path taken for switching off the
interaction potential to convert the system to an ideal gas of water
molecules.  In this study, we applied one of the most convenient
methods and integrated over the $\lambda^4$ path, where all
interaction parameters are scaled equally by this transformation
parameter.  This method has been shown to be reversible and provide
results in excellent agreement with other established
methods.\cite{Baez95b}

Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
Lennard-Jones interactions were gradually reduced by a cubic switching
function.  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, TIP3P, and SPC/E), simulations were
performed with longer cutoffs of 10.5, 12, 13.5, and 15~\AA\ to
compare with the 9~\AA\ cutoff results.  Finally, the effects of using
the Ewald summation were estimated for TIP3P and SPC/E by performing
single configuration Particle-Mesh Ewald (PME) calculations for each
of the ice polymorphs.\cite{Ponder87} 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{Initial Free Energy Results}

The calculated free energies of proton-ordered variants of three low
density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it i} or
Ice-$i^\prime$) and the stable higher density ice B are listed in
table \ref{tab:freeEnergy}.  Ice B was included because it has been
shown to be a minimum free energy structure for SPC/E at ambient
conditions.\cite{Baez95b} In addition to the free energies, the
relevant transition temperatures at standard pressure are also
displayed in table \ref{tab:freeEnergy}.  These free energy values
indicate that Ice-{\it i} is the most stable state for all of the
investigated water models.  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 results for the TIP3P and SSD/RF water models.  All
other models have similar structure, although the crossing points
between the phases move to different temperatures and pressures as
indicated from the transition temperatures in table
\ref{tab:freeEnergy}.  It is interesting to note that ice
I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear
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.} provides details on the phase
diagrams for SPC/E and TIP4P at higher pressures than those studied
here.\cite{Sanz04}
\begin{table}
\centering
\caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1 
ATMOSPHERE FOR SEVERAL WATER MODELS}

\footnotesize
\begin{tabular}{lccccccc}
\toprule
\toprule
Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\
\cmidrule(lr){2-6}
\cmidrule(l){7-8} 
& \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\
\midrule
TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\
TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\
TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\
SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\
SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\
SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\
\bottomrule
\end{tabular}
\label{tab:freeEnergy}
\end{table}

\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}

We note that all of the crystals investigated in this study are ideal
proton-ordered antiferroelectric structures. All of the structures
obey the Bernal-Fowler rules and should be able to form stable
proton-{\it disordered} crystals which have the traditional
$k_\textrm{B}$ln(3/2) residual entropy at 0~K.\cite{Bernal33,Pauling35}
Simulations of proton-disordered structures are relatively unstable
with all but the most expensive water models.\cite{Nada03} Our
simulations have therefore been performed with the ordered
antiferroelectric structures which do not require the residual entropy
term to be accounted for in the free energies. This may result in some
discrepancies when comparing our melting temperatures to the melting
temperatures that have been calculated via thermodynamic integrations
of the disordered structures.\cite{Sanz04}

Most of the water models have melting points that compare quite
favorably with the experimental value of 273~K.  The unfortunate
aspect of this result is that this phase change occurs between
Ice-{\it i} and the liquid state rather than ice I$_h$ and the liquid
state.  These results do not contradict other studies.  Studies of ice
I$_h$ using TIP4P predict a $T_m$ ranging from 191 to 238~K
(differences being attributed to choice of interaction truncation and
different ordered and disordered molecular
arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice
B and Ice-{\it i} were omitted, a $T_\textrm{m}$ value around 200~K
would be predicted from this work.  However, the $T_\textrm{m}$ from
Ice-{\it i} is calculated to be 262~K, indicating that these
simulation based structures ought to be included in studies probing
phase transitions with this model.  Also of interest in these results
is that SSD/E does not exhibit a melting point at 1 atm but does
sublime at 355~K.  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 resulted in the
spontaneous crystallization of Ice-{\it i} which led us to investigate
this structure.  These observations provide a warning that simulations
of SSD/E as a ``liquid'' near 300~K are actually metastable and run
the risk of spontaneous crystallization.  However, when a longer
cutoff radius is used, SSD/E prefers the liquid state under standard
temperature and pressure.

\section{Effects of Potential Truncation}

\begin{figure}
\includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
\caption{Free energy as a function of cutoff radius for SSD/E, TIP3P, 
SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
with an added Ewald correction term.  Error for the larger cutoff
points is equivalent to that observed at 9.0~\AA\ (see Table
\ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using
both 12 and 13.5~\AA\ cutoffs were omitted because the crystal was
prone to distortion and melting at 200~K.  Ice-$i^\prime$ is the
form of Ice-{\it i} used in the SPC/E simulations.}
\label{fig:incCutoff}
\end{figure}

For the more computationally efficient water models, we have also
investigated the effect of potential truncation on the computed free
energies as a function of the cutoff radius.  As seen in
Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with
water models lacking a long-range correction show significant cutoff
dependence.  In general, there is a narrowing of the free energy
differences while moving to greater cutoff radii.  As the free
energies for the polymorphs converge, the stability advantage that
Ice-{\it i} exhibits is reduced.  Adjacent to each of these plots are
results for systems with applied or estimated long-range corrections.
SSD/RF was parametrized for use with a reaction field, and the benefit
provided by this computationally inexpensive correction is apparent.
The free energies are largely independent of the size of the reaction
field cavity in this model, so small cutoff radii mimic bulk
calculations quite well under SSD/RF.
 
Although TIP3P was parametrized for use without the Ewald summation,
we have estimated the effect of this method for computing long-range
electrostatics for both TIP3P and SPC/E.  This was accomplished by
calculating the potential energy of identical crystals both with and
without particle mesh Ewald (PME).  Similar behavior to that observed
with reaction field is seen for both of these models.  The free
energies show reduced dependence on cutoff radius and span a narrower
range for the various polymorphs.  Like the dipolar water models,
TIP3P displays a relatively constant preference for the Ice-{\it i}
polymorph.  Crystal preference is much more difficult to determine for
SPC/E.  Without a long-range correction, each of the polymorphs
studied assumes the role of the preferred polymorph under different
cutoff radii.  The inclusion of the Ewald correction flattens and
narrows the gap in free energies such that the polymorphs are
isoenergetic within statistical uncertainty.  This suggests that other
conditions, such as the density in fixed-volume simulations, can
influence the polymorph expressed upon crystallization.

\section{Expanded Results Using Damped Shifted Force Electrostatics}


\section{Conclusions}

In this work, thermodynamic integration was used to determine the
absolute free energies of several ice polymorphs.  The new polymorph,
Ice-{\it i} was observed to be the stable crystalline state for {\it
all} the water models when using a 9.0~\AA\ cutoff.  However, the free
energy partially depends on simulation conditions (particularly on the
choice of long range correction method). Regardless, Ice-{\it i} was
still observed to be a stable polymorph for all of the studied water
models.

So what is the preferred solid polymorph for simulated water?  As
indicated above, the answer appears to be dependent both on the
conditions and the model used.  In the case of short cutoffs without a
long-range interaction correction, Ice-{\it i} and Ice-$i^\prime$ have
the lowest free energy of the studied polymorphs with all the models.
Ideally, crystallization of each model under constant pressure
conditions, as was done with SSD/E, would aid in the identification of
their respective preferred structures.  This work, however, helps
illustrate how studies involving one specific model can lead to
insight about important behavior of others.

We also note that none of the water models used in this study are
polarizable or flexible models.  It is entirely possible that the
polarizability of real water makes Ice-{\it i} substantially less
stable than ice I$_h$.  However, the calculations presented above seem
interesting enough to communicate before the role of polarizability
(or flexibility) has been thoroughly investigated. 

Finally, due to the stability of Ice-{\it i} in the investigated
simulation conditions, 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 to 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.  For the purpose of comparison with experimental
results, we have calculated the oxygen-oxygen pair correlation
function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$
for the two Ice-{\it i} variants (along with example ice
I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77~K, and they are shown in
figures \ref{fig:gofr} and \ref{fig:sofq} respectively.  It is
interesting to note that the structure factors for Ice-$i^\prime$ and
Ice-I$_c$ are quite similar.  The primary differences are small peaks
at 1.125, 2.29, and 2.53~\AA$^{-1}$, so particular attention to these
regions would be needed to identify the new $i^\prime$ variant from
the I$_\textrm{c}$ polymorph.


\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 truncation 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:	2981
Committed:	Tue Aug 29 23:34:48 2006 UTC (18 years ago) by chrisfen
Content type:	application/x-tex
File size:	25654 byte(s)
Log Message:	adding lots of tred stuff