--- trunk/iceiPaper/iceiPaper.tex 2004/09/21 20:17:12 1488 +++ trunk/iceiPaper/iceiPaper.tex 2004/12/02 14:23:11 1831 @@ -20,8 +20,8 @@ \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{Free Energy Analysis of Simulated Ice Polymorphs Using Simple +Dipolar and Charge Based Water Models} \author{Christopher J. Fennell and J. Daniel Gezelter \\ Department of Chemistry and Biochemistry\\ University of Notre Dame\\ @@ -84,23 +84,25 @@ constant pressure and temperature conditions. Crystall 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 +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 +rows of interlocking water tetramers. This crystal structure has a +limited resemblence to a recent two-dimensional ice tessellation +simulated on a silica surface.\cite{Yang04} 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 +(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 +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$. \begin{figure} @@ -116,8 +118,8 @@ down the (001) crystal face. The rows of water tetrame \begin{figure} \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 +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} \end{figure} @@ -128,7 +130,7 @@ contributions to the overall free energy. To address t 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 +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) @@ -139,7 +141,7 @@ was used in calculations involving SPC/E, TIP4P, and T 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 +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 @@ -154,7 +156,7 @@ propagated using the symplectic DLM integration method Canonical ensemble (NVT) molecular dynamics calculations were 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 +propagated using the symplectic DLM integration method. Details about the implementation of this technique can be found in a recent publication.\cite{Dullweber1997} @@ -167,13 +169,13 @@ generate phase diagrams. All simulations were carried 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. +generate phase diagrams. All simulations were carried out at +densities which correspond to a pressure of approximately 1 atm at +their respective temperatures. 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 +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} @@ -181,15 +183,15 @@ transformation parameter that scales the overall )}{\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 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. +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. For the thermodynamic integration of molecular crystals, the Einstein -crystal was chosen as the reference system. In an Einstein crystal, +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 @@ -201,9 +203,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} @@ -225,7 +231,7 @@ body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are \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} @@ -245,23 +251,22 @@ cubic switching between 100\% and 85\% of the cutoff v 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 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 -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, 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. +cubic switching between 100\% and 85\% of the cutoff value (9 \AA). +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 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, 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} @@ -273,35 +278,35 @@ as proton disordered or antiferroelectric variants of 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 +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. +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} -\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 } +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.} + +\begin{tabular}{lcccc} \hline Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\ \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)\\ +SPC/E & -12.87(2) & -13.05(2) & -13.26(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)\\ +SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\ \end{tabular} \label{freeEnergy} \end{center} @@ -310,25 +315,25 @@ models studied. With the calculated free energy at the 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 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, -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 +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, 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} \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} \end{figure} @@ -336,8 +341,8 @@ regime. Calculations producing these results were done \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 +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.} @@ -346,16 +351,17 @@ conservative charge based models.} \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 +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 } + +\begin{tabular}{lccccccc} \hline -Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\ +Equilibrium 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\\ @@ -367,92 +373,116 @@ calculated from this work. Surprisingly, most of these \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 +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 +liquid state. These results are actually not contrary to other +studies. Studies of ice $I_h$ using TIP4P predict a $T_m$ ranging +from 214 to 238 K (differences being attributed to choice of +interaction truncation and different ordered and disordered molecular arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be -predicted from this work. However, the $T_m$ from Ice-{\it i} is -calculated at 265 K, significantly higher in temperature than the -previous studies. Also of interest in these results is that SSD/E does +predicted from this work. However, the $T_m$ from Ice-{\it i} is +calculated to be 265 K, indicating that these simulation based +structures ought to be included in studies probing phase transitions +with this model. Also of interest in these results is that SSD/E does not exhibit a melting point at 1 atm, but it shows a sublimation point -at 355 K. This is due to the significant stability of Ice-{\it i} over -all other polymorphs for this particular model under these -conditions. While troubling, this behavior resulted in spontaneous +at 355 K. This is due to the significant stability of Ice-{\it i} +over all other polymorphs for this particular model under these +conditions. While troubling, this behavior resulted in spontaneous crystallization of Ice-{\it i} and led us to investigate this -structure. These observations provide a warning that simulations of +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 lessens when +risk of spontaneous crystallization. However, this risk lessens when applying a longer cutoff. \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. Calculations performed without a +long-range correction show noticable free energy dependence on the +cutoff radius and show some degree of converge at large cutoff radii. +Inclusion of a long-range correction reduces the cutoff radius +dependence of the free energy for all the models. 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 +involve potential truncation. As seen in Fig. \ref{incCutoff}, the +free energy of the ice polymorphs with water models lacking a +long-range correction show a significant cutoff radius dependence. In +general, there is a narrowing of the free energy differences while +moving to greater cutoff radii. As the free energies for the +polymorphs converge, the stability advantage that Ice-{\it i} exhibits +is reduced. Interestingly, increasing the cutoff radius a mere 1.5 +\AA\ with the SSD/E model destabilizes the Ice-{\it i} polymorph +enough 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. +greater than 9 \AA. Considering the stabilization of Ice-{\it i} with +smaller cutoffs, it is not surprising that crystallization was +observed with SSD/E. The choice of a 9 \AA\ cutoff in the previous +simulations gives the Ice-{\it i} polymorph a greater than 1 kcal/mol +lower free energy than the ice $I_\textrm{h}$ starting configurations. +Additionally, it should be noted that ice $I_c$ is not stable with +cutoff radii of 12 and 13.5 \AA\ using the TIP3P water model. These +simulations showed bulk distortions of the simulation cell that +rapidly deteriorated crystal integrity. -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 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. +Adjacent to each of these model plots is a system with an applied or +estimated long-range correction. SSD/RF was parametrized for use with +a reaction field, and the benefit provided by this computationally +inexpensive correction is apparent. Due to the relative independence +of the resultant free energies, calculations performed with a small +cutoff radius provide resultant properties similar to what one would +expect for the bulk material. In the cases of TIP3P and SPC/E, the +effect of an Ewald summation was estimated by applying the potential +energy difference do to its inclusion in systems in the presence and +absence of the correction. This was accomplished by calculation of +the potential energy of identical crystals both with and without +particle mesh Ewald (PME). Similar behavior to that observed with +reaction field is seen for both of these models. The free energies +show less dependence on cutoff radius and span a more narrowed range +for the various polymorphs. Like the dipolar water models, TIP3P +displays a relatively constant preference for the Ice-{\it i} +polymorph. Crystal preference is much more difficult to determine for +SPC/E. Without a long-range correction, each of the polymorphs +studied assumes the role of the preferred polymorph under different +cutoff conditions. The inclusion of the Ewald correction flattens and +narrows the sequences of free energies so much that they often overlap +within error (see Table \ref{spcecut}), indicating that other +conditions, such as cell volume in microcanonical simulations, can +influence the chosen polymorph upon crystallization. All of these +results support the finding that the Ice-{\it i} polymorph is a stable +crystal structure that should be considered when studying the phase +behavior of water models. \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 } + +\caption{The free energy versus cutoff radius for the studied ice +polymorphs using SPC/E after the inclusion of the PME long-range +interaction correction. Units are kcal/mol.} + +\begin{tabular}{ccccc} \hline -\ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\ +Cutoff (\AA) & $I_h$ & $I_c$ & 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)\\ +9.0 & -12.98(2) & -13.00(2) & -12.97(3) & -13.02(2) \\ +10.5 & -13.13(3) & -13.09(4) & -13.17(3) & -13.11(2) \\ +12.0 & -13.06(2) & -13.09(2) & -13.15(4) & -13.12(2) \\ +13.5 & -13.02(2) & -13.02(2) & -13.08(2) & -13.07(2) \\ +15.0 & -13.11(4) & -12.97(2) & -13.09(2) & -12.95(2) \\ \end{tabular} -\label{pmeShift} +\label{spcecut} \end{center} \end{minipage} \end{table*} @@ -462,19 +492,19 @@ via thermodynamic integration. All the water models st The free energy for proton ordered variants of hexagonal and cubic ice $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 +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 +\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 +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. +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 of the investigated simulation conditions, the question arises as to possible