ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceiPaper/iceiPaper.tex
(Generate patch)

Comparing trunk/iceiPaper/iceiPaper.tex (file contents):
Revision 1457 by chrisfen, Tue Sep 14 23:03:53 2004 UTC vs.
Revision 1468 by gezelter, Thu Sep 16 21:15:38 2004 UTC

# Line 1 | Line 1
1   %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 < \documentclass[preprint,aps,endfloats]{revtex4}
2 > \documentclass[11pt]{article}
3   %\documentclass[11pt]{article}
4 < %\usepackage{endfloat}
4 > \usepackage{endfloat}
5   \usepackage{amsmath}
6   \usepackage{epsf}
7   \usepackage{berkeley}
8 < %\usepackage{setspace}
9 < %\usepackage{tabularx}
8 > \usepackage{setspace}
9 > \usepackage{tabularx}
10   \usepackage{graphicx}
11 < %\usepackage[ref]{overcite}
12 < %\pagestyle{plain}
13 < %\pagenumbering{arabic}
14 < %\oddsidemargin 0.0cm \evensidemargin 0.0cm
15 < %\topmargin -21pt \headsep 10pt
16 < %\textheight 9.0in \textwidth 6.5in
17 < %\brokenpenalty=10000
11 > \usepackage[ref]{overcite}
12 > \pagestyle{plain}
13 > \pagenumbering{arabic}
14 > \oddsidemargin 0.0cm \evensidemargin 0.0cm
15 > \topmargin -21pt \headsep 10pt
16 > \textheight 9.0in \textwidth 6.5in
17 > \brokenpenalty=10000
18 > \renewcommand{\baselinestretch}{1.2}
19 > \renewcommand\citemid{\ } % no comma in optional reference note
20  
19 – %\renewcommand\citemid{\ } % no comma in optional reference note
20 –
21   \begin{document}
22  
23 < \title{A Free Energy Study of Low Temperature and Anomolous Ice}
23 > \title{Ice-{\it i}: a simulation-predicted ice polymorph which is more
24 > stable than Ice $I_h$ for point-charge and point-dipole water models}
25  
26 < \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
27 < \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
27 <
28 < \address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
26 > \author{Christopher J. Fennell and J. Daniel Gezelter \\
27 > Department of Chemistry and Biochemistry\\ University of Notre Dame\\
28   Notre Dame, Indiana 46556}
29  
30   \date{\today}
31  
32 < %\maketitle
32 > \maketitle
33   %\doublespacing
34  
35   \begin{abstract}
36 + The free energies of several ice polymorphs in the low pressure regime
37 + were calculated using thermodynamic integration.  These integrations
38 + were done for most of the common water models. Ice-{\it i}, a
39 + structure we recently observed to be stable in one of the single-point
40 + water models, was determined to be the stable crystalline state (at 1
41 + atm) for {\it all} the water models investigated.  Phase diagrams were
42 + generated, and phase coexistence lines were determined for all of the
43 + known low-pressure ice structures under all of the common water
44 + models.  Additionally, potential truncation was shown to have an
45 + effect on the calculated free energies, and can result in altered free
46 + energy landscapes.
47   \end{abstract}
48  
39 – \maketitle
40 –
41 – \newpage
42 –
49   %\narrowtext
50  
51   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 48 | Line 54 | Notre Dame, Indiana 46556}
54  
55   \section{Introduction}
56  
57 + Molecular dynamics is a valuable tool for studying the phase behavior
58 + of systems ranging from small or simple
59 + molecules\cite{Matsumoto02andOthers} to complex biological
60 + species.\cite{bigStuff} Many techniques have been developed to
61 + investigate the thermodynamic properites of model substances,
62 + providing both qualitative and quantitative comparisons between
63 + simulations and experiment.\cite{thermMethods} Investigation of these
64 + properties leads to the development of new and more accurate models,
65 + leading to better understanding and depiction of physical processes
66 + and intricate molecular systems.
67 +
68 + Water has proven to be a challenging substance to depict in
69 + simulations, and a variety of models have been developed to describe
70 + its behavior under varying simulation
71 + conditions.\cite{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
72 + These models have been used to investigate important physical
73 + phenomena like phase transitions and the hydrophobic
74 + effect.\cite{Yamada02} With the choice of models available, it
75 + is only natural to compare the models under interesting thermodynamic
76 + conditions in an attempt to clarify the limitations of each of the
77 + models.\cite{modelProps} Two important property to quantify are the
78 + Gibbs and Helmholtz free energies, particularly for the solid forms of
79 + water.  Difficulty in these types of studies typically arises from the
80 + assortment of possible crystalline polymorphs that water adopts over a
81 + wide range of pressures and temperatures. There are currently 13
82 + recognized forms of ice, and it is a challenging task to investigate
83 + the entire free energy landscape.\cite{Sanz04} Ideally, research is
84 + focused on the phases having the lowest free energy at a given state
85 + point, because these phases will dictate the true transition
86 + temperatures and pressures for their respective model.
87 +
88 + In this paper, standard reference state methods were applied to known
89 + crystalline water polymorphs in the low pressure regime. This work is
90 + unique in the fact that one of the crystal lattices was arrived at
91 + through crystallization of a computationally efficient water model
92 + under constant pressure and temperature conditions. Crystallization
93 + events are interesting in and of
94 + themselves;\cite{Matsumoto02,Yamada02} however, the crystal structure
95 + obtained in this case is different from any previously observed ice
96 + polymorphs in experiment or simulation.\cite{Fennell04} We have named
97 + this structure Ice-{\it i} to indicate its origin in computational
98 + simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
99 + water molecules that stack in rows of interlocking water
100 + tetramers. Proton ordering can be accomplished by orienting two of the
101 + molecules so that both of their donated hydrogen bonds are internal to
102 + their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
103 + constructed of water tetramers, the hydrogen bonds are not as linear
104 + as those observed in ice $I_h$, however the interlocking of these
105 + subunits appears to provide significant stabilization to the overall
106 + crystal. The arrangement of these tetramers results in surrounding
107 + open octagonal cavities that are typically greater than 6.3 \AA\ in
108 + diameter. This relatively open overall structure leads to crystals
109 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
110 +
111 + \begin{figure}
112 + \includegraphics[width=\linewidth]{unitCell.eps}
113 + \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
114 + elongated variant of Ice-{\it i}.  For Ice-{\it i}, the $a$ to $c$
115 + relation is given by $a = 2.1214c$, while for Ice-$i^\prime$, $a =
116 + 1.7850c$.}
117 + \label{iceiCell}
118 + \end{figure}
119 +
120 + \begin{figure}
121 + \includegraphics[width=\linewidth]{orderedIcei.eps}
122 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
123 + down the (001) crystal face. The rows of water tetramers surrounded by
124 + octagonal pores leads to a crystal structure that is significantly
125 + less dense than ice $I_h$.}
126 + \label{protOrder}
127 + \end{figure}
128 +
129 + Results from our previous study indicated that Ice-{\it i} is the
130 + minimum energy crystal structure for the single point water models we
131 + investigated (for discussions on these single point dipole models, see
132 + the previous work and related
133 + articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
134 + considered energetic stabilization and neglected entropic
135 + contributions to the overall free energy. To address this issue, the
136 + absolute free energy of this crystal was calculated using
137 + thermodynamic integration and compared to the free energies of cubic
138 + and hexagonal ice $I$ (the experimental low density ice polymorphs)
139 + and ice B (a higher density, but very stable crystal structure
140 + observed by B\`{a}ez and Clancy in free energy studies of
141 + SPC/E).\cite{Baez95b} This work includes results for the water model
142 + from which Ice-{\it i} was crystallized (SSD/E) in addition to several
143 + common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
144 + field parametrized single point dipole water model (SSD/RF). It should
145 + be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used
146 + in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of
147 + this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit
148 + it is extended in the direction of the (001) face and compressed along
149 + the other two faces.
150 +
151   \section{Methods}
152  
153   Canonical ensemble (NVT) molecular dynamics calculations were
154 < performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
155 < molecular mechanics package. All molecules were treated as rigid
156 < bodies, with orientational motion propogated using the symplectic DLM
157 < integration method. Details about the implementation of these
158 < techniques can be found in a recent publication.\cite{Meineke05}
154 > performed using the OOPSE molecular mechanics package.\cite{Meineke05}
155 > All molecules were treated as rigid bodies, with orientational motion
156 > propagated using the symplectic DLM integration method. Details about
157 > the implementation of these techniques can be found in a recent
158 > publication.\cite{Dullweber1997}
159  
160   Thermodynamic integration was utilized to calculate the free energy of
161   several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
162   SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
163   400 K for all of these water models were also determined using this
164 < same technique, in order to determine melting points and generate
165 < phase diagrams. All simulations were carried out at densities
166 < resulting in a pressure of approximately 1 atm at their respective
67 < temperatures.
164 > same technique in order to determine melting points and generate phase
165 > diagrams. All simulations were carried out at densities resulting in a
166 > pressure of approximately 1 atm at their respective temperatures.
167  
168 + A single thermodynamic integration involves a sequence of simulations
169 + over which the system of interest is converted into a reference system
170 + for which the free energy is known analytically. This transformation
171 + path is then integrated in order to determine the free energy
172 + difference between the two states:
173 + \begin{equation}
174 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
175 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
176 + \end{equation}
177 + where $V$ is the interaction potential and $\lambda$ is the
178 + transformation parameter that scales the overall
179 + potential. Simulations are distributed unevenly along this path in
180 + order to sufficiently sample the regions of greatest change in the
181 + potential. Typical integrations in this study consisted of $\sim$25
182 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
183 + (near the reference state) in length.
184 +
185   For the thermodynamic integration of molecular crystals, the Einstein
186 < Crystal is chosen as the reference state that the system is converted
71 < to over the course of the simulation. In an Einstein Crystal, the
186 > crystal was chosen as the reference state. In an Einstein crystal, the
187   molecules are harmonically restrained at their ideal lattice locations
188   and orientations. The partition function for a molecular crystal
189 < restrained in this fashion has been evaluated, and the Helmholtz Free
190 < Energy ({\it A}) is given by
189 > restrained in this fashion can be evaluated analytically, and the
190 > Helmholtz Free Energy ({\it A}) is given by
191   \begin{eqnarray}
192   A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
193   [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
# Line 92 | Line 207 | state.
207   minimum potential energy of the ideal crystal. In the case of
208   molecular liquids, the ideal vapor is chosen as the target reference
209   state.
210 +
211   \begin{figure}
212 < \includegraphics[scale=1.0]{rotSpring.eps}
212 > \includegraphics[width=\linewidth]{rotSpring.eps}
213   \caption{Possible orientational motions for a restrained molecule.
214   $\theta$ angles correspond to displacement from the body-frame {\it
215   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 104 | Line 220 | cubic switching between 100\% and 85\% of the cutoff v
220   \end{figure}
221  
222   Charge, dipole, and Lennard-Jones interactions were modified by a
223 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
224 < applying this function, these interactions are smoothly truncated,
225 < thereby avoiding poor energy conserving dynamics resulting from
226 < harsher truncation schemes. The effect of a long-range correction was
227 < also investigated on select model systems in a variety of manners. For
228 < the SSD/RF model, a reaction field with a fixed dielectric constant of
229 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
230 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
231 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
232 < compare with the 9 \AA\ cutoff results. Finally, results from the use
233 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
223 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
224 > ). By applying this function, these interactions are smoothly
225 > truncated, thereby avoiding the poor energy conservation which results
226 > from harsher truncation schemes. The effect of a long-range correction
227 > was also investigated on select model systems in a variety of
228 > manners. For the SSD/RF model, a reaction field with a fixed
229 > dielectric constant of 80 was applied in all
230 > simulations.\cite{Onsager36} For a series of the least computationally
231 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
232 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
233 > \AA\ cutoff results. Finally, results from the use of an Ewald
234 > summation were estimated for TIP3P and SPC/E by performing
235   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
236 < mechanics software package. TINKER was chosen because it can also
237 < propogate the motion of rigid-bodies, and provides the most direct
238 < comparison to the results from OOPSE. The calculated energy difference
239 < in the presence and absence of PME was applied to the previous results
123 < in order to predict changes in the free energy landscape.
236 > mechanics software package.\cite{Tinker} The calculated energy
237 > difference in the presence and absence of PME was applied to the
238 > previous results in order to predict changes to the free energy
239 > landscape.
240  
241   \section{Results and discussion}
242  
# Line 130 | Line 246 | Ice XI, the experimentally observed proton ordered var
246   as well as the higher density ice B, observed by B\`{a}ez and Clancy
247   and thought to be the minimum free energy structure for the SPC/E
248   model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
249 < Ice XI, the experimentally observed proton ordered variant of ice
250 < $I_h$, was investigated initially, but it was found not to be as
251 < stable as antiferroelectric variants of proton ordered or even proton
252 < disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
253 < ice $I_h$ used here is a simple antiferroelectric version that has an
254 < 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
255 < ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
256 < $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
249 > Ice XI, the experimentally-observed proton-ordered variant of ice
250 > $I_h$, was investigated initially, but was found to be not as stable
251 > as proton disordered or antiferroelectric variants of ice $I_h$. The
252 > proton ordered variant of ice $I_h$ used here is a simple
253 > antiferroelectric version that has an 8 molecule unit
254 > cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules
255 > for ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for
256 > ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
257   were necessary for simulations involving larger cutoff values.
258  
259   \begin{table*}
# Line 147 | Line 263 | kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF
263   \caption{Calculated free energies for several ice polymorphs with a
264   variety of common water models. All calculations used a cutoff radius
265   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
266 < kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
266 > kcal/mol. Calculated error of the final digits is in parentheses. *Ice
267 > $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
268   \begin{tabular}{ l  c  c  c  c }
269 < \hline \\[-7mm]
270 < \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
271 < \hline \\[-3mm]
272 < \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
273 < \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
274 < \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
275 < \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
276 < \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
277 < \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
269 > \hline
270 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
271 > \hline
272 > TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\
273 > TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\
274 > TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\
275 > SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\
276 > SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\
277 > SSD/RF & -11.51(4) & NA* & -12.08(5) & -12.29(4)\\
278   \end{tabular}
279   \label{freeEnergy}
280   \end{center}
# Line 167 | Line 284 | temperature and pressure dependence of the free energy
284   The free energy values computed for the studied polymorphs indicate
285   that Ice-{\it i} is the most stable state for all of the common water
286   models studied. With the free energy at these state points, the
287 < temperature and pressure dependence of the free energy was used to
288 < project to other state points and build phase diagrams. Figures
287 > Gibbs-Helmholtz equation was used to project to other state points and
288 > to build phase diagrams.  Figures
289   \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
290   from the free energy results. All other models have similar structure,
291 < only the crossing points between these phases exist at different
292 < temperatures and pressures. It is interesting to note that ice $I$
293 < does not exist in either cubic or hexagonal form in any of the phase
294 < diagrams for any of the models. For purposes of this study, ice B is
295 < representative of the dense ice polymorphs. A recent study by Sanz
296 < {\it et al.} goes into detail on the phase diagrams for SPC/E and
297 < TIP4P in the high pressure regime.\cite{Sanz04}
291 > although the crossing points between the phases exist at slightly
292 > different temperatures and pressures. It is interesting to note that
293 > ice $I$ does not exist in either cubic or hexagonal form in any of the
294 > phase diagrams for any of the models. For purposes of this study, ice
295 > B is representative of the dense ice polymorphs. A recent study by
296 > Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
297 > TIP4P in the high pressure regime.\cite{Sanz04}
298 >
299   \begin{figure}
300   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
301   \caption{Phase diagram for the TIP3P water model in the low pressure
# Line 187 | Line 305 | higher in energy and don't appear in the phase diagram
305   higher in energy and don't appear in the phase diagram.}
306   \label{tp3phasedia}
307   \end{figure}
308 +
309   \begin{figure}
310   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
311   \caption{Phase diagram for the SSD/RF water model in the low pressure
# Line 203 | Line 322 | temperatures of several common water models compared w
322   \renewcommand{\thefootnote}{\thempfootnote}
323   \begin{center}
324   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
325 < temperatures of several common water models compared with experiment.}
325 > temperatures at 1 atm for several common water models compared with
326 > experiment. The $T_m$ and $T_s$ values from simulation correspond to a
327 > transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
328 > liquid or gas state.}
329   \begin{tabular}{ l  c  c  c  c  c  c  c }
330 < \hline \\[-7mm]
331 < \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
332 < \hline \\[-3mm]
333 < \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
334 < \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
335 < \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
330 > \hline
331 > Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
332 > \hline
333 > $T_m$ (K)  & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\
334 > $T_b$ (K)  & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\
335 > $T_s$ (K)  & - & - & - & - & 355(3) & - & -\\
336   \end{tabular}
337   \label{meltandboil}
338   \end{center}
# Line 226 | Line 348 | ordered and disordered molecular arrangements). If the
348   studies in the literature. Earlier free energy studies of ice $I$
349   using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
350   being attributed to choice of interaction truncation and different
351 < ordered and disordered molecular arrangements). If the presence of ice
352 < B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
351 > ordered and disordered molecular
352 > arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
353 > Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
354   predicted from this work. However, the $T_m$ from Ice-{\it i} is
355   calculated at 265 K, significantly higher in temperature than the
356   previous studies. Also of interest in these results is that SSD/E does
# Line 235 | Line 358 | advantagious in that it facilitated the spontaneous cr
358   at 355 K. This is due to the significant stability of Ice-{\it i} over
359   all other polymorphs for this particular model under these
360   conditions. While troubling, this behavior turned out to be
361 < advantagious in that it facilitated the spontaneous crystallization of
361 > advantageous in that it facilitated the spontaneous crystallization of
362   Ice-{\it i}. These observations provide a warning that simulations of
363   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
364   risk of spontaneous crystallization. However, this risk changes when
365   applying a longer cutoff.
366  
367 + \begin{figure}
368 + \includegraphics[width=\linewidth]{cutoffChange.eps}
369 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
370 + TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
371 + \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
372 + \AA . These crystals are unstable at 200 K and rapidly convert into
373 + liquids. The connecting lines are qualitative visual aid.}
374 + \label{incCutoff}
375 + \end{figure}
376 +
377   Increasing the cutoff radius in simulations of the more
378   computationally efficient water models was done in order to evaluate
379   the trend in free energy values when moving to systems that do not
380   involve potential truncation. As seen in Fig. \ref{incCutoff}, the
381   free energy of all the ice polymorphs show a substantial dependence on
382   cutoff radius. In general, there is a narrowing of the free energy
383 < differences while moving to greater cutoff radius. This trend is much
384 < more subtle in the case of SSD/RF, indicating that the free energies
385 < calculated with a reaction field present provide a more accurate
386 < picture of the free energy landscape in the absence of potential
387 < truncation.
383 > differences while moving to greater cutoff radius. Interestingly, by
384 > increasing the cutoff radius, the free energy gap was narrowed enough
385 > in the SSD/E model that the liquid state is preferred under standard
386 > simulation conditions (298 K and 1 atm). Thus, it is recommended that
387 > simulations using this model choose interaction truncation radii
388 > greater than 9 \AA\. This narrowing trend is much more subtle in the
389 > case of SSD/RF, indicating that the free energies calculated with a
390 > reaction field present provide a more accurate picture of the free
391 > energy landscape in the absence of potential truncation.
392  
393   To further study the changes resulting to the inclusion of a
394   long-range interaction correction, the effect of an Ewald summation
# Line 263 | Line 400 | cutoff radius is observed in these results. Ice-{\it i
400   SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
401   are not fully supported in TINKER, so the results for these models
402   could not be estimated. The same trend pointed out through increase of
403 < cutoff radius is observed in these results. Ice-{\it i} is the
403 > cutoff radius is observed in these PME results. Ice-{\it i} is the
404   preferred polymorph at ambient conditions for both the TIP3P and SPC/E
405   water models; however, there is a narrowing of the free energy
406   differences between the various solid forms. In the case of SPC/E this
407 < narrowing is significant enough that it becomes less clear cut that
407 > narrowing is significant enough that it becomes less clear that
408   Ice-{\it i} is the most stable polymorph, and is possibly metastable
409   with respect to ice B and possibly ice $I_c$. However, these results
410   do not significantly alter the finding that the Ice-{\it i} polymorph
# Line 278 | Line 415 | the phase behavior of water models.
415   \begin{minipage}{\linewidth}
416   \renewcommand{\thefootnote}{\thempfootnote}
417   \begin{center}
418 < \caption{The free energy of the studied ice polymorphs after applying the energy difference attributed to the inclusion of the PME long-range interaction correction. Units are kcal/mol.}
418 > \caption{The free energy of the studied ice polymorphs after applying
419 > the energy difference attributed to the inclusion of the PME
420 > long-range interaction correction. Units are kcal/mol.}
421   \begin{tabular}{ l  c  c  c  c }
422 < \hline \\[-7mm]
422 > \hline
423   \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
424 < \hline \\[-3mm]
425 < \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
426 < \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
424 > \hline
425 > TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\
426 > SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\
427   \end{tabular}
428   \label{pmeShift}
429   \end{center}
# Line 293 | Line 432 | the phase behavior of water models.
432  
433   \section{Conclusions}
434  
435 + The free energy for proton ordered variants of hexagonal and cubic ice
436 + $I$, ice B, and recently discovered Ice-{\it i} were calculated under
437 + standard conditions for several common water models via thermodynamic
438 + integration. All the water models studied show Ice-{\it i} to be the
439 + minimum free energy crystal structure in the with a 9 \AA\ switching
440 + function cutoff. Calculated melting and boiling points show
441 + surprisingly good agreement with the experimental values; however, the
442 + solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
443 + interaction truncation was investigated through variation of the
444 + cutoff radius, use of a reaction field parameterized model, and
445 + estimation of the results in the presence of the Ewald
446 + summation. Interaction truncation has a significant effect on the
447 + computed free energy values, and may significantly alter the free
448 + energy landscape for the more complex multipoint water models. Despite
449 + these effects, these results show Ice-{\it i} to be an important ice
450 + polymorph that should be considered in simulation studies.
451 +
452 + Due to this relative stability of Ice-{\it i} in all manner of
453 + investigated simulation examples, the question arises as to possible
454 + experimental observation of this polymorph.  The rather extensive past
455 + and current experimental investigation of water in the low pressure
456 + regime makes us hesitant to ascribe any relevance of this work outside
457 + of the simulation community.  It is for this reason that we chose a
458 + name for this polymorph which involves an imaginary quantity.  That
459 + said, there are certain experimental conditions that would provide the
460 + most ideal situation for possible observation. These include the
461 + negative pressure or stretched solid regime, small clusters in vacuum
462 + deposition environments, and in clathrate structures involving small
463 + non-polar molecules.  Fig. \ref{fig:gofr} contains our predictions
464 + of both the pair distribution function ($g_{OO}(r)$) and the structure
465 + factor ($S(\vec{q})$ for this polymorph at a temperature of 77K.  We
466 + will leave it to our experimental colleagues to determine whether this
467 + ice polymorph should really be called Ice-{\it i} or if it should be
468 + promoted to Ice-0.
469 +
470 + \begin{figure}
471 + \includegraphics[width=\linewidth]{iceGofr.eps}
472 + \caption{Radial distribution functions of (A) Ice-{\it i} and (B) ice $I_c$ at 77 K from simulations of the SSD/RF water model.}
473 + \label{fig:gofr}
474 + \end{figure}
475 +
476   \section{Acknowledgments}
477   Support for this project was provided by the National Science
478   Foundation under grant CHE-0134881. Computation time was provided by
479 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
480 < DMR-0079647.
479 > the Notre Dame High Performance Computing Cluster and the Notre Dame
480 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
481  
482   \newpage
483  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines