trunk/iceiPaper/iceiPaper.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[11pt]{article}
%\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{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{Free Energy Analysis of Simulated Ice Polymorphs Using Simple 
Dipolar and Charge Based Water Models}

\author{Christopher J. Fennell and J. Daniel Gezelter \\
Department of Chemistry and Biochemistry\\ University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle
%\doublespacing

\begin{abstract}
The absolute free energies of several ice polymorphs which are stable
at low pressures were calculated using thermodynamic integration to a
reference system (the Einstein crystal).  These integrations were
performed for most of the common water models (SPC/E, TIP3P, TIP4P,
TIP5P, SSD/E, SSD/RF). Ice-{\it i}, a structure we recently observed
crystallizing at room temperature for one of the single-point water
models, was determined to be the stable crystalline state (at 1 atm)
for {\it all} the water models investigated.  Phase diagrams were
generated, and phase coexistence lines were determined for all of the
known low-pressure ice structures under all of these water models.
Additionally, potential truncation was shown to have an effect on the
calculated free energies, and can result in altered free energy
landscapes.  Structure factor predictions for the new crystal were
generated and we await experimental confirmation of the existence of
this new polymorph.
\end{abstract}

%\narrowtext 

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

\section{Introduction}

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,Berendsen98,Dill00,Mahoney00,Fennell04}
These models have been used to investigate important physical
phenomena like phase transitions, transport properties, 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,Clancy94,Mahoney01} Two
important properties to quantify are the Gibbs and Helmholtz free
energies, particularly for the solid forms of water.  Difficulty in
these types of studies typically arises from the assortment of
possible crystalline polymorphs 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 at a given state point, because
these phases will dictate the relevant transition temperatures and
pressures for the model.

