ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/fennellDissertation/iceChapter.tex
Revision: 2986
Committed: Wed Aug 30 22:14:37 2006 UTC (18 years, 8 months ago) by chrisfen
Content type: application/x-tex
File size: 26997 byte(s)
Log Message:
I may have messed some stuff up...

File Contents

# Content
1 \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