ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceiPaper/iceiPaper.tex
Revision: 1473
Committed: Fri Sep 17 14:54:49 2004 UTC (19 years, 9 months ago) by chrisfen
Content type: application/x-tex
File size: 26589 byte(s)
Log Message:
fixes in error values

File Contents

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