ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceiPaper/iceiPaper.tex
Revision: 1485
Committed: Tue Sep 21 18:27:12 2004 UTC (19 years, 9 months ago) by chrisfen
Content type: application/x-tex
File size: 26484 byte(s)
Log Message:
updated image and a little text change

File Contents

# Content
1 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 \documentclass[11pt]{article}
3 %\documentclass[11pt]{article}
4 \usepackage{endfloat}
5 \usepackage{amsmath}
6 \usepackage{epsf}
7 \usepackage{berkeley}
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
18 \renewcommand{\baselinestretch}{1.2}
19 \renewcommand\citemid{\ } % no comma in optional reference note
20
21 \begin{document}
22
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 \\
27 Department of Chemistry and Biochemistry\\ University of Notre Dame\\
28 Notre Dame, Indiana 46556}
29
30 \date{\today}
31
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 for the new crystal were generated and
49 we await experimental confirmation of the existence of this new
50 polymorph.
51 \end{abstract}
52
53 %\narrowtext
54
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 % BODY OF TEXT
57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 distibution 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. It is clear from Fig. \ref{waterSpring} that the values
205 of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from
206 $-\pi$ to $\pi$. The partition function for a molecular crystal
207 restrained in this fashion can be evaluated analytically, and the
208 Helmholtz Free Energy ({\it A}) is given by
209 \begin{eqnarray}
210 A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
211 [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
212 )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
213 )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
214 )^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
215 K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
216 (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
217 )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
218 \label{ecFreeEnergy}
219 \end{eqnarray}
220 where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
221 potential energy of the ideal crystal.\cite{Baez95a}
222
223 \begin{figure}
224 \includegraphics[width=\linewidth]{rotSpring.eps}
225 \caption{Possible orientational motions for a restrained molecule.
226 $\theta$ angles correspond to displacement from the body-frame {\it
227 z}-axis, while $\omega$ angles correspond to rotation about the
228 body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
229 constants for the harmonic springs restraining motion in the $\theta$
230 and $\omega$ directions.}
231 \label{waterSpring}
232 \end{figure}
233
234 In the case of molecular liquids, the ideal vapor is chosen as the
235 target reference state. There are several examples of liquid state
236 free energy calculations of water models present in the
237 literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
238 typically differ in regard to the path taken for switching off the
239 interaction potential to convert the system to an ideal gas of water
240 molecules. In this study, we applied of one of the most convenient
241 methods and integrated over the $\lambda^4$ path, where all
242 interaction parameters are scaled equally by this transformation
243 parameter. This method has been shown to be reversible and provide
244 results in excellent agreement with other established
245 methods.\cite{Baez95b}
246
247 Charge, dipole, and Lennard-Jones interactions were modified by a
248 cubic switching between 100\% and 85\% of the cutoff value (9 \AA
249 ). By applying this function, these interactions are smoothly
250 truncated, thereby avoiding the poor energy conservation which results
251 from harsher truncation schemes. The effect of a long-range correction
252 was also investigated on select model systems in a variety of
253 manners. For the SSD/RF model, a reaction field with a fixed
254 dielectric constant of 80 was applied in all
255 simulations.\cite{Onsager36} For a series of the least computationally
256 expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
257 performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
258 \AA\ cutoff results. Finally, the effects of utilizing an Ewald
259 summation were estimated for TIP3P and SPC/E by performing single
260 configuration calculations with Particle-Mesh Ewald (PME) in the
261 TINKER molecular mechanics software package.\cite{Tinker} The
262 calculated energy difference in the presence and absence of PME was
263 applied to the previous results in order to predict changes to the
264 free energy landscape.
265
266 \section{Results and discussion}
267
268 The free energy of proton-ordered Ice-{\it i} was calculated and
269 compared with the free energies of proton ordered variants of the
270 experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
271 as well as the higher density ice B, observed by B\`{a}ez and Clancy
272 and thought to be the minimum free energy structure for the SPC/E
273 model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
274 Ice XI, the experimentally-observed proton-ordered variant of ice
275 $I_h$, was investigated initially, but was found to be not as stable
276 as proton disordered or antiferroelectric variants of ice $I_h$. The
277 proton ordered variant of ice $I_h$ used here is a simple
278 antiferroelectric version that we devised, and it has an 8 molecule
279 unit cell similar to other predicted antiferroelectric $I_h$
280 crystals.\cite{Davidson84} The crystals contained 648 or 1728
281 molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
282 molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
283 crystal sizes were necessary for simulations involving larger cutoff
284 values.
285
286 \begin{table*}
287 \begin{minipage}{\linewidth}
288 \renewcommand{\thefootnote}{\thempfootnote}
289 \begin{center}
290 \caption{Calculated free energies for several ice polymorphs with a
291 variety of common water models. All calculations used a cutoff radius
292 of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
293 kcal/mol. Calculated error of the final digits is in parentheses. *Ice
294 $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
295 \begin{tabular}{ l c c c c }
296 \hline
297 Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
298 \hline
299 TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
300 TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
301 TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
302 SPC/E & -12.67(2) & -12.96(2) & -13.25(3) & -13.55(2)\\
303 SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
304 SSD/RF & -11.51(2) & NA* & -12.08(3) & -12.29(2)\\
305 \end{tabular}
306 \label{freeEnergy}
307 \end{center}
308 \end{minipage}
309 \end{table*}
310
311 The free energy values computed for the studied polymorphs indicate
312 that Ice-{\it i} is the most stable state for all of the common water
313 models studied. With the calculated free energy at these state points,
314 the Gibbs-Helmholtz equation was used to project to other state points
315 and to build phase diagrams. Figures
316 \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
317 from the free energy results. All other models have similar structure,
318 although the crossing points between the phases move to slightly
319 different temperatures and pressures. It is interesting to note that
320 ice $I$ does not exist in either cubic or hexagonal form in any of the
321 phase diagrams for any of the models. For purposes of this study, ice
322 B is representative of the dense ice polymorphs. A recent study by
323 Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
324 TIP4P at higher pressures than those studied here.\cite{Sanz04}
325
326 \begin{figure}
327 \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
328 \caption{Phase diagram for the TIP3P water model in the low pressure
329 regime. The displayed $T_m$ and $T_b$ values are good predictions of
330 the experimental values; however, the solid phases shown are not the
331 experimentally observed forms. Both cubic and hexagonal ice $I$ are
332 higher in energy and don't appear in the phase diagram.}
333 \label{tp3phasedia}
334 \end{figure}
335
336 \begin{figure}
337 \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
338 \caption{Phase diagram for the SSD/RF water model in the low pressure
339 regime. Calculations producing these results were done under an
340 applied reaction field. It is interesting to note that this
341 computationally efficient model (over 3 times more efficient than
342 TIP3P) exhibits phase behavior similar to the less computationally
343 conservative charge based models.}
344 \label{ssdrfphasedia}
345 \end{figure}
346
347 \begin{table*}
348 \begin{minipage}{\linewidth}
349 \renewcommand{\thefootnote}{\thempfootnote}
350 \begin{center}
351 \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
352 temperatures at 1 atm for several common water models compared with
353 experiment. The $T_m$ and $T_s$ values from simulation correspond to a
354 transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
355 liquid or gas state.}
356 \begin{tabular}{ l c c c c c c c }
357 \hline
358 Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
359 \hline
360 $T_m$ (K) & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
361 $T_b$ (K) & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
362 $T_s$ (K) & - & - & - & - & 355(2) & - & -\\
363 \end{tabular}
364 \label{meltandboil}
365 \end{center}
366 \end{minipage}
367 \end{table*}
368
369 Table \ref{meltandboil} lists the melting and boiling temperatures
370 calculated from this work. Surprisingly, most of these models have
371 melting points that compare quite favorably with experiment. The
372 unfortunate aspect of this result is that this phase change occurs
373 between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
374 liquid state. These results are actually not contrary to previous
375 studies in the literature. Earlier free energy studies of ice $I$
376 using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
377 being attributed to choice of interaction truncation and different
378 ordered and disordered molecular
379 arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
380 Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
381 predicted from this work. However, the $T_m$ from Ice-{\it i} is
382 calculated at 265 K, significantly higher in temperature than the
383 previous studies. Also of interest in these results is that SSD/E does
384 not exhibit a melting point at 1 atm, but it shows a sublimation point
385 at 355 K. This is due to the significant stability of Ice-{\it i} over
386 all other polymorphs for this particular model under these
387 conditions. While troubling, this behavior resulted in spontaneous
388 crystallization of Ice-{\it i} and led us to investigate this
389 structure. These observations provide a warning that simulations of
390 SSD/E as a ``liquid'' near 300 K are actually metastable and run the
391 risk of spontaneous crystallization. However, this risk lessens when
392 applying a longer cutoff.
393
394 \begin{figure}
395 \includegraphics[width=\linewidth]{cutoffChange.eps}
396 \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
397 TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
398 \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
399 \AA . These crystals are unstable at 200 K and rapidly convert into
400 liquids. The connecting lines are qualitative visual aid.}
401 \label{incCutoff}
402 \end{figure}
403
404 Increasing the cutoff radius in simulations of the more
405 computationally efficient water models was done in order to evaluate
406 the trend in free energy values when moving to systems that do not
407 involve potential truncation. As seen in Fig. \ref{incCutoff}, the
408 free energy of all the ice polymorphs show a substantial dependence on
409 cutoff radius. In general, there is a narrowing of the free energy
410 differences while moving to greater cutoff radius. Interestingly, by
411 increasing the cutoff radius, the free energy gap was narrowed enough
412 in the SSD/E model that the liquid state is preferred under standard
413 simulation conditions (298 K and 1 atm). Thus, it is recommended that
414 simulations using this model choose interaction truncation radii
415 greater than 9 \AA\ . This narrowing trend is much more subtle in the
416 case of SSD/RF, indicating that the free energies calculated with a
417 reaction field present provide a more accurate picture of the free
418 energy landscape in the absence of potential truncation.
419
420 To further study the changes resulting to the inclusion of a
421 long-range interaction correction, the effect of an Ewald summation
422 was estimated by applying the potential energy difference do to its
423 inclusion in systems in the presence and absence of the
424 correction. This was accomplished by calculation of the potential
425 energy of identical crystals both with and without PME. The free
426 energies for the investigated polymorphs using the TIP3P and SPC/E
427 water models are shown in Table \ref{pmeShift}. The same trend pointed
428 out through increase of cutoff radius is observed in these PME
429 results. Ice-{\it i} is the preferred polymorph at ambient conditions
430 for both the TIP3P and SPC/E water models; however, the narrowing of
431 the free energy differences between the various solid forms is
432 significant enough that it becomes less clear that it is the most
433 stable polymorph with the SPC/E model. The free energies of Ice-{\it
434 i} and ice B nearly overlap within error, with ice $I_c$ just outside
435 as well, indicating that Ice-{\it i} might be metastable with respect
436 to ice B and possibly ice $I_c$ with SPC/E. However, these results do
437 not significantly alter the finding that the Ice-{\it i} polymorph is
438 a stable crystal structure that should be considered when studying the
439 phase behavior of water models.
440
441 \begin{table*}
442 \begin{minipage}{\linewidth}
443 \renewcommand{\thefootnote}{\thempfootnote}
444 \begin{center}
445 \caption{The free energy of the studied ice polymorphs after applying
446 the energy difference attributed to the inclusion of the PME
447 long-range interaction correction. Units are kcal/mol.}
448 \begin{tabular}{ l c c c c }
449 \hline
450 \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
451 \hline
452 TIP3P & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3)\\
453 SPC/E & -12.77(2) & -12.92(2) & -12.96(3) & -13.02(2)\\
454 \end{tabular}
455 \label{pmeShift}
456 \end{center}
457 \end{minipage}
458 \end{table*}
459
460 \section{Conclusions}
461
462 The free energy for proton ordered variants of hexagonal and cubic ice
463 $I$, ice B, and our recently discovered Ice-{\it i} structure were
464 calculated under standard conditions for several common water models
465 via thermodynamic integration. All the water models studied show
466 Ice-{\it i} to be the minimum free energy crystal structure with a 9
467 \AA\ switching function cutoff. Calculated melting and boiling points
468 show surprisingly good agreement with the experimental values;
469 however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
470 effect of interaction truncation was investigated through variation of
471 the cutoff radius, use of a reaction field parameterized model, and
472 estimation of the results in the presence of the Ewald
473 summation. Interaction truncation has a significant effect on the
474 computed free energy values, and may significantly alter the free
475 energy landscape for the more complex multipoint water models. Despite
476 these effects, these results show Ice-{\it i} to be an important ice
477 polymorph that should be considered in simulation studies.
478
479 Due to this relative stability of Ice-{\it i} in all of the
480 investigated simulation conditions, the question arises as to possible
481 experimental observation of this polymorph. The rather extensive past
482 and current experimental investigation of water in the low pressure
483 regime makes us hesitant to ascribe any relevance of this work outside
484 of the simulation community. It is for this reason that we chose a
485 name for this polymorph which involves an imaginary quantity. That
486 said, there are certain experimental conditions that would provide the
487 most ideal situation for possible observation. These include the
488 negative pressure or stretched solid regime, small clusters in vacuum
489 deposition environments, and in clathrate structures involving small
490 non-polar molecules. Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
491 our predictions for both the pair distribution function ($g_{OO}(r)$)
492 and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
493 ice-{\it i} at a temperature of 77K. In studies of the high and low
494 density forms of amorphous ice, ``spurious'' diffraction peaks have
495 been observed experimentally.\cite{Bizid87} It is possible that a
496 variant of Ice-{\it i} could explain some of this behavior; however,
497 we will leave it to our experimental colleagues to make the final
498 determination on whether this ice polymorph is named appropriately
499 (i.e. with an imaginary number) or if it can be promoted to Ice-0.
500
501 \begin{figure}
502 \includegraphics[width=\linewidth]{iceGofr.eps}
503 \caption{Radial distribution functions of ice $I_h$, $I_c$,
504 Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
505 of the SSD/RF water model at 77 K.}
506 \label{fig:gofr}
507 \end{figure}
508
509 \begin{figure}
510 \includegraphics[width=\linewidth]{sofq.eps}
511 \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
512 and Ice-{\it i}$^\prime$ at 77 K. The raw structure factors have
513 been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
514 width) to compensate for the trunction effects in our finite size
515 simulations.}
516 \label{fig:sofq}
517 \end{figure}
518
519 \section{Acknowledgments}
520 Support for this project was provided by the National Science
521 Foundation under grant CHE-0134881. Computation time was provided by
522 the Notre Dame High Performance Computing Cluster and the Notre Dame
523 Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
524
525 \newpage
526
527 \bibliographystyle{jcp}
528 \bibliography{iceiPaper}
529
530
531 \end{document}