trunk/iceiPaper/iceiPaper.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[12pt]{article}
\usepackage{endfloat}
\usepackage{amsmath}
\usepackage{epsf}
\usepackage{times}
\usepackage{mathptm}
\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{Computational free energy studies of a new ice polymorph which
exhibits greater stability than Ice I$_h$}

\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 were calculated
using thermodynamic integration.  These polymorphs are predicted by
computer simulations using a variety of common water models to be
stable at low pressures.  A recently discovered ice polymorph that has
as yet {\it only} been observed in computer simulations (Ice-{\it i}),
was determined to be the stable crystalline state 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.  Additionally, potential truncation was shown to play
a role in the resulting shape of the free energy landscape.
\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.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two important
properties 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.  It is a challenging task to investigate the entire free
energy landscape\cite{Sanz04}; and 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 paper, 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 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-{\it i} and an axially-elongated
variant named Ice-{\it i}$^\prime$ both consist of eight water
molecules that stack in rows of interlocking water tetramers as
illustrated in figures \ref{unitcell}A and \ref{unitcell}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$_h$,
however the interlocking of these subunits appears to provide
significant stabilization to the overall crystal.  The arrangement of
these tetramers results in octagonal cavities that are typically
greater than 6.3 \AA\ in diameter (Fig. \ref{iCrystal}).  This open
structure leads to crystals that are typically 0.07 g/cm$^3$ less
dense than ice I$_h$.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{unitCell.eps}
\caption{(A) Unit cells for Ice-{\it i} and (B) Ice-{\it i}$^\prime$.  
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 2.1214 and 1.785 respectively.}
\label{unitcell}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{orderedIcei.eps}
\caption{A rendering of a proton ordered crystal of Ice-{\it i} looking
down the (001) crystal face.  The presence of large octagonal pores
leads to a polymorph that is less dense than ice I$_h$.}
\label{iCrystal}
\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
our previous work and related
articles).\cite{Fennell04,Liu96,Bratko85} Our earlier results
considered only 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 it to the free energies of ice
I$_c$ and ice I$_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-{\it i}$^\prime$, was used in calculations involving
SPC/E, TIP4P, and TIP5P.  The square tetramers in Ice-{\it i} distort
in Ice-{\it 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 a peak in the radial distribution function which corresponds
to diagonal sites in the tetramers.

\section{Methods}

Canonical ensemble (NVT) molecular dynamics calculations were
performed using the OOPSE molecular mechanics program.\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 was utilized to calculate the Helmholtz free
energies ($A$) of the listed water models at various state points
using the OOPSE molecular dynamics program.\cite{Meineke05}
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 during 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 (Einsten) crystal is
chosen as the reference state, while for liquid phases, the ideal gas
is taken as the reference state.  

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{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}
\centering
\includegraphics[width=4in]{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 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~\cite{Tinker} for each of the ice polymorphs.  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 calculated free energies of proton-ordered variants of three low
density polymorphs (I$_h$, I$_c$, and Ice-{\it i} or Ice-{\it
i}$^\prime$) and the stable higher density ice B are listed in Table
\ref{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{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.  Figure \ref{tp3PhaseDia} is
an example diagram built from the results for the TIP3P water model.
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{freeEnergy}.
It is interesting to note that ice I$_h$ (and ice I$_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*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Calculated free energies for several ice polymorphs along 
with the calculated melting (or sublimation) and boiling points for
the investigated water models.  All free energy calculations used a
cutoff radius of 9.0 \AA\ and were performed at 200 K and $\sim$1 atm.
Units of free energy are kcal/mol, while transition temperature are in
Kelvin.  Calculated error of the final digits is in parentheses.}
\begin{tabular}{lccccccc}
\hline 
Water Model & I$_h$ & I$_c$ & B & Ice-{\it i} & Ice-{\it i}$^\prime$ & $T_m$ (*$T_s$) & $T_b$\\
\hline 
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)\\
\end{tabular}
\label{freeEnergy}
\end{center}
\end{minipage}
\end{table*}

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

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_m$ value around 200 K would be
predicted from this work.  However, the $T_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.

\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.  Error for the larger cutoff
points is equivalent to that observed at 9.0\AA\ (see Table
\ref{freeEnergy}).  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}

For the more computationally efficient water models, we have also
investigated the effect of potential trunctaion on the computed free
energies as a function of the cutoff radius.  As seen in
Fig. \ref{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 paramaterized 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{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 observered 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-{\it
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_{OO}(r)$, and the structure factor, $S(\vec{q})$ for the
two Ice-{\it i} variants (along with example ice I$_h$ and I$_c$
plots) at 77K, 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-{\it 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 {\it i}$^\prime$ variant from the I$_{c}$ variant.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{iceGofr.eps}
\caption{Radial distribution functions of ice I$_h$, I$_c$, and
Ice-{\it i} calculated from from simulations of the SSD/RF water model
at 77 K.  The Ice-{\it i} distribution function was obtained from
simulations composed of TIP4P water.}
\label{fig:gofr}
\end{figure}

\begin{figure}
\centering
\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{achemso}
\bibliography{iceiPaper}


\end{document}
Revision:	2134
Committed:	Fri Mar 25 19:11:49 2005 UTC (19 years, 3 months ago) by chrisfen
Content type:	application/x-tex
File size:	24317 byte(s)
Log Message:	added a reference