ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/fennellDissertation/iceChapter.tex
(Generate patch)

Comparing trunk/fennellDissertation/iceChapter.tex (file contents):
Revision 2977 by chrisfen, Sun Aug 27 15:24:39 2006 UTC vs.
Revision 2978 by chrisfen, Sun Aug 27 19:34:53 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}
2  
3 – Molecular dynamics is a valuable tool for studying the phase behavior
4 – of systems ranging from small or simple
5 – molecules\cite{Matsumoto02,andOthers} to complex biological
6 – species.\cite{bigStuff} Many techniques have been developed to
7 – investigate the thermodynamic properites of model substances,
8 – providing both qualitative and quantitative comparisons between
9 – simulations and experiment.\cite{thermMethods} Investigation of these
10 – properties leads to the development of new and more accurate models,
11 – leading to better understanding and depiction of physical processes
12 – and intricate molecular systems.
13 –
3   Water has proven to be a challenging substance to depict in
4   simulations, and a variety of models have been developed to describe
5   its behavior under varying simulation
6 < conditions.\cite{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
6 > conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04}
7   These models have been used to investigate important physical
8   phenomena like phase transitions and the hydrophobic
9 < effect.\cite{Yamada02} With the choice of models available, it
10 < is only natural to compare the models under interesting thermodynamic
11 < conditions in an attempt to clarify the limitations of each of the
12 < models.\cite{modelProps} Two important property to quantify are the
13 < Gibbs and Helmholtz free energies, particularly for the solid forms of
14 < water.  Difficulty in these types of studies typically arises from the
15 < assortment of possible crystalline polymorphs that water adopts over a
16 < wide range of pressures and temperatures. There are currently 13
17 < recognized forms of ice, and it is a challenging task to investigate
18 < the entire free energy landscape.\cite{Sanz04} Ideally, research is
9 > effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models
10 > available, it is only natural to compare the models under interesting
11 > thermodynamic conditions in an attempt to clarify the limitations of
12 > each of the models.\cite{Jorgensen83,Jorgensen98,Baez94,Mahoney01} Two
13 > important property to quantify are the Gibbs and Helmholtz free
14 > energies, particularly for the solid forms of water, as these predict
15 > the thermodynamic stability of the various phases. Water has a
16 > particularly rich phase diagram and takes on a number of different and
17 > stable crystalline structures as the temperature and pressure are
18 > varied. This complexity makes it a challenging task to investigate the
19 > entire free energy landscape.\cite{Sanz04} Ideally, research is
20   focused on the phases having the lowest free energy at a given state
21 < point, because these phases will dictate the true transition
22 < temperatures and pressures for their respective model.
21 > point, because these phases will dictate the relevant transition
22 > temperatures and pressures for the model.
23  
24 < In this paper, standard reference state methods were applied to known
25 < crystalline water polymorphs in the low pressure regime. This work is
26 < unique in the fact that one of the crystal lattices was arrived at
27 < through crystallization of a computationally efficient water model
28 < under constant pressure and temperature conditions. Crystallization
29 < events are interesting in and of
30 < themselves;\cite{Matsumoto02,Yamada02} however, the crystal structure
31 < obtained in this case is different from any previously observed ice
32 < polymorphs in experiment or simulation.\cite{Fennell04} We have named
33 < this structure Ice-{\it i} to indicate its origin in computational
34 < simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
35 < water molecules that stack in rows of interlocking water
36 < tetramers. Proton ordering can be accomplished by orienting two of the
37 < molecules so that both of their donated hydrogen bonds are internal to
38 < their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
39 < constructed of water tetramers, the hydrogen bonds are not as linear
40 < as those observed in ice $I_h$, however the interlocking of these
41 < subunits appears to provide significant stabilization to the overall
42 < crystal. The arrangement of these tetramers results in surrounding
43 < open octagonal cavities that are typically greater than 6.3 \AA\ in
44 < diameter. This relatively open overall structure leads to crystals
45 < that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
24 > The high-pressure phases of water (ice II-ice X as well as ice XII)
25 > have been studied extensively both experimentally and
26 > computatuionally. In this chapter, standard reference state methods
27 > were applied in the {\it low} pressure regime to evaluate the free
28 > energies for a few known crystalline water polymorphs that might be
29 > stable at these pressures. This work is unique in the fact that one of
30 > the crystal lattices was arrived at through crystallization of a
31 > computationally efficient water model under constant pressure and
32 > temperature conditions. Crystallization events are interesting in and
33 > of themselves;\cite{Matsumoto02,Yamada02} however, the crystal
34 > structure obtained in this case is different from any previously
35 > observed ice polymorphs in experiment or simulation.\cite{Fennell04}
36 > We have named this structure Ice-{\it i} to indicate its origin in
37 > computational simulation. The unit cell of Ice-$i$ and an axially
38 > elongated variant named Ice-$i^\prime$ both consist of eight water
39 > molecules that stack in rows of interlocking water tetramers as
40 > illustrated in figure \ref{fig:unitCell}A,B. These tetramers form a
41 > crystal structure similar in appearance to a recent two-dimensional
42 > surface tessellation simulated on silica.\cite{Yang04} As expected in
43 > an ice crystal constructed of water tetramers, the hydrogen bonds are
44 > not as linear as those observed in ice I$_\textrm{h}$; however, the
45 > interlocking of these subunits appears to provide significant
46 > stabilization to the overall crystal. The arrangement of these
47 > tetramers results in open octagonal cavities that are typically
48 > greater than 6.3\AA\ in diameter (see figure
49 > \ref{fig:protOrder}). This open structure leads to crystals that are
50 > typically 0.07 g/cm$^3$ less dense than ice I$_\textrm{h}$.
51  
52   \begin{figure}
53   \includegraphics[width=\linewidth]{./figures/unitCell.pdf}
# Line 60 | Line 55 | relation is given by $a = 2.1214c$, while for Ice-$i^\
55   elongated variant of Ice-{\it i}.  For Ice-{\it i}, the $a$ to $c$
56   relation is given by $a = 2.1214c$, while for Ice-$i^\prime$, $a =
57   1.7850c$.}
58 < \label{iceiCell}
58 > \label{fig:iceiCell}
59   \end{figure}
60  
61   \begin{figure}
# Line 69 | Line 64 | octagonal pores leads to a crystal structure that is s
64   \caption{Image of a proton ordered crystal of Ice-{\it i} looking
65   down the (001) crystal face. The rows of water tetramers surrounded by
66   octagonal pores leads to a crystal structure that is significantly
67 < less dense than ice $I_h$.}
68 < \label{protOrder}
67 > less dense than ice I$_\textrm{h}$.}
68 > \label{fig:protOrder}
69   \end{figure}
70  
71   Results from our previous study indicated that Ice-{\it i} is the
72 < minimum energy crystal structure for the single point water models we
72 > minimum energy crystal structure for the single point water models
73   investigated (for discussions on these single point dipole models, see
74   the previous work and related
75 < articles\cite{Fennell04,Liu96,Bratko85}). Those results only
75 > articles\cite{Fennell04,Liu96,Bratko85}). Our earlier results only
76   considered energetic stabilization and neglected entropic
77 < contributions to the overall free energy. To address this issue, the
78 < absolute free energy of this crystal was calculated using
79 < thermodynamic integration and compared to the free energies of cubic
80 < and hexagonal ice $I$ (the experimental low density ice polymorphs)
81 < and ice B (a higher density, but very stable crystal structure
82 < observed by B\`{a}ez and Clancy in free energy studies of
77 > contributions to the overall free energy. To address this issue, we
78 > have calculated the absolute free energy of this crystal using
79 > thermodynamic integration and compared to the free energies of ice
80 > I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common low-density ice
81 > polymorphs) and ice B (a higher density, but very stable crystal
82 > structure observed by B\`{a}ez and Clancy in free energy studies of
83   SPC/E).\cite{Baez95b} This work includes results for the water model
84   from which Ice-{\it i} was crystallized (SSD/E) in addition to several
85   common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
86 < field parametrized single point dipole water model (SSD/RF). It should
87 < be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used
88 < in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of
89 < this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit
90 < it is extended in the direction of the (001) face and compressed along
91 < the other two faces.
86 > field parametrized single point dipole water model (SSD/RF). The
87 > axially elongated variant, Ice-$i^\prime$, was used in calculations
88 > involving SPC/E, TIP4P, and TIP5P. The square tetramer in Ice-$i$
89 > distorts in Ice-$i^\prime$ to form a rhombus with alternating 85 and
90 > 95 degree angles. Under SPC/E, TIP4P, and TIP5P, this geometry is
91 > better at forming favorable hydrogen bonds. The degree of rhomboid
92 > distortion depends on the water model used but is significant enough
93 > to split the peak in the radial distribution function which corresponds
94 > to diagonal sites in the tetramers.
95  
96   \section{Methods}
97  
98 < Canonical ensemble (NVT) molecular dynamics calculations were
98 > Canonical ensemble ({\it NVT}) molecular dynamics calculations were
99   performed using the OOPSE molecular mechanics package.\cite{Meineke05}
100 < All molecules were treated as rigid bodies, with orientational motion
101 < propagated using the symplectic DLM integration method described in
102 < section \ref{sec:IntroIntegration}.
100 > The densities chosen for the simulations were taken from
101 > isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and
102 > 200K. Each model (and each crystal structure) was allowed to relax for
103 > 300ps in the {\it NPT} ensemble before averaging the density to obtain
104 > the volumes for the {\it NVT} simulations.All molecules were treated
105 > as rigid bodies, with orientational motion propagated using the
106 > symplectic DLM integration method described in section
107 > \ref{sec:IntroIntegration}.
108  
109 < Thermodynamic integration was utilized to calculate the free energy of
110 < several ice crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E,
111 < SSD/RF, and SSD/E water models. Liquid state free energies at 300 and
112 < 400 K for all of these water models were also determined using this
113 < same technique in order to determine melting points and generate phase
114 < diagrams. All simulations were carried out at densities resulting in a
115 < pressure of approximately 1 atm at their respective temperatures.
116 <
117 < A single thermodynamic integration involves a sequence of simulations
118 < over which the system of interest is converted into a reference system
119 < for which the free energy is known analytically. This transformation
117 < path is then integrated in order to determine the free energy
118 < difference between the two states:
109 > We used thermodynamic integration to calculate the Helmholtz free
110 > energies ({\it A}) of the listed water models at various state
111 > points. Thermodynamic integration is an established technique that has
112 > been used extensively in the calculation of free energies for
113 > condensed phases of
114 > materials.\cite{Frenkel84,Hermans88,Meijer90,Baez95,Vlot99} This
115 > method uses a sequence of simulations during which the system of
116 > interest is converted into a reference system for which the free
117 > energy is known analytically ($A_0$). This transformation path is then
118 > integrated in order to determine the free energy difference between
119 > the two states:
120   \begin{equation}
121   \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
122   )}{\partial\lambda}\right\rangle_\lambda d\lambda,
# Line 125 | Line 126 | potential. Typical integrations in this study consiste
126   potential. Simulations are distributed unevenly along this path in
127   order to sufficiently sample the regions of greatest change in the
128   potential. Typical integrations in this study consisted of $\sim$25
129 < simulations ranging from 300 ps (for the unaltered system) to 75 ps
129 > simulations ranging from 300ps (for the unaltered system) to 75ps
130   (near the reference state) in length.
131  
132   For the thermodynamic integration of molecular crystals, the Einstein
# Line 143 | Line 144 | K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
144   K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
145   (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
146   )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
147 < \label{ecFreeEnergy}
147 > \label{eq:ecFreeEnergy}
148   \end{eqnarray}
149   where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
150 < \ref{ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
150 > \ref{eq:ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
151   $K_\mathrm{\omega}$ are the spring constants restraining translational
152   motion and deflection of and rotation around the principle axis of the
153   molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
# Line 167 | Line 168 | Charge, dipole, and Lennard-Jones interactions were mo
168   \end{figure}
169  
170   Charge, dipole, and Lennard-Jones interactions were modified by a
171 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA
172 < ). By applying this function, these interactions are smoothly
173 < truncated, thereby avoiding the poor energy conservation which results
174 < from harsher truncation schemes. The effect of a long-range correction
175 < was also investigated on select model systems in a variety of
176 < manners. For the SSD/RF model, a reaction field with a fixed
177 < dielectric constant of 80 was applied in all
178 < simulations.\cite{Onsager36} For a series of the least computationally
179 < expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
180 < performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
181 < \AA\ cutoff results. Finally, results from the use of an Ewald
181 < summation were estimated for TIP3P and SPC/E by performing
171 > cubic switching between 100\% and 85\% of the cutoff value (9\AA). By
172 > applying this function, these interactions are smoothly truncated,
173 > thereby avoiding the poor energy conservation which results from
174 > harsher truncation schemes. The effect of a long-range correction was
175 > also investigated on select model systems in a variety of manners. For
176 > the SSD/RF model, a reaction field with a fixed dielectric constant of
177 > 80 was applied in all simulations.\cite{Onsager36} For a series of the
178 > least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
179 > simulations were performed with longer cutoffs of 12 and 15\AA\ to
180 > compare with the 9\AA\ cutoff results. Finally, results from the use
181 > of an Ewald summation were estimated for TIP3P and SPC/E by performing
182   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
183   mechanics software package.\cite{Tinker} The calculated energy
184   difference in the presence and absence of PME was applied to the
# Line 189 | Line 189 | compared with the free energies of proton ordered vari
189  
190   The free energy of proton ordered Ice-{\it i} was calculated and
191   compared with the free energies of proton ordered variants of the
192 < experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
193 < as well as the higher density ice B, observed by B\`{a}ez and Clancy
194 < and thought to be the minimum free energy structure for the SPC/E
195 < model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
196 < Ice XI, the experimentally-observed proton-ordered variant of ice
197 < $I_h$, was investigated initially, but was found to be not as stable
198 < as proton disordered or antiferroelectric variants of ice $I_h$. The
199 < proton ordered variant of ice $I_h$ used here is a simple
200 < antiferroelectric version that has an 8 molecule unit
192 > experimentally observed low-density ice polymorphs, I$_\textrm{h}$ and
193 > I$_\textrm{c}$, as well as the higher density ice B, observed by
194 > B\`{a}ez and Clancy and thought to be the minimum free energy
195 > structure for the SPC/E model at ambient conditions.\cite{Baez95b} Ice
196 > XI, the experimentally-observed proton-ordered variant of ice
197 > I$_\textrm{h}$, was investigated initially, but was found to be not as
198 > stable as proton disordered or antiferroelectric variants of ice
199 > I$_\textrm{h}$. The proton ordered variant of ice I$_\textrm{h}$ used
200 > here is a simple antiferroelectric version that has an 8 molecule unit
201   cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules
202 < for ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for
203 < ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes
204 < were necessary for simulations involving larger cutoff values.
202 > for ice B, 1024 or 1280 molecules for ice I$_\textrm{h}$, 1000
203 > molecules for ice I$_\textrm{c}$, or 1024 molecules for Ice-{\it
204 > i}. The larger crystal sizes were necessary for simulations involving
205 > larger cutoff values.
206  
207 < \begin{table*}
208 < \begin{minipage}{\linewidth}
209 < \renewcommand{\thefootnote}{\thempfootnote}
210 < \begin{center}
210 < \caption{Calculated free energies for several ice polymorphs with a
211 < variety of common water models. All calculations used a cutoff radius
212 < of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
213 < kcal/mol. Calculated error of the final digits is in parentheses. *Ice
214 < $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.}
207 > \begin{table}
208 > \centering
209 > \caption{HELMHOLTZ FREE ENERGIES FOR SEVERAL ICE POLYMORPHS WITH A
210 > VARIETY OF COMMON WATER MODELS AT 200 KELVIN AND 1 ATMOSPHERE}
211   \begin{tabular}{ l  c  c  c  c }
212 < \hline
213 < Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\
214 < \hline
212 > \toprule
213 > \toprule
214 > Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or $i^\prime$) \\
215 > & kcal/mol & kcal/mol & kcal/mol & kcal/mol \\
216 > \midrule
217   TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\
218   TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\
219   TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\
220   SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\
221   SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\
222 < SSD/RF & -11.51(4) & NA* & -12.08(5) & -12.29(4)\\
222 > SSD/RF & -11.51(4) & NA & -12.08(5) & -12.29(4)\\
223 > \bottomrule
224   \end{tabular}
225 < \label{freeEnergy}
226 < \end{center}
228 < \end{minipage}
229 < \end{table*}
225 > \label{tab:freeEnergy}
226 > \end{table}
227  
228 < The free energy values computed for the studied polymorphs indicate
229 < that Ice-{\it i} is the most stable state for all of the common water
230 < models studied. With the free energy at these state points, the
231 < Gibbs-Helmholtz equation was used to project to other state points and
232 < to build phase diagrams.  Figures
233 < \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
234 < from the free energy results. All other models have similar structure,
238 < although the crossing points between the phases exist at slightly
239 < different temperatures and pressures. It is interesting to note that
240 < ice $I$ does not exist in either cubic or hexagonal form in any of the
241 < phase diagrams for any of the models. For purposes of this study, ice
242 < B is representative of the dense ice polymorphs. A recent study by
243 < Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
244 < TIP4P in the high pressure regime.\cite{Sanz04}
228 > Table \ref{tab:freeEnergy} shows the results of the free energy
229 > calculations with a cutoff radius of 9\AA. It should be noted that the
230 > ice I$_\textrm{c}$ crystal polymorph is not stable at 200K and 1 atm
231 > with the SSD/RF water model, hense omitted results for that cell. The
232 > free energy values displayed in this table, it is clear that Ice-{\it
233 > i} (or Ice-$i^\prime$ for TIP4P, TIP5P, and SPC/E) is the most stable
234 > state for all of the common water models studied.
235  
236 + With the free energy at these state points, the Gibbs-Helmholtz
237 + equation was used to project to other state points and to build phase
238 + diagrams.  Figures \ref{fig:tp3phasedia} and \ref{fig:ssdrfphasedia}
239 + are example diagrams built from the free energy results. All other
240 + models have similar structure, although the crossing points between
241 + the phases exist at slightly different temperatures and pressures. It
242 + is interesting to note that ice I does not exist in either cubic or
243 + hexagonal form in any of the phase diagrams for any of the models. For
244 + purposes of this study, ice B is representative of the dense ice
245 + polymorphs. A recent study by Sanz {\it et al.} goes into detail on
246 + the phase diagrams for SPC/E and TIP4P in the high pressure
247 + regime.\cite{Sanz04}
248 +
249   \begin{figure}
250   \centering
251   \includegraphics[width=\linewidth]{./figures/tp3PhaseDia.pdf}
# Line 251 | Line 254 | higher in energy and don't appear in the phase diagram
254   the experimental values; however, the solid phases shown are not the
255   experimentally observed forms. Both cubic and hexagonal ice $I$ are
256   higher in energy and don't appear in the phase diagram.}
257 < \label{tp3phasedia}
257 > \label{fig:tp3phasedia}
258   \end{figure}
259  
260   \begin{figure}
# Line 263 | Line 266 | conservative charge based models.}
266   computationally efficient model (over 3 times more efficient than
267   TIP3P) exhibits phase behavior similar to the less computationally
268   conservative charge based models.}
269 < \label{ssdrfphasedia}
269 > \label{fig:ssdrfphasedia}
270   \end{figure}
271  
272 < \begin{table*}
273 < \begin{minipage}{\linewidth}
271 < \renewcommand{\thefootnote}{\thempfootnote}
272 < \begin{center}
272 > \begin{table}
273 > \centering
274   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
275   temperatures at 1 atm for several common water models compared with
276   experiment. The $T_m$ and $T_s$ values from simulation correspond to a
277   transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
278   liquid or gas state.}
279   \begin{tabular}{ l  c  c  c  c  c  c  c }
280 < \hline
280 > \toprule
281 > \toprule
282   Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
283 < \hline
283 > \midrule
284   $T_m$ (K)  & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\
285   $T_b$ (K)  & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\
286   $T_s$ (K)  & - & - & - & - & 355(3) & - & -\\
287 + \bottomrule
288   \end{tabular}
289 < \label{meltandboil}
290 < \end{center}
288 < \end{minipage}
289 < \end{table*}
289 > \label{tab:meltandboil}
290 > \end{table}
291  
292 < Table \ref{meltandboil} lists the melting and boiling temperatures
292 > Table \ref{tab:meltandboil} lists the melting and boiling temperatures
293   calculated from this work. Surprisingly, most of these models have
294   melting points that compare quite favorably with experiment. The
295   unfortunate aspect of this result is that this phase change occurs
296 < between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
297 < liquid state. These results are actually not contrary to previous
298 < studies in the literature. Earlier free energy studies of ice $I$
299 < using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences
300 < being attributed to choice of interaction truncation and different
301 < ordered and disordered molecular
296 > between Ice-{\it i} and the liquid state rather than ice
297 > I$_\textrm{h}$ and the liquid state. These results are actually not
298 > contrary to previous studies in the literature. Earlier free energy
299 > studies of ice I using TIP4P predict a $T_m$ ranging from 214 to 238K
300 > (differences being attributed to choice of interaction truncation and
301 > different ordered and disordered molecular
302   arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
303 < Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
304 < predicted from this work. However, the $T_m$ from Ice-{\it i} is
305 < calculated at 265 K, significantly higher in temperature than the
306 < previous studies. Also of interest in these results is that SSD/E does
307 < not exhibit a melting point at 1 atm, but it shows a sublimation point
308 < at 355 K. This is due to the significant stability of Ice-{\it i} over
309 < all other polymorphs for this particular model under these
303 > Ice-{\it i} were omitted, a $T_m$ value around 210K would be predicted
304 > from this work. However, the $T_m$ from Ice-{\it i} is calculated at
305 > 265K, significantly higher in temperature than the previous
306 > studies. Also of interest in these results is that SSD/E does not
307 > exhibit a melting point at 1 atm, but it shows a sublimation point at
308 > 355K. This is due to the significant stability of Ice-{\it i} over all
309 > other polymorphs for this particular model under these
310   conditions. While troubling, this behavior turned out to be
311   advantageous in that it facilitated the spontaneous crystallization of
312   Ice-{\it i}. These observations provide a warning that simulations of
313 < SSD/E as a ``liquid'' near 300 K are actually metastable and run the
313 > SSD/E as a ``liquid'' near 300K are actually metastable and run the
314   risk of spontaneous crystallization. However, this risk changes when
315   applying a longer cutoff.
316  
317   \begin{figure}
318   \includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
319   \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
320 < TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
321 < \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
322 < \AA . These crystals are unstable at 200 K and rapidly convert into
323 < liquids. The connecting lines are qualitative visual aid.}
324 < \label{incCutoff}
320 > TIP3P, and (C) SSD/RF. Data points omitted include SSD/E:
321 > I$_\textrm{c}$ 12\AA\, TIP3P: I$_\textrm{c}$ 12\AA\ and B 12\AA\,
322 > and SSD/RF: I$_\textrm{c}$ 9\AA . These crystals are unstable at 200 K
323 > and rapidly convert into liquids. The connecting lines are qualitative
324 > visual aid.}
325 > \label{fig:incCutoff}
326   \end{figure}
327  
328   Increasing the cutoff radius in simulations of the more
329   computationally efficient water models was done in order to evaluate
330   the trend in free energy values when moving to systems that do not
331 < involve potential truncation. As seen in Fig. \ref{incCutoff}, the
332 < free energy of all the ice polymorphs show a substantial dependence on
333 < cutoff radius. In general, there is a narrowing of the free energy
334 < differences while moving to greater cutoff radius. Interestingly, by
335 < increasing the cutoff radius, the free energy gap was narrowed enough
336 < in the SSD/E model that the liquid state is preferred under standard
337 < simulation conditions (298 K and 1 atm). Thus, it is recommended that
338 < simulations using this model choose interaction truncation radii
339 < greater than 9 \AA\ . This narrowing trend is much more subtle in the
340 < case of SSD/RF, indicating that the free energies calculated with a
341 < reaction field present provide a more accurate picture of the free
342 < energy landscape in the absence of potential truncation.
331 > involve potential truncation. As seen in figure \ref{fig:incCutoff},
332 > the free energy of all the ice polymorphs show a substantial
333 > dependence on cutoff radius. In general, there is a narrowing of the
334 > free energy differences while moving to greater cutoff
335 > radius. Interestingly, by increasing the cutoff radius, the free
336 > energy gap was narrowed enough in the SSD/E model that the liquid
337 > state is preferred under standard simulation conditions (298K and 1
338 > atm). Thus, it is recommended that simulations using this model choose
339 > interaction truncation radii greater than 9\AA\ . This narrowing
340 > trend is much more subtle in the case of SSD/RF, indicating that the
341 > free energies calculated with a reaction field present provide a more
342 > accurate picture of the free energy landscape in the absence of
343 > potential truncation.
344  
345   To further study the changes resulting to the inclusion of a
346   long-range interaction correction, the effect of an Ewald summation
# Line 346 | Line 349 | free energies for the investigated polymorphs using th
349   correction. This was accomplished by calculation of the potential
350   energy of identical crystals with and without PME using TINKER. The
351   free energies for the investigated polymorphs using the TIP3P and
352 < SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
353 < are not fully supported in TINKER, so the results for these models
354 < could not be estimated. The same trend pointed out through increase of
355 < cutoff radius is observed in these PME results. Ice-{\it i} is the
356 < preferred polymorph at ambient conditions for both the TIP3P and SPC/E
357 < water models; however, there is a narrowing of the free energy
358 < differences between the various solid forms. In the case of SPC/E this
359 < narrowing is significant enough that it becomes less clear that
360 < Ice-{\it i} is the most stable polymorph, and is possibly metastable
361 < with respect to ice B and possibly ice $I_c$. However, these results
362 < do not significantly alter the finding that the Ice-{\it i} polymorph
363 < is a stable crystal structure that should be considered when studying
364 < the phase behavior of water models.
352 > SPC/E water models are shown in Table \ref{tab:pmeShift}. TIP4P and
353 > TIP5P are not fully supported in TINKER, so the results for these
354 > models could not be estimated. The same trend pointed out through
355 > increase of cutoff radius is observed in these PME results. Ice-{\it
356 > i} is the preferred polymorph at ambient conditions for both the TIP3P
357 > and SPC/E water models; however, there is a narrowing of the free
358 > energy differences between the various solid forms. In the case of
359 > SPC/E this narrowing is significant enough that it becomes less clear
360 > that Ice-{\it i} is the most stable polymorph, and is possibly
361 > metastable with respect to ice B and possibly ice
362 > I$_\textrm{c}$. However, these results do not significantly alter the
363 > finding that the Ice-{\it i} polymorph is a stable crystal structure
364 > that should be considered when studying the phase behavior of water
365 > models.
366  
367 < \begin{table*}
368 < \begin{minipage}{\linewidth}
365 < \renewcommand{\thefootnote}{\thempfootnote}
366 < \begin{center}
367 > \begin{table}
368 > \centering
369   \caption{The free energy of the studied ice polymorphs after applying
370   the energy difference attributed to the inclusion of the PME
371   long-range interaction correction. Units are kcal/mol.}
372   \begin{tabular}{ l  c  c  c  c }
373 < \hline
374 < \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
375 < \hline
373 > \toprule
374 > \toprule
375 > Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or Ice-$i^\prime$) \\
376 > \midrule
377   TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\
378   SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\
379 + \bottomrule
380   \end{tabular}
381 < \label{pmeShift}
382 < \end{center}
379 < \end{minipage}
380 < \end{table*}
381 > \label{tab:pmeShift}
382 > \end{table}
383  
384   \section{Conclusions}
385  
# Line 385 | Line 387 | integration. All the water models studied show Ice-{\i
387   $I$, ice B, and recently discovered Ice-{\it i} were calculated under
388   standard conditions for several common water models via thermodynamic
389   integration. All the water models studied show Ice-{\it i} to be the
390 < minimum free energy crystal structure in the with a 9 \AA\ switching
390 > minimum free energy crystal structure in the with a 9\AA\ switching
391   function cutoff. Calculated melting and boiling points show
392   surprisingly good agreement with the experimental values; however, the
393 < solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
394 < interaction truncation was investigated through variation of the
395 < cutoff radius, use of a reaction field parameterized model, and
393 > solid phase at 1 atm is Ice-{\it i}, not ice I$_\textrm{h}$. The
394 > effect of interaction truncation was investigated through variation of
395 > the cutoff radius, use of a reaction field parameterized model, and
396   estimation of the results in the presence of the Ewald
397   summation. Interaction truncation has a significant effect on the
398   computed free energy values, and may significantly alter the free
# Line 411 | Line 413 | our predictions for both the pair distribution functio
413   deposition environments, and in clathrate structures involving small
414   non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
415   our predictions for both the pair distribution function ($g_{OO}(r)$)
416 < and the structure factor ($S(\vec{q})$ for ice $I_c$ and for ice-{\it
417 < i} at a temperature of 77K.  In a quick comparison of the predicted
418 < S(q) for Ice-{\it i} and experimental studies of amorphous solid
419 < water, it is possible that some of the ``spurious'' peaks that could
420 < not be assigned in HDA could correspond to peaks labeled in this
416 > and the structure factor ($S(\vec{q})$ for ice I$_\textrm{c}$ and for
417 > ice-{\it i} at a temperature of 77K.  In a quick comparison of the
418 > predicted S(q) for Ice-{\it i} and experimental studies of amorphous
419 > solid water, it is possible that some of the ``spurious'' peaks that
420 > could not be assigned in HDA could correspond to peaks labeled in this
421   S(q).\cite{Bizid87} It should be noted that there is typically poor
422   agreement on crystal densities between simulation and experiment, so
423   such peak comparisons should be made with caution.  We will leave it
# Line 424 | Line 426 | is named appropriately or if it should be promoted to
426  
427   \begin{figure}
428   \includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
429 < \caption{Radial distribution functions of Ice-{\it i} and ice $I_c$
430 < calculated from from simulations of the SSD/RF water model at 77 K.}
429 > \caption{Radial distribution functions of Ice-{\it i} and ice
430 > I$_\textrm{c}$ calculated from from simulations of the SSD/RF water
431 > model at 77 K.}
432   \label{fig:gofr}
433   \end{figure}
434  
435   \begin{figure}
436   \includegraphics[width=\linewidth]{./figures/sofq.pdf}
437 < \caption{Predicted structure factors for Ice-{\it i} and ice $I_c$ at
438 < 77 K.  The raw structure factors have been convoluted with a gaussian
439 < instrument function (0.075 \AA$^{-1}$ width) to compensate for the
440 < trunction effects in our finite size simulations. The labeled peaks
441 < compared favorably with ``spurious'' peaks observed in experimental
442 < studies of amorphous solid water.\cite{Bizid87}}
437 > \caption{Predicted structure factors for Ice-{\it i} and ice
438 > I$_\textrm{c}$ at 77 K.  The raw structure factors have been
439 > convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
440 > width) to compensate for the trunction effects in our finite size
441 > simulations. The labeled peaks compared favorably with ``spurious''
442 > peaks observed in experimental studies of amorphous solid
443 > water.\cite{Bizid87}}
444   \label{fig:sofq}
445   \end{figure}
446  

Diff Legend

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