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

Comparing trunk/iceiPaper/iceiPaper.tex (file contents):
Revision 1457 by chrisfen, Tue Sep 14 23:03:53 2004 UTC vs.
Revision 1461 by chrisfen, Wed Sep 15 20:22:21 2004 UTC

# Line 1 | Line 1
1 +
2   %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
3   \documentclass[preprint,aps,endfloats]{revtex4}
4   %\documentclass[11pt]{article}
# Line 20 | Line 21
21  
22   \begin{document}
23  
24 < \title{A Free Energy Study of Low Temperature and Anomolous Ice}
24 > \title{A Free Energy Study of Low Temperature and Anomalous Ice}
25  
26   \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
27   \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
# Line 34 | Line 35 | Notre Dame, Indiana 46556}
35   %\doublespacing
36  
37   \begin{abstract}
38 + The free energies of several ice polymorphs in the low pressure regime
39 + were calculated using thermodynamic integration of systems consisting
40 + of a variety of common water models. Ice-{\it i}, a recent
41 + computationally observed solid structure, was determined to be the
42 + stable state with the lowest free energy for all the water models
43 + investigated. Phase diagrams were generated, and melting and boiling
44 + points for all the models were determined and show relatively good
45 + agreement with experiment, although the solid phase is different
46 + between simulation and experiment. In addition, potential truncation
47 + was shown to have an effect on the calculated free energies, and may
48 + result in altered free energy landscapes.
49   \end{abstract}
50  
51   \maketitle
# Line 48 | Line 60 | Notre Dame, Indiana 46556}
60  
61   \section{Introduction}
62  
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 +
74 + Water has proven to be a challenging substance to depict in
75 + simulations, and has resulted in a variety of models that attempt to
76 + describe its behavior under a varying simulation
77 + conditions.\cite{lotsOfWaterPapers} Many of these models have been
78 + used to investigate important physical phenomena like phase
79 + transitions and the hydrophobic effect.\cite{evenMorePapers} With the
80 + advent of numerous differing models, it is only natural that attention
81 + is placed on the properties of the models themselves in an attempt to
82 + clarify their benefits and limitations when applied to a system of
83 + interest.\cite{modelProps} One important but challenging property to
84 + quantify is the free energy, particularly of the solid forms of
85 + water. Difficulty in these types of studies typically arises from the
86 + assortment of possible crystalline polymorphs that water that water
87 + adopts over a wide range of pressures and temperatures. There are
88 + currently 13 recognized forms of ice, and it is a challenging task to
89 + investigate the entire free energy landscape.\cite{Sanz04} Ideally,
90 + research is focused on the phases having the lowest free energy,
91 + because these phases will dictate the true transition temperatures and
92 + pressures for their respective model.
93 +
94 + In this paper, standard reference state methods were applied to the
95 + study of crystalline water polymorphs in the low pressure regime. This
96 + work is unique in the fact that one of the crystal lattices was
97 + arrived at through crystallization of a computationally efficient
98 + water model under constant pressure and temperature
99 + conditions. Crystallization events are interesting in and of
100 + themselves\cite{nucleationStudies}; however, the crystal structure
101 + obtained in this case was different from any previously observed ice
102 + polymorphs, in experiment or simulation.\cite{Fennell04} This crystal
103 + was termed Ice-{\it i} in homage to its origin in computational
104 + simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
105 + water molecules that stack in rows of interlocking water
106 + tetramers. Proton ordering can be accomplished by orienting two of the
107 + waters so that both of their donated hydrogen bonds are internal to
108 + their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
109 + constructed of water tetramers, the hydrogen bonds are not as linear
110 + as those observed in ice $I_h$, however the interlocking of these
111 + subunits appears to provide significant stabilization to the overall
112 + crystal. The arrangement of these tetramers results in surrounding
113 + open octagonal cavities that are typically greater than 6.3 \AA\ in
114 + diameter. This relatively open overall structure leads to crystals
115 + that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
116 + \begin{figure}
117 + \includegraphics[scale=1.0]{unitCell.eps}
118 + \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$.}
119 + \label{iceiCell}
120 + \end{figure}
121 + \begin{figure}
122 + \includegraphics[scale=1.0]{orderedIcei.eps}
123 + \caption{Image of a proton ordered crystal of Ice-{\it i} looking
124 + down the (001) crystal face. The rows of water tetramers surrounded by
125 + octagonal pores leads to a crystal structure that is significantly
126 + less dense than ice $I_h$.}
127 + \label{protOrder}
128 + \end{figure}
129 +
130 + Results in the previous study indicated that Ice-{\it i} is the
131 + minimum energy crystal structure for the single point water models
132 + being studied (for discussions on these single point dipole models,
133 + see the previous work and related
134 + articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only
135 + consider energetic stabilization and neglect entropic contributions to
136 + the overall free energy. To address this issue, the absolute free
137 + energy of this crystal was calculated using thermodynamic integration
138 + and compared to the free energies of cubic and hexagonal ice $I$ (the
139 + experimental low density ice polymorphs) and ice B (a higher density,
140 + but very stable crystal structure observed by B\`{a}ez and Clancy in
141 + free energy studies of SPC/E).\cite{Baez95b} This work includes
142 + results for the water model from which Ice-{\it i} was crystallized
143 + (soft sticky dipole extended, SSD/E) in addition to several common
144 + water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field
145 + parametrized single point dipole water model (soft sticky dipole
146 + reaction field, SSD/RF). In should be noted that a second version of
147 + Ice-{\it i} (Ice-2{\it i}) was used in calculations involving SPC/E,
148 + TIP4P, and TIP5P. The unit cell of this crystal (Fig. \ref{iceiCell}B)
149 + is similar to the Ice-{\it i} unit it is extended in the direction of
150 + the (001) face and compressed along the other two faces.
151 +
152   \section{Methods}
153  
154   Canonical ensemble (NVT) molecular dynamics calculations were
155   performed using the OOPSE (Object-Oriented Parallel Simulation Engine)
156   molecular mechanics package. All molecules were treated as rigid
157 < bodies, with orientational motion propogated using the symplectic DLM
157 > bodies, with orientational motion propagated using the symplectic DLM
158   integration method. Details about the implementation of these
159   techniques can be found in a recent publication.\cite{Meineke05}
160  
# Line 66 | Line 167 | For the thermodynamic integration of molecular crystal
167   resulting in a pressure of approximately 1 atm at their respective
168   temperatures.
169  
170 + A single thermodynamic integration involves a sequence of simulations
171 + over which the system of interest is converted into a reference system
172 + for which the free energy is known. This transformation path is then
173 + integrated in order to determine the free energy difference between
174 + the two states:
175 + \begin{equation}
176 + \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
177 + )}{\partial\lambda}\right\rangle_\lambda d\lambda,
178 + \end{equation}
179 + where $V$ is the interaction potential and $\lambda$ is the
180 + transformation parameter that scales the overall
181 + potential. Simulations are distributed unevenly along this path in
182 + order to sufficiently sample the regions of greatest change in the
183 + potential. Typical integrations in this study consisted of $\sim$25
184 + simulations ranging from 300 ps (for the unaltered system) to 75 ps
185 + (near the reference state) in length.
186 +
187   For the thermodynamic integration of molecular crystals, the Einstein
188   Crystal is chosen as the reference state that the system is converted
189   to over the course of the simulation. In an Einstein Crystal, the
# Line 117 | Line 235 | propogate the motion of rigid-bodies, and provides the
235   of an Ewald summation were estimated for TIP3P and SPC/E by performing
236   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
237   mechanics software package. TINKER was chosen because it can also
238 < propogate the motion of rigid-bodies, and provides the most direct
238 > propagate the motion of rigid-bodies, and provides the most direct
239   comparison to the results from OOPSE. The calculated energy difference
240   in the presence and absence of PME was applied to the previous results
241   in order to predict changes in the free energy landscape.
# 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. This trend is much
379 < more subtle in the case of SSD/RF, indicating that the free energies
380 < calculated with a reaction field present provide a more accurate
381 < picture of the free energy landscape in the absence of potential
382 < truncation.
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
# Line 263 | Line 395 | cutoff radius is observed in these results. Ice-{\it i
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 results. Ice-{\it i} is the
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
# Line 278 | Line 410 | the phase behavior of water models.
410   \begin{minipage}{\linewidth}
411   \renewcommand{\thefootnote}{\thempfootnote}
412   \begin{center}
413 < \caption{The free energy of the studied ice polymorphs after applying the energy difference attributed to the inclusion of the PME long-range interaction correction. Units are kcal/mol.}
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 \\[-7mm]
418   \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\
# Line 293 | Line 427 | the phase behavior of water models.
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