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 1555 by chrisfen, Mon Oct 11 20:14:35 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 simulation-predicted ice polymorph which is more
24 > stable than Ice $I_h$ for point-charge and point-dipole water models}
25  
26 < \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
27 < \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
27 <
28 < \address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
26 > \author{Christopher J. Fennell and J. Daniel Gezelter \\
27 > Department of Chemistry and Biochemistry\\ University of Notre Dame\\
28   Notre Dame, Indiana 46556}
29  
30   \date{\today}
31  
32 < %\maketitle
32 > \maketitle
33   %\doublespacing
34  
35   \begin{abstract}
36 + The absolute free energies of several ice polymorphs which are stable
37 + at low pressures were calculated using thermodynamic integration to a
38 + reference system (the Einstein crystal).  These integrations were
39 + performed for most of the common water models (SPC/E, TIP3P, TIP4P,
40 + TIP5P, SSD/E, SSD/RF). Ice-{\it i}, a structure we recently observed
41 + crystallizing at room temperature for one of the single-point water
42 + models, was determined to be the stable crystalline state (at 1 atm)
43 + for {\it all} the water models investigated.  Phase diagrams were
44 + generated, and phase coexistence lines were determined for all of the
45 + known low-pressure ice structures under all of these water models.
46 + Additionally, potential truncation was shown to have an effect on the
47 + calculated free energies, and can result in altered free energy
48 + landscapes.  Structure factor predictions for the new crystal were
49 + generated and we await experimental confirmation of the existence of
50 + this new polymorph.
51   \end{abstract}
52  
39 – \maketitle
40 –
41 – \newpage
42 –
53   %\narrowtext
54  
55   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 48 | Line 58 | Notre Dame, Indiana 46556}
58  
59   \section{Introduction}
60  
61 + Water has proven to be a challenging substance to depict in
62 + simulations, and a variety of models have been developed to describe
63 + its behavior under varying simulation
64 + conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
65 + These models have been used to investigate important physical
66 + phenomena like phase transitions, transport properties, and the
67 + hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
68 + choice of models available, it is only natural to compare the models
69 + under interesting thermodynamic conditions in an attempt to clarify
70 + the limitations of each of the
71 + models.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two
72 + important properties to quantify are the Gibbs and Helmholtz free
73 + energies, particularly for the solid forms of water.  Difficulty in
74 + these types of studies typically arises from the assortment of
75 + possible crystalline polymorphs that water adopts over a wide range of
76 + pressures and temperatures.  There are currently 13 recognized forms
77 + of ice, and it is a challenging task to investigate the entire free
78 + energy landscape.\cite{Sanz04} Ideally, research is focused on the
79 + phases having the lowest free energy at a given state point, because
80 + these phases will dictate the relevant transition temperatures and
81 + pressures for the model.
82 +
83 + In this paper, standard reference state methods were applied to known
84 + crystalline water polymorphs in the low pressure regime.  This work is
85 + unique in that one of the crystal lattices was arrived at through
86 + crystallization of a computationally efficient water model under
87 + constant pressure and temperature conditions. Crystallization events
88 + are interesting in and of themselves;\cite{Matsumoto02,Yamada02}
89 + however, the crystal structure obtained in this case is different from
90 + any previously observed ice polymorphs in experiment or
91 + simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
92 + to indicate its origin in computational simulation. The unit cell
93 + (Fig. \ref{iceiCell}A) consists of eight water molecules that stack in
94 + rows of interlocking water tetramers. Proton ordering can be
95 + accomplished by orienting two of the molecules so that both of their
96 + donated hydrogen bonds are internal to their tetramer
97 + (Fig. \ref{protOrder}). As expected in an ice crystal constructed of
98 + water tetramers, the hydrogen bonds are not as linear as those
99 + observed in ice $I_h$, however the interlocking of these subunits
100 + appears to provide significant stabilization to the overall
101 + crystal. The arrangement of these tetramers results in surrounding
102 + open octagonal cavities that are typically greater than 6.3 \AA\ in
103 + diameter. This relatively open overall structure leads to crystals
104 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
105 +
106 + \begin{figure}
107 + \includegraphics[width=\linewidth]{unitCell.eps}
108 + \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it i}$^\prime$,
109 + the elongated variant of Ice-{\it i}.  The spheres represent the
110 + center-of-mass locations of the water molecules.  The $a$ to $c$
111 + ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
112 + $a:2.1214c$ and $a:1.7850c$ respectively.}
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 from our previous study indicated that Ice-{\it i} is the
126 + minimum energy crystal structure for the single point water models we
127 + had investigated (for discussions on these single point dipole models,
128 + see our previous work and related
129 + articles).\cite{Fennell04,Liu96,Bratko85} Those results only
130 + considered energetic stabilization and neglected entropic
131 + contributions to the overall free energy. To address this issue, we
132 + have calculated the absolute free energy of this crystal using
133 + thermodynamic integration and compared to the free energies of cubic
134 + and hexagonal ice $I$ (the experimental low density ice polymorphs)
135 + and ice B (a higher density, but very stable crystal structure
136 + observed by B\`{a}ez and Clancy in free energy studies of
137 + SPC/E).\cite{Baez95b} This work includes results for the water model
138 + from which Ice-{\it i} was crystallized (SSD/E) in addition to several
139 + common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
140 + field parametrized single point dipole water model (SSD/RF). It should
141 + be noted that a second version of Ice-{\it i} (Ice-{\it i}$^\prime$)
142 + was used in calculations involving SPC/E, TIP4P, and TIP5P. The unit
143 + cell of this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it
144 + i} unit it is extended in the direction of the (001) face and
145 + compressed along the other two faces.  There is typically a small
146 + distortion of proton ordered Ice-{\it i}$^\prime$ that converts the
147 + normally square tetramer into a rhombus with alternating approximately
148 + 85 and 95 degree angles.  The degree of this distortion is model
149 + dependent and significant enough to split the tetramer diagonal
150 + location peak in the radial distribution function.
151 +
152   \section{Methods}
153  
154 + Canonical ensemble (NVT) molecular dynamics calculations were
155 + performed using the OOPSE molecular mechanics package.\cite{Meineke05}
156 + All molecules were treated as rigid bodies, with orientational motion
157 + propagated using the symplectic DLM integration method. Details about
158 + the implementation of this technique can be found in a recent
159 + publication.\cite{Dullweber1997}
160 +
161 + Thermodynamic integration is an established technique for
162 + determination of free energies of condensed phases of
163 + materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
164 + method, implemented in the same manner illustrated by B\`{a}ez and
165 + Clancy, was utilized to calculate the free energy of several ice
166 + crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
167 + SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
168 + and 400 K for all of these water models were also determined using
169 + this same technique in order to determine melting points and to
170 + generate phase diagrams. All simulations were carried out at densities
171 + which correspond to a pressure of approximately 1 atm at their
172 + respective temperatures.
173 +
174 + Thermodynamic integration involves a sequence of simulations during
175 + which the system of interest is converted into a reference system for
176 + which the free energy is known analytically. This transformation path
177 + is then integrated in order to determine the free energy difference
178 + between the two states:
179 + \begin{equation}
180 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
181 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
182 + \end{equation}
183 + where $V$ is the interaction potential and $\lambda$ is the
184 + transformation parameter that scales the overall
185 + potential. Simulations are distributed strategically along this path
186 + in order to sufficiently sample the regions of greatest change in the
187 + potential. Typical integrations in this study consisted of $\sim$25
188 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
189 + (near the reference state) in length.
190 +
191 + For the thermodynamic integration of molecular crystals, the Einstein
192 + crystal was chosen as the reference system. In an Einstein crystal,
193 + the molecules are restrained at their ideal lattice locations and
194 + orientations. Using harmonic restraints, as applied by B\`{a}ez and
195 + Clancy, the total potential for this reference crystal
196 + ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
197 + \begin{equation}
198 + V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
199 + \frac{K_\omega\omega^2}{2},
200 + \end{equation}
201 + where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
202 + the spring constants restraining translational motion and deflection
203 + of and rotation around the principle axis of the molecule
204 + respectively.  These spring constants are typically calculated from
205 + the mean-square displacements of water molecules in an unrestrained
206 + ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
207 + mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
208 + 17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
209 + the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
210 + from $-\pi$ to $\pi$.  The partition function for a molecular crystal
211 + restrained in this fashion can be evaluated analytically, and the
212 + Helmholtz Free Energy ({\it A}) is given by
213 + \begin{eqnarray}
214 + A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
215 + [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
216 + )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
217 + )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
218 + )^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
219 + K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
220 + (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
221 + )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
222 + \label{ecFreeEnergy}
223 + \end{eqnarray}
224 + where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
225 + potential energy of the ideal crystal.\cite{Baez95a}
226 +
227 + \begin{figure}
228 + \includegraphics[width=\linewidth]{rotSpring.eps}
229 + \caption{Possible orientational motions for a restrained molecule.
230 + $\theta$ angles correspond to displacement from the body-frame {\it
231 + z}-axis, while $\omega$ angles correspond to rotation about the
232 + body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
233 + constants for the harmonic springs restraining motion in the $\theta$
234 + and $\omega$ directions.}
235 + \label{waterSpring}
236 + \end{figure}
237 +
238 + In the case of molecular liquids, the ideal vapor is chosen as the
239 + target reference state.  There are several examples of liquid state
240 + free energy calculations of water models present in the
241 + literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
242 + typically differ in regard to the path taken for switching off the
243 + interaction potential to convert the system to an ideal gas of water
244 + molecules.  In this study, we applied of one of the most convenient
245 + methods and integrated over the $\lambda^4$ path, where all
246 + interaction parameters are scaled equally by this transformation
247 + parameter.  This method has been shown to be reversible and provide
248 + results in excellent agreement with other established
249 + methods.\cite{Baez95b}
250 +
251 + Charge, dipole, and Lennard-Jones interactions were modified by a
252 + cubic switching between 100\% and 85\% of the cutoff value (9 \AA
253 + ). By applying this function, these interactions are smoothly
254 + truncated, thereby avoiding the poor energy conservation which results
255 + from harsher truncation schemes. The effect of a long-range correction
256 + was also investigated on select model systems in a variety of
257 + manners. For the SSD/RF model, a reaction field with a fixed
258 + dielectric constant of 80 was applied in all
259 + simulations.\cite{Onsager36} For a series of the least computationally
260 + expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
261 + performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
262 + \AA\ cutoff results. Finally, the effects of utilizing an Ewald
263 + summation were estimated for TIP3P and SPC/E by performing single
264 + configuration calculations with Particle-Mesh Ewald (PME) in the
265 + TINKER molecular mechanics software package.\cite{Tinker} The
266 + calculated energy difference in the presence and absence of PME was
267 + applied to the previous results in order to predict changes to the
268 + free energy landscape.
269 +
270   \section{Results and discussion}
271 +
272 + The free energy of proton-ordered Ice-{\it i} was calculated and
273 + compared with the free energies of proton ordered variants of the
274 + experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
275 + as well as the higher density ice B, observed by B\`{a}ez and Clancy
276 + and thought to be the minimum free energy structure for the SPC/E
277 + model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
278 + Ice XI, the experimentally-observed proton-ordered variant of ice
279 + $I_h$, was investigated initially, but was found to be not as stable
280 + as proton disordered or antiferroelectric variants of ice $I_h$. The
281 + proton ordered variant of ice $I_h$ used here is a simple
282 + antiferroelectric version that we devised, and it has an 8 molecule
283 + unit cell similar to other predicted antiferroelectric $I_h$
284 + crystals.\cite{Davidson84} The crystals contained 648 or 1728
285 + molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
286 + molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
287 + crystal sizes were necessary for simulations involving larger cutoff
288 + values.
289  
290 + \begin{table*}
291 + \begin{minipage}{\linewidth}
292 + \begin{center}
293 +
294 + \caption{Calculated free energies for several ice polymorphs with a
295 + variety of common water models. All calculations used a cutoff radius
296 + of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
297 + kcal/mol. Calculated error of the final digits is in parentheses.}
298 +
299 + \begin{tabular}{lcccc}
300 + \hline
301 + Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
302 + \hline
303 + TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
304 + TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
305 + TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
306 + SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & -13.55(2)\\
307 + SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
308 + SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\
309 + \end{tabular}
310 + \label{freeEnergy}
311 + \end{center}
312 + \end{minipage}
313 + \end{table*}
314 +
315 + The free energy values computed for the studied polymorphs indicate
316 + that Ice-{\it i} is the most stable state for all of the common water
317 + models studied. With the calculated free energy at these state points,
318 + the Gibbs-Helmholtz equation was used to project to other state points
319 + and to build phase diagrams.  Figures
320 + \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
321 + from the free energy results. All other models have similar structure,
322 + although the crossing points between the phases move to slightly
323 + different temperatures and pressures. It is interesting to note that
324 + ice $I$ does not exist in either cubic or hexagonal form in any of the
325 + phase diagrams for any of the models. For purposes of this study, ice
326 + B is representative of the dense ice polymorphs. A recent study by
327 + Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
328 + TIP4P at higher pressures than those studied here.\cite{Sanz04}
329 +
330 + \begin{figure}
331 + \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
332 + \caption{Phase diagram for the TIP3P water model in the low pressure
333 + regime. The displayed $T_m$ and $T_b$ values are good predictions of
334 + the experimental values; however, the solid phases shown are not the
335 + experimentally observed forms. Both cubic and hexagonal ice $I$ are
336 + higher in energy and don't appear in the phase diagram.}
337 + \label{tp3phasedia}
338 + \end{figure}
339 +
340 + \begin{figure}
341 + \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
342 + \caption{Phase diagram for the SSD/RF water model in the low pressure
343 + regime. Calculations producing these results were done under an
344 + applied reaction field. It is interesting to note that this
345 + computationally efficient model (over 3 times more efficient than
346 + TIP3P) exhibits phase behavior similar to the less computationally
347 + conservative charge based models.}
348 + \label{ssdrfphasedia}
349 + \end{figure}
350 +
351 + \begin{table*}
352 + \begin{minipage}{\linewidth}
353 + \begin{center}
354 +
355 + \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
356 + temperatures at 1 atm for several common water models compared with
357 + experiment. The $T_m$ and $T_s$ values from simulation correspond to a
358 + transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
359 + liquid or gas state.}
360 +
361 + \begin{tabular}{lccccccc}
362 + \hline
363 + Equilibrium Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
364 + \hline
365 + $T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
366 + $T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
367 + $T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
368 + \end{tabular}
369 + \label{meltandboil}
370 + \end{center}
371 + \end{minipage}
372 + \end{table*}
373 +
374 + Table \ref{meltandboil} lists the melting and boiling temperatures
375 + calculated from this work. Surprisingly, most of these models have
376 + melting points that compare quite favorably with experiment. The
377 + unfortunate aspect of this result is that this phase change occurs
378 + between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
379 + liquid state. These results are actually not contrary to previous
380 + studies in the literature. Earlier free energy studies of ice $I$
381 + using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
382 + being attributed to choice of interaction truncation and different
383 + ordered and disordered molecular
384 + arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
385 + Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
386 + predicted from this work. However, the $T_m$ from Ice-{\it i} is
387 + calculated at 265 K, significantly higher in temperature than the
388 + previous studies. Also of interest in these results is that SSD/E does
389 + not exhibit a melting point at 1 atm, but it shows a sublimation point
390 + at 355 K. This is due to the significant stability of Ice-{\it i} over
391 + all other polymorphs for this particular model under these
392 + conditions. While troubling, this behavior resulted in spontaneous
393 + crystallization of Ice-{\it i} and led us to investigate this
394 + structure. These observations provide a warning that simulations of
395 + SSD/E as a ``liquid'' near 300 K are actually metastable and run the
396 + risk of spontaneous crystallization. However, this risk lessens when
397 + applying a longer cutoff.
398 +
399 + \begin{figure}
400 + \includegraphics[width=\linewidth]{cutoffChange.eps}
401 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
402 + TIP3P, and (C) SSD/RF with a reaction field. Both SSD/E and TIP3P show
403 + significant cutoff radius dependence of the free energy and appear to
404 + converge when moving to cutoffs greater than 12 \AA. Use of a reaction
405 + field with SSD/RF results in free energies that exhibit minimal cutoff
406 + radius dependence.}
407 + \label{incCutoff}
408 + \end{figure}
409 +
410 + Increasing the cutoff radius in simulations of the more
411 + computationally efficient water models was done in order to evaluate
412 + the trend in free energy values when moving to systems that do not
413 + involve potential truncation. As seen in Fig. \ref{incCutoff}, the
414 + free energy of all the ice polymorphs for the SSD/E and TIP3P models
415 + show a substantial dependence on cutoff radius. In general, there is a
416 + narrowing of the free energy differences while moving to greater
417 + cutoff radii.  As the free energies for the polymorphs converge, the
418 + stability advantage that Ice-{\it i} exhibits is reduced; however, it
419 + remains the most stable polymorph for both of these models over the
420 + depicted range for both models. This narrowing trend is not
421 + significant in the case of SSD/RF, indicating that the free energies
422 + calculated with a reaction field present provide, at minimal
423 + computational cost, a more accurate picture of the free energy
424 + landscape in the absence of potential truncation.  Interestingly,
425 + increasing the cutoff radius a mere 1.5 \AA\ with the SSD/E model
426 + destabilizes the Ice-{\it i} polymorph enough that the liquid state is
427 + preferred under standard simulation conditions (298 K and 1
428 + atm). Thus, it is recommended that simulations using this model choose
429 + interaction truncation radii greater than 9 \AA. Considering this
430 + stabilization provided by smaller cutoffs, it is not surprising that
431 + crystallization into Ice-{\it i} was observed with SSD/E.  The choice
432 + of a 9 \AA\ cutoff in the previous simulations gives the Ice-{\it i}
433 + polymorph a greater than 1 kcal/mol lower free energy than the ice
434 + $I_\textrm{h}$ starting configurations.
435 +
436 + To further study the changes resulting to the inclusion of a
437 + long-range interaction correction, the effect of an Ewald summation
438 + was estimated by applying the potential energy difference do to its
439 + inclusion in systems in the presence and absence of the correction.
440 + This was accomplished by calculation of the potential energy of
441 + identical crystals both with and without PME.  The free energies for
442 + the investigated polymorphs using the TIP3P and SPC/E water models are
443 + shown in Table \ref{pmeShift}.  The same trend pointed out through
444 + increase of cutoff radius is observed in these PME results. Ice-{\it
445 + i} is the preferred polymorph at ambient conditions for both the TIP3P
446 + and SPC/E water models; however, the narrowing of the free energy
447 + differences between the various solid forms with the SPC/E model is
448 + significant enough that it becomes less clear that it is the most
449 + stable polymorph.  The free energies of Ice-{\it i} and $I_\textrm{c}$
450 + overlap within error, while ice B and $I_\textrm{h}$ are just outside
451 + at t slightly higher free energy.  This indicates that with SPC/E,
452 + Ice-{\it i} might be metastable with all the studied polymorphs,
453 + particularly ice $I_\textrm{c}$. However, these results do not
454 + significantly alter the finding that the Ice-{\it i} polymorph is a
455 + stable crystal structure that should be considered when studying the
456 + phase behavior of water models.
457 +
458 + \begin{table*}
459 + \begin{minipage}{\linewidth}
460 + \begin{center}
461 +
462 + \caption{The free energy of the studied ice polymorphs after applying
463 + the energy difference attributed to the inclusion of the PME
464 + long-range interaction correction. Units are kcal/mol.}
465 +
466 + \begin{tabular}{ccccc}
467 + \hline
468 + Water Model &  $I_h$ & $I_c$ &  B & Ice-{\it i} \\
469 + \hline
470 + TIP3P  & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3) \\
471 + SPC/E  & -12.97(2) & -13.00(2) & -12.96(3) & -13.02(2) \\
472 + \end{tabular}
473 + \label{pmeShift}
474 + \end{center}
475 + \end{minipage}
476 + \end{table*}
477 +
478   \section{Conclusions}
479  
480 + The free energy for proton ordered variants of hexagonal and cubic ice
481 + $I$, ice B, and our recently discovered Ice-{\it i} structure were
482 + calculated under standard conditions for several common water models
483 + via thermodynamic integration. All the water models studied show
484 + Ice-{\it i} to be the minimum free energy crystal structure with a 9
485 + \AA\ switching function cutoff. Calculated melting and boiling points
486 + show surprisingly good agreement with the experimental values;
487 + however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
488 + effect of interaction truncation was investigated through variation of
489 + the cutoff radius, use of a reaction field parameterized model, and
490 + estimation of the results in the presence of the Ewald
491 + summation. Interaction truncation has a significant effect on the
492 + computed free energy values, and may significantly alter the free
493 + energy landscape for the more complex multipoint water models. Despite
494 + these effects, these results show Ice-{\it i} to be an important ice
495 + polymorph that should be considered in simulation studies.
496 +
497 + Due to this relative stability of Ice-{\it i} in all of the
498 + investigated simulation conditions, the question arises as to possible
499 + experimental observation of this polymorph.  The rather extensive past
500 + and current experimental investigation of water in the low pressure
501 + regime makes us hesitant to ascribe any relevance of this work outside
502 + of the simulation community.  It is for this reason that we chose a
503 + name for this polymorph which involves an imaginary quantity.  That
504 + said, there are certain experimental conditions that would provide the
505 + most ideal situation for possible observation. These include the
506 + negative pressure or stretched solid regime, small clusters in vacuum
507 + deposition environments, and in clathrate structures involving small
508 + non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
509 + our predictions for both the pair distribution function ($g_{OO}(r)$)
510 + and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
511 + ice-{\it i} at a temperature of 77K.  In studies of the high and low
512 + density forms of amorphous ice, ``spurious'' diffraction peaks have
513 + been observed experimentally.\cite{Bizid87} It is possible that a
514 + variant of Ice-{\it i} could explain some of this behavior; however,
515 + we will leave it to our experimental colleagues to make the final
516 + determination on whether this ice polymorph is named appropriately
517 + (i.e. with an imaginary number) or if it can be promoted to Ice-0.
518 +
519 + \begin{figure}
520 + \includegraphics[width=\linewidth]{iceGofr.eps}
521 + \caption{Radial distribution functions of ice $I_h$, $I_c$,
522 + Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
523 + of the SSD/RF water model at 77 K.}
524 + \label{fig:gofr}
525 + \end{figure}
526 +
527 + \begin{figure}
528 + \includegraphics[width=\linewidth]{sofq.eps}
529 + \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
530 + and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
531 + been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
532 + width) to compensate for the trunction effects in our finite size
533 + simulations.}
534 + \label{fig:sofq}
535 + \end{figure}
536 +
537   \section{Acknowledgments}
538   Support for this project was provided by the National Science
539   Foundation under grant CHE-0134881. Computation time was provided by
540 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
541 < DMR-0079647.
540 > the Notre Dame High Performance Computing Cluster and the Notre Dame
541 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
542  
543   \newpage
544  

Diff Legend

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