ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chrisDissertation/Ice.tex
Revision: 2987
Committed: Wed Aug 30 22:36:06 2006 UTC (18 years ago) by chrisfen
Content type: application/x-tex
File size: 26997 byte(s)
Log Message:
Imported sources

File Contents

# User Rev Content
1 chrisfen 2987 \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}
2    
3     As discussed in the previous chapter, water has proven to be a
4     challenging substance to depict in simulations, and a variety of
5     models have been developed to describe its behavior under varying
6     simulation
7     conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04}
8     These models have been used to investigate important physical
9     phenomena like phase transitions and the hydrophobic
10     effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models
11     available, it is only natural to compare them under interesting
12     thermodynamic conditions in an attempt to clarify the limitations of
13     each.\cite{Jorgensen83,Jorgensen98b,Baez94,Mahoney01} Two important
14     property to quantify are the Gibbs and Helmholtz free energies,
15     particularly for the solid forms of water, as these predict the
16     thermodynamic stability of the various phases. Water has a
17     particularly rich phase diagram and takes on a number of different and
18     stable crystalline structures as the temperature and pressure are
19     varied. This complexity makes it a challenging task to investigate the
20     entire free energy landscape.\cite{Sanz04} Ideally, research is
21     focused on the phases having the lowest free energy at a given state
22     point, because these phases will dictate the relevant transition
23     temperatures and pressures for the model.
24    
25     The high-pressure phases of water (ice II-ice X as well as ice XII)
26     have been studied extensively both experimentally and
27     computationally. In this chapter, standard reference state methods
28     were applied in the {\it low} pressure regime to evaluate the free
29     energies for a few known crystalline water polymorphs that might be
30     stable at these pressures. This work is unique in the fact that one of
31     the crystal lattices was arrived at through crystallization of a
32     computationally efficient water model under constant pressure and
33     temperature conditions.
34    
35     While performing a series of melting simulations on an early iteration
36     of SSD/E, we observed several recrystallization events at a constant
37     pressure of 1 atm. After melting from ice I$_\textrm{h}$ at 235~K, two
38     of five systems recrystallized near 245~K. Crystallization events are
39     interesting in and of themselves;\cite{Matsumoto02,Yamada02} however,
40     the crystal structure extracted from these systems is different from
41     any previously observed ice polymorphs in experiment or
42     simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
43     to indicate its origin in computational simulation. The unit cell of
44     Ice-$i$ and an axially elongated variant named Ice-$i^\prime$ both
45     consist of eight water molecules that stack in rows of interlocking
46     water tetramers as illustrated in figure \ref{fig:iceiCell}A,B. These
47     tetramers form a crystal structure similar in appearance to a recent
48     two-dimensional surface tessellation simulated on silica.\cite{Yang04}
49     As expected in an ice crystal constructed of water tetramers, the
50     hydrogen bonds are not as linear as those observed in ice
51     I$_\textrm{h}$; however, the interlocking of these subunits appears to
52     provide significant stabilization to the overall crystal. The
53     arrangement of these tetramers results in open octagonal cavities that
54     are typically greater than 6.3~\AA\ in diameter (see figure
55     \ref{fig:protOrder}). This open structure leads to crystals that are
56     typically 0.07~g/cm$^3$ less dense than ice I$_\textrm{h}$.
57    
58     \begin{figure}
59     \includegraphics[width=\linewidth]{./figures/unitCell.pdf}
60     \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
61     elongated variant of Ice-{\it i}. For Ice-{\it i}, the $a$ to $c$
62     relation is given by $a = 2.1214c$, while for Ice-$i^\prime$, $a =
63     1.7850c$.}
64     \label{fig:iceiCell}
65     \end{figure}
66    
67     \begin{figure}
68     \centering
69     \includegraphics[width=3.5in]{./figures/orderedIcei.pdf}
70     \caption{Image of a proton ordered crystal of Ice-{\it i} looking
71     down the (001) crystal face. The rows of water tetramers surrounded by
72     octagonal pores leads to a crystal structure that is significantly
73     less dense than ice I$_\textrm{h}$.}
74     \label{fig:protOrder}
75     \end{figure}
76    
77     Results from our initial studies indicated that Ice-{\it i} is the
78     minimum energy crystal structure for the single point water models
79     investigated (for discussions on these single point dipole models, see
80     the previous work and related
81     articles\cite{Fennell04,Liu96,Bratko85}). These earlier results only
82     considered energetic stabilization and neglected entropic
83     contributions to the overall free energy. To address this issue, we
84     have calculated the absolute free energy of this crystal using
85     thermodynamic integration and compared to the free energies of ice
86     I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common low-density ice
87     polymorphs) and ice B (a higher density, but very stable crystal
88     structure observed by B\'{a}ez and Clancy in free energy studies of
89     SPC/E).\cite{Baez95b} This work includes results for the water model
90     from which Ice-{\it i} was crystallized (SSD/E) in addition to several
91     common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
92     field parametrized single point dipole water model (SSD/RF). The
93     axially elongated variant, Ice-$i^\prime$, was used in calculations
94     involving SPC/E, TIP4P, and TIP5P. The square tetramer in Ice-$i$
95     distorts in Ice-$i^\prime$ to form a rhombus with alternating 85 and
96     95 degree angles. Under SPC/E, TIP4P, and TIP5P, this geometry is
97     better at forming favorable hydrogen bonds. The degree of rhomboid
98     distortion depends on the water model used but is significant enough
99     to split the peak in the radial distribution function which corresponds
100     to diagonal sites in the tetramers.
101    
102     \section{Methods and Thermodynamic Integration}
103    
104     Canonical ensemble ({\it NVT}) molecular dynamics calculations were
105     performed using the OOPSE molecular mechanics package.\cite{Meineke05}
106     The densities chosen for the simulations were taken from
107     isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and
108     200~K. Each model (and each crystal structure) was allowed to relax for
109     300~ps in the {\it NPT} ensemble before averaging the density to obtain
110     the volumes for the {\it NVT} simulations.All molecules were treated
111     as rigid bodies, with orientational motion propagated using the
112     symplectic DLM integration method described in section
113     \ref{sec:IntroIntegrate}.
114    
115    
116     We used thermodynamic integration to calculate the Helmholtz free
117     energies ({\it A}) of the listed water models at various state
118     points. Thermodynamic integration is an established technique that has
119     been used extensively in the calculation of free energies for
120     condensed phases of
121     materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
122     method uses a sequence of simulations over which the system of
123     interest is converted into a reference system for which the free
124     energy is known analytically ($A_0$). The difference in potential
125     energy between the reference system and the system of interest
126     ($\Delta V$) is then integrated in order to determine the free energy
127     difference between the two states:
128     \begin{equation}
129     A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
130     \end{equation}
131     Here, $\lambda$ is the parameter that governs the transformation
132     between the reference system and the system of interest. For
133     crystalline phases, an harmonically-restrained (Einstein) crystal is
134     chosen as the reference state, while for liquid phases, the ideal gas
135     is taken as the reference state. Figure \ref{fig:integrationPath}
136     shows an example integration path for converting a crystalline system
137     to the Einstein crystal reference state.
138     \begin{figure}
139     \includegraphics[width=\linewidth]{./figures/integrationPath.pdf}
140     \caption{An example integration path to convert an unrestrained
141     crystal ($\lambda = 1$) to the Einstein crystal reference state
142     ($\lambda = 0$). Note the increase in samples at either end of the
143     path to improve the smoothness of the curve. For reversible processes,
144     conversion of the Einstein crystal back to the system of interest will
145     give an identical plot, thereby integrating to the same result.}
146     \label{fig:integrationPath}
147     \end{figure}
148    
149     In an Einstein crystal, the molecules are restrained at their ideal
150     lattice locations and orientations. Using harmonic restraints, as
151     applied by B\'{a}ez and Clancy, the total potential for this reference
152     crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
153     \begin{equation}
154     V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
155     \frac{K_\omega\omega^2}{2},
156     \end{equation}
157     where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
158     the spring constants restraining translational motion and deflection
159     of and rotation around the principle axis of the molecule
160     respectively. These spring constants are typically calculated from
161     the mean-square displacements of water molecules in an unrestrained
162     ice crystal at 200~K. For these studies, $K_\mathrm{v} =
163     4.29$~kcal~mol$^{-1}$~\AA$^{-2}$, $K_\theta\ =
164     13.88$~kcal~mol$^{-1}$~rad$^{-2}$, and $K_\omega\ =
165     17.75$~kcal~mol$^{-1}$~rad$^{-2}$. It is clear from
166     Fig. \ref{fig:waterSpring} that the values of $\theta$ range from $0$
167     to $\pi$, while $\omega$ ranges from $-\pi$ to $\pi$. The partition
168     function for a molecular crystal restrained in this fashion can be
169     evaluated analytically, and the Helmholtz Free Energy ({\it A}) is
170     given by
171     \begin{equation}
172     \begin{split}
173     A = E_m &- kT\ln\left(\frac{kT}{h\nu}\right)^3 \\
174     &- kT\ln\left[\pi^\frac{1}{2}\left(
175     \frac{8\pi^2I_\mathrm{A}kT}{h^2}\right)^\frac{1}{2}
176     \left(\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right)^\frac{1}{2}
177     \left(\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right)^\frac{1}{2}
178     \right] \\
179     &- kT\ln\left[\frac{kT}{2(\pi K_\omega K_\theta)^{\frac{1}{2}}}
180     \exp\left(-\frac{kT}{2K_\theta}\right)
181     \int_0^{\left(\frac{kT}{2K_\theta}\right)^\frac{1}{2}}
182     \exp(t^2)\mathrm{d}t\right],
183     \end{split}
184     \label{eq:ecFreeEnergy}
185     \end{equation}
186     where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
187     potential energy of the ideal crystal.\cite{Baez95a} The choice of an
188     Einstein crystal reference state is somewhat arbitrary. Any ideal
189     system for which the partition function is known exactly could be used
190     as a reference point as long as the system does not undergo a phase
191     transition during the integration path between the real and ideal
192     systems. Nada and van der Eerden have shown that the use of different
193     force constants in the Einstein crystal does not affect the total
194     free energy, and Gao {\it et al.} have shown that free energies
195     computed with the Debye crystal reference state differ from the
196     Einstein crystal by only a few tenths of a
197     kJ~mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can
198     lead to some uncertainty in the computed melting point of the solids.
199     \begin{figure}
200     \centering
201     \includegraphics[width=3.5in]{./figures/rotSpring.pdf}
202     \caption{Possible orientational motions for a restrained molecule.
203     $\theta$ angles correspond to displacement from the body-frame {\it
204     z}-axis, while $\omega$ angles correspond to rotation about the
205     body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
206     constants for the harmonic springs restraining motion in the $\theta$
207     and $\omega$ directions.}
208     \label{fig:waterSpring}
209     \end{figure}
210    
211     In the case of molecular liquids, the ideal vapor is chosen as the
212     target reference state. There are several examples of liquid state
213     free energy calculations of water models present in the
214     literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
215     typically differ in regard to the path taken for switching off the
216     interaction potential to convert the system to an ideal gas of water
217     molecules. In this study, we applied one of the most convenient
218     methods and integrated over the $\lambda^4$ path, where all
219     interaction parameters are scaled equally by this transformation
220     parameter. This method has been shown to be reversible and provide
221     results in excellent agreement with other established
222     methods.\cite{Baez95b}
223    
224     Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
225     Lennard-Jones interactions were gradually reduced by a cubic switching
226     function. By applying this function, these interactions are smoothly
227     truncated, thereby avoiding the poor energy conservation which results
228     from harsher truncation schemes. The effect of a long-range
229     correction was also investigated on select model systems in a variety
230     of manners. For the SSD/RF model, a reaction field with a fixed
231     dielectric constant of 80 was applied in all
232     simulations.\cite{Onsager36} For a series of the least computationally
233     expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
234     performed with longer cutoffs of 10.5, 12, 13.5, and 15~\AA\ to
235     compare with the 9~\AA\ cutoff results. Finally, the effects of using
236     the Ewald summation were estimated for TIP3P and SPC/E by performing
237     single configuration Particle-Mesh Ewald (PME) calculations for each
238     of the ice polymorphs.\cite{Ponder87} The calculated energy difference
239     in the presence and absence of PME was applied to the previous results
240     in order to predict changes to the free energy landscape.
241    
242     \section{Initial Free Energy Results}
243    
244     The calculated free energies of proton-ordered variants of three low
245     density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it i} or
246     Ice-$i^\prime$) and the stable higher density ice B are listed in
247     table \ref{tab:freeEnergy}. Ice B was included because it has been
248     shown to be a minimum free energy structure for SPC/E at ambient
249     conditions.\cite{Baez95b} In addition to the free energies, the
250     relevant transition temperatures at standard pressure are also
251     displayed in table \ref{tab:freeEnergy}. These free energy values
252     indicate that Ice-{\it i} is the most stable state for all of the
253     investigated water models. With the free energy at these state
254     points, the Gibbs-Helmholtz equation was used to project to other
255     state points and to build phase diagrams. Figures
256     \ref{fig:tp3PhaseDia} and \ref{fig:ssdrfPhaseDia} are example diagrams
257     built from the results for the TIP3P and SSD/RF water models. All
258     other models have similar structure, although the crossing points
259     between the phases move to different temperatures and pressures as
260     indicated from the transition temperatures in table
261     \ref{tab:freeEnergy}. It is interesting to note that ice
262     I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear
263     in any of the phase diagrams for any of the models. For purposes of
264     this study, ice B is representative of the dense ice polymorphs. A
265     recent study by Sanz {\it et al.} provides details on the phase
266     diagrams for SPC/E and TIP4P at higher pressures than those studied
267     here.\cite{Sanz04}
268     \begin{table}
269     \centering
270     \caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1
271     ATMOSPHERE FOR SEVERAL WATER MODELS}
272    
273     \footnotesize
274     \begin{tabular}{lccccccc}
275     \toprule
276     \toprule
277     Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\
278     \cmidrule(lr){2-6}
279     \cmidrule(l){7-8}
280     & \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\
281     \midrule
282     TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\
283     TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\
284     TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\
285     SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\
286     SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\
287     SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\
288     \bottomrule
289     \end{tabular}
290     \label{tab:freeEnergy}
291     \end{table}
292    
293     \begin{figure}
294     \centering
295     \includegraphics[width=\linewidth]{./figures/tp3PhaseDia.pdf}
296     \caption{Phase diagram for the TIP3P water model in the low pressure
297     regime. The displayed $T_m$ and $T_b$ values are good predictions of
298     the experimental values; however, the solid phases shown are not the
299     experimentally observed forms. Both cubic and hexagonal ice $I$ are
300     higher in energy and don't appear in the phase diagram.}
301     \label{fig:tp3PhaseDia}
302     \end{figure}
303    
304     \begin{figure}
305     \centering
306     \includegraphics[width=\linewidth]{./figures/ssdrfPhaseDia.pdf}
307     \caption{Phase diagram for the SSD/RF water model in the low pressure
308     regime. Calculations producing these results were done under an
309     applied reaction field. It is interesting to note that this
310     computationally efficient model (over 3 times more efficient than
311     TIP3P) exhibits phase behavior similar to the less computationally
312     conservative charge based models.}
313     \label{fig:ssdrfPhaseDia}
314     \end{figure}
315    
316     We note that all of the crystals investigated in this study are ideal
317     proton-ordered antiferroelectric structures. All of the structures
318     obey the Bernal-Fowler rules and should be able to form stable
319     proton-{\it disordered} crystals which have the traditional
320     $k_\textrm{B}$ln(3/2) residual entropy at 0~K.\cite{Bernal33,Pauling35}
321     Simulations of proton-disordered structures are relatively unstable
322     with all but the most expensive water models.\cite{Nada03} Our
323     simulations have therefore been performed with the ordered
324     antiferroelectric structures which do not require the residual entropy
325     term to be accounted for in the free energies. This may result in some
326     discrepancies when comparing our melting temperatures to the melting
327     temperatures that have been calculated via thermodynamic integrations
328     of the disordered structures.\cite{Sanz04}
329    
330     Most of the water models have melting points that compare quite
331     favorably with the experimental value of 273~K. The unfortunate
332     aspect of this result is that this phase change occurs between
333     Ice-{\it i} and the liquid state rather than ice I$_h$ and the liquid
334     state. These results do not contradict other studies. Studies of ice
335     I$_h$ using TIP4P predict a $T_m$ ranging from 191 to 238~K
336     (differences being attributed to choice of interaction truncation and
337     different ordered and disordered molecular
338     arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice
339     B and Ice-{\it i} were omitted, a $T_\textrm{m}$ value around 200~K
340     would be predicted from this work. However, the $T_\textrm{m}$ from
341     Ice-{\it i} is calculated to be 262~K, indicating that these
342     simulation based structures ought to be included in studies probing
343     phase transitions with this model. Also of interest in these results
344     is that SSD/E does not exhibit a melting point at 1 atm but does
345     sublime at 355~K. This is due to the significant stability of
346     Ice-{\it i} over all other polymorphs for this particular model under
347     these conditions. While troubling, this behavior resulted in the
348     spontaneous crystallization of Ice-{\it i} which led us to investigate
349     this structure. These observations provide a warning that simulations
350     of SSD/E as a ``liquid'' near 300~K are actually metastable and run
351     the risk of spontaneous crystallization. However, when a longer
352     cutoff radius is used, SSD/E prefers the liquid state under standard
353     temperature and pressure.
354    
355     \section{Effects of Potential Truncation}
356    
357     \begin{figure}
358     \includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
359     \caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
360     SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
361     with an added Ewald correction term. Error for the larger cutoff
362     points is equivalent to that observed at 9.0~\AA\ (see Table
363     \ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using
364     both 12 and 13.5~\AA\ cutoffs were omitted because the crystal was
365     prone to distortion and melting at 200~K. Ice-$i^\prime$ is the
366     form of Ice-{\it i} used in the SPC/E simulations.}
367     \label{fig:incCutoff}
368     \end{figure}
369    
370     For the more computationally efficient water models, we have also
371     investigated the effect of potential truncation on the computed free
372     energies as a function of the cutoff radius. As seen in
373     Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with
374     water models lacking a long-range correction show significant cutoff
375     dependence. In general, there is a narrowing of the free energy
376     differences while moving to greater cutoff radii. As the free
377     energies for the polymorphs converge, the stability advantage that
378     Ice-{\it i} exhibits is reduced. Adjacent to each of these plots are
379     results for systems with applied or estimated long-range corrections.
380     SSD/RF was parametrized for use with a reaction field, and the benefit
381     provided by this computationally inexpensive correction is apparent.
382     The free energies are largely independent of the size of the reaction
383     field cavity in this model, so small cutoff radii mimic bulk
384     calculations quite well under SSD/RF.
385    
386     Although TIP3P was parametrized for use without the Ewald summation,
387     we have estimated the effect of this method for computing long-range
388     electrostatics for both TIP3P and SPC/E. This was accomplished by
389     calculating the potential energy of identical crystals both with and
390     without particle mesh Ewald (PME). Similar behavior to that observed
391     with reaction field is seen for both of these models. The free
392     energies show reduced dependence on cutoff radius and span a narrower
393     range for the various polymorphs. Like the dipolar water models,
394     TIP3P displays a relatively constant preference for the Ice-{\it i}
395     polymorph. Crystal preference is much more difficult to determine for
396     SPC/E. Without a long-range correction, each of the polymorphs
397     studied assumes the role of the preferred polymorph under different
398     cutoff radii. The inclusion of the Ewald correction flattens and
399     narrows the gap in free energies such that the polymorphs are
400     isoenergetic within statistical uncertainty. This suggests that other
401     conditions, such as the density in fixed-volume simulations, can
402     influence the polymorph expressed upon crystallization.
403    
404     \section{Expanded Results Using Damped Shifted Force Electrostatics}
405    
406     In chapter \ref{chap:electrostatics}, we discussed in detail a
407     pairwise method for handling electrostatics (shifted force, {\sc sf})
408     that can be used as a simple and efficient replacement for the Ewald
409     summation. Answering the question of the free energies of these ice
410     polymorphs with varying water models would be an interesting
411     application of this technique. To this end, we set up thermodynamic
412     integrations of all of the previously discussed ice polymorphs using
413     the {\sc sf} technique with a cutoff radius of 12~\AA\ and an $\alpha$
414     of 0.2125~\AA . These calculations were performed on TIP5P-E and
415     TIP4P-Ew (variants of the root models optimized for the Ewald
416     summation) as well as SPC/E, SSD/RF, and TRED (see section
417     \ref{sec:tredWater}).
418    
419     \begin{table}
420     \centering
421     \caption{HELMHOLTZ FREE ENERGIES OF ICE POLYMORPHS USING THE DAMPED
422     SHIFTED FORCE CORRECTION}
423     \begin{tabular}{ lccccc }
424     \toprule
425     \toprule
426     Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\
427     \cmidrule(lr){2-6}
428     & \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\
429     \midrule
430     TIP5P-E & -10.76(4) & -10.72(4) & & - & -10.68(4) \\
431     TIP4P-Ew & & -11.77(3) & & - & -11.60(3) \\
432     SPC/E & -12.98(3) & -11.60(3) & & - & -12.93(3) \\
433     SSD/RF & -11.81(4) & -11.65(3) & & -12.41(4) & - \\
434     TRED & -12.58(3) & -12.44(3) & & -13.09(4) & - \\
435     \end{tabular}
436     \label{tab:dampedFreeEnergy}
437     \end{table}
438    
439    
440     \section{Conclusions}
441    
442     In this work, thermodynamic integration was used to determine the
443     absolute free energies of several ice polymorphs. The new polymorph,
444     Ice-{\it i} was observed to be the stable crystalline state for {\it
445     all} the water models when using a 9.0~\AA\ cutoff. However, the free
446     energy partially depends on simulation conditions (particularly on the
447     choice of long range correction method). Regardless, Ice-{\it i} was
448     still observed to be a stable polymorph for all of the studied water
449     models.
450    
451     So what is the preferred solid polymorph for simulated water? As
452     indicated above, the answer appears to be dependent both on the
453     conditions and the model used. In the case of short cutoffs without a
454     long-range interaction correction, Ice-{\it i} and Ice-$i^\prime$ have
455     the lowest free energy of the studied polymorphs with all the models.
456     Ideally, crystallization of each model under constant pressure
457     conditions, as was done with SSD/E, would aid in the identification of
458     their respective preferred structures. This work, however, helps
459     illustrate how studies involving one specific model can lead to
460     insight about important behavior of others.
461    
462     We also note that none of the water models used in this study are
463     polarizable or flexible models. It is entirely possible that the
464     polarizability of real water makes Ice-{\it i} substantially less
465     stable than ice I$_h$. However, the calculations presented above seem
466     interesting enough to communicate before the role of polarizability
467     (or flexibility) has been thoroughly investigated.
468    
469     Finally, due to the stability of Ice-{\it i} in the investigated
470     simulation conditions, the question arises as to possible experimental
471     observation of this polymorph. The rather extensive past and current
472     experimental investigation of water in the low pressure regime makes
473     us hesitant to ascribe any relevance to this work outside of the
474     simulation community. It is for this reason that we chose a name for
475     this polymorph which involves an imaginary quantity. That said, there
476     are certain experimental conditions that would provide the most ideal
477     situation for possible observation. These include the negative
478     pressure or stretched solid regime, small clusters in vacuum
479     deposition environments, and in clathrate structures involving small
480     non-polar molecules. For the purpose of comparison with experimental
481     results, we have calculated the oxygen-oxygen pair correlation
482     function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$
483     for the two Ice-{\it i} variants (along with example ice
484     I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77~K, and they are shown in
485     figures \ref{fig:gofr} and \ref{fig:sofq} respectively. It is
486     interesting to note that the structure factors for Ice-$i^\prime$ and
487     Ice-I$_c$ are quite similar. The primary differences are small peaks
488     at 1.125, 2.29, and 2.53~\AA$^{-1}$, so particular attention to these
489     regions would be needed to identify the new $i^\prime$ variant from
490     the I$_\textrm{c}$ polymorph.
491    
492    
493     \begin{figure}
494     \includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
495     \caption{Radial distribution functions of Ice-{\it i} and ice
496     I$_\textrm{c}$ calculated from from simulations of the SSD/RF water
497     model at 77~K.}
498     \label{fig:gofr}
499     \end{figure}
500    
501     \begin{figure}
502     \includegraphics[width=\linewidth]{./figures/sofq.pdf}
503     \caption{Predicted structure factors for Ice-{\it i} and ice
504     I$_\textrm{c}$ at 77~K. The raw structure factors have been
505     convoluted with a gaussian instrument function (0.075~\AA$^{-1}$
506     width) to compensate for the truncation effects in our finite size
507     simulations. The labeled peaks compared favorably with ``spurious''
508     peaks observed in experimental studies of amorphous solid
509     water.\cite{Bizid87}}
510     \label{fig:sofq}
511     \end{figure}
512