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 1812 by chrisfen, Wed Dec 1 17:35:14 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.  This crystal structure has a
95 + limited resemblence to a recent two-dimensional ice tessellation
96 + simulated on a silica surface.\cite{Yang04} Proton ordering can be
97 + accomplished by orienting two of the molecules so that both of their
98 + donated hydrogen bonds are internal to their tetramer
99 + (Fig. \ref{protOrder}).  As expected in an ice crystal constructed of
100 + water tetramers, the hydrogen bonds are not as linear as those
101 + observed in ice $I_h$, however the interlocking of these subunits
102 + appears to provide significant stabilization to the overall crystal.
103 + The arrangement of these tetramers results in surrounding open
104 + octagonal cavities that are typically greater than 6.3 \AA\ in
105 + diameter.  This relatively open overall structure leads to crystals
106 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
107 +
108 + \begin{figure}
109 + \includegraphics[width=\linewidth]{unitCell.eps}
110 + \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it i}$^\prime$,
111 + the elongated variant of Ice-{\it i}.  The spheres represent the
112 + center-of-mass locations of the water molecules.  The $a$ to $c$
113 + ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
114 + $a:2.1214c$ and $a:1.7850c$ respectively.}
115 + \label{iceiCell}
116 + \end{figure}
117 +
118 + \begin{figure}
119 + \includegraphics[width=\linewidth]{orderedIcei.eps}
120 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
121 + down the (001) crystal face. The rows of water tetramers surrounded by
122 + octagonal pores leads to a crystal structure that is significantly
123 + less dense than ice $I_h$.}
124 + \label{protOrder}
125 + \end{figure}
126 +
127 + Results from our previous study indicated that Ice-{\it i} is the
128 + minimum energy crystal structure for the single point water models we
129 + had investigated (for discussions on these single point dipole models,
130 + see our previous work and related
131 + articles).\cite{Fennell04,Liu96,Bratko85} Those results only
132 + considered energetic stabilization and neglected entropic
133 + contributions to the overall free energy.  To address this issue, we
134 + have calculated the absolute free energy of this crystal using
135 + thermodynamic integration and compared to the free energies of cubic
136 + and hexagonal ice $I$ (the experimental low density ice polymorphs)
137 + and ice B (a higher density, but very stable crystal structure
138 + observed by B\`{a}ez and Clancy in free energy studies of
139 + SPC/E).\cite{Baez95b} This work includes results for the water model
140 + from which Ice-{\it i} was crystallized (SSD/E) in addition to several
141 + common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
142 + field parametrized single point dipole water model (SSD/RF). It should
143 + be noted that a second version of Ice-{\it i} (Ice-{\it i}$^\prime$)
144 + was used in calculations involving SPC/E, TIP4P, and TIP5P.  The unit
145 + cell of this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it
146 + i} unit it is extended in the direction of the (001) face and
147 + compressed along the other two faces.  There is typically a small
148 + distortion of proton ordered Ice-{\it i}$^\prime$ that converts the
149 + normally square tetramer into a rhombus with alternating approximately
150 + 85 and 95 degree angles.  The degree of this distortion is model
151 + dependent and significant enough to split the tetramer diagonal
152 + location peak in the radial distribution function.
153 +
154   \section{Methods}
155  
156   Canonical ensemble (NVT) molecular dynamics calculations were
157 < performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
158 < molecular mechanics package. All molecules were treated as rigid
159 < bodies, with orientational motion propogated using the symplectic DLM
160 < integration method. Details about the implementation of these
161 < techniques can be found in a recent publication.\cite{Meineke05}
157 > performed using the OOPSE molecular mechanics package.\cite{Meineke05}
158 > All molecules were treated as rigid bodies, with orientational motion
159 > propagated using the symplectic DLM integration method.  Details about
160 > the implementation of this technique can be found in a recent
161 > publication.\cite{Dullweber1997}
162  
163 < Thermodynamic integration was utilized to calculate the free energy of
164 < several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
165 < SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
166 < 400 K for all of these water models were also determined using this
167 < same technique, in order to determine melting points and generate
168 < phase diagrams. All simulations were carried out at densities
169 < resulting in a pressure of approximately 1 atm at their respective
170 < temperatures.
163 > Thermodynamic integration is an established technique for
164 > determination of free energies of condensed phases of
165 > materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
166 > method, implemented in the same manner illustrated by B\`{a}ez and
167 > Clancy, was utilized to calculate the free energy of several ice
168 > crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
169 > SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
170 > and 400 K for all of these water models were also determined using
171 > this same technique in order to determine melting points and to
172 > generate phase diagrams.  All simulations were carried out at
173 > densities which correspond to a pressure of approximately 1 atm at
174 > their respective temperatures.
175  
176 + Thermodynamic integration involves a sequence of simulations during
177 + which the system of interest is converted into a reference system for
178 + which the free energy is known analytically.  This transformation path
179 + is then integrated in order to determine the free energy difference
180 + between the two states:
181 + \begin{equation}
182 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
183 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
184 + \end{equation}
185 + where $V$ is the interaction potential and $\lambda$ is the
186 + transformation parameter that scales the overall potential.
187 + Simulations are distributed strategically along this path in order to
188 + sufficiently sample the regions of greatest change in the potential.
189 + Typical integrations in this study consisted of $\sim$25 simulations
190 + ranging from 300 ps (for the unaltered system) to 75 ps (near the
191 + reference state) in length.
192 +
193   For the thermodynamic integration of molecular crystals, the Einstein
194 < Crystal is chosen as the reference state that the system is converted
195 < to over the course of the simulation. In an Einstein Crystal, the
196 < molecules are harmonically restrained at their ideal lattice locations
197 < and orientations. The partition function for a molecular crystal
198 < restrained in this fashion has been evaluated, and the Helmholtz Free
199 < Energy ({\it A}) is given by
194 > crystal was chosen as the reference system.  In an Einstein crystal,
195 > the molecules are restrained at their ideal lattice locations and
196 > orientations. Using harmonic restraints, as applied by B\`{a}ez and
197 > Clancy, the total potential for this reference crystal
198 > ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
199 > \begin{equation}
200 > V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
201 > \frac{K_\omega\omega^2}{2},
202 > \end{equation}
203 > where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
204 > the spring constants restraining translational motion and deflection
205 > of and rotation around the principle axis of the molecule
206 > respectively.  These spring constants are typically calculated from
207 > the mean-square displacements of water molecules in an unrestrained
208 > ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
209 > mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
210 > 17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
211 > the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
212 > from $-\pi$ to $\pi$.  The partition function for a molecular crystal
213 > restrained in this fashion can be evaluated analytically, and the
214 > Helmholtz Free Energy ({\it A}) is given by
215   \begin{eqnarray}
216   A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
217   [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
# Line 84 | Line 223 | where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a
223   )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
224   \label{ecFreeEnergy}
225   \end{eqnarray}
226 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
227 < \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
228 < $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.
226 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
227 > potential energy of the ideal crystal.\cite{Baez95a}
228 >
229   \begin{figure}
230 < \includegraphics[scale=1.0]{rotSpring.eps}
230 > \includegraphics[width=\linewidth]{rotSpring.eps}
231   \caption{Possible orientational motions for a restrained molecule.
232   $\theta$ angles correspond to displacement from the body-frame {\it
233   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 103 | Line 237 | Charge, dipole, and Lennard-Jones interactions were mo
237   \label{waterSpring}
238   \end{figure}
239  
240 + In the case of molecular liquids, the ideal vapor is chosen as the
241 + target reference state.  There are several examples of liquid state
242 + free energy calculations of water models present in the
243 + literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
244 + typically differ in regard to the path taken for switching off the
245 + interaction potential to convert the system to an ideal gas of water
246 + molecules.  In this study, we applied of one of the most convenient
247 + methods and integrated over the $\lambda^4$ path, where all
248 + interaction parameters are scaled equally by this transformation
249 + parameter.  This method has been shown to be reversible and provide
250 + results in excellent agreement with other established
251 + methods.\cite{Baez95b}
252 +
253   Charge, dipole, and Lennard-Jones interactions were modified by a
254 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
255 < applying this function, these interactions are smoothly truncated,
256 < thereby avoiding poor energy conserving dynamics resulting from
257 < harsher truncation schemes. The effect of a long-range correction was
258 < also investigated on select model systems in a variety of manners. For
259 < the SSD/RF model, a reaction field with a fixed dielectric constant of
260 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
261 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
262 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
263 < compare with the 9 \AA\ cutoff results. Finally, results from the use
264 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
265 < calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
266 < mechanics software package. TINKER was chosen because it can also
267 < propogate the motion of rigid-bodies, and provides the most direct
268 < comparison to the results from OOPSE. The calculated energy difference
269 < in the presence and absence of PME was applied to the previous results
123 < in order to predict changes in the free energy landscape.
254 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA).
255 > By applying this function, these interactions are smoothly truncated,
256 > thereby avoiding the poor energy conservation which results from
257 > harsher truncation schemes.  The effect of a long-range correction was
258 > also investigated on select model systems in a variety of manners.
259 > For the SSD/RF model, a reaction field with a fixed dielectric
260 > constant of 80 was applied in all simulations.\cite{Onsager36} For a
261 > series of the least computationally expensive models (SSD/E, SSD/RF,
262 > and TIP3P), simulations were performed with longer cutoffs of 12 and
263 > 15 \AA\ to compare with the 9 \AA\ cutoff results.  Finally, the
264 > effects of utilizing an Ewald summation were estimated for TIP3P and
265 > SPC/E by performing single configuration calculations with
266 > Particle-Mesh Ewald (PME) in the TINKER molecular mechanics software
267 > package.\cite{Tinker} The calculated energy difference in the presence
268 > and absence of PME was applied to the previous results in order to
269 > predict changes to the free energy landscape.
270  
271   \section{Results and discussion}
272  
273 < The free energy of proton ordered Ice-{\it i} was calculated and
273 > The free energy of proton-ordered Ice-{\it i} was calculated and
274   compared with the free energies of proton ordered variants of the
275   experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
276   as well as the higher density ice B, observed by B\`{a}ez and Clancy
277   and thought to be the minimum free energy structure for the SPC/E
278   model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
279 < Ice XI, the experimentally observed proton ordered variant of ice
280 < $I_h$, was investigated initially, but it was found not to be as
281 < stable as antiferroelectric variants of proton ordered or even proton
282 < disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
283 < ice $I_h$ used here is a simple antiferroelectric version that has an
284 < 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
285 < ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
286 < $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
287 < were necessary for simulations involving larger cutoff values.
279 > Ice XI, the experimentally-observed proton-ordered variant of ice
280 > $I_h$, was investigated initially, but was found to be not as stable
281 > as proton disordered or antiferroelectric variants of ice $I_h$.  The
282 > proton ordered variant of ice $I_h$ used here is a simple
283 > antiferroelectric version that we devised, and it has an 8 molecule
284 > unit cell similar to other predicted antiferroelectric $I_h$
285 > crystals.\cite{Davidson84} The crystals contained 648 or 1728
286 > molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
287 > molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}.  The
288 > larger crystal sizes were necessary for simulations involving larger
289 > cutoff values.
290  
291   \begin{table*}
292   \begin{minipage}{\linewidth}
145 – \renewcommand{\thefootnote}{\thempfootnote}
293   \begin{center}
294 +
295   \caption{Calculated free energies for several ice polymorphs with a
296   variety of common water models. All calculations used a cutoff radius
297   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
298 < kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
299 < \begin{tabular}{ l  c  c  c  c }
300 < \hline \\[-7mm]
301 < \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
302 < \hline \\[-3mm]
303 < \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
304 < \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
305 < \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
306 < \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
307 < \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
308 < \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
298 > kcal/mol. Calculated error of the final digits is in parentheses.}
299 >
300 > \begin{tabular}{lcccc}
301 > \hline
302 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
303 > \hline
304 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
305 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
306 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
307 > SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & -13.55(2)\\
308 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
309 > SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2)\\
310   \end{tabular}
311   \label{freeEnergy}
312   \end{center}
# Line 166 | Line 315 | models studied. With the free energy at these state po
315  
316   The free energy values computed for the studied polymorphs indicate
317   that Ice-{\it i} is the most stable state for all of the common water
318 < models studied. With the free energy at these state points, the
319 < temperature and pressure dependence of the free energy was used to
320 < project to other state points and build phase diagrams. Figures
321 < \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
322 < from the free energy results. All other models have similar structure,
323 < only the crossing points between these phases exist at different
324 < temperatures and pressures. It is interesting to note that ice $I$
318 > models studied.  With the calculated free energy at these state
319 > points, the Gibbs-Helmholtz equation was used to project to other
320 > state points and to build phase diagrams.  Figures \ref{tp3phasedia}
321 > and \ref{ssdrfphasedia} are example diagrams built from the free
322 > energy results.  All other models have similar structure, although the
323 > crossing points between the phases move to slightly different
324 > temperatures and pressures.  It is interesting to note that ice $I$
325   does not exist in either cubic or hexagonal form in any of the phase
326 < diagrams for any of the models. For purposes of this study, ice B is
327 < representative of the dense ice polymorphs. A recent study by Sanz
326 > diagrams for any of the models.  For purposes of this study, ice B is
327 > representative of the dense ice polymorphs.  A recent study by Sanz
328   {\it et al.} goes into detail on the phase diagrams for SPC/E and
329 < TIP4P in the high pressure regime.\cite{Sanz04}
329 > TIP4P at higher pressures than those studied here.\cite{Sanz04}
330 >
331   \begin{figure}
332   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
333   \caption{Phase diagram for the TIP3P water model in the low pressure
334 < regime. The displayed $T_m$ and $T_b$ values are good predictions of
334 > regime.  The displayed $T_m$ and $T_b$ values are good predictions of
335   the experimental values; however, the solid phases shown are not the
336 < experimentally observed forms. Both cubic and hexagonal ice $I$ are
336 > experimentally observed forms.  Both cubic and hexagonal ice $I$ are
337   higher in energy and don't appear in the phase diagram.}
338   \label{tp3phasedia}
339   \end{figure}
340 +
341   \begin{figure}
342   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
343   \caption{Phase diagram for the SSD/RF water model in the low pressure
344 < regime. Calculations producing these results were done under an
345 < applied reaction field. It is interesting to note that this
344 > regime.  Calculations producing these results were done under an
345 > applied reaction field.  It is interesting to note that this
346   computationally efficient model (over 3 times more efficient than
347   TIP3P) exhibits phase behavior similar to the less computationally
348   conservative charge based models.}
# Line 200 | Line 351 | conservative charge based models.}
351  
352   \begin{table*}
353   \begin{minipage}{\linewidth}
203 – \renewcommand{\thefootnote}{\thempfootnote}
354   \begin{center}
355 +
356   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
357 < temperatures of several common water models compared with experiment.}
358 < \begin{tabular}{ l  c  c  c  c  c  c  c }
359 < \hline \\[-7mm]
360 < \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
361 < \hline \\[-3mm]
362 < \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
363 < \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
364 < \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
357 > temperatures at 1 atm for several common water models compared with
358 > experiment.  The $T_m$ and $T_s$ values from simulation correspond to
359 > a transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
360 > liquid or gas state.}
361 >
362 > \begin{tabular}{lccccccc}
363 > \hline
364 > Equilibrium Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
365 > \hline
366 > $T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
367 > $T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
368 > $T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
369   \end{tabular}
370   \label{meltandboil}
371   \end{center}
# Line 218 | Line 373 | calculated from this work. Surprisingly, most of these
373   \end{table*}
374  
375   Table \ref{meltandboil} lists the melting and boiling temperatures
376 < calculated from this work. Surprisingly, most of these models have
377 < melting points that compare quite favorably with experiment. The
376 > calculated from this work.  Surprisingly, most of these models have
377 > melting points that compare quite favorably with experiment.  The
378   unfortunate aspect of this result is that this phase change occurs
379   between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
380 < liquid state. These results are actually not contrary to previous
381 < studies in the literature. Earlier free energy studies of ice $I$
382 < using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
383 < being attributed to choice of interaction truncation and different
384 < ordered and disordered molecular arrangements). If the presence of ice
385 < B and 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
380 > liquid state.  These results are actually not contrary to other
381 > studies.  Studies of ice $I_h$ using TIP4P predict a $T_m$ ranging
382 > from 214 to 238 K (differences being attributed to choice of
383 > interaction truncation and different 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 to be 265 K, indicating that these simulation based
388 > structures ought to be included in studies probing phase transitions
389 > with this model.  Also of interest in these results is that SSD/E does
390   not exhibit a melting point at 1 atm, but it shows a sublimation point
391 < at 355 K. This is due to the significant stability of Ice-{\it i} over
392 < all other polymorphs for this particular model under these
393 < conditions. While troubling, this behavior turned out to be
394 < advantagious in that it facilitated the spontaneous crystallization of
395 < Ice-{\it i}. These observations provide a warning that simulations of
391 > at 355 K.  This is due to the significant stability of Ice-{\it i}
392 > over all other polymorphs for this particular model under these
393 > conditions.  While troubling, this behavior resulted in spontaneous
394 > crystallization of Ice-{\it i} and led us to investigate this
395 > structure.  These observations provide a warning that simulations of
396   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
397 < risk of spontaneous crystallization. However, this risk changes when
397 > risk of spontaneous crystallization.  However, this risk lessens when
398   applying a longer cutoff.
399  
400 + \begin{figure}
401 + \includegraphics[width=\linewidth]{cutoffChange.eps}
402 + \caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
403 + SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
404 + with an added Ewald correction term.  Calculations performed without a
405 + long-range correction show noticable free energy dependence on the
406 + cutoff radius and show some degree of converge at large cutoff radii.
407 + Inclusion of a long-range correction reduces the cutoff radius
408 + dependence of the free energy for all the models.  Data for ice I$_c$
409 + with TIP3P using 12 and 13.5 \AA\ cutoff radii were omitted being that
410 + the crystal was prone to distortion and melting at 200 K.}
411 + \label{incCutoff}
412 + \end{figure}
413 +
414   Increasing the cutoff radius in simulations of the more
415   computationally efficient water models was done in order to evaluate
416   the trend in free energy values when moving to systems that do not
417 < involve potential truncation. As seen in Fig. \ref{incCutoff}, the
418 < free energy of all the ice polymorphs show a substantial dependence on
419 < cutoff radius. In general, there is a narrowing of the free energy
420 < differences while moving to greater cutoff radius. This trend is much
421 < more subtle in the case of SSD/RF, indicating that the free energies
422 < calculated with a reaction field present provide a more accurate
423 < picture of the free energy landscape in the absence of potential
424 < truncation.
417 > involve potential truncation.  As seen in Fig. \ref{incCutoff}, the
418 > free energy of the ice polymorphs with water models lacking a
419 > long-range correction show a significant cutoff radius dependence.  In
420 > general, there is a narrowing of the free energy differences while
421 > moving to greater cutoff radii.  As the free energies for the
422 > polymorphs converge, the stability advantage that Ice-{\it i} exhibits
423 > is reduced.  Interestingly, increasing the cutoff radius a mere 1.5
424 > \AA\ with the SSD/E model destabilizes the Ice-{\it i} polymorph
425 > enough that the liquid state is preferred under standard simulation
426 > conditions (298 K and 1 atm).  Thus, it is recommended that
427 > simulations using this model choose interaction truncation radii
428 > greater than 9 \AA.  Considering the stabilization of Ice-{\it i} with
429 > smaller cutoffs, it is not surprising that crystallization was
430 > observed with SSD/E.  The choice of a 9 \AA\ cutoff in the previous
431 > simulations gives the Ice-{\it i} polymorph a greater than 1 kcal/mol
432 > lower free energy than the ice $I_\textrm{h}$ starting configurations.
433 > Additionally, it should be noted that ice $I_c$ is not stable with
434 > cutoff radii of 12 and 13.5 \AA\ using the TIP3P water model.  These
435 > simulations showed bulk distortions of the simulation cell that
436 > rapidly deteriorated crystal integrity.
437  
438 < To further study the changes resulting to the inclusion of a
439 < long-range interaction correction, the effect of an Ewald summation
440 < was estimated by applying the potential energy difference do to its
441 < inclusion in systems in the presence and absence of the
442 < correction. This was accomplished by calculation of the potential
443 < energy of identical crystals with and without PME using TINKER. The
444 < free energies for the investigated polymorphs using the TIP3P and
445 < SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
446 < are not fully supported in TINKER, so the results for these models
447 < could not be estimated. The same trend pointed out through increase of
448 < cutoff radius is observed in these results. Ice-{\it i} is the
449 < preferred polymorph at ambient conditions for both the TIP3P and SPC/E
450 < water models; however, there is a narrowing of the free energy
451 < differences between the various solid forms. In the case of SPC/E this
452 < narrowing is significant enough that it becomes less clear cut that
453 < Ice-{\it i} is the most stable polymorph, and is possibly metastable
454 < with respect to ice B and possibly ice $I_c$. However, these results
455 < do not significantly alter the finding that the Ice-{\it i} polymorph
456 < is a stable crystal structure that should be considered when studying
457 < the phase behavior of water models.
438 > Adjacent to each of these model plots is a system with an applied or
439 > estimated long-range correction.  SSD/RF was parametrized for use with
440 > a reaction field, and the benefit provided by this computationally
441 > inexpensive correction is apparent.  Due to the relative independence
442 > of the resultant free energies, calculations performed with a small
443 > cutoff radius provide resultant properties similar to what one would
444 > expect for the bulk material.  In the cases of TIP3P and SPC/E, the
445 > effect of an Ewald summation was estimated by applying the potential
446 > energy difference do to its inclusion in systems in the presence and
447 > absence of the correction.  This was accomplished by calculation of
448 > the potential energy of identical crystals both with and without
449 > particle mesh Ewald (PME).  Similar behavior to that observed with
450 > reaction field is seen for both of these models.  The free energies
451 > show less dependence on cutoff radius and span a more narrowed range
452 > for the various polymorphs.  Like the dipolar water models, TIP3P
453 > displays a relatively constant preference for the Ice-{\it i}
454 > polymorph.  Crystal preference is much more difficult to determine for
455 > SPC/E.  Without a long-range correction, each of the polymorphs
456 > studied assumes the role of the preferred polymorph under different
457 > cutoff conditions.  The inclusion of the Ewald correction flattens and
458 > narrows the sequences of free energies so much that they often overlap
459 > within error, indicating that other conditions, such as cell volume in
460 > microcanonical simulations, can influence the chosen polymorph upon
461 > crystallization.  All of these results support the finding that the
462 > Ice-{\it i} polymorph is a stable crystal structure that should be
463 > considered when studying the phase behavior of water models.
464  
465   \begin{table*}
466   \begin{minipage}{\linewidth}
279 – \renewcommand{\thefootnote}{\thempfootnote}
467   \begin{center}
468 < \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.}
469 < \begin{tabular}{ l  c  c  c  c }
470 < \hline \\[-7mm]
471 < \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
472 < \hline \\[-3mm]
473 < \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
474 < \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
468 >
469 > \caption{The free energy versus cutoff radius for the studied ice
470 > polymorphs using SPC/E after the inclusion of the PME long-range
471 > interaction correction. Units are kcal/mol.}
472 >
473 > \begin{tabular}{ccccc}
474 > \hline
475 > Cutoff (\AA) &  $I_h$ & $I_c$ &  B & Ice-{\it i} \\
476 > \hline
477 > 9.0   & -12.98(2) & -13.00(2) & -12.97(3) & -13.02(2) \\
478 > 10.5  & -13.13(3) & -13.09(4) & -13.17(3) & -13.11(2) \\
479 > 12.0  & -13.06(2) & -13.09(2) & -13.15(4) & -13.12(2) \\
480 > 13.5  & -13.02(2) & -13.02(2) & -13.08(2) & -13.07(2) \\
481 > 15.0  & -13.11(4) & -12.97(2) & -13.09(2) & -12.95(2) \\
482   \end{tabular}
483   \label{pmeShift}
484   \end{center}
# Line 293 | Line 487 | the phase behavior of water models.
487  
488   \section{Conclusions}
489  
490 + The free energy for proton ordered variants of hexagonal and cubic ice
491 + $I$, ice B, and our recently discovered Ice-{\it i} structure were
492 + calculated under standard conditions for several common water models
493 + via thermodynamic integration.  All the water models studied show
494 + Ice-{\it i} to be the minimum free energy crystal structure with a 9
495 + \AA\ switching function cutoff.  Calculated melting and boiling points
496 + show surprisingly good agreement with the experimental values;
497 + however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$.  The
498 + effect of interaction truncation was investigated through variation of
499 + the cutoff radius, use of a reaction field parameterized model, and
500 + estimation of the results in the presence of the Ewald summation.
501 + Interaction truncation has a significant effect on the computed free
502 + energy values, and may significantly alter the free energy landscape
503 + for the more complex multipoint water models.  Despite these effects,
504 + these results show Ice-{\it i} to be an important ice polymorph that
505 + should be considered in simulation studies.
506 +
507 + Due to this relative stability of Ice-{\it i} in all of the
508 + investigated simulation conditions, the question arises as to possible
509 + experimental observation of this polymorph.  The rather extensive past
510 + and current experimental investigation of water in the low pressure
511 + regime makes us hesitant to ascribe any relevance of this work outside
512 + of the simulation community.  It is for this reason that we chose a
513 + name for this polymorph which involves an imaginary quantity.  That
514 + said, there are certain experimental conditions that would provide the
515 + most ideal situation for possible observation. These include the
516 + negative pressure or stretched solid regime, small clusters in vacuum
517 + deposition environments, and in clathrate structures involving small
518 + non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
519 + our predictions for both the pair distribution function ($g_{OO}(r)$)
520 + and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
521 + ice-{\it i} at a temperature of 77K.  In studies of the high and low
522 + density forms of amorphous ice, ``spurious'' diffraction peaks have
523 + been observed experimentally.\cite{Bizid87} It is possible that a
524 + variant of Ice-{\it i} could explain some of this behavior; however,
525 + we will leave it to our experimental colleagues to make the final
526 + determination on whether this ice polymorph is named appropriately
527 + (i.e. with an imaginary number) or if it can be promoted to Ice-0.
528 +
529 + \begin{figure}
530 + \includegraphics[width=\linewidth]{iceGofr.eps}
531 + \caption{Radial distribution functions of ice $I_h$, $I_c$,
532 + Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
533 + of the SSD/RF water model at 77 K.}
534 + \label{fig:gofr}
535 + \end{figure}
536 +
537 + \begin{figure}
538 + \includegraphics[width=\linewidth]{sofq.eps}
539 + \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
540 + and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
541 + been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
542 + width) to compensate for the trunction effects in our finite size
543 + simulations.}
544 + \label{fig:sofq}
545 + \end{figure}
546 +
547   \section{Acknowledgments}
548   Support for this project was provided by the National Science
549   Foundation under grant CHE-0134881. Computation time was provided by
550 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
551 < DMR-0079647.
550 > the Notre Dame High Performance Computing Cluster and the Notre Dame
551 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
552  
553   \newpage
554  

Diff Legend

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