trunk/iceiPaper/iceiPaper.tex


%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[preprint,aps,endfloats]{revtex4}
%\documentclass[11pt]{article}
%\usepackage{endfloat}
\usepackage{amsmath}
\usepackage{epsf}
\usepackage{berkeley}
%\usepackage{setspace}
%\usepackage{tabularx}
\usepackage{graphicx}
%\usepackage[ref]{overcite}
%\pagestyle{plain}
%\pagenumbering{arabic}
%\oddsidemargin 0.0cm \evensidemargin 0.0cm
%\topmargin -21pt \headsep 10pt
%\textheight 9.0in \textwidth 6.5in
%\brokenpenalty=10000

%\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{A Free Energy Study of Low Temperature and Anomalous Ice}

\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}

\address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

%\maketitle
%\doublespacing

\begin{abstract}
The free energies of several ice polymorphs in the low pressure regime
were calculated using thermodynamic integration of systems consisting
of a variety of common water models. Ice-{\it i}, a recent
computationally observed solid structure, was determined to be the
stable state with the lowest free energy for all the water models
investigated. Phase diagrams were generated, and melting and boiling
points for all the models were determined and show relatively good
agreement with experiment, although the solid phase is different
between simulation and experiment. In addition, potential truncation
was shown to have an effect on the calculated free energies, and may
result in altered free energy landscapes.
\end{abstract}

\maketitle

\newpage

%\narrowtext 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

Molecular dynamics has developed into a valuable tool for studying the
phase behavior of systems ranging from small or simple
molecules\cite{Matsumoto02andOthers} to complex biological
species.\cite{bigStuff} Many techniques have been developed in order
to investigate the thermodynamic properites of model substances,
providing both qualitative and quantitative comparisons between
simulations and experiment.\cite{thermMethods} Investigation of these
properties leads to the development of new and more accurate models,
leading to better understanding and depiction of physical processes
and intricate molecular systems.

Water has proven to be a challenging substance to depict in
simulations, and has resulted in a variety of models that attempt to
describe its behavior under a varying simulation
conditions.\cite{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
Many of these models have been used to investigate important physical
phenomena like phase transitions and the hydrophobic
effect.\cite{evenMorePapers} With the advent of numerous differing
models, it is only natural that attention is placed on the properties
of the models themselves in an attempt to clarify their benefits and
limitations when applied to a system of interest.\cite{modelProps} One
important but challenging property to quantify is the free energy,
particularly of the solid forms of water. Difficulty in these types of
studies typically arises from the assortment of possible crystalline
polymorphs that water that water adopts over a wide range of pressures
and temperatures. There are currently 13 recognized forms of ice, and
it is a challenging task to investigate the entire free energy
landscape.\cite{Sanz04} Ideally, research is focused on the phases
having the lowest free energy, because these phases will dictate the
true transition temperatures and pressures for their respective model.

In this paper, standard reference state methods were applied to the
study of crystalline water polymorphs in the low pressure regime. 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 was different from any previously observed ice
polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
was termed Ice-{\it i} in homage to its origin in computational
simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
water molecules that stack in rows of interlocking water
tetramers. Proton ordering can be accomplished by orienting two of the
waters so that both of their donated hydrogen bonds are internal to
their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
constructed of water tetramers, the hydrogen bonds are not as linear
as those observed in ice $I_h$, however the interlocking of these
subunits appears to provide significant stabilization to the overall
crystal. The arrangement of these tetramers results in surrounding
open octagonal cavities that are typically greater than 6.3 \AA\ in
diameter. This relatively open overall structure leads to crystals
that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
\begin{figure}
\includegraphics[scale=1.0]{unitCell.eps}
\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-2{\it i}, 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-2{\it i}, $a = 1.7850c$.}
\label{iceiCell}
\end{figure}
\begin{figure}
\includegraphics[scale=1.0]{orderedIcei.eps}
\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_h$.}
\label{protOrder}
\end{figure}

