ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceiPaper/iceiPaper.tex
Revision: 1528
Committed: Wed Oct 6 14:45:47 2004 UTC (19 years, 9 months ago) by chrisfen
Content type: application/x-tex
File size: 26886 byte(s)
Log Message:
preparing for better cutoff plots

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 gezelter 1475 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 gezelter 1463 generated, and phase coexistence lines were determined for all of the
45 gezelter 1475 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 chrisfen 1487 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 chrisfen 1453 \end{abstract}
52    
53     %\narrowtext
54    
55     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56     % BODY OF TEXT
57     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58    
59     \section{Introduction}
60    
61 chrisfen 1459 Water has proven to be a challenging substance to depict in
62 gezelter 1463 simulations, and a variety of models have been developed to describe
63     its behavior under varying simulation
64 gezelter 1477 conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
65 gezelter 1463 These models have been used to investigate important physical
66 gezelter 1475 phenomena like phase transitions, transport properties, and the
67 chrisfen 1471 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 gezelter 1475 important properties to quantify are the Gibbs and Helmholtz free
73 chrisfen 1471 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 gezelter 1475 these phases will dictate the relevant transition temperatures and
81 chrisfen 1471 pressures for the model.
82 chrisfen 1459
83 gezelter 1465 In this paper, standard reference state methods were applied to known
84 chrisfen 1471 crystalline water polymorphs in the low pressure regime. This work is
85 gezelter 1475 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 chrisfen 1459 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 gezelter 1463
106 chrisfen 1460 \begin{figure}
107 gezelter 1463 \includegraphics[width=\linewidth]{unitCell.eps}
108 chrisfen 1483 \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 chrisfen 1471 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 chrisfen 1460 \label{iceiCell}
114     \end{figure}
115 gezelter 1463
116 chrisfen 1460 \begin{figure}
117 gezelter 1463 \includegraphics[width=\linewidth]{orderedIcei.eps}
118 chrisfen 1460 \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 chrisfen 1459
125 gezelter 1465 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 gezelter 1475 had investigated (for discussions on these single point dipole models,
128     see our previous work and related
129 chrisfen 1471 articles).\cite{Fennell04,Liu96,Bratko85} Those results only
130 gezelter 1465 considered energetic stabilization and neglected entropic
131 gezelter 1475 contributions to the overall free energy. To address this issue, we
132     have calculated the absolute free energy of this crystal using
133 gezelter 1465 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 chrisfen 1483 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 chrisfen 1484 compressed along the other two faces. There is typically a small
146 chrisfen 1485 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 chrisfen 1488 location peak in the radial distribution function.
151 chrisfen 1459
152 chrisfen 1453 \section{Methods}
153    
154 chrisfen 1454 Canonical ensemble (NVT) molecular dynamics calculations were
155 gezelter 1465 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 chrisfen 1471 the implementation of this technique can be found in a recent
159 gezelter 1468 publication.\cite{Dullweber1997}
160 chrisfen 1454
161 chrisfen 1471 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 gezelter 1475 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 chrisfen 1454
174 gezelter 1475 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 chrisfen 1458 \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 chrisfen 1459 transformation parameter that scales the overall
185 chrisfen 1471 potential. Simulations are distributed strategically along this path
186     in order to sufficiently sample the regions of greatest change in the
187 chrisfen 1459 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 chrisfen 1458
191 chrisfen 1454 For the thermodynamic integration of molecular crystals, the Einstein
192 chrisfen 1471 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 gezelter 1465 restrained in this fashion can be evaluated analytically, and the
208     Helmholtz Free Energy ({\it A}) is given by
209 chrisfen 1454 \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 chrisfen 1471 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 gezelter 1463
223 chrisfen 1456 \begin{figure}
224 gezelter 1463 \includegraphics[width=\linewidth]{rotSpring.eps}
225 chrisfen 1456 \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 chrisfen 1454
234 chrisfen 1471 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 gezelter 1475 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 chrisfen 1471
247 chrisfen 1456 Charge, dipole, and Lennard-Jones interactions were modified by a
248 chrisfen 1462 cubic switching between 100\% and 85\% of the cutoff value (9 \AA
249     ). By applying this function, these interactions are smoothly
250 gezelter 1465 truncated, thereby avoiding the poor energy conservation which results
251 chrisfen 1462 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 gezelter 1475 \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 chrisfen 1454
266 chrisfen 1456 \section{Results and discussion}
267 chrisfen 1454
268 gezelter 1475 The free energy of proton-ordered Ice-{\it i} was calculated and
269 chrisfen 1456 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 gezelter 1465 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 gezelter 1475 antiferroelectric version that we devised, and it has an 8 molecule
279 chrisfen 1473 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 chrisfen 1454
286 chrisfen 1456 \begin{table*}
287     \begin{minipage}{\linewidth}
288     \begin{center}
289 gezelter 1489
290 chrisfen 1456 \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 chrisfen 1528 kcal/mol. Calculated error of the final digits is in parentheses.}
294 gezelter 1489
295     \begin{tabular}{lcccc}
296 gezelter 1463 \hline
297 chrisfen 1466 Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
298 gezelter 1463 \hline
299 chrisfen 1473 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 chrisfen 1528 SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\
305 chrisfen 1456 \end{tabular}
306     \label{freeEnergy}
307     \end{center}
308     \end{minipage}
309     \end{table*}
310 chrisfen 1453
311 chrisfen 1456 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 gezelter 1475 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 chrisfen 1456 \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
317     from the free energy results. All other models have similar structure,
318 gezelter 1475 although the crossing points between the phases move to slightly
319 gezelter 1465 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 gezelter 1475 TIP4P at higher pressures than those studied here.\cite{Sanz04}
325 gezelter 1463
326 chrisfen 1456 \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 gezelter 1463
336 chrisfen 1456 \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     \begin{center}
350 gezelter 1489
351 chrisfen 1456 \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
352 chrisfen 1466 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 gezelter 1489
357     \begin{tabular}{lccccccc}
358 gezelter 1463 \hline
359 gezelter 1489 Equilibrium 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 gezelter 1475 conditions. While troubling, this behavior resulted in spontaneous
389     crystallization of Ice-{\it i} and led us to investigate this
390     structure. These observations provide a warning that simulations of
391 chrisfen 1456 SSD/E as a ``liquid'' near 300 K are actually metastable and run the
392 gezelter 1475 risk of spontaneous crystallization. However, this risk lessens when
393 chrisfen 1456 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 chrisfen 1528 TIP3P, and (C) SSD/RF with a reaction field. Both SSD/E and TIP3P show
399     significant cutoff radius dependence of the free energy and appear to
400     converge when moving to cutoffs greater than 12 \AA. Use of a reaction
401     field with SSD/RF results in free energies that exhibit minimal cutoff
402     radius dependence.}
403 chrisfen 1458 \label{incCutoff}
404     \end{figure}
405    
406 chrisfen 1457 Increasing the cutoff radius in simulations of the more
407     computationally efficient water models was done in order to evaluate
408     the trend in free energy values when moving to systems that do not
409     involve potential truncation. As seen in Fig. \ref{incCutoff}, the
410 chrisfen 1528 free energy of all the ice polymorphs for the SSD/E and TIP3P models
411     show a substantial dependence on cutoff radius. In general, there is a
412     narrowing of the free energy differences while moving to greater
413     cutoff radii. As the free energies for the polymorphs converge, the
414     stability advantage that Ice-{\it i} exhibits is reduced; however, it
415     remains the most stable polymorph for both of these models over the
416     depicted range for both models. This narrowing trend is not
417     significant in the case of SSD/RF, indicating that the free energies
418     calculated with a reaction field present provide, at minimal
419     computational cost, a more accurate picture of the free energy
420     landscape in the absence of potential truncation. Interestingly,
421     increasing the cutoff radius a mere 1.5 \AA\ with the SSD/E model
422     destabilizes the Ice-{\it i} polymorph enough that the liquid state is
423     preferred under standard simulation conditions (298 K and 1
424     atm). Thus, it is recommended that simulations using this model choose
425     interaction truncation radii greater than 9 \AA. Considering this
426     stabilization provided by smaller cutoffs, it is not surprising that
427     crystallization into Ice-{\it i} was observed with SSD/E. The choice
428     of a 9 \AA\ cutoff in the previous simulations gives the Ice-{\it i}
429     polymorph a greater than 1 kcal/mol lower free energy than the ice
430     $I_\textrm{h}$ starting configurations.
431 chrisfen 1456
432 chrisfen 1457 To further study the changes resulting to the inclusion of a
433     long-range interaction correction, the effect of an Ewald summation
434     was estimated by applying the potential energy difference do to its
435     inclusion in systems in the presence and absence of the
436     correction. This was accomplished by calculation of the potential
437 gezelter 1475 energy of identical crystals both with and without PME. The free
438     energies for the investigated polymorphs using the TIP3P and SPC/E
439     water models are shown in Table \ref{pmeShift}. The same trend pointed
440     out through increase of cutoff radius is observed in these PME
441 chrisfen 1471 results. Ice-{\it i} is the preferred polymorph at ambient conditions
442     for both the TIP3P and SPC/E water models; however, the narrowing of
443     the free energy differences between the various solid forms is
444     significant enough that it becomes less clear that it is the most
445 chrisfen 1474 stable polymorph with the SPC/E model. The free energies of Ice-{\it
446     i} and ice B nearly overlap within error, with ice $I_c$ just outside
447     as well, indicating that Ice-{\it i} might be metastable with respect
448     to ice B and possibly ice $I_c$ with SPC/E. However, these results do
449     not significantly alter the finding that the Ice-{\it i} polymorph is
450     a stable crystal structure that should be considered when studying the
451     phase behavior of water models.
452 chrisfen 1456
453 chrisfen 1457 \begin{table*}
454     \begin{minipage}{\linewidth}
455     \begin{center}
456 gezelter 1489
457 chrisfen 1458 \caption{The free energy of the studied ice polymorphs after applying
458     the energy difference attributed to the inclusion of the PME
459     long-range interaction correction. Units are kcal/mol.}
460 gezelter 1489
461     \begin{tabular}{ccccc}
462 gezelter 1463 \hline
463 gezelter 1489 Water Model & $I_h$ & $I_c$ & B & Ice-{\it i} \\
464 gezelter 1463 \hline
465 gezelter 1489 TIP3P & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3) \\
466     SPC/E & -12.77(2) & -12.92(2) & -12.96(3) & -13.02(2) \\
467 chrisfen 1457 \end{tabular}
468     \label{pmeShift}
469     \end{center}
470     \end{minipage}
471     \end{table*}
472    
473 chrisfen 1453 \section{Conclusions}
474    
475 chrisfen 1458 The free energy for proton ordered variants of hexagonal and cubic ice
476 gezelter 1475 $I$, ice B, and our recently discovered Ice-{\it i} structure were
477     calculated under standard conditions for several common water models
478     via thermodynamic integration. All the water models studied show
479     Ice-{\it i} to be the minimum free energy crystal structure with a 9
480     \AA\ switching function cutoff. Calculated melting and boiling points
481     show surprisingly good agreement with the experimental values;
482     however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
483     effect of interaction truncation was investigated through variation of
484     the cutoff radius, use of a reaction field parameterized model, and
485 gezelter 1465 estimation of the results in the presence of the Ewald
486     summation. Interaction truncation has a significant effect on the
487 chrisfen 1459 computed free energy values, and may significantly alter the free
488     energy landscape for the more complex multipoint water models. Despite
489     these effects, these results show Ice-{\it i} to be an important ice
490     polymorph that should be considered in simulation studies.
491 chrisfen 1458
492 gezelter 1475 Due to this relative stability of Ice-{\it i} in all of the
493     investigated simulation conditions, the question arises as to possible
494 gezelter 1465 experimental observation of this polymorph. The rather extensive past
495 chrisfen 1459 and current experimental investigation of water in the low pressure
496 gezelter 1465 regime makes us hesitant to ascribe any relevance of this work outside
497     of the simulation community. It is for this reason that we chose a
498     name for this polymorph which involves an imaginary quantity. That
499     said, there are certain experimental conditions that would provide the
500     most ideal situation for possible observation. These include the
501     negative pressure or stretched solid regime, small clusters in vacuum
502     deposition environments, and in clathrate structures involving small
503 gezelter 1469 non-polar molecules. Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
504     our predictions for both the pair distribution function ($g_{OO}(r)$)
505 chrisfen 1479 and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
506     ice-{\it i} at a temperature of 77K. In studies of the high and low
507     density forms of amorphous ice, ``spurious'' diffraction peaks have
508     been observed experimentally.\cite{Bizid87} It is possible that a
509     variant of Ice-{\it i} could explain some of this behavior; however,
510     we will leave it to our experimental colleagues to make the final
511     determination on whether this ice polymorph is named appropriately
512     (i.e. with an imaginary number) or if it can be promoted to Ice-0.
513 chrisfen 1459
514 chrisfen 1467 \begin{figure}
515     \includegraphics[width=\linewidth]{iceGofr.eps}
516 chrisfen 1479 \caption{Radial distribution functions of ice $I_h$, $I_c$,
517     Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
518     of the SSD/RF water model at 77 K.}
519 chrisfen 1467 \label{fig:gofr}
520     \end{figure}
521    
522 gezelter 1469 \begin{figure}
523     \includegraphics[width=\linewidth]{sofq.eps}
524 chrisfen 1479 \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
525     and Ice-{\it i}$^\prime$ at 77 K. The raw structure factors have
526     been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
527     width) to compensate for the trunction effects in our finite size
528     simulations.}
529 gezelter 1469 \label{fig:sofq}
530     \end{figure}
531    
532 chrisfen 1453 \section{Acknowledgments}
533     Support for this project was provided by the National Science
534     Foundation under grant CHE-0134881. Computation time was provided by
535 chrisfen 1458 the Notre Dame High Performance Computing Cluster and the Notre Dame
536     Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
537 chrisfen 1453
538     \newpage
539    
540     \bibliographystyle{jcp}
541     \bibliography{iceiPaper}
542    
543    
544     \end{document}