--- trunk/iceiPaper/iceiPaper.tex	2004/09/17 14:56:05	1474
+++ trunk/iceiPaper/iceiPaper.tex	2005/01/06 20:04:03	1905
@@ -1,6 +1,5 @@
 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
 \documentclass[11pt]{article}
-%\documentclass[11pt]{article}
 \usepackage{endfloat}
 \usepackage{amsmath}
 \usepackage{epsf}
@@ -20,11 +19,12 @@
 
 \begin{document}
 
-\title{Ice-{\it i}: a simulation-predicted ice polymorph which is more
-stable than Ice $I_h$ for point-charge and point-dipole water models}
+\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\\
+Department of Chemistry and Biochemistry\\ 
+University of Notre Dame\\
 Notre Dame, Indiana 46556}
 
 \date{\today}
@@ -33,17 +33,16 @@ The free energies of several ice polymorphs in the low
 %\doublespacing
 
 \begin{abstract}
-The free energies of several ice polymorphs in the low pressure regime
-were calculated using thermodynamic integration.  These integrations
-were done for most of the common water models. Ice-{\it i}, a
-structure we recently observed to be stable in one of the single-point
-water models, was determined to be the stable crystalline state (at 1
-atm) for {\it all} the water models investigated.  Phase diagrams were
-generated, and phase coexistence lines were determined for all of the
-known low-pressure ice structures under all of the common water
-models.  Additionally, potential truncation was shown to have an
-effect on the calculated free energies, and can result in altered free
-energy landscapes.
+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 
@@ -54,148 +53,135 @@ Computer simulations are a valuable tool for studying 
 
 \section{Introduction}
 
-Computer simulations are a valuable tool for studying the phase
-behavior of systems ranging from small or simple molecules to complex
-biological species.\cite{Matsumoto02,Li96,Marrink01} Useful techniques
-have been developed to investigate the thermodynamic properites of
-model substances, providing both qualitative and quantitative
-comparisons between simulations and
-experiment.\cite{Widom63,Frenkel84} Investigation of these properties
-leads to the development of new and more accurate models, leading to
-better understanding and depiction of physical processes and intricate
-molecular systems.
-
 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{Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Berendsen98,Mahoney00,Fennell04}
+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, molecule transport, and the
+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 property to quantify are the Gibbs and Helmholtz free
-energies, particularly for the solid forms of water.  Difficulty in
-these types of studies typically arises from the assortment of
-possible crystalline polymorphs that water adopts over a wide range of
-pressures and temperatures.  There are currently 13 recognized forms
-of ice, and it is a challenging task to investigate the entire free
-energy landscape.\cite{Sanz04} Ideally, research is focused on the
+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 true transition temperatures and
-pressures for the model.
+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 the fact 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
+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 (Fig. \ref{iceiCell}A) consists of eight
-water molecules that stack in rows of interlocking water
-tetramers. Proton ordering can be accomplished by orienting two of the
-molecules so that both of their donated hydrogen bonds are internal to
-their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
-constructed of water tetramers, the hydrogen bonds are not as linear
-as those observed in ice $I_h$, however the interlocking of these
-subunits appears to provide significant stabilization to the overall
-crystal. The arrangement of these tetramers results in surrounding
-open octagonal cavities that are typically greater than 6.3 \AA\ in
-diameter. This relatively open overall structure leads to crystals
-that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
+simulation. The unit cell of Ice-{\it i} and an axially-elongated
+variant named Ice-{\it i}$^\prime$ both consist of eight water
+molecules that stack in rows of interlocking water tetramers as
+illustrated in figures \ref{unitcell}A and \ref{unitcell}B.  These
+tetramers form a crystal structure similar in appearance to a recent
+two-dimensional surface tessellation simulated on silica.\cite{Yang04}
+As expected in an ice crystal constructed of water tetramers, the
+hydrogen bonds are not as linear as those observed in ice $I_h$,
+however the interlocking of these subunits appears to provide
+significant stabilization to the overall crystal.  The arrangement of
+these tetramers results in surrounding open octagonal cavities that
+are typically greater than 6.3 \AA\ in diameter
+(Fig. \ref{iCrystal}).  This open structure leads to crystals that
+are typically 0.07 g/cm$^3$ less dense than ice $I_h$.
 
 \begin{figure}