Results in the previous study indicated that Ice-{\it i} is the
minimum energy crystal structure for the single point water models
being studied (for discussions on these single point dipole models,
see the previous work and related
articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
consider energetic stabilization and neglect entropic contributions to
the overall free energy. To address this issue, the absolute free
energy of this crystal was calculated using thermodynamic integration
and compared to the free energies of cubic and hexagonal ice $I$ (the
experimental 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
(soft sticky dipole extended, 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 (soft sticky dipole
reaction field, SSD/RF). In should be noted that a second version of
Ice-{\it i} (Ice-2{\it i}) was used in calculations involving SPC/E,
TIP4P, and TIP5P. The unit cell of this crystal (Fig. \ref{iceiCell}B)
is similar to the Ice-{\it i} unit it is extended in the direction of
the (001) face and compressed along the other two faces.

\section{Methods}

Canonical ensemble (NVT) molecular dynamics calculations were
performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
molecular mechanics package. All molecules were treated as rigid
bodies, with orientational motion propagated using the symplectic DLM
integration method. Details about the implementation of these
techniques can be found in a recent publication.\cite{Meineke05} 

Thermodynamic integration was utilized to calculate the free energy of
several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
400 K for all of these water models were also determined using this
same technique, in order to determine melting points and generate
phase diagrams. All simulations were carried out at densities
resulting in a pressure of approximately 1 atm at their respective
temperatures.

A single thermodynamic integration involves a sequence of simulations
over which the system of interest is converted into a reference system
for which the free energy is known. 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 300 ps (for the unaltered system) to 75 ps
(near the reference state) in length.

For the thermodynamic integration of molecular crystals, the Einstein
Crystal is chosen as the reference state that the system is converted
to over the course of the simulation. 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 has been evaluated, 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{ecFreeEnergy}
\end{eqnarray}
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
\ref{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}
\includegraphics[scale=1.0]{rotSpring.eps}
\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 poor energy conserving dynamics resulting
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} TINKER was chosen because it
can also propagate the motion of rigid-bodies, and provides the most
direct comparison to the results from OOPSE. The calculated energy
difference in the presence and absence of PME was applied to the
previous results in order to predict changes in 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_h$ and $I_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 (Table \ref{freeEnergy}).\cite{Baez95b}
Ice XI, the experimentally observed proton ordered variant of ice
$I_h$, was investigated initially, but it was found not to be as
stable as antiferroelectric variants of proton ordered or even proton
disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
ice $I_h$ used here is a simple antiferroelectric version that has an
8 molecule unit cell. The crystals contained 648 or 1728 molecules for
ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
$I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
were necessary for simulations involving larger cutoff values.

\begin{table*}
\begin{minipage}{\linewidth}
\renewcommand{\thefootnote}{\thempfootnote}
\begin{center}
\caption{Calculated free energies for several ice polymorphs with a 
variety of common water models. All calculations used a cutoff radius
of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
\begin{tabular}{ l  c  c  c  c }
\hline \\[-7mm]
\ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
\hline \\[-3mm]
\ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
\ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
\ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
\ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
\ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
\ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
\end{tabular}
\label{freeEnergy}
\end{center}
\end{minipage}
\end{table*}

The free energy values computed for the studied polymorphs indicate
that Ice-{\it i} is the most stable state for all of the common water
models studied. With the free energy at these state points, the
temperature and pressure dependence of the free energy was used to
project to other state points and build phase diagrams. Figures
\ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
from the free energy results. All other models have similar structure,
only the crossing points between these phases exist at 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}
\includegraphics[width=\linewidth]{tp3PhaseDia.eps}
\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{tp3phasedia}
\end{figure}
\begin{figure}
\includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
\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{ssdrfphasedia}
\end{figure}

\begin{table*}
\begin{minipage}{\linewidth}
\renewcommand{\thefootnote}{\thempfootnote}
\begin{center}
\caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) 
temperatures of several common water models compared with experiment.}
\begin{tabular}{ l  c  c  c  c  c  c  c }
\hline \\[-7mm]
\ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
\hline \\[-3mm]
\ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
\ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
\ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
\end{tabular}
\label{meltandboil}
\end{center}
\end{minipage}
\end{table*}

Table \ref{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_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 238 K (differences
being attributed to choice of interaction truncation and different
ordered and disordered molecular arrangements). If the presence of ice
B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
predicted from this work. However, the $T_m$ from Ice-{\it i} is
calculated at 265 K, 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 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 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 300 K are actually metastable and run the
risk of spontaneous crystallization. However, this risk changes when
applying a longer cutoff.

\begin{figure}
\includegraphics[width=\linewidth]{cutoffChange.eps}
\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_c$ 12
\AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
\AA\. These crystals are unstable at 200 K and rapidly convert into a
liquid. The connecting lines are qualitative visual aid.}
\label{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 Fig. \ref{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 (298 K 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{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 cut that
Ice-{\it i} is the most stable polymorph, and is possibly metastable
with respect to ice B and possibly ice $I_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*}
\begin{minipage}{\linewidth}
\renewcommand{\thefootnote}{\thempfootnote}
\begin{center}
\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 }
\hline \\[-7mm]
\ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
\hline \\[-3mm]
\ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
\ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
\end{tabular}
\label{pmeShift}
\end{center}
\end{minipage}
\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} where 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_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
correction. 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 leads the authors to be hesitant in ascribing relevance outside
of computational models, hence the descriptive name presented. That
being 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.

\section{Acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0134881. Computation time was provided by
the Notre Dame High Performance Computing Cluster and the Notre Dame
Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).

\newpage

\bibliographystyle{jcp}
\bibliography{iceiPaper}


\end{document}
Revision:	1462
Committed:	Wed Sep 15 20:47:13 2004 UTC (19 years, 9 months ago) by chrisfen
Content type:	application/x-tex
File size:	24180 byte(s)
Log Message:	Added references