ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceiPaper/iceiPaper.tex
(Generate patch)

Comparing trunk/iceiPaper/iceiPaper.tex (file contents):
Revision 1453 by chrisfen, Mon Sep 13 18:39:53 2004 UTC vs.
Revision 1463 by gezelter, Wed Sep 15 21:11:19 2004 UTC

# Line 1 | Line 1
1   %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 < \documentclass[preprint,aps,endfloats]{revtex4}
2 > \documentclass[11pt]{article}
3   %\documentclass[11pt]{article}
4 < %\usepackage{endfloat}
4 > \usepackage{endfloat}
5   \usepackage{amsmath}
6   \usepackage{epsf}
7   \usepackage{berkeley}
8 < %\usepackage{setspace}
9 < %\usepackage{tabularx}
8 > \usepackage{setspace}
9 > \usepackage{tabularx}
10   \usepackage{graphicx}
11 < %\usepackage[ref]{overcite}
12 < %\pagestyle{plain}
13 < %\pagenumbering{arabic}
14 < %\oddsidemargin 0.0cm \evensidemargin 0.0cm
15 < %\topmargin -21pt \headsep 10pt
16 < %\textheight 9.0in \textwidth 6.5in
17 < %\brokenpenalty=10000
11 > \usepackage[ref]{overcite}
12 > \pagestyle{plain}
13 > \pagenumbering{arabic}
14 > \oddsidemargin 0.0cm \evensidemargin 0.0cm
15 > \topmargin -21pt \headsep 10pt
16 > \textheight 9.0in \textwidth 6.5in
17 > \brokenpenalty=10000
18 > \renewcommand{\baselinestretch}{1.2}
19 > \renewcommand\citemid{\ } % no comma in optional reference note
20  
19 – %\renewcommand\citemid{\ } % no comma in optional reference note
20 –
21   \begin{document}
22  
23 < \title{A Free Energy Study of Low Temperature and Anomolous Ice}
23 > \title{Ice-{\it i}: a novel ice polymorph predicted via computer simulation}
24  
25 < \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
26 < \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
27 <
28 < \address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
25 > \author{Christopher J. Fennell and J. Daniel Gezelter \\
26 > Department of Chemistry and Biochemistry\\ University of Notre Dame\\
27   Notre Dame, Indiana 46556}
28  
29   \date{\today}
30  
31 < %\maketitle
31 > \maketitle
32   %\doublespacing
33  
34   \begin{abstract}
35 + The free energies of several ice polymorphs in the low pressure regime
36 + were calculated using thermodynamic integration.  These integrations
37 + were done for most of the common water models. Ice-{\it i}, a
38 + structure we recently observed to be stable in one of the single-point
39 + water models, was determined to be the stable crystalline state (at 1
40 + atm) for {\it all} the water models investigated.  Phase diagrams were
41 + generated, and phase coexistence lines were determined for all of the
42 + known low-pressure ice structures under all of the common water
43 + models.  Additionally, potential truncation was shown to have an
44 + effect on the calculated free energies, and can result in altered free
45 + energy landscapes.
46   \end{abstract}
47  
39 – \maketitle
40 –
41 – \newpage
42 –
48   %\narrowtext
49  
50   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 48 | Line 53 | Notre Dame, Indiana 46556}
53  
54   \section{Introduction}
55  
56 + Molecular dynamics is a valuable tool for studying the phase behavior
57 + of systems ranging from small or simple
58 + molecules\cite{Matsumoto02andOthers} to complex biological
59 + species.\cite{bigStuff} Many techniques have been developed to
60 + investigate the thermodynamic properites of model substances,
61 + providing both qualitative and quantitative comparisons between
62 + simulations and experiment.\cite{thermMethods} Investigation of these
63 + properties leads to the development of new and more accurate models,
64 + leading to better understanding and depiction of physical processes
65 + and intricate molecular systems.
66 +
67 + Water has proven to be a challenging substance to depict in
68 + simulations, and a variety of models have been developed to describe
69 + its behavior under varying simulation
70 + conditions.\cite{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
71 + These models have been used to investigate important physical
72 + phenomena like phase transitions and the hydrophobic
73 + effect.\cite{evenMorePapers} With the choice of models available, it
74 + is only natural to compare the models under interesting thermodynamic
75 + conditions in an attempt to clarify the limitations of each of the
76 + models.\cite{modelProps} Two important property to quantify are the
77 + Gibbs and Helmholtz free energies, particularly for the solid forms of
78 + water.  Difficulty in these types of studies typically arises from the
79 + assortment of possible crystalline polymorphs that water adopts over a
80 + wide range of pressures and temperatures. There are currently 13
81 + recognized forms of ice, and it is a challenging task to investigate
82 + the entire free energy landscape.\cite{Sanz04} Ideally, research is
83 + focused on the phases having the lowest free energy at a given state
84 + point, because these phases will dictate the true transition
85 + temperatures and pressures for their respective model.
86 +
87 + In this paper, standard reference state methods were applied to the
88 + study of crystalline water polymorphs in the low pressure regime. This
89 + work is unique in the fact that one of the crystal lattices was
90 + arrived at through crystallization of a computationally efficient
91 + water model under constant pressure and temperature
92 + conditions. Crystallization events are interesting in and of
93 + themselves\cite{Matsumoto02,Yamada02}; however, the crystal structure
94 + obtained in this case was different from any previously observed ice
95 + polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
96 + was termed Ice-{\it i} in homage to its origin in computational
97 + simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
98 + water molecules that stack in rows of interlocking water
99 + tetramers. Proton ordering can be accomplished by orienting two of the
100 + waters so that both of their donated hydrogen bonds are internal to
101 + their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
102 + constructed of water tetramers, the hydrogen bonds are not as linear
103 + as those observed in ice $I_h$, however the interlocking of these
104 + subunits appears to provide significant stabilization to the overall
105 + crystal. The arrangement of these tetramers results in surrounding
106 + open octagonal cavities that are typically greater than 6.3 \AA\ in
107 + diameter. This relatively open overall structure leads to crystals
108 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
109 +
110 + \begin{figure}
111 + \includegraphics[width=\linewidth]{unitCell.eps}
112 + \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$.}
113 + \label{iceiCell}
114 + \end{figure}
115 +
116 + \begin{figure}
117 + \includegraphics[width=\linewidth]{orderedIcei.eps}
118 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
119 + down the (001) crystal face. The rows of water tetramers surrounded by
120 + octagonal pores leads to a crystal structure that is significantly
121 + less dense than ice $I_h$.}
122 + \label{protOrder}
123 + \end{figure}
124 +
125 + Results in the previous study indicated that Ice-{\it i} is the
126 + minimum energy crystal structure for the single point water models
127 + being studied (for discussions on these single point dipole models,
128 + see the previous work and related
129 + articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
130 + consider energetic stabilization and neglect entropic contributions to
131 + the overall free energy. To address this issue, the absolute free
132 + energy of this crystal was calculated using thermodynamic integration
133 + and compared to the free energies of cubic and hexagonal ice $I$ (the
134 + experimental low density ice polymorphs) and ice B (a higher density,
135 + but very stable crystal structure observed by B\`{a}ez and Clancy in
136 + free energy studies of SPC/E).\cite{Baez95b} This work includes
137 + results for the water model from which Ice-{\it i} was crystallized
138 + (soft sticky dipole extended, SSD/E) in addition to several common
139 + water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field
140 + parametrized single point dipole water model (soft sticky dipole
141 + reaction field, SSD/RF). In should be noted that a second version of
142 + Ice-{\it i} (Ice-2{\it i}) was used in calculations involving SPC/E,
143 + TIP4P, and TIP5P. The unit cell of this crystal (Fig. \ref{iceiCell}B)
144 + is similar to the Ice-{\it i} unit it is extended in the direction of
145 + the (001) face and compressed along the other two faces.
146 +
147   \section{Methods}
148 +
149 + Canonical ensemble (NVT) molecular dynamics calculations were
150 + performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
151 + molecular mechanics package. All molecules were treated as rigid
152 + bodies, with orientational motion propagated using the symplectic DLM
153 + integration method. Details about the implementation of these
154 + techniques can be found in a recent publication.\cite{Meineke05}
155 +
156 + Thermodynamic integration was utilized to calculate the free energy of
157 + several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
158 + SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
159 + 400 K for all of these water models were also determined using this
160 + same technique, in order to determine melting points and generate
161 + phase diagrams. All simulations were carried out at densities
162 + resulting in a pressure of approximately 1 atm at their respective
163 + temperatures.
164 +
165 + A single thermodynamic integration involves a sequence of simulations
166 + over which the system of interest is converted into a reference system
167 + for which the free energy is known. This transformation path is then
168 + integrated in order to determine the free energy difference between
169 + the two states:
170 + \begin{equation}
171 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
172 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
173 + \end{equation}
174 + where $V$ is the interaction potential and $\lambda$ is the
175 + transformation parameter that scales the overall
176 + potential. Simulations are distributed unevenly along this path in
177 + order to sufficiently sample the regions of greatest change in the
178 + potential. Typical integrations in this study consisted of $\sim$25
179 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
180 + (near the reference state) in length.
181 +
182 + For the thermodynamic integration of molecular crystals, the Einstein
183 + Crystal is chosen as the reference state that the system is converted
184 + to over the course of the simulation. In an Einstein Crystal, the
185 + molecules are harmonically restrained at their ideal lattice locations
186 + and orientations. The partition function for a molecular crystal
187 + restrained in this fashion has been evaluated, and the Helmholtz Free
188 + Energy ({\it A}) is given by
189 + \begin{eqnarray}
190 + A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
191 + [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
192 + )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
193 + )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
194 + )^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
195 + K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
196 + (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
197 + )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
198 + \label{ecFreeEnergy}
199 + \end{eqnarray}
200 + where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
201 + \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
202 + $K_\mathrm{\omega}$ are the spring constants restraining translational
203 + motion and deflection of and rotation around the principle axis of the
204 + molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
205 + minimum potential energy of the ideal crystal. In the case of
206 + molecular liquids, the ideal vapor is chosen as the target reference
207 + state.
208 +
209 + \begin{figure}
210 + \includegraphics[width=\linewidth]{rotSpring.eps}
211 + \caption{Possible orientational motions for a restrained molecule.
212 + $\theta$ angles correspond to displacement from the body-frame {\it
213 + z}-axis, while $\omega$ angles correspond to rotation about the
214 + body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
215 + constants for the harmonic springs restraining motion in the $\theta$
216 + and $\omega$ directions.}
217 + \label{waterSpring}
218 + \end{figure}
219 +
220 + Charge, dipole, and Lennard-Jones interactions were modified by a
221 + cubic switching between 100\% and 85\% of the cutoff value (9 \AA
222 + ). By applying this function, these interactions are smoothly
223 + truncated, thereby avoiding poor energy conserving dynamics resulting
224 + from harsher truncation schemes. The effect of a long-range correction
225 + was also investigated on select model systems in a variety of
226 + manners. For the SSD/RF model, a reaction field with a fixed
227 + dielectric constant of 80 was applied in all
228 + simulations.\cite{Onsager36} For a series of the least computationally
229 + expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
230 + performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
231 + \AA\ cutoff results. Finally, results from the use of an Ewald
232 + summation were estimated for TIP3P and SPC/E by performing
233 + calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
234 + mechanics software package.\cite{Tinker} TINKER was chosen because it
235 + can also propagate the motion of rigid-bodies, and provides the most
236 + direct comparison to the results from OOPSE. The calculated energy
237 + difference in the presence and absence of PME was applied to the
238 + previous results in order to predict changes in the free energy
239 + landscape.
240  
241   \section{Results and discussion}
242  
243 + The free energy of proton ordered Ice-{\it i} was calculated and
244 + compared with the free energies of proton ordered variants of the
245 + experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
246 + as well as the higher density ice B, observed by B\`{a}ez and Clancy
247 + and thought to be the minimum free energy structure for the SPC/E
248 + model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
249 + Ice XI, the experimentally observed proton ordered variant of ice
250 + $I_h$, was investigated initially, but it was found not to be as
251 + stable as antiferroelectric variants of proton ordered or even proton
252 + disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
253 + ice $I_h$ used here is a simple antiferroelectric version that has an
254 + 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
255 + ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
256 + $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
257 + were necessary for simulations involving larger cutoff values.
258 +
259 + \begin{table*}
260 + \begin{minipage}{\linewidth}
261 + \renewcommand{\thefootnote}{\thempfootnote}
262 + \begin{center}
263 + \caption{Calculated free energies for several ice polymorphs with a
264 + variety of common water models. All calculations used a cutoff radius
265 + of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
266 + kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
267 + \begin{tabular}{ l  c  c  c  c }
268 + \hline
269 + \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
270 + \hline
271 + \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
272 + \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
273 + \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
274 + \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
275 + \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
276 + \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
277 + \end{tabular}
278 + \label{freeEnergy}
279 + \end{center}
280 + \end{minipage}
281 + \end{table*}
282 +
283 + The free energy values computed for the studied polymorphs indicate
284 + that Ice-{\it i} is the most stable state for all of the common water
285 + models studied. With the free energy at these state points, the
286 + temperature and pressure dependence of the free energy was used to
287 + project to other state points and build phase diagrams. Figures
288 + \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
289 + from the free energy results. All other models have similar structure,
290 + only the crossing points between these phases exist at different
291 + temperatures and pressures. It is interesting to note that ice $I$
292 + does not exist in either cubic or hexagonal form in any of the phase
293 + diagrams for any of the models. For purposes of this study, ice B is
294 + representative of the dense ice polymorphs. A recent study by Sanz
295 + {\it et al.} goes into detail on the phase diagrams for SPC/E and
296 + TIP4P in the high pressure regime.\cite{Sanz04}
297 +
298 + \begin{figure}
299 + \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
300 + \caption{Phase diagram for the TIP3P water model in the low pressure
301 + regime. The displayed $T_m$ and $T_b$ values are good predictions of
302 + the experimental values; however, the solid phases shown are not the
303 + experimentally observed forms. Both cubic and hexagonal ice $I$ are
304 + higher in energy and don't appear in the phase diagram.}
305 + \label{tp3phasedia}
306 + \end{figure}
307 +
308 + \begin{figure}
309 + \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
310 + \caption{Phase diagram for the SSD/RF water model in the low pressure
311 + regime. Calculations producing these results were done under an
312 + applied reaction field. It is interesting to note that this
313 + computationally efficient model (over 3 times more efficient than
314 + TIP3P) exhibits phase behavior similar to the less computationally
315 + conservative charge based models.}
316 + \label{ssdrfphasedia}
317 + \end{figure}
318 +
319 + \begin{table*}
320 + \begin{minipage}{\linewidth}
321 + \renewcommand{\thefootnote}{\thempfootnote}
322 + \begin{center}
323 + \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
324 + temperatures of several common water models compared with experiment.}
325 + \begin{tabular}{ l  c  c  c  c  c  c  c }
326 + \hline
327 + \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
328 + \hline
329 + \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
330 + \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
331 + \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
332 + \end{tabular}
333 + \label{meltandboil}
334 + \end{center}
335 + \end{minipage}
336 + \end{table*}
337 +
338 + Table \ref{meltandboil} lists the melting and boiling temperatures
339 + calculated from this work. Surprisingly, most of these models have
340 + melting points that compare quite favorably with experiment. The
341 + unfortunate aspect of this result is that this phase change occurs
342 + between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
343 + liquid state. These results are actually not contrary to previous
344 + studies in the literature. Earlier free energy studies of ice $I$
345 + using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
346 + being attributed to choice of interaction truncation and different
347 + ordered and disordered molecular arrangements). If the presence of ice
348 + B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
349 + predicted from this work. However, the $T_m$ from Ice-{\it i} is
350 + calculated at 265 K, significantly higher in temperature than the
351 + previous studies. Also of interest in these results is that SSD/E does
352 + not exhibit a melting point at 1 atm, but it shows a sublimation point
353 + at 355 K. This is due to the significant stability of Ice-{\it i} over
354 + all other polymorphs for this particular model under these
355 + conditions. While troubling, this behavior turned out to be
356 + advantageous in that it facilitated the spontaneous crystallization of
357 + Ice-{\it i}. These observations provide a warning that simulations of
358 + SSD/E as a ``liquid'' near 300 K are actually metastable and run the
359 + risk of spontaneous crystallization. However, this risk changes when
360 + applying a longer cutoff.
361 +
362 + \begin{figure}
363 + \includegraphics[width=\linewidth]{cutoffChange.eps}
364 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
365 + TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
366 + \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
367 + \AA\. These crystals are unstable at 200 K and rapidly convert into a
368 + liquid. The connecting lines are qualitative visual aid.}
369 + \label{incCutoff}
370 + \end{figure}
371 +
372 + Increasing the cutoff radius in simulations of the more
373 + computationally efficient water models was done in order to evaluate
374 + the trend in free energy values when moving to systems that do not
375 + involve potential truncation. As seen in Fig. \ref{incCutoff}, the
376 + free energy of all the ice polymorphs show a substantial dependence on
377 + cutoff radius. In general, there is a narrowing of the free energy
378 + differences while moving to greater cutoff radius. Interestingly, by
379 + increasing the cutoff radius, the free energy gap was narrowed enough
380 + in the SSD/E model that the liquid state is preferred under standard
381 + simulation conditions (298 K and 1 atm). Thus, it is recommended that
382 + simulations using this model choose interaction truncation radii
383 + greater than 9 \AA\. This narrowing trend is much more subtle in the
384 + case of SSD/RF, indicating that the free energies calculated with a
385 + reaction field present provide a more accurate picture of the free
386 + energy landscape in the absence of potential truncation.
387 +
388 + To further study the changes resulting to the inclusion of a
389 + long-range interaction correction, the effect of an Ewald summation
390 + was estimated by applying the potential energy difference do to its
391 + inclusion in systems in the presence and absence of the
392 + correction. This was accomplished by calculation of the potential
393 + energy of identical crystals with and without PME using TINKER. The
394 + free energies for the investigated polymorphs using the TIP3P and
395 + SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
396 + are not fully supported in TINKER, so the results for these models
397 + could not be estimated. The same trend pointed out through increase of
398 + cutoff radius is observed in these PME results. Ice-{\it i} is the
399 + preferred polymorph at ambient conditions for both the TIP3P and SPC/E
400 + water models; however, there is a narrowing of the free energy
401 + differences between the various solid forms. In the case of SPC/E this
402 + narrowing is significant enough that it becomes less clear cut that
403 + Ice-{\it i} is the most stable polymorph, and is possibly metastable
404 + with respect to ice B and possibly ice $I_c$. However, these results
405 + do not significantly alter the finding that the Ice-{\it i} polymorph
406 + is a stable crystal structure that should be considered when studying
407 + the phase behavior of water models.
408 +
409 + \begin{table*}
410 + \begin{minipage}{\linewidth}
411 + \renewcommand{\thefootnote}{\thempfootnote}
412 + \begin{center}
413 + \caption{The free energy of the studied ice polymorphs after applying
414 + the energy difference attributed to the inclusion of the PME
415 + long-range interaction correction. Units are kcal/mol.}
416 + \begin{tabular}{ l  c  c  c  c }
417 + \hline
418 + \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
419 + \hline
420 + \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
421 + \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
422 + \end{tabular}
423 + \label{pmeShift}
424 + \end{center}
425 + \end{minipage}
426 + \end{table*}
427 +
428   \section{Conclusions}
429  
430 + The free energy for proton ordered variants of hexagonal and cubic ice
431 + $I$, ice B, and recently discovered Ice-{\it i} where calculated under
432 + standard conditions for several common water models via thermodynamic
433 + integration. All the water models studied show Ice-{\it i} to be the
434 + minimum free energy crystal structure in the with a 9 \AA\ switching
435 + function cutoff. Calculated melting and boiling points show
436 + surprisingly good agreement with the experimental values; however, the
437 + solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
438 + interaction truncation was investigated through variation of the
439 + cutoff radius, use of a reaction field parameterized model, and
440 + estimation of the results in the presence of the Ewald summation
441 + correction. Interaction truncation has a significant effect on the
442 + computed free energy values, and may significantly alter the free
443 + energy landscape for the more complex multipoint water models. Despite
444 + these effects, these results show Ice-{\it i} to be an important ice
445 + polymorph that should be considered in simulation studies.
446 +
447 + Due to this relative stability of Ice-{\it i} in all manner of
448 + investigated simulation examples, the question arises as to possible
449 + experimental observation of this polymorph. The rather extensive past
450 + and current experimental investigation of water in the low pressure
451 + regime leads the authors to be hesitant in ascribing relevance outside
452 + of computational models, hence the descriptive name presented. That
453 + being said, there are certain experimental conditions that would
454 + provide the most ideal situation for possible observation. These
455 + include the negative pressure or stretched solid regime, small
456 + clusters in vacuum deposition environments, and in clathrate
457 + structures involving small non-polar molecules.
458 +
459   \section{Acknowledgments}
460   Support for this project was provided by the National Science
461   Foundation under grant CHE-0134881. Computation time was provided by
462 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
463 < DMR-0079647.
462 > the Notre Dame High Performance Computing Cluster and the Notre Dame
463 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
464  
465   \newpage
466  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines