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 1464 by chrisfen, Wed Sep 15 21:44:27 2004 UTC

# Line 1 | Line 1
1   %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 < \documentclass[preprint,aps,endfloats]{revtex4}
2 > \documentclass[11pt]{article}
3   %\documentclass[11pt]{article}
4 < %\usepackage{endfloat}
4 > \usepackage{endfloat}
5   \usepackage{amsmath}
6   \usepackage{epsf}
7   \usepackage{berkeley}
8 < %\usepackage{setspace}
9 < %\usepackage{tabularx}
8 > \usepackage{setspace}
9 > \usepackage{tabularx}
10   \usepackage{graphicx}
11 < %\usepackage[ref]{overcite}
12 < %\pagestyle{plain}
13 < %\pagenumbering{arabic}
14 < %\oddsidemargin 0.0cm \evensidemargin 0.0cm
15 < %\topmargin -21pt \headsep 10pt
16 < %\textheight 9.0in \textwidth 6.5in
17 < %\brokenpenalty=10000
11 > \usepackage[ref]{overcite}
12 > \pagestyle{plain}
13 > \pagenumbering{arabic}
14 > \oddsidemargin 0.0cm \evensidemargin 0.0cm
15 > \topmargin -21pt \headsep 10pt
16 > \textheight 9.0in \textwidth 6.5in
17 > \brokenpenalty=10000
18 > \renewcommand{\baselinestretch}{1.2}
19 > \renewcommand\citemid{\ } % no comma in optional reference note
20  
19 – %\renewcommand\citemid{\ } % no comma in optional reference note
20 –
21   \begin{document}
22  
23 < \title{A Free Energy Study of Low Temperature and Anomolous Ice}
23 > \title{Ice-{\it i}: a novel ice polymorph predicted via computer simulation}
24  
25 < \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
26 < \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
27 <
28 < \address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
25 > \author{Christopher J. Fennell and J. Daniel Gezelter \\
26 > Department of Chemistry and Biochemistry\\ University of Notre Dame\\
27   Notre Dame, Indiana 46556}
28  
29   \date{\today}
30  
31 < %\maketitle
31 > \maketitle
32   %\doublespacing
33  
34   \begin{abstract}
35 + The free energies of several ice polymorphs in the low pressure regime
36 + were calculated using thermodynamic integration.  These integrations
37 + were done for most of the common water models. Ice-{\it i}, a
38 + structure we recently observed to be stable in one of the single-point
39 + water models, was determined to be the stable crystalline state (at 1
40 + atm) for {\it all} the water models investigated.  Phase diagrams were
41 + generated, and phase coexistence lines were determined for all of the
42 + known low-pressure ice structures under all of the common water
43 + models.  Additionally, potential truncation was shown to have an
44 + effect on the calculated free energies, and can result in altered free
45 + energy landscapes.
46   \end{abstract}
47  
39 – \maketitle
40 –
41 – \newpage
42 –
48   %\narrowtext
49  
50   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 48 | Line 53 | Notre Dame, Indiana 46556}
53  
54   \section{Introduction}
55  
56 + Molecular dynamics is a valuable tool for studying the phase behavior
57 + of systems ranging from small or simple
58 + molecules\cite{Matsumoto02andOthers} to complex biological
59 + species.\cite{bigStuff} Many techniques have been developed to
60 + investigate the thermodynamic properites of model substances,
61 + providing both qualitative and quantitative comparisons between
62 + simulations and experiment.\cite{thermMethods} Investigation of these
63 + properties leads to the development of new and more accurate models,
64 + leading to better understanding and depiction of physical processes
65 + and intricate molecular systems.
66 +
67 + Water has proven to be a challenging substance to depict in
68 + simulations, and a variety of models have been developed to describe
69 + its behavior under varying simulation
70 + conditions.\cite{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
71 + These models have been used to investigate important physical
72 + phenomena like phase transitions and the hydrophobic
73 + effect.\cite{Yamada02} With the choice of models available, it
74 + is only natural to compare the models under interesting thermodynamic
75 + conditions in an attempt to clarify the limitations of each of the
76 + models.\cite{modelProps} Two important property to quantify are the
77 + Gibbs and Helmholtz free energies, particularly for the solid forms of
78 + water.  Difficulty in these types of studies typically arises from the
79 + assortment of possible crystalline polymorphs that water adopts over a
80 + wide range of pressures and temperatures. There are currently 13
81 + recognized forms of ice, and it is a challenging task to investigate
82 + the entire free energy landscape.\cite{Sanz04} Ideally, research is
83 + focused on the phases having the lowest free energy at a given state
84 + point, because these phases will dictate the true transition
85 + temperatures and pressures for their respective model.
86 +
87 + In this paper, standard reference state methods were applied to the
88 + study of crystalline water polymorphs in the low pressure regime. This
89 + work is unique in the fact that one of the crystal lattices was
90 + arrived at through crystallization of a computationally efficient
91 + water model under constant pressure and temperature
92 + conditions. Crystallization events are interesting in and of
93 + themselves\cite{Matsumoto02,Yamada02}; however, the crystal structure
94 + obtained in this case was different from any previously observed ice
95 + polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
96 + was termed Ice-{\it i} in homage to its origin in computational
97 + simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
98 + water molecules that stack in rows of interlocking water
99 + tetramers. Proton ordering can be accomplished by orienting two of the
100 + waters so that both of their donated hydrogen bonds are internal to
101 + their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
102 + constructed of water tetramers, the hydrogen bonds are not as linear
103 + as those observed in ice $I_h$, however the interlocking of these
104 + subunits appears to provide significant stabilization to the overall
105 + crystal. The arrangement of these tetramers results in surrounding
106 + open octagonal cavities that are typically greater than 6.3 \AA\ in
107 + diameter. This relatively open overall structure leads to crystals
108 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
109 +
110 + \begin{figure}
111 + \includegraphics[width=\linewidth]{unitCell.eps}
112 + \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-2{\it i}, the elongated variant of Ice-{\it i}.  For Ice-{\it i}, the $a$ to $c$ relation is given by $a = 2.1214c$, while for Ice-2{\it i}, $a = 1.7850c$.}
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 in the previous study indicated that Ice-{\it i} is the
126 + minimum energy crystal structure for the single point water models
127 + being studied (for discussions on these single point dipole models,
128 + see the previous work and related
129 + articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
130 + consider energetic stabilization and neglect entropic contributions to
131 + the overall free energy. To address this issue, the absolute free
132 + energy of this crystal was calculated using thermodynamic integration
133 + and compared to the free energies of cubic and hexagonal ice $I$ (the
134 + experimental low density ice polymorphs) and ice B (a higher density,
135 + but very stable crystal structure observed by B\`{a}ez and Clancy in
136 + free energy studies of SPC/E).\cite{Baez95b} This work includes
137 + results for the water model from which Ice-{\it i} was crystallized
138 + (soft sticky dipole extended, SSD/E) in addition to several common
139 + water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field
140 + parametrized single point dipole water model (soft sticky dipole
141 + reaction field, SSD/RF). In should be noted that a second version of
142 + Ice-{\it i} (Ice-2{\it i}) was used in calculations involving SPC/E,
143 + TIP4P, and TIP5P. The unit cell of this crystal (Fig. \ref{iceiCell}B)
144 + is similar to the Ice-{\it i} unit it is extended in the direction of
145 + the (001) face and compressed 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 propogated using the symplectic DLM
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}
155  
# Line 66 | Line 162 | For the thermodynamic integration of molecular crystal
162   resulting in a pressure of approximately 1 atm at their respective
163   temperatures.
164  
165 + A single thermodynamic integration involves a sequence of simulations
166 + over which the system of interest is converted into a reference system
167 + for which the free energy is known. This transformation path is then
168 + integrated in order to determine the free energy difference between
169 + the two states:
170 + \begin{equation}
171 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
172 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
173 + \end{equation}
174 + where $V$ is the interaction potential and $\lambda$ is the
175 + transformation parameter that scales the overall
176 + potential. Simulations are distributed unevenly along this path in
177 + order to sufficiently sample the regions of greatest change in the
178 + potential. Typical integrations in this study consisted of $\sim$25
179 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
180 + (near the reference state) in length.
181 +
182   For the thermodynamic integration of molecular crystals, the Einstein
183   Crystal is chosen as the reference state that the system is converted
184   to over the course of the simulation. In an Einstein Crystal, the
# Line 92 | Line 205 | state.
205   minimum potential energy of the ideal crystal. In the case of
206   molecular liquids, the ideal vapor is chosen as the target reference
207   state.
208 +
209   \begin{figure}
210 < \includegraphics[scale=1.0]{rotSpring.eps}
210 > \includegraphics[width=\linewidth]{rotSpring.eps}
211   \caption{Possible orientational motions for a restrained molecule.
212   $\theta$ angles correspond to displacement from the body-frame {\it
213   z}-axis, while $\omega$ angles correspond to rotation about the
# Line 104 | Line 218 | cubic switching between 100\% and 85\% of the cutoff v
218   \end{figure}
219  
220   Charge, dipole, and Lennard-Jones interactions were modified by a
221 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
222 < applying this function, these interactions are smoothly truncated,
223 < thereby avoiding poor energy conserving dynamics resulting from
224 < harsher truncation schemes. The effect of a long-range correction was
225 < also investigated on select model systems in a variety of manners. For
226 < the SSD/RF model, a reaction field with a fixed dielectric constant of
227 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
228 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
229 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
230 < compare with the 9 \AA\ cutoff results. Finally, results from the use
231 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
221 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
222 > ). By applying this function, these interactions are smoothly
223 > truncated, thereby avoiding poor energy conserving dynamics resulting
224 > from harsher truncation schemes. The effect of a long-range correction
225 > was also investigated on select model systems in a variety of
226 > manners. For the SSD/RF model, a reaction field with a fixed
227 > dielectric constant of 80 was applied in all
228 > simulations.\cite{Onsager36} For a series of the least computationally
229 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
230 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
231 > \AA\ cutoff results. Finally, results from the use of an Ewald
232 > summation were estimated for TIP3P and SPC/E by performing
233   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
234 < mechanics software package. TINKER was chosen because it can also
235 < propogate the motion of rigid-bodies, and provides the most direct
236 < comparison to the results from OOPSE. The calculated energy difference
237 < in the presence and absence of PME was applied to the previous results
238 < in order to predict changes in the free energy landscape.
234 > mechanics software package.\cite{Tinker} TINKER was chosen because it
235 > can also propagate the motion of rigid-bodies, and provides the most
236 > direct comparison to the results from OOPSE. The calculated energy
237 > difference in the presence and absence of PME was applied to the
238 > previous results in order to predict changes in the free energy
239 > landscape.
240  
241   \section{Results and discussion}
242  
# Line 149 | Line 265 | kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF
265   of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are
266   kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.}
267   \begin{tabular}{ l  c  c  c  c }
268 < \hline \\[-7mm]
268 > \hline
269   \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \  & \ \quad \ \ \ \ B \ \  & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\
270 < \hline \\[-3mm]
270 > \hline
271   \ \quad \ TIP3P  & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\
272   \ \quad \ TIP4P  & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\
273   \ \quad \ TIP5P  & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\
# Line 178 | Line 294 | TIP4P in the high pressure regime.\cite{Sanz04}
294   representative of the dense ice polymorphs. A recent study by Sanz
295   {\it et al.} goes into detail on the phase diagrams for SPC/E and
296   TIP4P in the high pressure regime.\cite{Sanz04}
297 +
298   \begin{figure}
299   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
300   \caption{Phase diagram for the TIP3P water model in the low pressure
# Line 187 | Line 304 | higher in energy and don't appear in the phase diagram
304   higher in energy and don't appear in the phase diagram.}
305   \label{tp3phasedia}
306   \end{figure}
307 +
308   \begin{figure}
309   \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
310   \caption{Phase diagram for the SSD/RF water model in the low pressure
# Line 205 | Line 323 | temperatures of several common water models compared w
323   \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
324   temperatures of several common water models compared with experiment.}
325   \begin{tabular}{ l  c  c  c  c  c  c  c }
326 < \hline \\[-7mm]
326 > \hline
327   \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\
328 < \hline \\[-3mm]
328 > \hline
329   \ \ $T_m$ (K)  & \ \ 269 & \ \ 265 & \ \ 271 &  297 & \ \ - & \ \ 278 & \ \ 273\\
330   \ \ $T_b$ (K)  & \ \ 357 & \ \ 354 & \ \ 337 &  396 & \ \ - & \ \ 349 & \ \ 373\\
331   \ \ $T_s$ (K)  & \ \ - & \ \ - & \ \ - &  - & \ \ 355 & \ \ - & \ \ -\\
# Line 235 | Line 353 | advantagious in that it facilitated the spontaneous cr
353   at 355 K. This is due to the significant stability of Ice-{\it i} over
354   all other polymorphs for this particular model under these
355   conditions. While troubling, this behavior turned out to be
356 < advantagious in that it facilitated the spontaneous crystallization of
356 > advantageous in that it facilitated the spontaneous crystallization of
357   Ice-{\it i}. These observations provide a warning that simulations of
358   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
359   risk of spontaneous crystallization. However, this risk changes when
360   applying a longer cutoff.
361  
362 + \begin{figure}
363 + \includegraphics[width=\linewidth]{cutoffChange.eps}
364 + \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
365 + TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12
366 + \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9
367 + \AA\. These crystals are unstable at 200 K and rapidly convert into a
368 + liquid. The connecting lines are qualitative visual aid.}
369 + \label{incCutoff}
370 + \end{figure}
371  
372 + Increasing the cutoff radius in simulations of the more
373 + computationally efficient water models was done in order to evaluate
374 + the trend in free energy values when moving to systems that do not
375 + involve potential truncation. As seen in Fig. \ref{incCutoff}, the
376 + free energy of all the ice polymorphs show a substantial dependence on
377 + cutoff radius. In general, there is a narrowing of the free energy
378 + differences while moving to greater cutoff radius. Interestingly, by
379 + increasing the cutoff radius, the free energy gap was narrowed enough
380 + in the SSD/E model that the liquid state is preferred under standard
381 + simulation conditions (298 K and 1 atm). Thus, it is recommended that
382 + simulations using this model choose interaction truncation radii
383 + greater than 9 \AA\. This narrowing trend is much more subtle in the
384 + case of SSD/RF, indicating that the free energies calculated with a
385 + reaction field present provide a more accurate picture of the free
386 + energy landscape in the absence of potential truncation.
387  
388 + To further study the changes resulting to the inclusion of a
389 + long-range interaction correction, the effect of an Ewald summation
390 + was estimated by applying the potential energy difference do to its
391 + inclusion in systems in the presence and absence of the
392 + correction. This was accomplished by calculation of the potential
393 + energy of identical crystals with and without PME using TINKER. The
394 + free energies for the investigated polymorphs using the TIP3P and
395 + SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
396 + are not fully supported in TINKER, so the results for these models
397 + could not be estimated. The same trend pointed out through increase of
398 + cutoff radius is observed in these PME results. Ice-{\it i} is the
399 + preferred polymorph at ambient conditions for both the TIP3P and SPC/E
400 + water models; however, there is a narrowing of the free energy
401 + differences between the various solid forms. In the case of SPC/E this
402 + narrowing is significant enough that it becomes less clear cut that
403 + Ice-{\it i} is the most stable polymorph, and is possibly metastable
404 + with respect to ice B and possibly ice $I_c$. However, these results
405 + do not significantly alter the finding that the Ice-{\it i} polymorph
406 + is a stable crystal structure that should be considered when studying
407 + the phase behavior of water models.
408 +
409 + \begin{table*}
410 + \begin{minipage}{\linewidth}
411 + \renewcommand{\thefootnote}{\thempfootnote}
412 + \begin{center}
413 + \caption{The free energy of the studied ice polymorphs after applying
414 + the energy difference attributed to the inclusion of the PME
415 + long-range interaction correction. Units are kcal/mol.}
416 + \begin{tabular}{ l  c  c  c  c }
417 + \hline
418 + \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
419 + \hline
420 + \ \ TIP3P  & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\
421 + \ \ SPC/E  & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\
422 + \end{tabular}
423 + \label{pmeShift}
424 + \end{center}
425 + \end{minipage}
426 + \end{table*}
427 +
428   \section{Conclusions}
429  
430 + The free energy for proton ordered variants of hexagonal and cubic ice
431 + $I$, ice B, and recently discovered Ice-{\it i} where calculated under
432 + standard conditions for several common water models via thermodynamic
433 + integration. All the water models studied show Ice-{\it i} to be the
434 + minimum free energy crystal structure in the with a 9 \AA\ switching
435 + function cutoff. Calculated melting and boiling points show
436 + surprisingly good agreement with the experimental values; however, the
437 + solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
438 + interaction truncation was investigated through variation of the
439 + cutoff radius, use of a reaction field parameterized model, and
440 + estimation of the results in the presence of the Ewald summation
441 + correction. Interaction truncation has a significant effect on the
442 + computed free energy values, and may significantly alter the free
443 + energy landscape for the more complex multipoint water models. Despite
444 + these effects, these results show Ice-{\it i} to be an important ice
445 + polymorph that should be considered in simulation studies.
446 +
447 + Due to this relative stability of Ice-{\it i} in all manner of
448 + investigated simulation examples, the question arises as to possible
449 + experimental observation of this polymorph. The rather extensive past
450 + and current experimental investigation of water in the low pressure
451 + regime leads the authors to be hesitant in ascribing relevance outside
452 + of computational models, hence the descriptive name presented. That
453 + being said, there are certain experimental conditions that would
454 + provide the most ideal situation for possible observation. These
455 + include the negative pressure or stretched solid regime, small
456 + clusters in vacuum deposition environments, and in clathrate
457 + structures involving small non-polar molecules.
458 +
459   \section{Acknowledgments}
460   Support for this project was provided by the National Science
461   Foundation under grant CHE-0134881. Computation time was provided by
462 < the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
463 < DMR-0079647.
462 > the Notre Dame High Performance Computing Cluster and the Notre Dame
463 > Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
464  
465   \newpage
466  

Diff Legend

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