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 1459 by chrisfen, Wed Sep 15 19:06:05 2004 UTC vs.
Revision 1479 by chrisfen, Mon Sep 20 20:57:46 2004 UTC

# Line 1 | Line 1
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  
20 %\renewcommand\citemid{\ } % no comma in optional reference note
21
21   \begin{document}
22  
23 < \title{A Free Energy Study of Low Temperature and Anomalous 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}}
28 <
29 < \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 of systems consisting
38 < of a variety of common water models. Ice-{\it i}, a recent
39 < computationally observed solid structure, was determined to be the
40 < stable state with the lowest free energy for all the water models
41 < investigated. Phase diagrams were generated, and melting and boiling
42 < points for all the models were determined and show relatively good
43 < agreement with experiment, although the solid phase is different
44 < between simulation and experiment. In addition, potential truncation
45 < was shown to have an effect on the calculated free energies, and may
46 < result in altered free energy landscapes.
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 for the new crystal were generated and
49 > we await experimental confirmation of the existence of this new
50 > polymorph.  
51   \end{abstract}
52  
51 \maketitle
52
53 \newpage
54
53   %\narrowtext
54  
55   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 60 | Line 58 | Molecular dynamics has developed into a valuable tool
58  
59   \section{Introduction}
60  
63 Molecular dynamics has developed into a valuable tool for studying the
64 phase behavior of systems ranging from small or simple
65 molecules\cite{smallStuff} to complex biological
66 species.\cite{bigStuff} Many techniques have been developed in order
67 to investigate the thermodynamic properites of model substances,
68 providing both qualitative and quantitative comparisons between
69 simulations and experiment.\cite{thermMethods} Investigation of these
70 properties leads to the development of new and more accurate models,
71 leading to better understanding and depiction of physical processes
72 and intricate molecular systems.
73
61   Water has proven to be a challenging substance to depict in
62 < simulations, and has resulted in a variety of models that attempt to
63 < describe its behavior under a varying simulation
64 < conditions.\cite{lotsOfWaterPapers} Many of these models have been
65 < used to investigate important physical phenomena like phase
66 < transitions and the hydrophobic effect.\cite{evenMorePapers} With the
67 < advent of numerous differing models, it is only natural that attention
68 < is placed on the properties of the models themselves in an attempt to
69 < clarify their benefits and limitations when applied to a system of
70 < interest.\cite{modelProps} One important but challenging property to
71 < quantify is the free energy, particularly of the solid forms of
72 < water. Difficulty in these types of studies typically arises from the
73 < assortment of possible crystalline polymorphs that water that water
74 < adopts over a wide range of pressures and temperatures. There are
75 < currently 13 recognized forms of ice, and it is a challenging task to
76 < investigate the entire free energy landscape.\cite{Sanz04} Ideally,
77 < research is focused on the phases having the lowest free energy,
78 < because these phases will dictate the true transition temperatures and
79 < pressures for their respective model.
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 the
84 < study of crystalline water polymorphs in the low pressure regime. This
85 < work is unique in the fact that one of the crystal lattices was
86 < arrived at through crystallization of a computationally efficient
87 < water model under constant pressure and temperature
88 < conditions. Crystallization events are interesting in and of
89 < themselves\cite{nucleationStudies}; however, the crystal structure
90 < obtained in this case was different from any previously observed ice
91 < polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
92 < was termed Ice-{\it i} in homage to its origin in computational
93 < simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
94 < water molecules that stack in rows of interlocking water
95 < tetramers. Proton ordering can be accomplished by orienting two of the
96 < waters so that both of their donated hydrogen bonds are internal to
97 < their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
98 < constructed of water tetramers, the hydrogen bonds are not as linear
99 < as those observed in ice $I_h$, however the interlocking of these
100 < subunits appears to provide significant stabilization to the overall
83 > In this paper, standard reference state methods were applied to known
84 > crystalline water polymorphs in the low pressure regime.  This work is
85 > unique in that one of the crystal lattices was arrived at through
86 > crystallization of a computationally efficient water model under
87 > constant pressure and temperature conditions. Crystallization events
88 > are interesting in and of themselves;\cite{Matsumoto02,Yamada02}
89 > however, the crystal structure obtained in this case is different from
90 > any previously observed ice polymorphs in experiment or
91 > simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
92 > to indicate its origin in computational simulation. The unit cell
93 > (Fig. \ref{iceiCell}A) consists of eight water molecules that stack in
94 > rows of interlocking water tetramers. Proton ordering can be
95 > accomplished by orienting two of the molecules so that both of their
96 > donated hydrogen bonds are internal to their tetramer
97 > (Fig. \ref{protOrder}). As expected in an ice crystal constructed of
98 > water tetramers, the hydrogen bonds are not as linear as those
99 > observed in ice $I_h$, however the interlocking of these subunits
100 > appears to provide significant stabilization to the overall
101   crystal. The arrangement of these tetramers results in surrounding
102   open octagonal cavities that are typically greater than 6.3 \AA\ in
103   diameter. This relatively open overall structure leads to crystals
104   that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
105  
106 < Results in the previous study indicated that Ice-{\it i} is the
107 < minimum energy crystal structure for the single point water models
108 < being studied (for discussions on these single point dipole models,
109 < see the previous work and related
110 < articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
111 < consider energetic stabilization and neglect entropic contributions to
112 < the overall free energy. To address this issue, the absolute free
113 < energy of this crystal was calculated using thermodynamic integration
114 < and compared to the free energies of cubic and hexagonal ice $I$ (the
126 < experimental low density ice polymorphs) and ice B (a higher density,
127 < but very stable crystal structure observed by B\`{a}ez and Clancy in
128 < free energy studies of SPC/E).\cite{Baez95b} This work includes
129 < results for the water model from which Ice-{\it i} was crystallized
130 < (soft sticky dipole extended, SSD/E) in addition to several common
131 < water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field
132 < parametrized single point dipole water model (soft sticky dipole
133 < reaction field, SSD/RF). In should be noted that a second version of
134 < Ice-{\it i} (Ice-2{\it i}) was used in calculations involving SPC/E,
135 < TIP4P, and TIP5P. The unit cell of this crystal (Fig. \ref{iceiCell}B)
136 < is similar to the Ice-{\it i} unit it is extended in the direction of
137 < the (001) face and compressed along the other two faces.
106 > \begin{figure}
107 > \includegraphics[width=\linewidth]{unitCell.eps}
108 > \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
109 > elongated variant of Ice-{\it i}.  The spheres represent the
110 > center-of-mass locations of the water molecules.  The $a$ to $c$
111 > ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
112 > $a:2.1214c$ and $a:1.7850c$ respectively.}
113 > \label{iceiCell}
114 > \end{figure}
115  
116 + \begin{figure}
117 + \includegraphics[width=\linewidth]{orderedIcei.eps}
118 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
119 + down the (001) crystal face. The rows of water tetramers surrounded by
120 + octagonal pores leads to a crystal structure that is significantly
121 + less dense than ice $I_h$.}
122 + \label{protOrder}
123 + \end{figure}
124 +
125 + Results from our previous study indicated that Ice-{\it i} is the
126 + minimum energy crystal structure for the single point water models we
127 + had investigated (for discussions on these single point dipole models,
128 + see our previous work and related
129 + articles).\cite{Fennell04,Liu96,Bratko85} Those results only
130 + considered energetic stabilization and neglected entropic
131 + contributions to the overall free energy. To address this issue, we
132 + have calculated the absolute free energy of this crystal using
133 + thermodynamic integration and compared to the free energies of cubic
134 + and hexagonal ice $I$ (the experimental low density ice polymorphs)
135 + and ice B (a higher density, but very stable crystal structure
136 + observed by B\`{a}ez and Clancy in free energy studies of
137 + SPC/E).\cite{Baez95b} This work includes results for the water model
138 + from which Ice-{\it i} was crystallized (SSD/E) in addition to several
139 + common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
140 + field parametrized single point dipole water model (SSD/RF). It should
141 + be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was
142 + used in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell
143 + of this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i}
144 + unit it is extended in the direction of the (001) face and compressed
145 + along the other two faces.
146 +
147   \section{Methods}
148  
149   Canonical ensemble (NVT) molecular dynamics calculations were
150 < performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
151 < molecular mechanics package. All molecules were treated as rigid
152 < bodies, with orientational motion propagated using the symplectic DLM
153 < integration method. Details about the implementation of these
154 < techniques can be found in a recent publication.\cite{Meineke05}
150 > performed using the OOPSE molecular mechanics package.\cite{Meineke05}
151 > All molecules were treated as rigid bodies, with orientational motion
152 > propagated using the symplectic DLM integration method. Details about
153 > the implementation of this technique can be found in a recent
154 > publication.\cite{Dullweber1997}
155  
156 < Thermodynamic integration was utilized to calculate the free energy of
157 < several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
158 < SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
159 < 400 K for all of these water models were also determined using this
160 < same technique, in order to determine melting points and generate
161 < phase diagrams. All simulations were carried out at densities
162 < resulting in a pressure of approximately 1 atm at their respective
163 < temperatures.
156 > Thermodynamic integration is an established technique for
157 > determination of free energies of condensed phases of
158 > materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
159 > method, implemented in the same manner illustrated by B\`{a}ez and
160 > Clancy, was utilized to calculate the free energy of several ice
161 > crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
162 > SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
163 > and 400 K for all of these water models were also determined using
164 > this same technique in order to determine melting points and to
165 > generate phase diagrams. All simulations were carried out at densities
166 > which correspond to a pressure of approximately 1 atm at their
167 > respective temperatures.
168  
169 < A single thermodynamic integration involves a sequence of simulations
170 < over which the system of interest is converted into a reference system
171 < for which the free energy is known. This transformation path is then
172 < integrated in order to determine the free energy difference between
173 < the two states:
169 > Thermodynamic integration involves a sequence of simulations during
170 > which the system of interest is converted into a reference system for
171 > which the free energy is known analytically. This transformation path
172 > is then integrated in order to determine the free energy difference
173 > between the two states:
174   \begin{equation}
175   \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
176   )}{\partial\lambda}\right\rangle_\lambda d\lambda,
177   \end{equation}
178   where $V$ is the interaction potential and $\lambda$ is the
179   transformation parameter that scales the overall
180 < potential. Simulations are distributed unevenly along this path in
181 < order to sufficiently sample the regions of greatest change in the
180 > potential. Simulations are distributed strategically along this path
181 > in order to sufficiently sample the regions of greatest change in the
182   potential. Typical integrations in this study consisted of $\sim$25
183   simulations ranging from 300 ps (for the unaltered system) to 75 ps
184   (near the reference state) in length.
185  
186   For the thermodynamic integration of molecular crystals, the Einstein
187 < Crystal is chosen as the reference state that the system is converted
188 < to over the course of the simulation. In an Einstein Crystal, the
189 < molecules are harmonically restrained at their ideal lattice locations
190 < and orientations. The partition function for a molecular crystal
191 < restrained in this fashion has been evaluated, and the Helmholtz Free
192 < Energy ({\it A}) is given by
187 > crystal was chosen as the reference system. In an Einstein crystal,
188 > the molecules are restrained at their ideal lattice locations and
189 > orientations. Using harmonic restraints, as applied by B\`{a}ez and
190 > Clancy, the total potential for this reference crystal
191 > ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
192 > \begin{equation}
193 > V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
194 > \frac{K_\omega\omega^2}{2},
195 > \end{equation}
196 > where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
197 > the spring constants restraining translational motion and deflection
198 > of and rotation around the principle axis of the molecule
199 > respectively.  It is clear from Fig. \ref{waterSpring} that the values
200 > of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from
201 > $-\pi$ to $\pi$.  The partition function for a molecular crystal
202 > restrained in this fashion can be evaluated analytically, and the
203 > Helmholtz Free Energy ({\it A}) is given by
204   \begin{eqnarray}
205   A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
206   [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
# Line 189 | Line 212 | where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a
212   )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
213   \label{ecFreeEnergy}
214   \end{eqnarray}
215 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
216 < \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
217 < $K_\mathrm{\omega}$ are the spring constants restraining translational
195 < motion and deflection of and rotation around the principle axis of the
196 < molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
197 < minimum potential energy of the ideal crystal. In the case of
198 < molecular liquids, the ideal vapor is chosen as the target reference
199 < state.
215 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
216 > potential energy of the ideal crystal.\cite{Baez95a}
217 >
218   \begin{figure}
219 < \includegraphics[scale=1.0]{rotSpring.eps}
219 > \includegraphics[width=\linewidth]{rotSpring.eps}
220   \caption{Possible orientational motions for a restrained molecule.
221   $\theta$ angles correspond to displacement from the body-frame {\it
222   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 208 | Line 226 | Charge, dipole, and Lennard-Jones interactions were mo
226   \label{waterSpring}
227   \end{figure}
228  
229 + In the case of molecular liquids, the ideal vapor is chosen as the
230 + target reference state.  There are several examples of liquid state
231 + free energy calculations of water models present in the
232 + literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
233 + typically differ in regard to the path taken for switching off the
234 + interaction potential to convert the system to an ideal gas of water
235 + molecules.  In this study, we applied of one of the most convenient
236 + methods and integrated over the $\lambda^4$ path, where all
237 + interaction parameters are scaled equally by this transformation
238 + parameter.  This method has been shown to be reversible and provide
239 + results in excellent agreement with other established
240 + methods.\cite{Baez95b}
241 +
242   Charge, dipole, and Lennard-Jones interactions were modified by a
243 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
244 < applying this function, these interactions are smoothly truncated,
245 < thereby avoiding poor energy conserving dynamics resulting from
246 < harsher truncation schemes. The effect of a long-range correction was
247 < also investigated on select model systems in a variety of manners. For
248 < the SSD/RF model, a reaction field with a fixed dielectric constant of
249 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
250 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
251 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
252 < compare with the 9 \AA\ cutoff results. Finally, results from the use
253 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
254 < calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
255 < mechanics software package. TINKER was chosen because it can also
256 < propagate the motion of rigid-bodies, and provides the most direct
257 < comparison to the results from OOPSE. The calculated energy difference
258 < in the presence and absence of PME was applied to the previous results
259 < in order to predict changes in the free energy landscape.
243 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
244 > ). By applying this function, these interactions are smoothly
245 > truncated, thereby avoiding the poor energy conservation which results
246 > from harsher truncation schemes. The effect of a long-range correction
247 > was also investigated on select model systems in a variety of
248 > manners. For the SSD/RF model, a reaction field with a fixed
249 > dielectric constant of 80 was applied in all
250 > simulations.\cite{Onsager36} For a series of the least computationally
251 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
252 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
253 > \AA\ cutoff results. Finally, the effects of utilizing an Ewald
254 > summation were estimated for TIP3P and SPC/E by performing single
255 > configuration calculations with Particle-Mesh Ewald (PME) in the
256 > TINKER molecular mechanics software package.\cite{Tinker} The
257 > calculated energy difference in the presence and absence of PME was
258 > applied to the previous results in order to predict changes to the
259 > free energy landscape.
260  
261   \section{Results and discussion}
262  
263 < The free energy of proton ordered Ice-{\it i} was calculated and
263 > The free energy of proton-ordered Ice-{\it i} was calculated and
264   compared with the free energies of proton ordered variants of the
265   experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
266   as well as the higher density ice B, observed by B\`{a}ez and Clancy
267   and thought to be the minimum free energy structure for the SPC/E
268   model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
269 < Ice XI, the experimentally observed proton ordered variant of ice
270 < $I_h$, was investigated initially, but it was found not to be as
271 < stable as antiferroelectric variants of proton ordered or even proton
272 < disordered ice$I_h$.\cite{Davidson84} The proton ordered variant of
273 < ice $I_h$ used here is a simple antiferroelectric version that has an
274 < 8 molecule unit cell. The crystals contained 648 or 1728 molecules for
275 < ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice
276 < $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
277 < were necessary for simulations involving larger cutoff values.
269 > Ice XI, the experimentally-observed proton-ordered variant of ice
270 > $I_h$, was investigated initially, but was found to be not as stable
271 > as proton disordered or antiferroelectric variants of ice $I_h$. The
272 > proton ordered variant of ice $I_h$ used here is a simple
273 > antiferroelectric version that we devised, and it has an 8 molecule
274 > unit cell similar to other predicted antiferroelectric $I_h$
275 > crystals.\cite{Davidson84} The crystals contained 648 or 1728
276 > molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
277 > molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger
278 > crystal sizes were necessary for simulations involving larger cutoff
279 > values.
280  
281   \begin{table*}
282   \begin{minipage}{\linewidth}
# Line 252 | Line 285 | kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF
285   \caption{Calculated free energies for several ice polymorphs with a
286   variety of common water models. All calculations used a cutoff radius
287   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
288 < kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
288 > kcal/mol. Calculated error of the final digits is in parentheses. *Ice
289 > $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
290   \begin{tabular}{ l  c  c  c  c }
291 < \hline \\[-7mm]
292 < \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
293 < \hline \\[-3mm]
294 < \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
295 < \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
296 < \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
297 < \ \quad \ SPC/E  & \ \quad \ -12.67 & \ \quad \ -12.96 & \ \quad \ -13.25 & \quad -13.55\\
298 < \ \quad \ SSD/E  & \ \quad \ -11.27 & \ \quad \ -11.19 & \ \quad \ -12.09 & \quad -12.54\\
299 < \ \quad \ SSD/RF & \ \quad \ -11.51 & \ \quad \ NA* & \ \quad \ -12.08 & \quad -12.29\\
291 > \hline
292 > Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
293 > \hline
294 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3)\\
295 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & -12.33(3)\\
296 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & -12.29(2)\\
297 > SPC/E & -12.67(2) & -12.96(2) & -13.25(3) & -13.55(2)\\
298 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2)\\
299 > SSD/RF & -11.51(2) & NA* & -12.08(3) & -12.29(2)\\
300   \end{tabular}
301   \label{freeEnergy}
302   \end{center}
# Line 271 | Line 305 | models studied. With the free energy at these state po
305  
306   The free energy values computed for the studied polymorphs indicate
307   that Ice-{\it i} is the most stable state for all of the common water
308 < models studied. With the free energy at these state points, the
309 < temperature and pressure dependence of the free energy was used to
310 < project to other state points and build phase diagrams. Figures
308 > models studied. With the calculated free energy at these state points,
309 > the Gibbs-Helmholtz equation was used to project to other state points
310 > and to build phase diagrams.  Figures
311   \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
312   from the free energy results. All other models have similar structure,
313 < only the crossing points between these phases exist at different
314 < temperatures and pressures. It is interesting to note that ice $I$
315 < does not exist in either cubic or hexagonal form in any of the phase
316 < diagrams for any of the models. For purposes of this study, ice B is
317 < representative of the dense ice polymorphs. A recent study by Sanz
318 < {\it et al.} goes into detail on the phase diagrams for SPC/E and
319 < TIP4P in the high pressure regime.\cite{Sanz04}
313 > although the crossing points between the phases move to slightly
314 > different temperatures and pressures. It is interesting to note that
315 > ice $I$ does not exist in either cubic or hexagonal form in any of the
316 > phase diagrams for any of the models. For purposes of this study, ice
317 > B is representative of the dense ice polymorphs. A recent study by
318 > Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
319 > TIP4P at higher pressures than those studied here.\cite{Sanz04}
320 >
321   \begin{figure}
322   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
323   \caption{Phase diagram for the TIP3P water model in the low pressure
# Line 292 | Line 327 | higher in energy and don't appear in the phase diagram
327   higher in energy and don't appear in the phase diagram.}
328   \label{tp3phasedia}
329   \end{figure}
330 +
331   \begin{figure}
332   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
333   \caption{Phase diagram for the SSD/RF water model in the low pressure
# Line 308 | Line 344 | temperatures of several common water models compared w
344   \renewcommand{\thefootnote}{\thempfootnote}
345   \begin{center}
346   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
347 < temperatures of several common water models compared with experiment.}
347 > temperatures at 1 atm for several common water models compared with
348 > experiment. The $T_m$ and $T_s$ values from simulation correspond to a
349 > transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
350 > liquid or gas state.}
351   \begin{tabular}{ l  c  c  c  c  c  c  c }
352 < \hline \\[-7mm]
353 < \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
354 < \hline \\[-3mm]
355 < \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
356 < \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
357 < \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
352 > \hline
353 > Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
354 > \hline
355 > $T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
356 > $T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
357 > $T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
358   \end{tabular}
359   \label{meltandboil}
360   \end{center}
# Line 331 | Line 370 | ordered and disordered molecular arrangements). If the
370   studies in the literature. Earlier free energy studies of ice $I$
371   using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
372   being attributed to choice of interaction truncation and different
373 < ordered and disordered molecular arrangements). If the presence of ice
374 < B and Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
373 > ordered and disordered molecular
374 > arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
375 > Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
376   predicted from this work. However, the $T_m$ from Ice-{\it i} is
377   calculated at 265 K, significantly higher in temperature than the
378   previous studies. Also of interest in these results is that SSD/E does
379   not exhibit a melting point at 1 atm, but it shows a sublimation point
380   at 355 K. This is due to the significant stability of Ice-{\it i} over
381   all other polymorphs for this particular model under these
382 < conditions. While troubling, this behavior turned out to be
383 < advantageous in that it facilitated the spontaneous crystallization of
384 < Ice-{\it i}. These observations provide a warning that simulations of
382 > conditions. While troubling, this behavior resulted in spontaneous
383 > crystallization of Ice-{\it i} and led us to investigate this
384 > structure. These observations provide a warning that simulations of
385   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
386 < risk of spontaneous crystallization. However, this risk changes when
386 > risk of spontaneous crystallization. However, this risk lessens when
387   applying a longer cutoff.
388  
389   \begin{figure}
# Line 351 | Line 391 | TIP3P, and (C) SSD/RF. Data points omitted include SSD
391   \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
392   TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
393   \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
394 < \AA\. These crystals are unstable at 200 K and rapidly convert into a
395 < liquid. The connecting lines are qualitative visual aid.}
394 > \AA . These crystals are unstable at 200 K and rapidly convert into
395 > liquids. The connecting lines are qualitative visual aid.}
396   \label{incCutoff}
397   \end{figure}
398  
# Line 367 | Line 407 | greater than 9 \AA\. This narrowing trend is much more
407   in the SSD/E model that the liquid state is preferred under standard
408   simulation conditions (298 K and 1 atm). Thus, it is recommended that
409   simulations using this model choose interaction truncation radii
410 < greater than 9 \AA\. This narrowing trend is much more subtle in the
410 > greater than 9 \AA\ . This narrowing trend is much more subtle in the
411   case of SSD/RF, indicating that the free energies calculated with a
412   reaction field present provide a more accurate picture of the free
413   energy landscape in the absence of potential truncation.
# Line 377 | Line 417 | energy of identical crystals with and without PME usin
417   was estimated by applying the potential energy difference do to its
418   inclusion in systems in the presence and absence of the
419   correction. This was accomplished by calculation of the potential
420 < energy of identical crystals with and without PME using TINKER. The
421 < free energies for the investigated polymorphs using the TIP3P and
422 < SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
423 < are not fully supported in TINKER, so the results for these models
424 < could not be estimated. The same trend pointed out through increase of
425 < cutoff radius is observed in these PME results. Ice-{\it i} is the
426 < preferred polymorph at ambient conditions for both the TIP3P and SPC/E
427 < water models; however, there is a narrowing of the free energy
428 < differences between the various solid forms. In the case of SPC/E this
429 < narrowing is significant enough that it becomes less clear cut that
430 < Ice-{\it i} is the most stable polymorph, and is possibly metastable
431 < with respect to ice B and possibly ice $I_c$. However, these results
432 < do not significantly alter the finding that the Ice-{\it i} polymorph
433 < is a stable crystal structure that should be considered when studying
434 < the phase behavior of water models.
420 > energy of identical crystals both with and without PME. The free
421 > energies for the investigated polymorphs using the TIP3P and SPC/E
422 > water models are shown in Table \ref{pmeShift}. The same trend pointed
423 > out through increase of cutoff radius is observed in these PME
424 > results. Ice-{\it i} is the preferred polymorph at ambient conditions
425 > for both the TIP3P and SPC/E water models; however, the narrowing of
426 > the free energy differences between the various solid forms is
427 > significant enough that it becomes less clear that it is the most
428 > stable polymorph with the SPC/E model.  The free energies of Ice-{\it
429 > i} and ice B nearly overlap within error, with ice $I_c$ just outside
430 > as well, indicating that Ice-{\it i} might be metastable with respect
431 > to ice B and possibly ice $I_c$ with SPC/E. However, these results do
432 > not significantly alter the finding that the Ice-{\it i} polymorph is
433 > a stable crystal structure that should be considered when studying the
434 > phase behavior of water models.
435  
436   \begin{table*}
437   \begin{minipage}{\linewidth}
# Line 401 | Line 441 | long-range interaction correction. Units are kcal/mol.
441   the energy difference attributed to the inclusion of the PME
442   long-range interaction correction. Units are kcal/mol.}
443   \begin{tabular}{ l  c  c  c  c }
444 < \hline \\[-7mm]
444 > \hline
445   \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
446 < \hline \\[-3mm]
447 < \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
448 < \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
446 > \hline
447 > TIP3P  & -11.53(2) & -11.24(3) & -11.51(3) & -11.67(3)\\
448 > SPC/E  & -12.77(2) & -12.92(2) & -12.96(3) & -13.02(2)\\
449   \end{tabular}
450   \label{pmeShift}
451   \end{center}
# Line 415 | Line 455 | $I$, ice B, and recently discovered Ice-{\it i} where
455   \section{Conclusions}
456  
457   The free energy for proton ordered variants of hexagonal and cubic ice
458 < $I$, ice B, and recently discovered Ice-{\it i} where calculated under
459 < standard conditions for several common water models via thermodynamic
460 < integration. All the water models studied show Ice-{\it i} to be the
461 < minimum free energy crystal structure in the with a 9 \AA\ switching
462 < function cutoff. Calculated melting and boiling points show
463 < surprisingly good agreement with the experimental values; however, the
464 < solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
465 < interaction truncation was investigated through variation of the
466 < cutoff radius, use of a reaction field parameterized model, and
467 < estimation of the results in the presence of the Ewald summation
468 < correction. Interaction truncation has a significant effect on the
458 > $I$, ice B, and our recently discovered Ice-{\it i} structure were
459 > calculated under standard conditions for several common water models
460 > via thermodynamic integration. All the water models studied show
461 > Ice-{\it i} to be the minimum free energy crystal structure with a 9
462 > \AA\ switching function cutoff. Calculated melting and boiling points
463 > show surprisingly good agreement with the experimental values;
464 > however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
465 > effect of interaction truncation was investigated through variation of
466 > the cutoff radius, use of a reaction field parameterized model, and
467 > estimation of the results in the presence of the Ewald
468 > summation. Interaction truncation has a significant effect on the
469   computed free energy values, and may significantly alter the free
470   energy landscape for the more complex multipoint water models. Despite
471   these effects, these results show Ice-{\it i} to be an important ice
472   polymorph that should be considered in simulation studies.
473  
474 < Due to this relative stability of Ice-{\it i} in all manner of
475 < investigated simulation examples, the question arises as to possible
476 < experimental observation of this polymorph. The rather extensive past
474 > Due to this relative stability of Ice-{\it i} in all of the
475 > investigated simulation conditions, the question arises as to possible
476 > experimental observation of this polymorph.  The rather extensive past
477   and current experimental investigation of water in the low pressure
478 < regime leads the authors to be hesitant in ascribing relevance outside
479 < of computational models, hence the descriptive name presented. That
480 < being said, there are certain experimental conditions that would
481 < provide the most ideal situation for possible observation. These
482 < include the negative pressure or stretched solid regime, small
483 < clusters in vacuum deposition environments, and in clathrate
484 < structures involving small non-polar molecules.
478 > regime makes us hesitant to ascribe any relevance of this work outside
479 > of the simulation community.  It is for this reason that we chose a
480 > name for this polymorph which involves an imaginary quantity.  That
481 > said, there are certain experimental conditions that would provide the
482 > most ideal situation for possible observation. These include the
483 > negative pressure or stretched solid regime, small clusters in vacuum
484 > deposition environments, and in clathrate structures involving small
485 > non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
486 > our predictions for both the pair distribution function ($g_{OO}(r)$)
487 > and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
488 > ice-{\it i} at a temperature of 77K.  In studies of the high and low
489 > density forms of amorphous ice, ``spurious'' diffraction peaks have
490 > been observed experimentally.\cite{Bizid87} It is possible that a
491 > variant of Ice-{\it i} could explain some of this behavior; however,
492 > we will leave it to our experimental colleagues to make the final
493 > determination on whether this ice polymorph is named appropriately
494 > (i.e. with an imaginary number) or if it can be promoted to Ice-0.
495  
496 + \begin{figure}
497 + \includegraphics[width=\linewidth]{iceGofr.eps}
498 + \caption{Radial distribution functions of ice $I_h$, $I_c$,
499 + Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
500 + of the SSD/RF water model at 77 K.}
501 + \label{fig:gofr}
502 + \end{figure}
503 +
504 + \begin{figure}
505 + \includegraphics[width=\linewidth]{sofq.eps}
506 + \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
507 + and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
508 + been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
509 + width) to compensate for the trunction effects in our finite size
510 + simulations.}
511 + \label{fig:sofq}
512 + \end{figure}
513 +
514   \section{Acknowledgments}
515   Support for this project was provided by the National Science
516   Foundation under grant CHE-0134881. Computation time was provided by

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines