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 1474 by chrisfen, Fri Sep 17 14:56:05 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. This trend is much
412 < more subtle in the case of SSD/RF, indicating that the free energies
413 < calculated with a reaction field present provide a more accurate
414 < picture of the free energy landscape in the absence of potential
415 < truncation.
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
# Line 260 | Line 425 | SPC/E water models are shown in Table \ref{pmeShift}.
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}. TIP4P and TIP5P
429 < are not fully supported in TINKER, so the results for these models
430 < could not be estimated. The same trend pointed out through increase of
431 < cutoff radius is observed in these results. Ice-{\it i} is the
432 < preferred polymorph at ambient conditions for both the TIP3P and SPC/E
433 < water models; however, there is a narrowing of the free energy
434 < differences between the various solid forms. In the case of SPC/E this
435 < narrowing is significant enough that it becomes less clear cut that
436 < Ice-{\it i} is the most stable polymorph, and is possibly metastable
437 < with respect to ice B and possibly ice $I_c$. However, these results
438 < do not significantly alter the finding that the Ice-{\it i} polymorph
439 < is a stable crystal structure that should be considered when studying
440 < the phase behavior of water models.
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 with the SPC/E model.  The free energies of Ice-{\it
435 > i} and ice B nearly overlap within error, with ice $I_c$ just outside
436 > as well, indicating that Ice-{\it i} might be metastable with respect
437 > to ice B and possibly ice $I_c$ with SPC/E. However, these results do
438 > not significantly alter the finding that the Ice-{\it i} polymorph is
439 > a stable crystal structure that should be considered when studying the
440 > 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 the energy difference attributed to the inclusion of the PME long-range interaction correction. Units are kcal/mol.}
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 \\[-7mm]
450 > \hline
451   \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
452 < \hline \\[-3mm]
453 < \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
454 < \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
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}
# Line 293 | Line 460 | the phase behavior of water models.
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