--- trunk/iceiPaper/iceiPaper.tex	2004/09/15 21:11:19	1463
+++ trunk/iceiPaper/iceiPaper.tex	2004/09/21 20:17:12	1488
@@ -20,7 +20,8 @@
 
 \begin{document}
 
-\title{Ice-{\it i}: a novel ice polymorph predicted via computer simulation}
+\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}
 
 \author{Christopher J. Fennell and J. Daniel Gezelter \\
 Department of Chemistry and Biochemistry\\ University of Notre Dame\\
@@ -32,17 +33,21 @@ 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
+The absolute free energies of several ice polymorphs which are stable
+at low pressures were calculated using thermodynamic integration to a
+reference system (the Einstein crystal).  These integrations were
+performed for most of the common water models (SPC/E, TIP3P, TIP4P,
+TIP5P, SSD/E, SSD/RF). Ice-{\it i}, a structure we recently observed
+crystallizing at room temperature for 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.
+known low-pressure ice structures under all of these water models.
+Additionally, potential truncation was shown to have an effect on the
+calculated free energies, and can result in altered free energy
+landscapes.  Structure factor predictions for the new crystal were
+generated and we await experimental confirmation of the existence of
+this new polymorph.
 \end{abstract}
 
 %\narrowtext 
@@ -53,55 +58,46 @@ Molecular dynamics is a valuable tool for studying the
 
 \section{Introduction}
 
-Molecular dynamics is a valuable tool for studying the phase behavior
-of systems ranging from small or simple
-molecules\cite{Matsumoto02andOthers} to complex biological
-species.\cite{bigStuff} Many techniques have been developed to
-investigate the thermodynamic properites of model substances,
-providing both qualitative and quantitative comparisons between
-simulations and experiment.\cite{thermMethods} 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{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,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 and the hydrophobic
-effect.\cite{evenMorePapers} 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{modelProps} 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 phases having the lowest free energy at a given state
-point, because these phases will dictate the true transition
-temperatures and pressures for their respective model.
+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.  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
+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 the
-study of 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
-obtained in this case was different from any previously observed ice
-polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
-was termed Ice-{\it i} in homage to 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
-waters 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
+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
+(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
@@ -109,7 +105,11 @@ that are 0.07 g/cm$^3$ less dense on average than ice 
 
 \begin{figure}
 \includegraphics[width=\linewidth]{unitCell.eps}
-\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-2{\it i}, the elongated variant of Ice-{\it i}.  For Ice-{\it i}, the $a$ to $c$ relation is given by $a = 2.1214c$, while for Ice-2{\it i}, $a = 1.7850c$.}
+\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it 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}
 \end{figure}
 
@@ -122,70 +122,90 @@ Results in the previous study indicated that Ice-{\it 
 \label{protOrder}
 \end{figure}
 
-Results in the previous study indicated that Ice-{\it i} is the
-minimum energy crystal structure for the single point water models
-being studied (for discussions on these single point dipole models,
-see the previous work and related
-articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
-consider energetic stabilization and neglect 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
-(soft sticky dipole extended, 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 (soft sticky dipole
-reaction field, SSD/RF). In should be noted that a second version of
-Ice-{\it i} (Ice-2{\it i}) 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.
+Results from our previous study indicated that Ice-{\it i} is the
+minimum energy crystal structure for the single point water models we
+had 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). It should
+be noted that a second version of Ice-{\it i} (Ice-{\it 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.  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.
 
 \section{Methods}
 
 Canonical ensemble (NVT) molecular dynamics calculations were
-performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
-molecular mechanics package. All molecules were treated as rigid
-bodies, with orientational motion propagated using the symplectic DLM
-integration method. Details about the implementation of these
-techniques can be found in a recent publication.\cite{Meineke05} 
+performed using the OOPSE molecular mechanics package.\cite{Meineke05}
+All molecules were treated as rigid bodies, with orientational motion
+propagated using the symplectic DLM integration method. Details about
+the implementation of this technique can be found in a recent
+publication.\cite{Dullweber1997}
 
-Thermodynamic integration was utilized to calculate the free energy of
-several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
-SSD/RF, and SSD/E water models. 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.
+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 to
+generate phase diagrams. All simulations were carried out at densities
+which correspond to 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. This transformation path is then
-integrated in order to determine the free energy difference between
-the two states:
+Thermodynamic integration 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. Simulations are distributed unevenly along this path in
-order to sufficiently sample the regions of greatest change in the
+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.
 
 For the thermodynamic integration of molecular crystals, the Einstein
-Crystal is chosen as the reference state that the system is converted
-to over the course of the simulation. In an Einstein Crystal, the
-molecules are harmonically restrained at their ideal lattice locations
-and orientations. The partition function for a molecular crystal
-restrained in this fashion has been evaluated, and the Helmholtz Free
-Energy ({\it A}) is given by
+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,
+\begin{equation}
+V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
+\frac{K_\omega\omega^2}{2},
+\end{equation}
+where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
+the spring constants restraining translational motion and deflection
+of and rotation around the principle axis of the molecule
+respectively.  It is clear from Fig. \ref{waterSpring} that the values
+of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from
+$-\pi$ to $\pi$.  The partition function for a molecular crystal
+restrained in this fashion can be evaluated analytically, and the
+Helmholtz Free Energy ({\it A}) is given by
 \begin{eqnarray}
 A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
 [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
@@ -197,14 +217,8 @@ where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a
 )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
 \label{ecFreeEnergy}
 \end{eqnarray}
-where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
-\ref{ecFreeEnergy}, $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 (Fig. \ref{waterSpring}), and $E_m$ is the
-minimum potential energy of the ideal crystal. In the case of
-molecular liquids, the ideal vapor is chosen as the target reference
-state.
+where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
+potential energy of the ideal crystal.\cite{Baez95a}
 
 \begin{figure}
 \includegraphics[width=\linewidth]{rotSpring.eps}
@@ -217,10 +231,23 @@ Charge, dipole, and Lennard-Jones interactions were mo
 \label{waterSpring}
 \end{figure}
 
+In the case of molecular liquids, the ideal vapor is chosen as the
+target reference state.  There are several examples of liquid state
+free energy calculations of water models present in the
+literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
+typically differ in regard to the path taken for switching off the
+interaction potential to convert the system to an ideal gas of water
+molecules.  In this study, we applied of 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
-truncated, thereby avoiding poor energy conserving dynamics resulting
+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
@@ -228,33 +255,33 @@ performed with longer cutoffs of 12 and 15 \AA\ to com
 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} TINKER was chosen because it
-can also propagate the motion of rigid-bodies, and provides the most
-direct comparison to the results from OOPSE. The calculated energy
-difference in the presence and absence of PME was applied to the
-previous results in order to predict changes in the free energy
-landscape.
+\AA\ cutoff results. Finally, the effects of utilizing an Ewald
+summation were estimated for TIP3P and SPC/E by performing single
+configuration 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.
 
 \section{Results and discussion}
 
-The free energy of proton ordered Ice-{\it i} was calculated and
+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 it was found not to be as
-stable as antiferroelectric variants of proton ordered or even proton
-disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
-ice $I_h$ used here is a simple antiferroelectric version that has an
-8 molecule unit cell. 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.
+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 devised, 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.
 
 \begin{table*}
 \begin{minipage}{\linewidth}
@@ -263,17 +290,18 @@ kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF
 \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. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
+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 }
 \hline 
-\ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
+Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
 \hline 
-\ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
-\ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
-\ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
-\ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
-\ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
-\ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
+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)\\
 \end{tabular}
 \label{freeEnergy}
 \end{center}
@@ -282,18 +310,18 @@ models studied. With the free energy at these state po
 
 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
-temperature and pressure dependence of the free energy was used to
-project to other state points and build phase diagrams. Figures
+models studied. With the calculated 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,
-only the crossing points between these phases exist at 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} 
+although the crossing points between the phases move to 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 at higher pressures than those studied here.\cite{Sanz04}
 
 \begin{figure}
 \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
@@ -321,14 +349,17 @@ temperatures of several common water models compared w
 \renewcommand{\thefootnote}{\thempfootnote}
 \begin{center}
 \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) 
