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 1458 by chrisfen, Wed Sep 15 06:34:49 2004 UTC vs.
Revision 1462 by chrisfen, Wed Sep 15 20:47:13 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 47 | Line 59 | Notre Dame, Indiana 46556}
59   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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{Matsumoto02andOthers} 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{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04}
78 + Many of these models have been used to investigate important physical
79 + phenomena like phase transitions and the hydrophobic
80 + effect.\cite{evenMorePapers} With the advent of numerous differing
81 + models, it is only natural that attention is placed on the properties
82 + of the models themselves in an attempt to clarify their benefits and
83 + limitations when applied to a system of interest.\cite{modelProps} One
84 + important but challenging property to quantify is the free energy,
85 + particularly of the solid forms of water. Difficulty in these types of
86 + studies typically arises from the assortment of possible crystalline
87 + polymorphs that water that water adopts over a wide range of pressures
88 + and temperatures. There are currently 13 recognized forms of ice, and
89 + it is a challenging task to investigate the entire free energy
90 + landscape.\cite{Sanz04} Ideally, research is focused on the phases
91 + having the lowest free energy, because these phases will dictate the
92 + true transition temperatures and 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{Matsumoto02,Yamada02}; 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 72 | Line 173 | the two states:
173   integrated in order to determine the free energy difference between
174   the two states:
175   \begin{equation}
75 – \begin{center}
176   \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
177   )}{\partial\lambda}\right\rangle_\lambda d\lambda,
78 – \end{center}
178   \end{equation}
179   where $V$ is the interaction potential and $\lambda$ is the
180 < transformation parameter. Simulations are distributed unevenly along
181 < this path in order to sufficiently sample the regions of greatest
182 < change in the potential. Typical integrations in this study consisted
183 < of $\sim$25 simulations ranging from 300 ps (for the unaltered system)
184 < to 75 ps (near the reference state) in length.
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
# Line 122 | Line 222 | cubic switching between 100\% and 85\% of the cutoff v
222   \end{figure}
223  
224   Charge, dipole, and Lennard-Jones interactions were modified by a
225 < cubic switching between 100\% and 85\% of the cutoff value (9 \AA ). By
226 < applying this function, these interactions are smoothly truncated,
227 < thereby avoiding poor energy conserving dynamics resulting from
228 < harsher truncation schemes. The effect of a long-range correction was
229 < also investigated on select model systems in a variety of manners. For
230 < the SSD/RF model, a reaction field with a fixed dielectric constant of
231 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
232 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
233 < simulations were performed with longer cutoffs of 12 and 15 \AA\ to
234 < compare with the 9 \AA\ cutoff results. Finally, results from the use
235 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
225 > cubic switching between 100\% and 85\% of the cutoff value (9 \AA
226 > ). By applying this function, these interactions are smoothly
227 > truncated, thereby avoiding poor energy conserving dynamics resulting
228 > from harsher truncation schemes. The effect of a long-range correction
229 > was also investigated on select model systems in a variety of
230 > manners. For the SSD/RF model, a reaction field with a fixed
231 > dielectric constant of 80 was applied in all
232 > simulations.\cite{Onsager36} For a series of the least computationally
233 > expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
234 > performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
235 > \AA\ cutoff results. Finally, results from the use of an Ewald
236 > summation were estimated for TIP3P and SPC/E by performing
237   calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
238 < mechanics software package. TINKER was chosen because it can also
239 < propogate the motion of rigid-bodies, and provides the most direct
240 < comparison to the results from OOPSE. The calculated energy difference
241 < in the presence and absence of PME was applied to the previous results
242 < in order to predict changes in the free energy landscape.
238 > mechanics software package.\cite{Tinker} TINKER was chosen because it
239 > can also propagate the motion of rigid-bodies, and provides the most
240 > direct comparison to the results from OOPSE. The calculated energy
241 > difference in the presence and absence of PME was applied to the
242 > previous results in order to predict changes in the free energy
243 > landscape.
244  
245   \section{Results and discussion}
246  
# Line 253 | Line 355 | advantagious in that it facilitated the spontaneous cr
355   at 355 K. This is due to the significant stability of Ice-{\it i} over
356   all other polymorphs for this particular model under these
357   conditions. While troubling, this behavior turned out to be
358 < advantagious in that it facilitated the spontaneous crystallization of
358 > advantageous in that it facilitated the spontaneous crystallization of
359   Ice-{\it i}. These observations provide a warning that simulations of
360   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
361   risk of spontaneous crystallization. However, this risk changes when
# Line 275 | Line 377 | differences while moving to greater cutoff radius. Thi
377   involve potential truncation. As seen in Fig. \ref{incCutoff}, the
378   free energy of all the ice polymorphs show a substantial dependence on
379   cutoff radius. In general, there is a narrowing of the free energy
380 < differences while moving to greater cutoff radius. This trend is much
381 < more subtle in the case of SSD/RF, indicating that the free energies
382 < calculated with a reaction field present provide a more accurate
383 < picture of the free energy landscape in the absence of potential
384 < truncation.
380 > differences while moving to greater cutoff radius. Interestingly, by
381 > increasing the cutoff radius, the free energy gap was narrowed enough
382 > in the SSD/E model that the liquid state is preferred under standard
383 > simulation conditions (298 K and 1 atm). Thus, it is recommended that
384 > simulations using this model choose interaction truncation radii
385 > greater than 9 \AA\. This narrowing trend is much more subtle in the
386 > case of SSD/RF, indicating that the free energies calculated with a
387 > reaction field present provide a more accurate picture of the free
388 > energy landscape in the absence of potential truncation.
389  
390   To further study the changes resulting to the inclusion of a
391   long-range interaction correction, the effect of an Ewald summation
# Line 291 | Line 397 | cutoff radius is observed in these results. Ice-{\it i
397   SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P
398   are not fully supported in TINKER, so the results for these models
399   could not be estimated. The same trend pointed out through increase of
400 < cutoff radius is observed in these results. Ice-{\it i} is the
400 > cutoff radius is observed in these PME results. Ice-{\it i} is the
401   preferred polymorph at ambient conditions for both the TIP3P and SPC/E
402   water models; however, there is a narrowing of the free energy
403   differences between the various solid forms. In the case of SPC/E this
# Line 335 | Line 441 | computed free energy values, but Ice-{\it i} is still
441   cutoff radius, use of a reaction field parameterized model, and
442   estimation of the results in the presence of the Ewald summation
443   correction. Interaction truncation has a significant effect on the
444 < computed free energy values, but Ice-{\it i} is still observed to be a
445 < relavent ice polymorph in simulation studies.
444 > computed free energy values, and may significantly alter the free
445 > energy landscape for the more complex multipoint water models. Despite
446 > these effects, these results show Ice-{\it i} to be an important ice
447 > polymorph that should be considered in simulation studies.
448  
449 + Due to this relative stability of Ice-{\it i} in all manner of
450 + investigated simulation examples, the question arises as to possible
451 + experimental observation of this polymorph. The rather extensive past
452 + and current experimental investigation of water in the low pressure
453 + regime leads the authors to be hesitant in ascribing relevance outside
454 + of computational models, hence the descriptive name presented. That
455 + being said, there are certain experimental conditions that would
456 + provide the most ideal situation for possible observation. These
457 + include the negative pressure or stretched solid regime, small
458 + clusters in vacuum deposition environments, and in clathrate
459 + structures involving small non-polar molecules.
460 +
461   \section{Acknowledgments}
462   Support for this project was provided by the National Science
463   Foundation under grant CHE-0134881. Computation time was provided by

Diff Legend

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