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 1555 by chrisfen, Mon Oct 11 20:14:35 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 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 + generated, and phase coexistence lines were determined for all of the
45 + 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 + 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   \end{abstract}
52  
39 – \maketitle
40 –
41 – \newpage
42 –
53   %\narrowtext
54  
55   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 48 | Line 58 | Notre Dame, Indiana 46556}
58  
59   \section{Introduction}
60  
61 + Water has proven to be a challenging substance to depict in
62 + simulations, and a variety of models have been developed to describe
63 + its behavior under varying simulation
64 + conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
65 + These models have been used to investigate important physical
66 + phenomena like phase transitions, transport properties, and the
67 + 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 + important properties to quantify are the Gibbs and Helmholtz free
73 + 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 + these phases will dictate the relevant transition temperatures and
81 + pressures for the model.
82 +
83 + In this paper, standard reference state methods were applied to known
84 + crystalline water polymorphs in the low pressure regime.  This work is
85 + 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 + 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 +
106 + \begin{figure}
107 + \includegraphics[width=\linewidth]{unitCell.eps}
108 + \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 + 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 + \label{iceiCell}
114 + \end{figure}
115 +
116 + \begin{figure}
117 + \includegraphics[width=\linewidth]{orderedIcei.eps}
118 + \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 +
125 + 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 + had investigated (for discussions on these single point dipole models,
128 + see our previous work and related
129 + articles).\cite{Fennell04,Liu96,Bratko85} Those results only
130 + considered energetic stabilization and neglected entropic
131 + contributions to the overall free energy. To address this issue, we
132 + have calculated the absolute free energy of this crystal using
133 + 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 + 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 + compressed along the other two faces.  There is typically a small
146 + 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 + location peak in the radial distribution function.
151 +
152   \section{Methods}
153  
154   Canonical ensemble (NVT) molecular dynamics calculations were
155 < performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
156 < molecular mechanics package. All molecules were treated as rigid
157 < bodies, with orientational motion propogated using the symplectic DLM
158 < integration method. Details about the implementation of these
159 < techniques can be found in a recent publication.\cite{Meineke05}
155 > 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 > the implementation of this technique can be found in a recent
159 > publication.\cite{Dullweber1997}
160  
161 < Thermodynamic integration was utilized to calculate the free energy of
162 < several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
163 < SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
164 < 400 K for all of these water models were also determined using this
165 < same technique, in order to determine melting points and generate
166 < phase diagrams. All simulations were carried out at densities
167 < resulting in a pressure of approximately 1 atm at their respective
168 < temperatures.
161 > 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 > 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  
174 + 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 + \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 + transformation parameter that scales the overall
185 + potential. Simulations are distributed strategically along this path
186 + in order to sufficiently sample the regions of greatest change in the
187 + 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 +
191   For the thermodynamic integration of molecular crystals, the Einstein
192 < Crystal is chosen as the reference state that the system is converted
193 < to over the course of the simulation. In an Einstein Crystal, the
194 < molecules are harmonically restrained at their ideal lattice locations
195 < and orientations. The partition function for a molecular crystal
196 < restrained in this fashion has been evaluated, and the Helmholtz Free
197 < Energy ({\it A}) is given by
192 > 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.  These spring constants are typically calculated from
205 > the mean-square displacements of water molecules in an unrestrained
206 > ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
207 > mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
208 > 17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
209 > the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
210 > from $-\pi$ to $\pi$.  The partition function for a molecular crystal
211 > restrained in this fashion can be evaluated analytically, and the
212 > Helmholtz Free Energy ({\it A}) is given by
213   \begin{eqnarray}
214   A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
215   [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
# Line 84 | Line 221 | where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a
221   )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
222   \label{ecFreeEnergy}
223   \end{eqnarray}
224 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
225 < \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
226 < $K_\mathrm{\omega}$ are the spring constants restraining translational
90 < motion and deflection of and rotation around the principle axis of the
91 < molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
92 < minimum potential energy of the ideal crystal. In the case of
93 < molecular liquids, the ideal vapor is chosen as the target reference
94 < state.
224 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
225 > potential energy of the ideal crystal.\cite{Baez95a}
226 >
227   \begin{figure}
228 < \includegraphics[scale=1.0]{rotSpring.eps}
228 > \includegraphics[width=\linewidth]{rotSpring.eps}
229   \caption{Possible orientational motions for a restrained molecule.
230   $\theta$ angles correspond to displacement from the body-frame {\it
231   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 103 | Line 235 | Charge, dipole, and Lennard-Jones interactions were mo
235   \label{waterSpring}
236   \end{figure}
237  
238 + In the case of molecular liquids, the ideal vapor is chosen as the
239 + target reference state.  There are several examples of liquid state
240 + free energy calculations of water models present in the
241 + literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
242 + typically differ in regard to the path taken for switching off the
243 + interaction potential to convert the system to an ideal gas of water
244 + molecules.  In this study, we applied of one of the most convenient
245 + methods and integrated over the $\lambda^4$ path, where all
246 + interaction parameters are scaled equally by this transformation
247 + parameter.  This method has been shown to be reversible and provide
248 + results in excellent agreement with other established
249 + methods.\cite{Baez95b}
250 +
251   Charge, dipole, and Lennard-Jones interactions were modified by a
252 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
253 < applying this function, these interactions are smoothly truncated,
254 < thereby avoiding poor energy conserving dynamics resulting from
255 < harsher truncation schemes. The effect of a long-range correction was
256 < also investigated on select model systems in a variety of manners. For
257 < the SSD/RF model, a reaction field with a fixed dielectric constant of
258 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
259 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
260 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
261 < compare with the 9 \AA\ cutoff results. Finally, results from the use
262 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
263 < calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
264 < mechanics software package. TINKER was chosen because it can also
265 < propogate the motion of rigid-bodies, and provides the most direct
266 < comparison to the results from OOPSE. The calculated energy difference
267 < in the presence and absence of PME was applied to the previous results
268 < in order to predict changes in the free energy landscape.
252 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
253 > ). By applying this function, these interactions are smoothly
254 > truncated, thereby avoiding the poor energy conservation which results
255 > from harsher truncation schemes. The effect of a long-range correction
256 > was also investigated on select model systems in a variety of
257 > manners. For the SSD/RF model, a reaction field with a fixed
258 > dielectric constant of 80 was applied in all
259 > simulations.\cite{Onsager36} For a series of the least computationally
260 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
261 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
262 > \AA\ cutoff results. Finally, the effects of utilizing an Ewald
263 > summation were estimated for TIP3P and SPC/E by performing single
264 > configuration calculations with Particle-Mesh Ewald (PME) in the
265 > TINKER molecular mechanics software package.\cite{Tinker} The
266 > calculated energy difference in the presence and absence of PME was
267 > applied to the previous results in order to predict changes to the
268 > free energy landscape.
269  
270   \section{Results and discussion}
271  
272 < The free energy of proton ordered Ice-{\it i} was calculated and
272 > The free energy of proton-ordered Ice-{\it i} was calculated and
273   compared with the free energies of proton ordered variants of the
274   experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
275   as well as the higher density ice B, observed by B\`{a}ez and Clancy
276   and thought to be the minimum free energy structure for the SPC/E
277   model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
278 < Ice XI, the experimentally observed proton ordered variant of ice
279 < $I_h$, was investigated initially, but it was found not to be as
280 < stable as antiferroelectric variants of proton ordered or even proton
281 < disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
282 < ice $I_h$ used here is a simple antiferroelectric version that has an
283 < 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
284 < ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
285 < $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
286 < were necessary for simulations involving larger cutoff values.
278 > Ice XI, the experimentally-observed proton-ordered variant of ice
279 > $I_h$, was investigated initially, but was found to be not as stable
280 > as proton disordered or antiferroelectric variants of ice $I_h$. The
281 > proton ordered variant of ice $I_h$ used here is a simple
282 > antiferroelectric version that we devised, and it has an 8 molecule
283 > unit cell similar to other predicted antiferroelectric $I_h$
284 > crystals.\cite{Davidson84} The crystals contained 648 or 1728
285 > molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
286 > molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
287 > crystal sizes were necessary for simulations involving larger cutoff
288 > values.
289  
290   \begin{table*}
291   \begin{minipage}{\linewidth}
145 – \renewcommand{\thefootnote}{\thempfootnote}
292   \begin{center}
293 +
294   \caption{Calculated free energies for several ice polymorphs with a
295   variety of common water models. All calculations used a cutoff radius
296   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
297 < kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
298 < \begin{tabular}{ l  c  c  c  c }
299 < \hline \\[-7mm]
300 < \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
301 < \hline \\[-3mm]
302 < \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
303 < \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
304 < \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
305 < \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
306 < \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
307 < \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
297 > kcal/mol. Calculated error of the final digits is in parentheses.}
298 >
299 > \begin{tabular}{lcccc}
300 > \hline
301 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
302 > \hline
303 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
304 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
305 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
306 > SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & -13.55(2)\\
307 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
308 > SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\
309   \end{tabular}
310   \label{freeEnergy}
311   \end{center}
# Line 166 | Line 314 | models studied. With the free energy at these state po
314  
315   The free energy values computed for the studied polymorphs indicate
316   that Ice-{\it i} is the most stable state for all of the common water
317 < models studied. With the free energy at these state points, the
318 < temperature and pressure dependence of the free energy was used to
319 < project to other state points and build phase diagrams. Figures
317 > models studied. With the calculated free energy at these state points,
318 > the Gibbs-Helmholtz equation was used to project to other state points
319 > and to build phase diagrams.  Figures
320   \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
321   from the free energy results. All other models have similar structure,
322 < only the crossing points between these phases exist at different
323 < temperatures and pressures. It is interesting to note that ice $I$
324 < does not exist in either cubic or hexagonal form in any of the phase
325 < diagrams for any of the models. For purposes of this study, ice B is
326 < representative of the dense ice polymorphs. A recent study by Sanz
327 < {\it et al.} goes into detail on the phase diagrams for SPC/E and
328 < TIP4P in the high pressure regime.\cite{Sanz04}
322 > although the crossing points between the phases move to slightly
323 > different temperatures and pressures. It is interesting to note that
324 > ice $I$ does not exist in either cubic or hexagonal form in any of the
325 > phase diagrams for any of the models. For purposes of this study, ice
326 > B is representative of the dense ice polymorphs. A recent study by
327 > Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
328 > TIP4P at higher pressures than those studied here.\cite{Sanz04}
329 >
330   \begin{figure}
331   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
332   \caption{Phase diagram for the TIP3P water model in the low pressure
# Line 187 | Line 336 | higher in energy and don't appear in the phase diagram
336   higher in energy and don't appear in the phase diagram.}
337   \label{tp3phasedia}
338   \end{figure}
339 +
340   \begin{figure}
341   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
342   \caption{Phase diagram for the SSD/RF water model in the low pressure
# Line 200 | Line 350 | conservative charge based models.}
350  
351   \begin{table*}
352   \begin{minipage}{\linewidth}
203 – \renewcommand{\thefootnote}{\thempfootnote}
353   \begin{center}
354 +
355   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
356 < temperatures of several common water models compared with experiment.}
357 < \begin{tabular}{ l  c  c  c  c  c  c  c }
358 < \hline \\[-7mm]
359 < \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
360 < \hline \\[-3mm]
361 < \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
362 < \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
363 < \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
356 > temperatures at 1 atm for several common water models compared with
357 > experiment. The $T_m$ and $T_s$ values from simulation correspond to a
358 > transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
359 > liquid or gas state.}
360 >
361 > \begin{tabular}{lccccccc}
362 > \hline
363 > Equilibrium Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
364 > \hline
365 > $T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
366 > $T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
367 > $T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
368   \end{tabular}
369   \label{meltandboil}
370   \end{center}
# Line 226 | Line 380 | ordered and disordered molecular arrangements). If the
380   studies in the literature. Earlier free energy studies of ice $I$
381   using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
382   being attributed to choice of interaction truncation and different
383 < ordered and disordered molecular arrangements). If the presence of ice
384 < B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
383 > ordered and disordered molecular
384 > arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
385 > Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
386   predicted from this work. However, the $T_m$ from Ice-{\it i} is
387   calculated at 265 K, significantly higher in temperature than the
388   previous studies. Also of interest in these results is that SSD/E does
389   not exhibit a melting point at 1 atm, but it shows a sublimation point
390   at 355 K. This is due to the significant stability of Ice-{\it i} over
391   all other polymorphs for this particular model under these
392 < conditions. While troubling, this behavior turned out to be
393 < advantagious in that it facilitated the spontaneous crystallization of
394 < Ice-{\it i}. These observations provide a warning that simulations of
392 > conditions. While troubling, this behavior resulted in spontaneous
393 > crystallization of Ice-{\it i} and led us to investigate this
394 > structure. These observations provide a warning that simulations of
395   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
396 < risk of spontaneous crystallization. However, this risk changes when
396 > risk of spontaneous crystallization. However, this risk lessens when
397   applying a longer cutoff.
398  
399 + \begin{figure}
400 + \includegraphics[width=\linewidth]{cutoffChange.eps}
401 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
402 + TIP3P, and (C) SSD/RF with a reaction field. Both SSD/E and TIP3P show
403 + significant cutoff radius dependence of the free energy and appear to
404 + converge when moving to cutoffs greater than 12 \AA. Use of a reaction
405 + field with SSD/RF results in free energies that exhibit minimal cutoff
406 + radius dependence.}
407 + \label{incCutoff}
408 + \end{figure}
409 +
410   Increasing the cutoff radius in simulations of the more
411   computationally efficient water models was done in order to evaluate
412   the trend in free energy values when moving to systems that do not
413   involve potential truncation. As seen in Fig. \ref{incCutoff}, the
414 < free energy of all the ice polymorphs show a substantial dependence on
415 < cutoff radius. In general, there is a narrowing of the free energy
416 < differences while moving to greater cutoff radius. This trend is much
417 < more subtle in the case of SSD/RF, indicating that the free energies
418 < calculated with a reaction field present provide a more accurate
419 < picture of the free energy landscape in the absence of potential
420 < truncation.
414 > free energy of all the ice polymorphs for the SSD/E and TIP3P models
415 > show a substantial dependence on cutoff radius. In general, there is a
416 > narrowing of the free energy differences while moving to greater
417 > cutoff radii.  As the free energies for the polymorphs converge, the
418 > stability advantage that Ice-{\it i} exhibits is reduced; however, it
419 > remains the most stable polymorph for both of these models over the
420 > depicted range for both models. This narrowing trend is not
421 > significant in the case of SSD/RF, indicating that the free energies
422 > calculated with a reaction field present provide, at minimal
423 > computational cost, a more accurate picture of the free energy
424 > landscape in the absence of potential truncation.  Interestingly,
425 > increasing the cutoff radius a mere 1.5 \AA\ with the SSD/E model
426 > destabilizes the Ice-{\it i} polymorph enough that the liquid state is
427 > preferred under standard simulation conditions (298 K and 1
428 > atm). Thus, it is recommended that simulations using this model choose
429 > interaction truncation radii greater than 9 \AA. Considering this
430 > stabilization provided by smaller cutoffs, it is not surprising that
431 > crystallization into Ice-{\it i} was observed with SSD/E.  The choice
432 > of a 9 \AA\ cutoff in the previous simulations gives the Ice-{\it i}
433 > polymorph a greater than 1 kcal/mol lower free energy than the ice
434 > $I_\textrm{h}$ starting configurations.
435  
436   To further study the changes resulting to the inclusion of a
437   long-range interaction correction, the effect of an Ewald summation
438   was estimated by applying the potential energy difference do to its
439 < inclusion in systems in the presence and absence of the
440 < correction. This was accomplished by calculation of the potential
441 < energy of identical crystals with and without PME using TINKER. The
442 < free energies for the investigated polymorphs using the TIP3P and
443 < SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
444 < are not fully supported in TINKER, so the results for these models
445 < could not be estimated. The same trend pointed out through increase of
446 < cutoff radius is observed in these results. Ice-{\it i} is the
447 < preferred polymorph at ambient conditions for both the TIP3P and SPC/E
448 < water models; however, there is a narrowing of the free energy
449 < differences between the various solid forms. In the case of SPC/E this
450 < narrowing is significant enough that it becomes less clear cut that
451 < Ice-{\it i} is the most stable polymorph, and is possibly metastable
452 < with respect to ice B and possibly ice $I_c$. However, these results
453 < do not significantly alter the finding that the Ice-{\it i} polymorph
454 < is a stable crystal structure that should be considered when studying
455 < the phase behavior of water models.
439 > inclusion in systems in the presence and absence of the correction.
440 > This was accomplished by calculation of the potential energy of
441 > identical crystals both with and without PME.  The free energies for
442 > the investigated polymorphs using the TIP3P and SPC/E water models are
443 > shown in Table \ref{pmeShift}.  The same trend pointed out through
444 > increase of cutoff radius is observed in these PME results. Ice-{\it
445 > i} is the preferred polymorph at ambient conditions for both the TIP3P
446 > and SPC/E water models; however, the narrowing of the free energy
447 > differences between the various solid forms with the SPC/E model is
448 > significant enough that it becomes less clear that it is the most
449 > stable polymorph.  The free energies of Ice-{\it i} and $I_\textrm{c}$
450 > overlap within error, while ice B and $I_\textrm{h}$ are just outside
451 > at t slightly higher free energy.  This indicates that with SPC/E,
452 > Ice-{\it i} might be metastable with all the studied polymorphs,
453 > particularly ice $I_\textrm{c}$. However, these results do not
454 > significantly alter the finding that the Ice-{\it i} polymorph is a
455 > stable crystal structure that should be considered when studying the
456 > phase behavior of water models.
457  
458   \begin{table*}
459   \begin{minipage}{\linewidth}
279 – \renewcommand{\thefootnote}{\thempfootnote}
460   \begin{center}
461 < \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.}
462 < \begin{tabular}{ l  c  c  c  c }
463 < \hline \\[-7mm]
464 < \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
465 < \hline \\[-3mm]
466 < \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
467 < \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
461 >
462 > \caption{The free energy of the studied ice polymorphs after applying
463 > the energy difference attributed to the inclusion of the PME
464 > long-range interaction correction. Units are kcal/mol.}
465 >
466 > \begin{tabular}{ccccc}
467 > \hline
468 > Water Model &  $I_h$ & $I_c$ &  B & Ice-{\it i} \\
469 > \hline
470 > TIP3P  & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3) \\
471 > SPC/E  & -12.97(2) & -13.00(2) & -12.96(3) & -13.02(2) \\
472   \end{tabular}
473   \label{pmeShift}
474   \end{center}
# Line 293 | Line 477 | the phase behavior of water models.
477  
478   \section{Conclusions}
479  
480 + The free energy for proton ordered variants of hexagonal and cubic ice
481 + $I$, ice B, and our recently discovered Ice-{\it i} structure were
482 + calculated under standard conditions for several common water models
483 + via thermodynamic integration. All the water models studied show
484 + Ice-{\it i} to be the minimum free energy crystal structure with a 9
485 + \AA\ switching function cutoff. Calculated melting and boiling points
486 + show surprisingly good agreement with the experimental values;
487 + however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
488 + effect of interaction truncation was investigated through variation of
489 + the cutoff radius, use of a reaction field parameterized model, and
490 + estimation of the results in the presence of the Ewald
491 + summation. Interaction truncation has a significant effect on the
492 + computed free energy values, and may significantly alter the free
493 + energy landscape for the more complex multipoint water models. Despite
494 + these effects, these results show Ice-{\it i} to be an important ice
495 + polymorph that should be considered in simulation studies.
496 +
497 + Due to this relative stability of Ice-{\it i} in all of the
498 + investigated simulation conditions, the question arises as to possible
499 + experimental observation of this polymorph.  The rather extensive past
500 + and current experimental investigation of water in the low pressure
501 + regime makes us hesitant to ascribe any relevance of this work outside
502 + of the simulation community.  It is for this reason that we chose a
503 + name for this polymorph which involves an imaginary quantity.  That
504 + said, there are certain experimental conditions that would provide the
505 + most ideal situation for possible observation. These include the
506 + negative pressure or stretched solid regime, small clusters in vacuum
507 + deposition environments, and in clathrate structures involving small
508 + non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
509 + our predictions for both the pair distribution function ($g_{OO}(r)$)
510 + and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
511 + ice-{\it i} at a temperature of 77K.  In studies of the high and low
512 + density forms of amorphous ice, ``spurious'' diffraction peaks have
513 + been observed experimentally.\cite{Bizid87} It is possible that a
514 + variant of Ice-{\it i} could explain some of this behavior; however,
515 + we will leave it to our experimental colleagues to make the final
516 + determination on whether this ice polymorph is named appropriately
517 + (i.e. with an imaginary number) or if it can be promoted to Ice-0.
518 +
519 + \begin{figure}
520 + \includegraphics[width=\linewidth]{iceGofr.eps}
521 + \caption{Radial distribution functions of ice $I_h$, $I_c$,
522 + Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
523 + of the SSD/RF water model at 77 K.}
524 + \label{fig:gofr}
525 + \end{figure}
526 +
527 + \begin{figure}
528 + \includegraphics[width=\linewidth]{sofq.eps}
529 + \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
530 + and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
531 + been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
532 + width) to compensate for the trunction effects in our finite size
533 + simulations.}
534 + \label{fig:sofq}
535 + \end{figure}
536 +
537   \section{Acknowledgments}
538   Support for this project was provided by the National Science
539   Foundation under grant CHE-0134881. Computation time was provided by
540 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
541 < DMR-0079647.
540 > the Notre Dame High Performance Computing Cluster and the Notre Dame
541 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
542  
543   \newpage
544  

Diff Legend

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