-temperatures of several common water models compared with experiment.}
+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 \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
+Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
 \hline 
-\ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
-\ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
-\ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
+$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}
@@ -344,19 +375,20 @@ ordered and disordered molecular arrangements). If the
 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 arrangements). If the presence of ice
-B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
+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
+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, this risk changes when
+risk of spontaneous crystallization. However, this risk lessens when
 applying a longer cutoff.
 
 \begin{figure}
@@ -364,8 +396,8 @@ TIP3P, and (C) SSD/RF. Data points omitted include SSD
 \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 a
-liquid. The connecting lines are qualitative visual aid.}
+\AA . These crystals are unstable at 200 K and rapidly convert into
+liquids. The connecting lines are qualitative visual aid.}
 \label{incCutoff}
 \end{figure}
 
@@ -380,7 +412,7 @@ greater than 9 \AA\. This narrowing trend is much more
 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
+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.
@@ -390,21 +422,21 @@ energy of identical crystals with and without PME usin
 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}. TIP4P and TIP5P
-are not fully supported in TINKER, so the results for these models
-could not be estimated. 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, there is a narrowing of the free energy
-differences between the various solid forms. In the case of SPC/E this
-narrowing is significant enough that it becomes less clear cut that
-Ice-{\it i} is the most stable polymorph, and is possibly metastable
-with respect to ice B and possibly ice $I_c$. 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.
+energy of identical crystals both with and without PME. 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.
 
 \begin{table*}
 \begin{minipage}{\linewidth}