In this paper, standard reference state methods were applied to known
crystalline water polymorphs in the low pressure regime.  This work is
unique in 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
(Fig. \ref{iceiCell}A) consists of eight water molecules that stack in
rows of interlocking water tetramers.  This crystal structure has a
limited resemblence to a recent two-dimensional ice tessellation
simulated on a silica surface.\cite{Yang04} Proton ordering can be
accomplished by orienting two of the molecules 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[width=\linewidth]{unitCell.eps}
\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it i}$^\prime$, 
the elongated variant of Ice-{\it i}.  The spheres represent the
center-of-mass locations of the water molecules.  The $a$ to $c$
ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
$a:2.1214c$ and $a:1.7850c$ respectively.}
\label{iceiCell}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{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 from our previous study indicated that Ice-{\it i} is the
minimum energy crystal structure for the single point water models we
had investigated (for discussions on these single point dipole models,
see our previous work and related
articles).\cite{Fennell04,Liu96,Bratko85} Those 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 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 (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). It should
be noted that a second version of Ice-{\it i} (Ice-{\it i}$^\prime$)
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.  There is typically a small
distortion of proton ordered Ice-{\it i}$^\prime$ that converts the
normally square tetramer into a rhombus with alternating approximately
85 and 95 degree angles.  The degree of this distortion is model
dependent and significant enough to split the tetramer diagonal
location peak in the radial distribution function.

\section{Methods}

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

Thermodynamic integration is an established technique for
determination of free energies of condensed phases of
materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
method, implemented in the same manner illustrated by B\`{a}ez and
Clancy, 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.\cite{Baez95a} 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 to
generate phase diagrams.  All simulations were carried out at
densities which correspond to a pressure of approximately 1 atm at
their respective temperatures.

Thermodynamic integration involves a sequence of simulations during
which the system of interest is converted into a reference system for
which the free energy is known analytically.  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 strategically 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 was chosen as the reference system.  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{r} = 4.29$ kcal
mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{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{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}$, and $E_m$ is the minimum
potential energy of the ideal crystal.\cite{Baez95a}

\begin{figure}
\includegraphics[width=\linewidth]{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}

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

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, the
effects of utilizing an Ewald summation were estimated for TIP3P and
SPC/E by performing single configuration 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_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 was found to be not as stable
as proton disordered or antiferroelectric variants of ice $I_h$.  The
proton ordered variant of ice $I_h$ used here is a simple
antiferroelectric version that we devised, and it has an 8 molecule
unit cell similar to other predicted antiferroelectric $I_h$
crystals.\cite{Davidson84} 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}
\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.  Calculated error of the final digits is in parentheses.}

\begin{tabular}{lcccc}
\hline 
Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
\hline 
TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & -13.55(2)\\
SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\
\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 calculated 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{tp3phasedia}
and \ref{ssdrfphasedia} are example diagrams built from the free
energy results.  All other models have similar structure, although the
crossing points between the phases move to 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 at higher pressures than those studied here.\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}
\begin{center}

\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}{lccccccc}
\hline 
Equilibrium Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
\hline 
$T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
$T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
$T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
\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 other
studies.  Studies of ice $I_h$ 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).\cite{Vlot99,Gao00,Sanz04} 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 to be 265 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 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 resulted in spontaneous
crystallization of Ice-{\it i} and 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, this risk lessens when
applying a longer cutoff.

\begin{figure}
\includegraphics[width=\linewidth]{cutoffChange.eps}
\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.  Calculations performed without a
long-range correction show noticable free energy dependence on the
cutoff radius and show some degree of converge at large cutoff radii.
Inclusion of a long-range correction reduces the cutoff radius
dependence of the free energy for all the models.  Data for ice I$_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-{\it
i}$^\prime$ is the form of Ice-{\it i} used in the SPC/E simulations.}
\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 the ice polymorphs with water models lacking a
long-range correction show a significant cutoff radius 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.  Interestingly, increasing the cutoff radius a mere 1.5
\AA\ with the SSD/E model destabilizes the Ice-{\it i} polymorph
enough 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.  Considering the stabilization of Ice-{\it i} with
smaller cutoffs, it is not surprising that crystallization was
observed with SSD/E.  The choice of a 9 \AA\ cutoff in the previous
simulations gives the Ice-{\it i} polymorph a greater than 1 kcal/mol
lower free energy than the ice $I_\textrm{h}$ starting configurations.
Additionally, it should be noted that ice $I_c$ is not stable with
cutoff radii of 12 and 13.5 \AA\ using the TIP3P water model.  These
simulations showed bulk distortions of the simulation cell that
rapidly deteriorated crystal integrity.

Adjacent to each of these model plots is a system with an applied or
estimated long-range correction.  SSD/RF was parametrized for use with
a reaction field, and the benefit provided by this computationally
inexpensive correction is apparent.  Due to the relative independence
of the resultant free energies, calculations performed with a small
cutoff radius provide resultant properties similar to what one would
expect for the bulk material.  In the cases of TIP3P and SPC/E, 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 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 less dependence on cutoff radius and span a more narrowed 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 conditions.  The inclusion of the Ewald correction flattens and
narrows the sequences of free energies so much that they often overlap
within error (see Table \ref{spcecut}), indicating that other
conditions, such as cell volume in microcanonical simulations, can
influence the chosen polymorph upon crystallization.  All of these
results support 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}
\begin{center}

\caption{The free energy versus cutoff radius for the studied ice 
polymorphs using SPC/E after the inclusion of the PME long-range
interaction correction. Units are kcal/mol.}

\begin{tabular}{ccccc}
\hline 
Cutoff (\AA) & $I_h$ & $I_c$ & B & Ice-{\it i} \\
\hline 
9.0   & -12.98(2) & -13.00(2) & -12.97(3) & -13.02(2) \\
10.5  & -13.13(3) & -13.09(4) & -13.17(3) & -13.11(2) \\
12.0  & -13.06(2) & -13.09(2) & -13.15(4) & -13.12(2) \\
13.5  & -13.02(2) & -13.02(2) & -13.08(2) & -13.07(2) \\
15.0  & -13.11(4) & -12.97(2) & -13.09(2) & -12.95(2) \\
\end{tabular}
\label{spcecut}
\end{center}
\end{minipage}
\end{table*}

\section{Conclusions}

The free energy for proton ordered variants of hexagonal and cubic ice
$I$, ice B, and our recently discovered Ice-{\it i} structure 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 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.
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 of 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 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_h$, $I_c$, and for
ice-{\it i} at a temperature of 77K.  In studies of the high and low
density forms of amorphous ice, ``spurious'' diffraction peaks have
been observed experimentally.\cite{Bizid87} It is possible that a
variant of Ice-{\it i} could explain some of this behavior; however,
we will leave it to our experimental colleagues to make the final
determination on whether this ice polymorph is named appropriately
(i.e. with an imaginary number) or if it can be promoted to Ice-0.

\begin{figure}
\includegraphics[width=\linewidth]{iceGofr.eps}
\caption{Radial distribution functions of ice $I_h$, $I_c$, 
Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
of the SSD/RF water model at 77 K.}
\label{fig:gofr}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{sofq.eps}
\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
 and Ice-{\it i}$^\prime$ 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.}
\label{fig:sofq}
\end{figure}

\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:	1831
Committed:	Thu Dec 2 14:23:11 2004 UTC (19 years, 7 months ago) by chrisfen
Content type:	application/x-tex
File size:	28157 byte(s)
Log Message:	new title and text adjustments