+\centering
 \includegraphics[width=\linewidth]{unitCell.eps}
-\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
-elongated variant of Ice-{\it i}.  The spheres represent the
-center-of-mass locations of the water molecules.  The $a$ to $c$
-ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
-$a:2.1214c$ and $a:1.7850c$ respectively.}
-\label{iceiCell}
+\caption{(A) Unit cells for Ice-{\it i} and (B) Ice-{\it i}$^\prime$.  
+The spheres represent the center-of-mass locations of the water
+molecules.  The $a$ to $c$ ratios for Ice-{\it i} and Ice-{\it
+i}$^\prime$ are given by 2.1214 and 1.785 respectively.}
+\label{unitcell}
 \end{figure}
 
 \begin{figure}
+\centering
 \includegraphics[width=\linewidth]{orderedIcei.eps}
-\caption{Image of a proton ordered crystal of Ice-{\it i} looking 
-down the (001) crystal face. The rows of water tetramers surrounded by
-octagonal pores leads to a crystal structure that is significantly
-less dense than ice $I_h$.}
-\label{protOrder}
+\caption{A rendering of a proton ordered crystal of Ice-{\it i} looking
+down the (001) crystal face.  The presence of large octagonal pores
+leads to a polymorph that is less dense than ice $I_h$.}
+\label{iCrystal}
 \end{figure}
 
 Results from our previous study indicated that Ice-{\it i} is the
