start/scienceIcei/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}

\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 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 show 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 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.  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 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 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).  The
extruded variant, Ice-{\it i}$^\prime$, was used in calculations
involving SPC/E, TIP4P, and TIP5P due to its enhanced stability with
these models.  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.

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 varients 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}.  The reason for inclusion of ice B was that it was 
shown to be a minimum free energy structure for SPC/E at ambient
conditions.\cite{Baez95b} In addition to the free energies, the
relavent 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, and 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$ 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{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
so much that they often overlap within error, 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.

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.  Regardless of possible experimental observation,
the presence of these stable ice polymorphs has implications in the
understanding and depiction of phase changes involving the common
water models used in simulations.

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