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 1456 by chrisfen, Tue Sep 14 21:55:24 2004 UTC vs.
Revision 1473 by chrisfen, Fri Sep 17 14:54:49 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 + Computer simulations are a valuable tool for studying the phase
58 + behavior of systems ranging from small or simple molecules to complex
59 + biological species.\cite{Matsumoto02,Li96,Marrink01} Useful techniques
60 + have been developed to investigate the thermodynamic properites of
61 + model substances, providing both qualitative and quantitative
62 + comparisons between simulations and
63 + experiment.\cite{Widom63,Frenkel84} Investigation of these properties
64 + leads to the development of new and more accurate models, leading to
65 + better understanding and depiction of physical processes and intricate
66 + 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{Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Berendsen98,Mahoney00,Fennell04}
72 + These models have been used to investigate important physical
73 + phenomena like phase transitions, molecule transport, and the
74 + hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
75 + choice of models available, it is only natural to compare the models
76 + under interesting thermodynamic conditions in an attempt to clarify
77 + the limitations of each of the
78 + models.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two
79 + important property to quantify are the Gibbs and Helmholtz free
80 + energies, particularly for the solid forms of water.  Difficulty in
81 + these types of studies typically arises from the assortment of
82 + possible crystalline polymorphs that water adopts over a wide range of
83 + pressures and temperatures.  There are currently 13 recognized forms
84 + of ice, and it is a challenging task to investigate the entire free
85 + energy landscape.\cite{Sanz04} Ideally, research is focused on the
86 + phases having the lowest free energy at a given state point, because
87 + these phases will dictate the true transition temperatures and
88 + pressures for the model.
89 +
90 + In this paper, standard reference state methods were applied to known
91 + crystalline water polymorphs in the low pressure regime.  This work is
92 + unique in the fact that one of the crystal lattices was arrived at
93 + through crystallization of a computationally efficient water model
94 + under constant pressure and temperature conditions. Crystallization
95 + events are interesting in and of
96 + themselves;\cite{Matsumoto02,Yamada02} however, the crystal structure
97 + obtained in this case is different from any previously observed ice
98 + polymorphs in experiment or simulation.\cite{Fennell04} We have named
99 + this structure Ice-{\it i} to indicate its origin in computational
100 + simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
101 + water molecules that stack in rows of interlocking water
102 + tetramers. Proton ordering can be accomplished by orienting two of the
103 + molecules so that both of their donated hydrogen bonds are internal to
104 + their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
105 + constructed of water tetramers, the hydrogen bonds are not as linear
106 + as those observed in ice $I_h$, however the interlocking of these
107 + subunits appears to provide significant stabilization to the overall
108 + crystal. The arrangement of these tetramers results in surrounding
109 + open octagonal cavities that are typically greater than 6.3 \AA\ in
110 + diameter. This relatively open overall structure leads to crystals
111 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
112 +
113 + \begin{figure}
114 + \includegraphics[width=\linewidth]{unitCell.eps}
115 + \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
116 + elongated variant of Ice-{\it i}.  The spheres represent the
117 + center-of-mass locations of the water molecules.  The $a$ to $c$
118 + ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
119 + $a:2.1214c$ and $a:1.7850c$ respectively.}
120 + \label{iceiCell}
121 + \end{figure}
122 +
123 + \begin{figure}
124 + \includegraphics[width=\linewidth]{orderedIcei.eps}
125 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
126 + down the (001) crystal face. The rows of water tetramers surrounded by
127 + octagonal pores leads to a crystal structure that is significantly
128 + less dense than ice $I_h$.}
129 + \label{protOrder}
130 + \end{figure}
131 +
132 + Results from our previous study indicated that Ice-{\it i} is the
133 + minimum energy crystal structure for the single point water models we
134 + investigated (for discussions on these single point dipole models, see
135 + our previous work and related
136 + articles).\cite{Fennell04,Liu96,Bratko85} Those results only
137 + considered energetic stabilization and neglected entropic
138 + contributions to the overall free energy. To address this issue, the
139 + absolute free energy of this crystal was calculated using
140 + thermodynamic integration and compared to the free energies of cubic
141 + and hexagonal ice $I$ (the experimental low density ice polymorphs)
142 + and ice B (a higher density, but very stable crystal structure
143 + observed by B\`{a}ez and Clancy in free energy studies of
144 + SPC/E).\cite{Baez95b} This work includes results for the water model
145 + from which Ice-{\it i} was crystallized (SSD/E) in addition to several
146 + common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
147 + field parametrized single point dipole water model (SSD/RF). It should
148 + be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used
149 + in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of
150 + this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit
151 + it is extended in the direction of the (001) face and compressed along
152 + the other two faces.
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
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 generate
172   phase diagrams. All simulations were carried out at densities
173   resulting in a pressure of approximately 1 atm at their respective
174   temperatures.
175  
176 + A single thermodynamic integration involves a sequence of simulations
177 + over which the system of interest is converted into a reference system
178 + for which the free energy is known analytically. This transformation
179 + path is then integrated in order to determine the free energy
180 + difference 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
187 + potential. Simulations are distributed strategically along this path
188 + in order to sufficiently sample the regions of greatest change in the
189 + potential. Typical integrations in this study consisted of $\sim$25
190 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
191 + (near the 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.  It is clear from Fig. \ref{waterSpring} that the values
207 > of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from
208 > $-\pi$ to $\pi$.  The partition function for a molecular crystal
209 > restrained in this fashion can be evaluated analytically, and the
210 > Helmholtz Free Energy ({\it A}) is given by
211   \begin{eqnarray}
212   A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
213   [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
# Line 84 | Line 219 | where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a
219   )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
220   \label{ecFreeEnergy}
221   \end{eqnarray}
222 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
223 < \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
224 < $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.
222 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
223 > potential energy of the ideal crystal.\cite{Baez95a}
224 >
225   \begin{figure}
226 < \includegraphics[scale=1.0]{rotSpring.eps}
226 > \includegraphics[width=\linewidth]{rotSpring.eps}
227   \caption{Possible orientational motions for a restrained molecule.
228   $\theta$ angles correspond to displacement from the body-frame {\it
229   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 103 | Line 233 | Charge, dipole, and Lennard-Jones interactions were mo
233   \label{waterSpring}
234   \end{figure}
235  
236 < Charge, dipole, and Lennard-Jones interactions were modified by a
237 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
238 < applying this function, these interactions are smoothly truncated,
239 < thereby avoiding poor energy conserving dynamics resulting from
240 < harsher truncation schemes. The effect of a long-range correction was
241 < also investigated on select model systems in a variety of manners. For
242 < the SSD/RF model, a reaction field with a fixed dielectric constant of
243 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
244 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
245 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
246 < compare with the 9 \AA\ cutoff results. Finally, results from the use
247 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
236 > In the case of molecular liquids, the ideal vapor is chosen as the
237 > target reference state.  There are several examples of liquid state
238 > free energy calculations of water models present in the
239 > literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
240 > typically differ in regard to the path taken for switching off the
241 > interaction potential to convert the system to an ideal gas of water
242 > molecules.  In this study, we apply of one of the most convenient
243 > methods and integrate over the $\lambda^4$ path, where all interaction
244 > parameters are scaled equally by this transformation parameter.  This
245 > method has been shown to be reversible and provide results in
246 > excellent agreement with other established methods.\cite{Baez95b}
247 >
248 > Charge, dipole, and Lennard-Jones interactions were modified by a
249 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
250 > ). By applying this function, these interactions are smoothly
251 > truncated, thereby avoiding the poor energy conservation which results
252 > from harsher truncation schemes. The effect of a long-range correction
253 > was also investigated on select model systems in a variety of
254 > manners. For the SSD/RF model, a reaction field with a fixed
255 > dielectric constant of 80 was applied in all
256 > simulations.\cite{Onsager36} For a series of the least computationally
257 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
258 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
259 > \AA\ cutoff results. Finally, results from the use of an Ewald
260 > summation were estimated for TIP3P and SPC/E by performing
261   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
262 < mechanics software package. TINKER was chosen because it can also
263 < propogate the motion of rigid-bodies, and provides the most direct
264 < comparison to the results from OOPSE. The calculated energy difference
265 < in the presence and absence of PME was applied to the previous results
123 < in order to predict changes in the free energy landscape.
262 > mechanics software package.\cite{Tinker} The calculated energy
263 > difference in the presence and absence of PME was applied to the
264 > previous results in order to predict changes to the free energy
265 > landscape.
266  
267   \section{Results and discussion}
268  
# Line 130 | Line 272 | Ice XI, the experimentally observed proton ordered var
272   as well as the higher density ice B, observed by B\`{a}ez and Clancy
273   and thought to be the minimum free energy structure for the SPC/E
274   model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
275 < Ice XI, the experimentally observed proton ordered variant of ice
276 < $I_h$, was investigated initially, but it was found not to be as
277 < stable as antiferroelectric variants of proton ordered or even proton
278 < disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
279 < ice $I_h$ used here is a simple antiferroelectric version that has an
280 < 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
281 < ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
282 < $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
283 < were necessary for simulations involving larger cutoff values.
275 > Ice XI, the experimentally-observed proton-ordered variant of ice
276 > $I_h$, was investigated initially, but was found to be not as stable
277 > as proton disordered or antiferroelectric variants of ice $I_h$. The
278 > proton ordered variant of ice $I_h$ used here is a simple
279 > antiferroelectric version that we divised, and it has an 8 molecule
280 > unit cell similar to other predicted antiferroelectric $I_h$
281 > crystals.\cite{Davidson84} The crystals contained 648 or 1728
282 > molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
283 > molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
284 > crystal sizes were necessary for simulations involving larger cutoff
285 > values.
286  
287   \begin{table*}
288   \begin{minipage}{\linewidth}
# Line 147 | Line 291 | kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF
291   \caption{Calculated free energies for several ice polymorphs with a
292   variety of common water models. All calculations used a cutoff radius
293   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
294 < kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
294 > kcal/mol. Calculated error of the final digits is in parentheses. *Ice
295 > $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
296   \begin{tabular}{ l  c  c  c  c }
297 < \hline \\[-7mm]
298 < \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
299 < \hline \\[-3mm]
300 < \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
301 < \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
302 < \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
303 < \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
304 < \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
305 < \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
297 > \hline
298 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
299 > \hline
300 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
301 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
302 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
303 > SPC/E & -12.67(2) & -12.96(2) & -13.25(3) & -13.55(2)\\
304 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
305 > SSD/RF & -11.51(2) & NA* & -12.08(3) & -12.29(2)\\
306   \end{tabular}
307   \label{freeEnergy}
308   \end{center}
# Line 167 | Line 312 | temperature and pressure dependence of the free energy
312   The free energy values computed for the studied polymorphs indicate
313   that Ice-{\it i} is the most stable state for all of the common water
314   models studied. With the free energy at these state points, the
315 < temperature and pressure dependence of the free energy was used to
316 < project to other state points and build phase diagrams. Figures
315 > Gibbs-Helmholtz equation was used to project to other state points and
316 > to build phase diagrams.  Figures
317   \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
318   from the free energy results. All other models have similar structure,
319 < only the crossing points between these phases exist at different
320 < temperatures and pressures. It is interesting to note that ice $I$
321 < does not exist in either cubic or hexagonal form in any of the phase
322 < diagrams for any of the models. For purposes of this study, ice B is
323 < representative of the dense ice polymorphs. A recent study by Sanz
324 < {\it et al.} goes into detail on the phase diagrams for SPC/E and
325 < TIP4P in the high pressure regime.\cite{Sanz04}
319 > although the crossing points between the phases exist at slightly
320 > different temperatures and pressures. It is interesting to note that
321 > ice $I$ does not exist in either cubic or hexagonal form in any of the
322 > phase diagrams for any of the models. For purposes of this study, ice
323 > B is representative of the dense ice polymorphs. A recent study by
324 > Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
325 > TIP4P in the high pressure regime.\cite{Sanz04}
326 >
327   \begin{figure}
328   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
329   \caption{Phase diagram for the TIP3P water model in the low pressure
# Line 187 | Line 333 | higher in energy and don't appear in the phase diagram
333   higher in energy and don't appear in the phase diagram.}
334   \label{tp3phasedia}
335   \end{figure}
336 +
337   \begin{figure}
338   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
339   \caption{Phase diagram for the SSD/RF water model in the low pressure
# Line 203 | Line 350 | temperatures of several common water models compared w
350   \renewcommand{\thefootnote}{\thempfootnote}
351   \begin{center}
352   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
353 < temperatures of several common water models compared with experiment.}
353 > temperatures at 1 atm for several common water models compared with
354 > experiment. The $T_m$ and $T_s$ values from simulation correspond to a
355 > transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
356 > liquid or gas state.}
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 & \ \ - & \ \ -\\
358 > \hline
359 > Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
360 > \hline
361 > $T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
362 > $T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
363 > $T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
364   \end{tabular}
365   \label{meltandboil}
366   \end{center}
# Line 226 | Line 376 | ordered and disordered molecular arrangements). If the
376   studies in the literature. Earlier free energy studies of ice $I$
377   using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
378   being attributed to choice of interaction truncation and different
379 < ordered and disordered molecular arrangements). If the presence of ice
380 < B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
379 > ordered and disordered molecular
380 > arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
381 > Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
382   predicted from this work. However, the $T_m$ from Ice-{\it i} is
383   calculated at 265 K, significantly higher in temperature than the
384   previous studies. Also of interest in these results is that SSD/E does
# Line 235 | Line 386 | advantagious in that it facilitated the spontaneous cr
386   at 355 K. This is due to the significant stability of Ice-{\it i} over
387   all other polymorphs for this particular model under these
388   conditions. While troubling, this behavior turned out to be
389 < advantagious in that it facilitated the spontaneous crystallization of
389 > advantageous in that it facilitated the spontaneous crystallization of
390   Ice-{\it i}. These observations provide a warning that simulations of
391   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
392   risk of spontaneous crystallization. However, this risk changes when
393   applying a longer cutoff.
394  
395 + \begin{figure}
396 + \includegraphics[width=\linewidth]{cutoffChange.eps}
397 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
398 + TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
399 + \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
400 + \AA . These crystals are unstable at 200 K and rapidly convert into
401 + liquids. The connecting lines are qualitative visual aid.}
402 + \label{incCutoff}
403 + \end{figure}
404  
405 + Increasing the cutoff radius in simulations of the more
406 + computationally efficient water models was done in order to evaluate
407 + the trend in free energy values when moving to systems that do not
408 + involve potential truncation. As seen in Fig. \ref{incCutoff}, the
409 + free energy of all the ice polymorphs show a substantial dependence on
410 + cutoff radius. In general, there is a narrowing of the free energy
411 + differences while moving to greater cutoff radius. Interestingly, by
412 + increasing the cutoff radius, the free energy gap was narrowed enough
413 + in the SSD/E model that the liquid state is preferred under standard
414 + simulation conditions (298 K and 1 atm). Thus, it is recommended that
415 + simulations using this model choose interaction truncation radii
416 + greater than 9 \AA\ . This narrowing trend is much more subtle in the
417 + case of SSD/RF, indicating that the free energies calculated with a
418 + reaction field present provide a more accurate picture of the free
419 + energy landscape in the absence of potential truncation.
420  
421 + To further study the changes resulting to the inclusion of a
422 + long-range interaction correction, the effect of an Ewald summation
423 + was estimated by applying the potential energy difference do to its
424 + inclusion in systems in the presence and absence of the
425 + correction. This was accomplished by calculation of the potential
426 + energy of identical crystals with and without PME using TINKER. The
427 + free energies for the investigated polymorphs using the TIP3P and
428 + SPC/E water models are shown in Table \ref{pmeShift}. The same trend
429 + pointed out through increase of cutoff radius is observed in these PME
430 + results. Ice-{\it i} is the preferred polymorph at ambient conditions
431 + for both the TIP3P and SPC/E water models; however, the narrowing of
432 + the free energy differences between the various solid forms is
433 + significant enough that it becomes less clear that it is the most
434 + stable polymorph in the SPC/E.  The free energies of Ice-{\it i} and
435 + ice B nearly overlap within error, with ice $I_c$ just outside as
436 + well, indicating that Ice-{\it i} might be metastable with respect to
437 + ice B and possibly ice $I_c$ in the SPC/E water model. However, these
438 + results do not significantly alter the finding that the Ice-{\it i}
439 + polymorph is a stable crystal structure that should be considered when
440 + studying the phase behavior of water models.
441 +
442 + \begin{table*}
443 + \begin{minipage}{\linewidth}
444 + \renewcommand{\thefootnote}{\thempfootnote}
445 + \begin{center}
446 + \caption{The free energy of the studied ice polymorphs after applying
447 + the energy difference attributed to the inclusion of the PME
448 + long-range interaction correction. Units are kcal/mol.}
449 + \begin{tabular}{ l  c  c  c  c }
450 + \hline
451 + \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
452 + \hline
453 + TIP3P  & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3)\\
454 + SPC/E  & -12.77(2) & -12.92(2) & -12.96(3) & -13.02(2)\\
455 + \end{tabular}
456 + \label{pmeShift}
457 + \end{center}
458 + \end{minipage}
459 + \end{table*}
460 +
461   \section{Conclusions}
462  
463 + The free energy for proton ordered variants of hexagonal and cubic ice
464 + $I$, ice B, and recently discovered Ice-{\it i} were calculated under
465 + standard conditions for several common water models via thermodynamic
466 + integration. All the water models studied show Ice-{\it i} to be the
467 + minimum free energy crystal structure in the with a 9 \AA\ switching
468 + function cutoff. Calculated melting and boiling points show
469 + surprisingly good agreement with the experimental values; however, the
470 + solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
471 + interaction truncation was investigated through variation of the
472 + cutoff radius, use of a reaction field parameterized model, and
473 + estimation of the results in the presence of the Ewald
474 + summation. Interaction truncation has a significant effect on the
475 + computed free energy values, and may significantly alter the free
476 + energy landscape for the more complex multipoint water models. Despite
477 + these effects, these results show Ice-{\it i} to be an important ice
478 + polymorph that should be considered in simulation studies.
479 +
480 + Due to this relative stability of Ice-{\it i} in all manner of
481 + investigated simulation examples, the question arises as to possible
482 + experimental observation of this polymorph.  The rather extensive past
483 + and current experimental investigation of water in the low pressure
484 + regime makes us hesitant to ascribe any relevance of this work outside
485 + of the simulation community.  It is for this reason that we chose a
486 + name for this polymorph which involves an imaginary quantity.  That
487 + said, there are certain experimental conditions that would provide the
488 + most ideal situation for possible observation. These include the
489 + negative pressure or stretched solid regime, small clusters in vacuum
490 + deposition environments, and in clathrate structures involving small
491 + non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
492 + our predictions for both the pair distribution function ($g_{OO}(r)$)
493 + and the structure factor ($S(\vec{q})$ for ice $I_c$ and for ice-{\it
494 + i} at a temperature of 77K.  In a quick comparison of the predicted
495 + S(q) for Ice-{\it i} and experimental studies of amorphous solid
496 + water, it is possible that some of the ``spurious'' peaks that could
497 + not be assigned in HDA could correspond to peaks labeled in this
498 + S(q).\cite{Bizid87} It should be noted that there is typically poor
499 + agreement on crystal densities between simulation and experiment, so
500 + such peak comparisons should be made with caution.  We will leave it
501 + to our experimental colleagues to determine whether this ice polymorph
502 + is named appropriately or if it should be promoted to Ice-0.
503 +
504 + \begin{figure}
505 + \includegraphics[width=\linewidth]{iceGofr.eps}
506 + \caption{Radial distribution functions of Ice-{\it i} and ice $I_c$
507 + calculated from from simulations of the SSD/RF water model at 77 K.}
508 + \label{fig:gofr}
509 + \end{figure}
510 +
511 + \begin{figure}
512 + \includegraphics[width=\linewidth]{sofq.eps}
513 + \caption{Predicted structure factors for Ice-{\it i} and ice $I_c$ at
514 + 77 K.  The raw structure factors have been convoluted with a gaussian
515 + instrument function (0.075 \AA$^{-1}$ width) to compensate for the
516 + trunction effects in our finite size simulations. The labeled peaks
517 + compared favorably with ``spurious'' peaks observed in experimental
518 + studies of amorphous solid water.\cite{Bizid87}}
519 + \label{fig:sofq}
520 + \end{figure}
521 +
522   \section{Acknowledgments}
523   Support for this project was provided by the National Science
524   Foundation under grant CHE-0134881. Computation time was provided by
525 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
526 < DMR-0079647.
525 > the Notre Dame High Performance Computing Cluster and the Notre Dame
526 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
527  
528   \newpage
529  

Diff Legend

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