-minimum energy crystal structure for the single point water models we
+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, the
-absolute free energy of this crystal was calculated using
-thermodynamic integration and compared to the free energies of cubic
-and hexagonal ice $I$ (the experimental low density ice polymorphs)
-and ice B (a higher density, but very stable crystal structure
-observed by B\`{a}ez and Clancy in free energy studies of
-SPC/E).\cite{Baez95b} This work includes results for the water model
-from which Ice-{\it i} was crystallized (SSD/E) in addition to several
-common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
-field parametrized single point dipole water model (SSD/RF). It should
-be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used
-in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of
-this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit
-it is extended in the direction of the (001) face and compressed along
-the other two faces.
+articles).\cite{Fennell04,Liu96,Bratko85} Our earlier results
+considered only energetic stabilization and neglected entropic
+contributions to the overall free energy.  To address this issue, we
+have calculated the absolute free energy of this crystal using
+thermodynamic integration and compared it to the free energies of ice
+$I_c$ and ice $I_h$ (the common low density ice polymorphs) and ice B
+(a higher density, but very stable crystal structure observed by
+B\`{a}ez and Clancy in free energy studies of SPC/E).\cite{Baez95b}
+This work includes results for the water model from which Ice-{\it i}
+was crystallized (SSD/E) in addition to several common water models
+(TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field parametrized
+single point dipole water model (SSD/RF).  The axially-elongated
+variant, Ice-{\it i}$^\prime$, was used in calculations involving
+SPC/E, TIP4P, and TIP5P.  The square tetramers in Ice-{\it i} distort
+in Ice-{\it i}$^\prime$ to form a rhombus with alternating 85 and 95
+degree angles.  Under SPC/E, TIP4P, and TIP5P, this geometry is better
+at forming favorable hydrogen bonds.  The degree of rhomboid
+distortion depends on the water model used, but is significant enough
+to split a peak in the radial distribution function which corresponds
+to diagonal sites in the tetramers.
 
 \section{Methods}
 
 Canonical ensemble (NVT) molecular dynamics calculations were
-performed using the OOPSE molecular mechanics package.\cite{Meineke05}
+performed using the OOPSE molecular mechanics program.\cite{Meineke05}
 All molecules were treated as rigid bodies, with orientational motion
-propagated using the symplectic DLM integration method. Details about
+propagated using the symplectic DLM integration method.  Details about
 the implementation of this technique can be found in a recent
 publication.\cite{Dullweber1997}
 
-Thermodynamic integration is an established technique for
-determination of free energies of condensed phases of
-materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
-method, implemented in the same manner illustrated by B\`{a}ez and
-Clancy, was utilized to calculate the free energy of several ice
-crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
-SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
-and 400 K for all of these water models were also determined using
-this same technique in order to determine melting points and generate
-phase diagrams. All simulations were carried out at densities
-resulting in a pressure of approximately 1 atm at their respective
-temperatures.
-
-A single thermodynamic integration involves a sequence of simulations
-over 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
+Thermodynamic integration was utilized to calculate the Helmholtz free
+energies ($A$) of the listed water models at various state points
+using the OOPSE molecular dynamics program.\cite{Meineke05}
+Thermodynamic integration is an established technique that has been
+used extensively in the calculation of free energies for condensed
+phases of
+materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}.  This
+method uses a sequence of simulations during which the system of
+interest is converted into a reference system for which the free
+energy is known analytically ($A_0$).  The difference in potential
+energy between the reference system and the system of interest
+($\Delta V$) is then integrated in order to determine the free energy
 difference between the two states:
 \begin{equation}
-\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda 
-)}{\partial\lambda}\right\rangle_\lambda d\lambda,
+ A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
 \end{equation}
-where $V$ is the interaction potential and $\lambda$ is the
-transformation parameter that scales the overall
-potential. Simulations are distributed strategically along this path
-in order to sufficiently sample the regions of greatest change in the
-potential. Typical integrations in this study consisted of $\sim$25
-simulations ranging from 300 ps (for the unaltered system) to 75 ps
-(near the reference state) in length.
+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.  
 
-For the thermodynamic integration of molecular crystals, the Einstein
-crystal was chosen as the reference system. In an Einstein crystal,
-the molecules are restrained at their ideal lattice locations and
-orientations. Using harmonic restraints, as applied by B\`{a}ez and
-Clancy, the total potential for this reference crystal
-($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
+In an Einstein crystal, the molecules are restrained at their ideal
+lattice locations and orientations. Using harmonic restraints, as
+applied by B\`{a}ez and Clancy, the total potential for this reference
+crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
 \begin{equation}
 V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
 \frac{K_\omega\omega^2}{2},
@@ -203,9 +189,13 @@ respectively.  It is clear from Fig. \ref{waterSpring}
 where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
 the spring constants restraining translational motion and deflection
 of and rotation around the principle axis of the molecule
-respectively.  It is clear from Fig. \ref{waterSpring} that the values
-of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from
-$-\pi$ to $\pi$.  The partition function for a molecular crystal
+respectively.  These spring constants are typically calculated from
+the mean-square displacements of water molecules in an unrestrained
+ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
+mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
+17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
+the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
+from $-\pi$ to $\pi$.  The partition function for a molecular crystal
 restrained in this fashion can be evaluated analytically, and the
 Helmholtz Free Energy ({\it A}) is given by
 \begin{eqnarray}
@@ -223,11 +213,12 @@ potential energy of the ideal crystal.\cite{Baez95a}
 potential energy of the ideal crystal.\cite{Baez95a}
 
 \begin{figure}
-\includegraphics[width=\linewidth]{rotSpring.eps}
+\centering
+\includegraphics[width=4in]{rotSpring.eps}
 \caption{Possible orientational motions for a restrained molecule. 
 $\theta$ angles correspond to displacement from the body-frame {\it
 z}-axis, while $\omega$ angles correspond to rotation about the
-body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
+body-frame {\it z}-axis.  $K_\theta$ and $K_\omega$ are spring
 constants for the harmonic springs restraining motion in the $\theta$
 and $\omega$ directions.}
 \label{waterSpring}
@@ -239,283 +230,220 @@ molecules.  In this study, we apply of one of the most
 literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
 typically differ in regard to the path taken for switching off the
 interaction potential to convert the system to an ideal gas of water
-molecules.  In this study, we apply of one of the most convenient
-methods and integrate over the $\lambda^4$ path, where all interaction
-parameters are scaled equally by this transformation parameter.  This
-method has been shown to be reversible and provide results in
-excellent agreement with other established methods.\cite{Baez95b}
+molecules.  In this study, we applied one of the most convenient
+methods and integrated over the $\lambda^4$ path, where all
+interaction parameters are scaled equally by this transformation
+parameter.  This method has been shown to be reversible and provide
+results in excellent agreement with other established
+methods.\cite{Baez95b}
 
-Charge, dipole, and Lennard-Jones interactions were modified by a
-cubic switching between 100\% and 85\% of the cutoff value (9 \AA
-). By applying this function, these interactions are smoothly
+Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
+Lennard-Jones interactions were gradually reduced by a cubic switching
+function.  By applying this function, these interactions are smoothly
 truncated, thereby avoiding the poor energy conservation which results
-from harsher truncation schemes. The effect of a long-range correction
-was also investigated on select model systems in a variety of
-manners. For the SSD/RF model, a reaction field with a fixed
+from harsher truncation schemes.  The effect of a long-range
+correction was also investigated on select model systems in a variety
+of manners.  For the SSD/RF model, a reaction field with a fixed
 dielectric constant of 80 was applied in all
 simulations.\cite{Onsager36} For a series of the least computationally
-expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
-performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
-\AA\ cutoff results. Finally, results from the use of an Ewald
-summation were estimated for TIP3P and SPC/E by performing
-calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
-mechanics software package.\cite{Tinker} The calculated energy
-difference in the presence and absence of PME was applied to the
-previous results in order to predict changes to the free energy
-landscape.
+expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
+performed with longer cutoffs of 10.5, 12, 13.5, and 15 \AA\ to
+compare with the 9 \AA\ cutoff results.  Finally, the effects of using
+the Ewald summation were estimated for TIP3P and SPC/E by performing
+single configuration Particle-Mesh Ewald (PME)
+calculations~\cite{Tinker} for each of the ice polymorphs.  The
+calculated energy difference in the presence and absence of PME was
+applied to the previous results in order to predict changes to the
+free energy landscape.
 
-\section{Results and discussion}
+\section{Results and Discussion}
 
-The free energy of proton ordered Ice-{\it i} was calculated and
-compared with the free energies of proton ordered variants of the
-experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
-as well as the higher density ice B, observed by B\`{a}ez and Clancy
-and thought to be the minimum free energy structure for the SPC/E
-model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
-Ice XI, the experimentally-observed proton-ordered variant of ice
-$I_h$, was investigated initially, but was found to be not as stable
-as proton disordered or antiferroelectric variants of ice $I_h$. The
-proton ordered variant of ice $I_h$ used here is a simple
-antiferroelectric version that we divised, and it has an 8 molecule
-unit cell similar to other predicted antiferroelectric $I_h$
-crystals.\cite{Davidson84} The crystals contained 648 or 1728
-molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
-molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
-crystal sizes were necessary for simulations involving larger cutoff
-values.
+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}
-\renewcommand{\thefootnote}{\thempfootnote}
 \begin{center}
-\caption{Calculated free energies for several ice polymorphs with a 
-variety of common water models. All calculations used a cutoff radius
-of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
-kcal/mol. Calculated error of the final digits is in parentheses. *Ice
-$I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
-\begin{tabular}{ l  c  c  c  c }
+\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}\\
+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)\\
-TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
-TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
-SPC/E & -12.67(2) & -12.96(2) & -13.25(3) & -13.55(2)\\
-SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
-SSD/RF & -11.51(2) & NA* & -12.08(3) & -12.29(2)\\
+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*}
 
-The free energy values computed for the studied polymorphs indicate
-that Ice-{\it i} is the most stable state for all of the common water
-models studied. With the free energy at these state points, the
-Gibbs-Helmholtz equation was used to project to other state points and
-to build phase diagrams.  Figures
-\ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
-from the free energy results. All other models have similar structure,
-although the crossing points between the phases exist at slightly
-different temperatures and pressures. It is interesting to note that
-ice $I$ does not exist in either cubic or hexagonal form in any of the
-phase diagrams for any of the models. For purposes of this study, ice
-B is representative of the dense ice polymorphs. A recent study by
-Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
-TIP4P in the high pressure regime.\cite{Sanz04}
-
 \begin{figure}
 \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
 \caption{Phase diagram for the TIP3P water model in the low pressure 
-regime. The displayed $T_m$ and $T_b$ values are good predictions of
+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
+experimentally observed forms.  Both cubic and hexagonal ice $I$ are
 higher in energy and don't appear in the phase diagram.}
-\label{tp3phasedia}
+\label{tp3PhaseDia}
 \end{figure}
 
-\begin{figure}
-\includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
-\caption{Phase diagram for the SSD/RF water model in the low pressure 
-regime. Calculations producing these results were done under an
-applied reaction field. It is interesting to note that this
-computationally efficient model (over 3 times more efficient than
-TIP3P) exhibits phase behavior similar to the less computationally
-conservative charge based models.}
-\label{ssdrfphasedia}
-\end{figure}
-
-\begin{table*}
-\begin{minipage}{\linewidth}
-\renewcommand{\thefootnote}{\thempfootnote}
-\begin{center}
-\caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) 
-temperatures at 1 atm for several common water models compared with
-experiment. The $T_m$ and $T_s$ values from simulation correspond to a
-transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
-liquid or gas state.}
-\begin{tabular}{ l  c  c  c  c  c  c  c }
-\hline 
-Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
-\hline 
-$T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
-$T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
-$T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
-\end{tabular}
-\label{meltandboil}
-\end{center}
-\end{minipage}
-\end{table*}
-
-Table \ref{meltandboil} lists the melting and boiling temperatures
-calculated from this work. Surprisingly, most of these models have
-melting points that compare quite favorably with experiment. The
-unfortunate aspect of this result is that this phase change occurs
-between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
-liquid state. These results are actually not contrary to previous
-studies in the literature. Earlier free energy studies of ice $I$
-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
+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 at 265 K, significantly higher in temperature than the
-previous studies. Also of interest in these results is that SSD/E does
-not exhibit a melting point at 1 atm, but it shows a sublimation point
-at 355 K. This is due to the significant stability of Ice-{\it i} over
-all other polymorphs for this particular model under these
-conditions. While troubling, this behavior turned out to be
-advantageous in that it facilitated the spontaneous crystallization of
-Ice-{\it i}. These observations provide a warning that simulations of
-SSD/E as a ``liquid'' near 300 K are actually metastable and run the
-risk of spontaneous crystallization. However, this risk changes when
-applying a longer cutoff.
+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 (A) SSD/E, (B) 
-TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
-\AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
-\AA . These crystals are unstable at 200 K and rapidly convert into
-liquids. The connecting lines are qualitative visual aid.}
+\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 all the ice polymorphs show a substantial dependence on
-cutoff radius. In general, there is a narrowing of the free energy
-differences while moving to greater cutoff radius. Interestingly, by
-increasing the cutoff radius, the free energy gap was narrowed enough
-in the SSD/E model that the liquid state is preferred under standard
-simulation conditions (298 K and 1 atm). Thus, it is recommended that
-simulations using this model choose interaction truncation radii
-greater than 9 \AA\ . This narrowing trend is much more subtle in the
-case of SSD/RF, indicating that the free energies calculated with a
-reaction field present provide a more accurate picture of the free
-energy landscape in the absence of potential truncation.
+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.
 
-To further study the changes resulting to the inclusion of a
-long-range interaction correction, 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 with and without PME using TINKER. The
-free energies for the investigated polymorphs using the TIP3P and
-SPC/E water models are shown in Table \ref{pmeShift}. The same trend
-pointed out through increase of cutoff radius is observed in these PME
-results. Ice-{\it i} is the preferred polymorph at ambient conditions
-for both the TIP3P and SPC/E water models; however, the narrowing of
-the free energy differences between the various solid forms is
-significant enough that it becomes less clear that it is the most
-stable polymorph with the SPC/E model.  The free energies of Ice-{\it
-i} and ice B nearly overlap within error, with ice $I_c$ just outside
-as well, indicating that Ice-{\it i} might be metastable with respect
-to ice B and possibly ice $I_c$ with SPC/E. However, these results do
-not significantly alter 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.
+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.
 
-\begin{table*}
-\begin{minipage}{\linewidth}
-\renewcommand{\thefootnote}{\thempfootnote}
-\begin{center}
-\caption{The free energy of the studied ice polymorphs after applying 
-the energy difference attributed to the inclusion of the PME
-long-range interaction correction. Units are kcal/mol.}
-\begin{tabular}{ l  c  c  c  c }
-\hline 
-\ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
-\hline 
-TIP3P  & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3)\\
-SPC/E  & -12.77(2) & -12.92(2) & -12.96(3) & -13.02(2)\\
-\end{tabular}
-\label{pmeShift}
-\end{center}
-\end{minipage}
-\end{table*}
+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. 
 
-\section{Conclusions}
-
-The free energy for proton ordered variants of hexagonal and cubic ice
-$I$, ice B, and recently discovered Ice-{\it i} were calculated under
-standard conditions for several common water models via thermodynamic
-integration. All the water models studied show Ice-{\it i} to be the
-minimum free energy crystal structure in the with a 9 \AA\ switching
-function cutoff. Calculated melting and boiling points show
-surprisingly good agreement with the experimental values; however, the
-solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
-interaction truncation was investigated through variation of the
-cutoff radius, use of a reaction field parameterized model, and
-estimation of the results in the presence of the Ewald
-summation. Interaction truncation has a significant effect on the
-computed free energy values, and may significantly alter the free
-energy landscape for the more complex multipoint water models. Despite
-these effects, these results show Ice-{\it i} to be an important ice
-polymorph that should be considered in simulation studies.
-
-Due to this relative stability of Ice-{\it i} in all manner of
-investigated simulation examples, 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
+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.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
-our predictions for both the pair distribution function ($g_{OO}(r)$)
-and the structure factor ($S(\vec{q})$ for ice $I_c$ and for ice-{\it
-i} at a temperature of 77K.  In a quick comparison of the predicted
-S(q) for Ice-{\it i} and experimental studies of amorphous solid
-water, it is possible that some of the ``spurious'' peaks that could
-not be assigned in HDA could correspond to peaks labeled in this
-S(q).\cite{Bizid87} It should be noted that there is typically poor
-agreement on crystal densities between simulation and experiment, so
-such peak comparisons should be made with caution.  We will leave it
-to our experimental colleagues to determine whether this ice polymorph
-is named appropriately or if it should be promoted to Ice-0.
+non-polar molecules.  For experimental comparison purposes, example
+$g_{OO}(r)$ and $S(\vec{q})$ plots were generated for the two Ice-{\it
+i} variants (along with example ice $I_h$ and $I_c$ plots) at 77K, and
+they are shown in figures \ref{fig:gofr} and \ref{fig:sofq}
+respectively.
 
 \begin{figure}
+\centering
 \includegraphics[width=\linewidth]{iceGofr.eps}
-\caption{Radial distribution functions of Ice-{\it i} and ice $I_c$ 
-calculated from from simulations of the SSD/RF water model at 77 K.}
+\caption{Radial distribution functions of ice $I_h$, $I_c$, and
+Ice-{\it i} calculated from from simulations of the SSD/RF water model
+at 77 K.  The Ice-{\it i} distribution function was obtained from
+simulations composed of TIP4P water.}
 \label{fig:gofr}
 \end{figure}
 
 \begin{figure}
+\centering
 \includegraphics[width=\linewidth]{sofq.eps}
-\caption{Predicted structure factors for Ice-{\it i} and ice $I_c$ at
-77 K.  The raw structure factors have been convoluted with a gaussian
-instrument function (0.075 \AA$^{-1}$ width) to compensate for the
-trunction effects in our finite size simulations. The labeled peaks
-compared favorably with ``spurious'' peaks observed in experimental
-studies of amorphous solid water.\cite{Bizid87}}
+\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
+ and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
+ been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
+ width) to compensate for the trunction effects in our finite size
+ simulations.}
 \label{fig:sofq}
 \end{figure}