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 1834 by chrisfen, Thu Dec 2 18:58:25 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{Free Energy Analysis of Simulated Ice Polymorphs Using Simple
24 > Dipolar and Charge Based 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
122 + by 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
234 < body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
234 > body-frame {\it z}-axis.  $K_\theta$ and $K_\omega$ are spring
235   constants for the harmonic springs restraining motion in the $\theta$
236   and $\omega$ directions.}
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\\
296 > variety of common water models.  All calculations used a cutoff radius
297 > of 9.0 \AA\ and were performed at 200 K and $\sim$1 atm.  Units are
298 > kcal/mol.  Calculated error of the final digits is in parentheses.}
299 >
300 > \begin{tabular}{lccccc}
301 > \hline
302 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i} & Ice-{\it i}$^\prime$\\
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.  Error for the
409 + larger cutoff points is equivalent to that observed at 9.0 \AA\ (see
410 + Table \ref{freeEnergy}).  Data for ice I$_c$ with TIP3P using both 12
411 + and 13.5 \AA\ cutoffs were omitted because the crystal was prone to
412 + distortion and melting at 200 K.  Ice-{\it i}$^\prime$ is the form of
413 + Ice-{\it i} used in the SPC/E simulations.}
414 + \label{incCutoff}
415 + \end{figure}
416  
417 + Increasing the cutoff radius in simulations of the more
418 + computationally efficient water models was done in order to evaluate
419 + the trend in free energy values when moving to systems that do not
420 + involve potential truncation.  As seen in Fig. \ref{incCutoff}, the
421 + free energy of the ice polymorphs with water models lacking a
422 + long-range correction show a significant cutoff radius dependence.  In
423 + general, there is a narrowing of the free energy differences while
424 + moving to greater cutoff radii.  As the free energies for the
425 + polymorphs converge, the stability advantage that Ice-{\it i} exhibits
426 + is reduced.  Interestingly, increasing the cutoff radius a mere 1.5
427 + \AA\ with the SSD/E model destabilizes the Ice-{\it i} polymorph
428 + enough that the liquid state is preferred under standard simulation
429 + conditions (298 K and 1 atm).  Thus, it is recommended that
430 + simulations using this model choose interaction truncation radii
431 + greater than 9 \AA.  Considering the stabilization of Ice-{\it i} with
432 + smaller cutoffs, it is not surprising that crystallization was
433 + observed with SSD/E.  The choice of a 9 \AA\ cutoff in the previous
434 + simulations gives the Ice-{\it i} polymorph a greater than 1 kcal/mol
435 + lower free energy than the ice $I_\textrm{h}$ starting configurations.
436 + Additionally, it should be noted that ice $I_c$ is not stable with
437 + cutoff radii of 12 and 13.5 \AA\ using the TIP3P water model.  These
438 + simulations showed bulk distortions of the simulation cell that
439 + rapidly deteriorated crystal integrity.
440  
441 + Adjacent to each of these model plots is a system with an applied or
442 + estimated long-range correction.  SSD/RF was parametrized for use with
443 + a reaction field, and the benefit provided by this computationally
444 + inexpensive correction is apparent.  Due to the relative independence
445 + of the resultant free energies, calculations performed with a small
446 + cutoff radius provide resultant properties similar to what one would
447 + expect for the bulk material.  In the cases of TIP3P and SPC/E, the
448 + effect of an Ewald summation was estimated by applying the potential
449 + energy difference do to its inclusion in systems in the presence and
450 + absence of the correction.  This was accomplished by calculation of
451 + the potential energy of identical crystals both with and without
452 + particle mesh Ewald (PME).  Similar behavior to that observed with
453 + reaction field is seen for both of these models.  The free energies
454 + show less dependence on cutoff radius and span a more narrowed range
455 + for the various polymorphs.  Like the dipolar water models, TIP3P
456 + displays a relatively constant preference for the Ice-{\it i}
457 + polymorph.  Crystal preference is much more difficult to determine for
458 + SPC/E.  Without a long-range correction, each of the polymorphs
459 + studied assumes the role of the preferred polymorph under different
460 + cutoff conditions.  The inclusion of the Ewald correction flattens and
461 + narrows the sequences of free energies so much that they often overlap
462 + within error, indicating that other conditions, such as cell volume in
463 + microcanonical simulations, can influence the chosen polymorph upon
464 + crystallization.  All of these results support the finding that the
465 + Ice-{\it i} polymorph is a stable crystal structure that should be
466 + considered when studying the phase behavior of water models.
467 +
468   \section{Conclusions}
469  
470 + The free energy for proton ordered variants of hexagonal and cubic ice
471 + $I$, ice B, and our recently discovered Ice-{\it i} structure were
472 + calculated under standard conditions for several common water models
473 + via thermodynamic integration.  All the water models studied show
474 + Ice-{\it i} to be the minimum free energy crystal structure with a 9
475 + \AA\ switching function cutoff.  Calculated melting and boiling points
476 + show surprisingly good agreement with the experimental values;
477 + however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$.  The
478 + effect of interaction truncation was investigated through variation of
479 + the cutoff radius, use of a reaction field parameterized model, and
480 + estimation of the results in the presence of the Ewald summation.
481 + Interaction truncation has a significant effect on the computed free
482 + energy values, and may significantly alter the free energy landscape
483 + for the more complex multipoint water models.  Despite these effects,
484 + these results show Ice-{\it i} to be an important ice polymorph that
485 + should be considered in simulation studies.
486 +
487 + Due to this relative stability of Ice-{\it i} in all of the
488 + investigated simulation conditions, the question arises as to possible
489 + experimental observation of this polymorph.  The rather extensive past
490 + and current experimental investigation of water in the low pressure
491 + regime makes us hesitant to ascribe any relevance of this work outside
492 + of the simulation community.  It is for this reason that we chose a
493 + name for this polymorph which involves an imaginary quantity.  That
494 + said, there are certain experimental conditions that would provide the
495 + most ideal situation for possible observation. These include the
496 + negative pressure or stretched solid regime, small clusters in vacuum
497 + deposition environments, and in clathrate structures involving small
498 + non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
499 + our predictions for both the pair distribution function ($g_{OO}(r)$)
500 + and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
501 + ice-{\it i} at a temperature of 77K.  In studies of the high and low
502 + density forms of amorphous ice, ``spurious'' diffraction peaks have
503 + been observed experimentally.\cite{Bizid87} It is possible that a
504 + variant of Ice-{\it i} could explain some of this behavior; however,
505 + we will leave it to our experimental colleagues to make the final
506 + determination on whether this ice polymorph is named appropriately
507 + (i.e. with an imaginary number) or if it can be promoted to Ice-0.
508 +
509 + \begin{figure}
510 + \includegraphics[width=\linewidth]{iceGofr.eps}
511 + \caption{Radial distribution functions of ice $I_h$, $I_c$,
512 + Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
513 + of the SSD/RF water model at 77 K.}
514 + \label{fig:gofr}
515 + \end{figure}
516 +
517 + \begin{figure}
518 + \includegraphics[width=\linewidth]{sofq.eps}
519 + \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
520 + and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
521 + been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
522 + width) to compensate for the trunction effects in our finite size
523 + simulations.}
524 + \label{fig:sofq}
525 + \end{figure}
526 +
527   \section{Acknowledgments}
528   Support for this project was provided by the National Science
529   Foundation under grant CHE-0134881. Computation time was provided by
530 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
531 < DMR-0079647.
530 > the Notre Dame High Performance Computing Cluster and the Notre Dame
531 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
532  
533   \newpage
534  

Diff Legend

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