trunk/scienceIcei/iceiPaper.tex

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

\begin{document}

\title{Free Energy Analysis of Simulated Ice Polymorphs}

\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 with
a variety of common water models.  A recently discovered ice polymorph
that has as yet 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.  Difficulties in studies
addressing these thermodynamic quantities typically arise from the
assortment of possible crystalline polymorphs that water adopts over a
wide range of pressures and temperatures.  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.

In this paper, standard reference state methods were applied to known
crystalline water polymorphs to evaluate their free energy 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 of Ice-{\it i} and an extruded 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{iCrystal}A and
\ref{iCrystal}B.  These tetramers make the crystal structure similar
in appearance to a recent two-dimensional ice tessellation simulated
on a silica surface.\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 surrounding
open octagonal cavities that are typically greater than 6.3 \AA\ in
diameter (Fig. \ref{iCrystal}C).  This open structure leads to
crystals that are typically 0.07 g/cm$^3$ less dense than ice $I_h$.

\begin{figure}
\includegraphics[width=4in]{iCrystal.eps}
\caption{(A) Unit cell for Ice-{\it i}, (B) Ice-{\it i}$^\prime$, 
and (C) a rendering of a proton ordered crystal of Ice-{\it i} looking
down the (001) crystal face.  In the unit cells, 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.785c$ respectively. The presence of large
octagonal pores in both crystal forms lead 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} 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 it to the free energies of ice
$I_c$ and ice $I_h$ (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).  The extruded variant,
Ice-{\it i}$^\prime$, was used in calculations involving SPC/E, TIP4P,
and TIP5P.  These models exhibit enhanced stability with Ice-{\it
i}$^\prime$ because of their more tetrahedrally arranged internal
charge distributions.  Additionally, there is often 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.

Thermodynamic integration was utilized to calculate the free energies
of the listed water models at various state points using a modified
form of the OOPSE molecular dynamics package.\cite{Meineke05}
This calculation method 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.  For
liquid and solid phases, the ideal gas and harmonically restrained
crystal are chosen as the reference states respectively. 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}. 

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(4) & 357(2)\\
TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 266(5) & 354(2)\\
TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 271(4) & 337(2)\\
SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 296(3) & 396(2)\\
SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(2) & -\\
SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2) & - & 278(4) & 349(2)\\
\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.  Surprisingly, these results are 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 rather 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, when applying a
longer cutoff, the liquid state is preferred under standard
conditions.

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

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.  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
such that they often overlap within error, indicating that other
conditions, such as the density in fixed volume simulations, can
influence the chosen polymorph upon crystallization.  

So what is the preferred solid polymorph for simulated water?  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.  In general, the above 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.

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.