@@ -417,8 +449,8 @@ long-range interaction correction. Units are kcal/mol.
 \hline 
 \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
 \hline 
-\ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
-\ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
+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}
@@ -428,34 +460,62 @@ $I$, ice B, and recently discovered Ice-{\it i} where 
 \section{Conclusions}
 
 The free energy for proton ordered variants of hexagonal and cubic ice
-$I$, ice B, and recently discovered Ice-{\it i} where 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
-correction. Interaction truncation has a significant effect on the
+$I$, ice B, and our recently discovered Ice-{\it i} structure 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 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
+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 leads the authors to be hesitant in ascribing relevance outside
-of computational models, hence the descriptive name presented. That
-being 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.
+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.  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_h$, $I_c$, and for
+ice-{\it i} at a temperature of 77K.  In studies of the high and low
+density forms of amorphous ice, ``spurious'' diffraction peaks have
+been observed experimentally.\cite{Bizid87} It is possible that a
+variant of Ice-{\it i} could explain some of this behavior; however,
+we will leave it to our experimental colleagues to make the final
+determination on whether this ice polymorph is named appropriately
+(i.e. with an imaginary number) or if it can be promoted to Ice-0.
 
+\begin{figure}
+\includegraphics[width=\linewidth]{iceGofr.eps}
+\caption{Radial distribution functions of ice $I_h$, $I_c$, 
+Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
+of the SSD/RF water model at 77 K.}
+\label{fig:gofr}
+\end{figure}
+
+\begin{figure}
+\includegraphics[width=\linewidth]{sofq.eps}
+\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
+ and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
+ been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
+ width) to compensate for the trunction effects in our finite size
+ simulations.}
+\label{fig:sofq}
+\end{figure}
+
 \section{Acknowledgments}
 Support for this project was provided by the National Science
 Foundation under grant CHE-0134881. Computation time was provided by
