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{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{iCrystal}A and
\ref{iCrystal}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 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}
\centering
\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 2.1214
and 1.785 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} 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.

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} 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.  There is some subtlety in the path
taken as $\lambda$ goes from $0$ to $1$.  However, 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.  These results do not contradict 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 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.

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.

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.

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