\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:	1888
Committed:	Wed Dec 15 18:47:55 2004 UTC (20 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	17208 byte(s)
Log Message:	fixed grammer and typos
#	User	Rev	Content
1	chrisfen	1845	%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2			\documentclass[11pt]{article}
3			\usepackage{endfloat}
4			\usepackage{amsmath}
5			\usepackage{epsf}
6			\usepackage{berkeley}
7			\usepackage{setspace}
8			\usepackage{tabularx}
9			\usepackage{graphicx}
10			\usepackage[ref]{overcite}
11			\pagestyle{plain}
12			\pagenumbering{arabic}
13			\oddsidemargin 0.0cm \evensidemargin 0.0cm
14			\topmargin -21pt \headsep 10pt
15			\textheight 9.0in \textwidth 6.5in
16			\brokenpenalty=10000
17			\renewcommand{\baselinestretch}{1.2}
18			\renewcommand\citemid{\ } % no comma in optional reference note
19
20			\begin{document}
21
22			\title{Free Energy Analysis of Simulated Ice Polymorphs}
23
24			\author{Christopher J. Fennell and J. Daniel Gezelter \\
25			Department of Chemistry and Biochemistry\\ University of Notre Dame\\
26			Notre Dame, Indiana 46556}
27
28			\date{\today}
29
30			\maketitle
31			%\doublespacing
32
33			\begin{abstract}
34			The absolute free energies of several ice polymorphs which are stable
35			at low pressures were calculated using thermodynamic integration with
36			a variety of common water models. A recently discovered ice polymorph
37	chrisfen	1888	that has as yet only been observed in computer simulations (Ice-{\it
38			i}), was determined to be the stable crystalline state for {\it all}
39			the water models investigated. Phase diagrams were generated, and
40			phase coexistence lines were determined for all of the known
41			low-pressure ice structures. Additionally, potential truncation was
42			shown to play a role in the resulting shape of the free energy
43			landscape.
44	chrisfen	1845	\end{abstract}
45
46			%\narrowtext
47
48			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49			% BODY OF TEXT
50			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52			\section{Introduction}
53
54			Water has proven to be a challenging substance to depict in
55			simulations, and a variety of models have been developed to describe
56			its behavior under varying simulation
57			conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
58			These models have been used to investigate important physical
59			phenomena like phase transitions, transport properties, and the
60			hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
61			choice of models available, it is only natural to compare the models
62			under interesting thermodynamic conditions in an attempt to clarify
63	chrisfen	1888	the limitations of
64			each.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two important
65			properties to quantify are the Gibbs and Helmholtz free energies,
66			particularly for the solid forms of water. Difficulties in studies
67			addressing these thermodynamic quantities typically arise from the
68			assortment of possible crystalline polymorphs that water adopts over a
69			wide range of pressures and temperatures. It is a challenging task to
70			investigate the entire free energy landscape\cite{Sanz04}; and
71			ideally, research is focused on the phases having the lowest free
72	chrisfen	1845	energy at a given state point, because these phases will dictate the
73			relevant transition temperatures and pressures for the model.
74
75			In this paper, standard reference state methods were applied to known
76	chrisfen	1860	crystalline water polymorphs to evaluate their free energy in the low
77			pressure regime. This work is unique in that one of the crystal
78			lattices was arrived at through crystallization of a computationally
79			efficient water model under constant pressure and temperature
80			conditions. Crystallization events are interesting in and of
81			themselves\cite{Matsumoto02,Yamada02}; however, the crystal structure
82			obtained in this case is different from any previously observed ice
83			polymorphs in experiment or simulation.\cite{Fennell04} We have named
84			this structure Ice-{\it i} to indicate its origin in computational
85			simulation. The unit cell of Ice-{\it i} and an extruded variant named
86			Ice-{\it i}$^\prime$ both consist of eight water molecules that stack
87			in rows of interlocking water tetramers as illustrated in figures
88			\ref{iCrystal}A and
89	chrisfen	1845	\ref{iCrystal}B. These tetramers make the crystal structure similar
90			in appearance to a recent two-dimensional ice tessellation simulated
91			on a silica surface.\cite{Yang04} As expected in an ice crystal
92			constructed of water tetramers, the hydrogen bonds are not as linear
93			as those observed in ice $I_h$, however the interlocking of these
94			subunits appears to provide significant stabilization to the overall
95			crystal. The arrangement of these tetramers results in surrounding
96			open octagonal cavities that are typically greater than 6.3 \AA\ in
97			diameter (Fig. \ref{iCrystal}C). This open structure leads to
98			crystals that are typically 0.07 g/cm$^3$ less dense than ice $I_h$.
99
100			\begin{figure}
101			\includegraphics[width=4in]{iCrystal.eps}
102			\caption{(A) Unit cell for Ice-{\it i}, (B) Ice-{\it i}$^\prime$,
103			and (C) a rendering of a proton ordered crystal of Ice-{\it i} looking
104			down the (001) crystal face. In the unit cells, the spheres represent
105			the center-of-mass locations of the water molecules. The $a$ to $c$
106			ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
107			$a:2.1214c$ and $a:1.785c$ respectively. The presence of large
108			octagonal pores in both crystal forms lead to a polymorph that is less
109			dense than ice $I_h$.}
110			\label{iCrystal}
111			\end{figure}
112
113			Results from our previous study indicated that Ice-{\it i} is the
114			minimum energy crystal structure for the single point water models
115			investigated (for discussions on these single point dipole models, see
116			our previous work and related
117			articles).\cite{Fennell04,Liu96,Bratko85} Those results only
118			considered energetic stabilization and neglected entropic
119			contributions to the overall free energy. To address this issue, we
120			have calculated the absolute free energy of this crystal using
121	chrisfen	1888	thermodynamic integration and compared it to the free energies of ice
122			$I_c$ and ice $I_h$ (the experimental low density ice polymorphs) and
123			ice B (a higher density, but very stable crystal structure observed by
124			B\`{a}ez and Clancy in free energy studies of SPC/E).\cite{Baez95b}
125			This work includes results for the water model from which Ice-{\it i}
126			was crystallized (SSD/E) in addition to several common water models
127			(TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field parametrized
128			single point dipole water model (SSD/RF). The extruded variant,
129			Ice-{\it i}$^\prime$, was used in calculations involving SPC/E, TIP4P,
130			and TIP5P. These models exhibit enhanced stability with Ice-{\it
131			i}$^\prime$ because of their more tetrahedrally arranged internal
132			charge distributions. Additionally, there is often a small distortion
133			of proton ordered Ice-{\it i}$^\prime$ that converts the normally
134			square tetramer into a rhombus with alternating approximately 85 and
135			95 degree angles. The degree of this distortion is model dependent
136			and significant enough to split the tetramer diagonal location peak in
137			the radial distribution function.
138	chrisfen	1845
139			Thermodynamic integration was utilized to calculate the free energies
140			of the listed water models at various state points using a modified
141			form of the OOPSE molecular dynamics package.\cite{Meineke05}
142			This calculation method involves a sequence of simulations during
143			which the system of interest is converted into a reference system for
144			which the free energy is known analytically. This transformation path
145	chrisfen	1872	is then integrated, in order to determine the free energy difference
146	chrisfen	1845	between the two states:
147			\begin{equation}
148			\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
149			)}{\partial\lambda}\right\rangle_\lambda d\lambda,
150			\end{equation}
151			where $V$ is the interaction potential and $\lambda$ is the
152			transformation parameter that scales the overall potential. For
153			liquid and solid phases, the ideal gas and harmonically restrained
154			crystal are chosen as the reference states respectively. Thermodynamic
155			integration is an established technique that has been used extensively
156			in the calculation of free energies for condensed phases of
157			materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}.
158
159	chrisfen	1888	The calculated free energies of proton-ordered variants of three low
160	chrisfen	1845	density polymorphs ($I_h$, $I_c$, and Ice-{\it i} or Ice-{\it
161			i}$^\prime$) and the stable higher density ice B are listed in Table
162	chrisfen	1888	\ref{freeEnergy}. Ice B was included because it has been
163	chrisfen	1845	shown to be a minimum free energy structure for SPC/E at ambient
164			conditions.\cite{Baez95b} In addition to the free energies, the
165	chrisfen	1888	relevant transition temperatures at standard pressure are also
166	chrisfen	1845	displayed in Table \ref{freeEnergy}. These free energy values
167			indicate that Ice-{\it i} is the most stable state for all of the
168			investigated water models. With the free energy at these state
169			points, the Gibbs-Helmholtz equation was used to project to other
170	chrisfen	1888	state points and to build phase diagrams. Figure \ref{tp3PhaseDia} is
171			an example diagram built from the results for the TIP3P water model.
172			All other models have similar structure, although the crossing points
173			between the phases move to different temperatures and pressures as
174			indicated from the transition temperatures in Table \ref{freeEnergy}.
175			It is interesting to note that ice $I_h$ (and ice $I_c$ for that
176			matter) do not appear in any of the phase diagrams for any of the
177			models. For purposes of this study, ice B is representative of the
178			dense ice polymorphs. A recent study by Sanz {\it et al.} provides
179			details on the phase diagrams for SPC/E and TIP4P at higher pressures
180			than those studied here.\cite{Sanz04}
181	chrisfen	1845
182			\begin{table*}
183			\begin{minipage}{\linewidth}
184			\begin{center}
185			\caption{Calculated free energies for several ice polymorphs along
186			with the calculated melting (or sublimation) and boiling points for
187			the investigated water models. All free energy calculations used a
188			cutoff radius of 9.0 \AA\ and were performed at 200 K and $\sim$1 atm.
189			Units of free energy are kcal/mol, while transition temperature are in
190			Kelvin. Calculated error of the final digits is in parentheses.}
191			\begin{tabular}{lccccccc}
192			\hline
193			Water Model & $I_h$ & $I_c$ & B & Ice-{\it i} & Ice-{\it i}$^\prime$ & $T_m$ (*$T_s$) & $T_b$\\
194			\hline
195			TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(4) & 357(2)\\
196			TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 266(5) & 354(2)\\
197			TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 271(4) & 337(2)\\
198			SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 296(3) & 396(2)\\
199			SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(2) & -\\
200			SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2) & - & 278(4) & 349(2)\\
201			\end{tabular}
202			\label{freeEnergy}
203			\end{center}
204			\end{minipage}
205			\end{table*}
206
207			\begin{figure}
208			\includegraphics[width=\linewidth]{tp3PhaseDia.eps}
209			\caption{Phase diagram for the TIP3P water model in the low pressure
210			regime. The displayed $T_m$ and $T_b$ values are good predictions of
211			the experimental values; however, the solid phases shown are not the
212			experimentally observed forms. Both cubic and hexagonal ice $I$ are
213			higher in energy and don't appear in the phase diagram.}
214			\label{tp3PhaseDia}
215			\end{figure}
216
217			Most of the water models have melting points that compare quite
218			favorably with the experimental value of 273 K. The unfortunate
219			aspect of this result is that this phase change occurs between
220			Ice-{\it i} and the liquid state rather than ice $I_h$ and the liquid
221			state. Surprisingly, these results are not contrary to other studies.
222			Studies of ice $I_h$ using TIP4P predict a $T_m$ ranging from 214 to
223			238 K (differences being attributed to choice of interaction
224			truncation and different ordered and disordered molecular
225			arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
226			Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
227			predicted from this work. However, the $T_m$ from Ice-{\it i} is
228			calculated to be 265 K, indicating that these simulation based
229			structures ought to be included in studies probing phase transitions
230			with this model. Also of interest in these results is that SSD/E does
231			not exhibit a melting point at 1 atm, but it rather shows a
232			sublimation point at 355 K. This is due to the significant stability
233			of Ice-{\it i} over all other polymorphs for this particular model
234			under these conditions. While troubling, this behavior resulted in
235			spontaneous crystallization of Ice-{\it i} and led us to investigate
236			this structure. These observations provide a warning that simulations
237			of SSD/E as a ``liquid'' near 300 K are actually metastable and run
238			the risk of spontaneous crystallization. However, when applying a
239			longer cutoff, the liquid state is preferred under standard
240			conditions.
241
242			\begin{figure}
243			\includegraphics[width=\linewidth]{cutoffChange.eps}
244			\caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
245			SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
246			with an added Ewald correction term. Error for the larger cutoff
247			points is equivalent to that observed at 9.0\AA\ (see Table
248			\ref{freeEnergy}). Data for ice I$_c$ with TIP3P using both 12 and
249			13.5 \AA\ cutoffs were omitted because the crystal was prone to
250			distortion and melting at 200 K. Ice-{\it i}$^\prime$ is the form of
251			Ice-{\it i} used in the SPC/E simulations.}
252			\label{incCutoff}
253			\end{figure}
254
255			Increasing the cutoff radius in simulations of the more
256			computationally efficient water models was done in order to evaluate
257			the trend in free energy values when moving to systems that do not
258			involve potential truncation. As seen in Fig. \ref{incCutoff}, the
259			free energy of the ice polymorphs with water models lacking a
260			long-range correction show a significant cutoff radius dependence. In
261			general, there is a narrowing of the free energy differences while
262			moving to greater cutoff radii. As the free energies for the
263			polymorphs converge, the stability advantage that Ice-{\it i} exhibits
264			is reduced. Adjacent to each of these model plots is a system with an
265			applied or estimated long-range correction. SSD/RF was parametrized
266			for use with a reaction field, and the benefit provided by this
267			computationally inexpensive correction is apparent. Due to the
268			relative independence of the resultant free energies, calculations
269			performed with a small cutoff radius provide resultant properties
270			similar to what one would expect for the bulk material. In the cases
271			of TIP3P and SPC/E, the effect of an Ewald summation was estimated by
272			applying the potential energy difference do to its inclusion in
273			systems in the presence and absence of the correction. This was
274			accomplished by calculation of the potential energy of identical
275			crystals both with and without particle mesh Ewald (PME). Similar
276			behavior to that observed with reaction field is seen for both of
277			these models. The free energies show less dependence on cutoff radius
278			and span a more narrowed range for the various polymorphs. Like the
279			dipolar water models, TIP3P displays a relatively constant preference
280			for the Ice-{\it i} polymorph. Crystal preference is much more
281			difficult to determine for SPC/E. Without a long-range correction,
282			each of the polymorphs studied assumes the role of the preferred
283			polymorph under different cutoff conditions. The inclusion of the
284			Ewald correction flattens and narrows the sequences of free energies
285	chrisfen	1872	such that they often overlap within error, indicating that other
286			conditions, such as the density in fixed volume simulations, can
287	chrisfen	1860	influence the chosen polymorph upon crystallization.
288	chrisfen	1845
289	chrisfen	1860	So what is the preferred solid polymorph for simulated water? The
290	chrisfen	1888	answer appears to be dependent both on the conditions and the model
291	chrisfen	1860	used. In the case of short cutoffs without a long-range interaction
292			correction, Ice-{\it i} and Ice-{\it i}$^\prime$ have the lowest free
293			energy of the studied polymorphs with all the models. Ideally,
294			crystallization of each model under constant pressure conditions, as
295			was done with SSD/E, would aid in the identification of their
296			respective preferred structures. This work, however, helps illustrate
297			how studies involving one specific model can lead to insight about
298			important behavior of others. In general, the above results support
299			the finding that the Ice-{\it i} polymorph is a stable crystal
300			structure that should be considered when studying the phase behavior
301			of water models.
302
303			Finally, due to the stability of Ice-{\it i} in the investigated
304			simulation conditions, the question arises as to possible experimental
305			observation of this polymorph. The rather extensive past and current
306			experimental investigation of water in the low pressure regime makes
307	chrisfen	1888	us hesitant to ascribe any relevance to this work outside of the
308	chrisfen	1860	simulation community. It is for this reason that we chose a name for
309			this polymorph which involves an imaginary quantity. That said, there
310			are certain experimental conditions that would provide the most ideal
311			situation for possible observation. These include the negative
312			pressure or stretched solid regime, small clusters in vacuum
313	chrisfen	1845	deposition environments, and in clathrate structures involving small
314	chrisfen	1860	non-polar molecules.
315	chrisfen	1845
316			\section{Acknowledgments}
317			Support for this project was provided by the National Science
318			Foundation under grant CHE-0134881. Computation time was provided by
319			the Notre Dame High Performance Computing Cluster and the Notre Dame
320			Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
321
322			\newpage
323
324			\bibliographystyle{jcp}
325			\bibliography{iceiPaper}
326
327
328			\end{document}