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 1473 by chrisfen, Fri Sep 17 14:54:49 2004 UTC vs.
Revision 1488 by chrisfen, Tue Sep 21 20:17:12 2004 UTC

# Line 33 | Line 33 | The free energies of several ice polymorphs in the low
33   %\doublespacing
34  
35   \begin{abstract}
36 < The free energies of several ice polymorphs in the low pressure regime
37 < were calculated using thermodynamic integration.  These integrations
38 < were done for most of the common water models. Ice-{\it i}, a
39 < structure we recently observed to be stable in one of the single-point
40 < water models, was determined to be the stable crystalline state (at 1
41 < atm) for {\it all} the water models investigated.  Phase diagrams were
36 > The absolute free energies of several ice polymorphs which are stable
37 > at low pressures were calculated using thermodynamic integration to a
38 > reference system (the Einstein crystal).  These integrations were
39 > performed for most of the common water models (SPC/E, TIP3P, TIP4P,
40 > TIP5P, SSD/E, SSD/RF). Ice-{\it i}, a structure we recently observed
41 > crystallizing at room temperature for one of the single-point water
42 > models, was determined to be the stable crystalline state (at 1 atm)
43 > for {\it all} the water models investigated.  Phase diagrams were
44   generated, and phase coexistence lines were determined for all of the
45 < known low-pressure ice structures under all of the common water
46 < models.  Additionally, potential truncation was shown to have an
47 < effect on the calculated free energies, and can result in altered free
48 < energy landscapes.
45 > known low-pressure ice structures under all of these water models.
46 > Additionally, potential truncation was shown to have an effect on the
47 > calculated free energies, and can result in altered free energy
48 > landscapes.  Structure factor predictions for the new crystal were
49 > generated and we await experimental confirmation of the existence of
50 > this new polymorph.
51   \end{abstract}
52  
53   %\narrowtext
# Line 54 | Line 58 | Computer simulations are a valuable tool for studying
58  
59   \section{Introduction}
60  
57 Computer simulations are a valuable tool for studying the phase
58 behavior of systems ranging from small or simple molecules to complex
59 biological species.\cite{Matsumoto02,Li96,Marrink01} Useful techniques
60 have been developed to investigate the thermodynamic properites of
61 model substances, providing both qualitative and quantitative
62 comparisons between simulations and
63 experiment.\cite{Widom63,Frenkel84} Investigation of these properties
64 leads to the development of new and more accurate models, leading to
65 better understanding and depiction of physical processes and intricate
66 molecular systems.
67
61   Water has proven to be a challenging substance to depict in
62   simulations, and a variety of models have been developed to describe
63   its behavior under varying simulation
64 < conditions.\cite{Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Berendsen98,Mahoney00,Fennell04}
64 > conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
65   These models have been used to investigate important physical
66 < phenomena like phase transitions, molecule transport, and the
66 > phenomena like phase transitions, transport properties, and the
67   hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
68   choice of models available, it is only natural to compare the models
69   under interesting thermodynamic conditions in an attempt to clarify
70   the limitations of each of the
71   models.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two
72 < important property to quantify are the Gibbs and Helmholtz free
72 > important properties to quantify are the Gibbs and Helmholtz free
73   energies, particularly for the solid forms of water.  Difficulty in
74   these types of studies typically arises from the assortment of
75   possible crystalline polymorphs that water adopts over a wide range of
# Line 84 | Line 77 | these phases will dictate the true transition temperat
77   of ice, and it is a challenging task to investigate the entire free
78   energy landscape.\cite{Sanz04} Ideally, research is focused on the
79   phases having the lowest free energy at a given state point, because
80 < these phases will dictate the true transition temperatures and
80 > these phases will dictate the relevant transition temperatures and
81   pressures for the model.
82  
83   In this paper, standard reference state methods were applied to known
84   crystalline water polymorphs in the low pressure regime.  This work is
85 < unique in the fact that one of the crystal lattices was arrived at
86 < through crystallization of a computationally efficient water model
87 < under constant pressure and temperature conditions. Crystallization
88 < events are interesting in and of
89 < themselves;\cite{Matsumoto02,Yamada02} however, the crystal structure
90 < obtained in this case is different from any previously observed ice
91 < polymorphs in experiment or simulation.\cite{Fennell04} We have named
92 < this structure Ice-{\it i} to indicate its origin in computational
93 < simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight
94 < water molecules that stack in rows of interlocking water
95 < tetramers. Proton ordering can be accomplished by orienting two of the
96 < molecules so that both of their donated hydrogen bonds are internal to
97 < their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal
98 < constructed of water tetramers, the hydrogen bonds are not as linear
99 < as those observed in ice $I_h$, however the interlocking of these
100 < subunits appears to provide significant stabilization to the overall
85 > unique in that one of the crystal lattices was arrived at through
86 > crystallization of a computationally efficient water model under
87 > constant pressure and temperature conditions. Crystallization events
88 > are interesting in and of themselves;\cite{Matsumoto02,Yamada02}
89 > however, the crystal structure obtained in this case is different from
90 > any previously observed ice polymorphs in experiment or
91 > simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
92 > to indicate its origin in computational simulation. The unit cell
93 > (Fig. \ref{iceiCell}A) consists of eight water molecules that stack in
94 > rows of interlocking water tetramers. Proton ordering can be
95 > accomplished by orienting two of the molecules so that both of their
96 > donated hydrogen bonds are internal to their tetramer
97 > (Fig. \ref{protOrder}). As expected in an ice crystal constructed of
98 > water tetramers, the hydrogen bonds are not as linear as those
99 > observed in ice $I_h$, however the interlocking of these subunits
100 > appears to provide significant stabilization to the overall
101   crystal. The arrangement of these tetramers results in surrounding
102   open octagonal cavities that are typically greater than 6.3 \AA\ in
103   diameter. This relatively open overall structure leads to crystals
# Line 112 | Line 105 | that are 0.07 g/cm$^3$ less dense on average than ice
105  
106   \begin{figure}
107   \includegraphics[width=\linewidth]{unitCell.eps}
108 < \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the
109 < elongated variant of Ice-{\it i}.  The spheres represent the
108 > \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it i}$^\prime$,
109 > the elongated variant of Ice-{\it i}.  The spheres represent the
110   center-of-mass locations of the water molecules.  The $a$ to $c$
111   ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
112   $a:2.1214c$ and $a:1.7850c$ respectively.}
# Line 131 | Line 124 | investigated (for discussions on these single point di
124  
125   Results from our previous study indicated that Ice-{\it i} is the
126   minimum energy crystal structure for the single point water models we
127 < investigated (for discussions on these single point dipole models, see
128 < our previous work and related
127 > had investigated (for discussions on these single point dipole models,
128 > see our previous work and related
129   articles).\cite{Fennell04,Liu96,Bratko85} Those results only
130   considered energetic stabilization and neglected entropic
131 < contributions to the overall free energy. To address this issue, the
132 < absolute free energy of this crystal was calculated using
131 > contributions to the overall free energy. To address this issue, we
132 > have calculated the absolute free energy of this crystal using
133   thermodynamic integration and compared to the free energies of cubic
134   and hexagonal ice $I$ (the experimental low density ice polymorphs)
135   and ice B (a higher density, but very stable crystal structure
# Line 145 | Line 138 | be noted that a second version of Ice-{\it i} (Ice-$i^
138   from which Ice-{\it i} was crystallized (SSD/E) in addition to several
139   common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
140   field parametrized single point dipole water model (SSD/RF). It should
141 < be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used
142 < in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of
143 < this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit
144 < it is extended in the direction of the (001) face and compressed along
145 < the other two faces.
141 > be noted that a second version of Ice-{\it i} (Ice-{\it i}$^\prime$)
142 > was used in calculations involving SPC/E, TIP4P, and TIP5P. The unit
143 > cell of this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it
144 > i} unit it is extended in the direction of the (001) face and
145 > compressed along the other two faces.  There is typically a small
146 > distortion of proton ordered Ice-{\it i}$^\prime$ that converts the
147 > normally square tetramer into a rhombus with alternating approximately
148 > 85 and 95 degree angles.  The degree of this distortion is model
149 > dependent and significant enough to split the tetramer diagonal
150 > location peak in the radial distribution function.
151  
152   \section{Methods}
153  
# Line 168 | Line 166 | this same technique in order to determine melting poin
166   crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
167   SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
168   and 400 K for all of these water models were also determined using
169 < this same technique in order to determine melting points and generate
170 < phase diagrams. All simulations were carried out at densities
171 < resulting in a pressure of approximately 1 atm at their respective
172 < temperatures.
169 > this same technique in order to determine melting points and to
170 > generate phase diagrams. All simulations were carried out at densities
171 > which correspond to a pressure of approximately 1 atm at their
172 > respective temperatures.
173  
174 < A single thermodynamic integration involves a sequence of simulations
175 < over which the system of interest is converted into a reference system
176 < for which the free energy is known analytically. This transformation
177 < path is then integrated in order to determine the free energy
178 < difference between the two states:
174 > Thermodynamic integration involves a sequence of simulations during
175 > which the system of interest is converted into a reference system for
176 > which the free energy is known analytically. This transformation path
177 > is then integrated in order to determine the free energy difference
178 > between the two states:
179   \begin{equation}
180   \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
181   )}{\partial\lambda}\right\rangle_\lambda d\lambda,
# Line 239 | Line 237 | molecules.  In this study, we apply of one of the most
237   literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
238   typically differ in regard to the path taken for switching off the
239   interaction potential to convert the system to an ideal gas of water
240 < molecules.  In this study, we apply of one of the most convenient
241 < methods and integrate over the $\lambda^4$ path, where all interaction
242 < parameters are scaled equally by this transformation parameter.  This
243 < method has been shown to be reversible and provide results in
244 < excellent agreement with other established methods.\cite{Baez95b}
240 > molecules.  In this study, we applied of one of the most convenient
241 > methods and integrated over the $\lambda^4$ path, where all
242 > interaction parameters are scaled equally by this transformation
243 > parameter.  This method has been shown to be reversible and provide
244 > results in excellent agreement with other established
245 > methods.\cite{Baez95b}
246  
247   Charge, dipole, and Lennard-Jones interactions were modified by a
248   cubic switching between 100\% and 85\% of the cutoff value (9 \AA
# Line 256 | Line 255 | performed with longer cutoffs of 12 and 15 \AA\ to com
255   simulations.\cite{Onsager36} For a series of the least computationally
256   expensive models (SSD/E, SSD/RF, and TIP3P), simulations were
257   performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9
258 < \AA\ cutoff results. Finally, results from the use of an Ewald
259 < summation were estimated for TIP3P and SPC/E by performing
260 < calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
261 < mechanics software package.\cite{Tinker} The calculated energy
262 < difference in the presence and absence of PME was applied to the
263 < previous results in order to predict changes to the free energy
264 < landscape.
258 > \AA\ cutoff results. Finally, the effects of utilizing an Ewald
259 > summation were estimated for TIP3P and SPC/E by performing single
260 > configuration calculations with Particle-Mesh Ewald (PME) in the
261 > TINKER molecular mechanics software package.\cite{Tinker} The
262 > calculated energy difference in the presence and absence of PME was
263 > applied to the previous results in order to predict changes to the
264 > free energy landscape.
265  
266   \section{Results and discussion}
267  
268 < The free energy of proton ordered Ice-{\it i} was calculated and
268 > The free energy of proton-ordered Ice-{\it i} was calculated and
269   compared with the free energies of proton ordered variants of the
270   experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
271   as well as the higher density ice B, observed by B\`{a}ez and Clancy
# Line 276 | Line 275 | antiferroelectric version that we divised, and it has
275   $I_h$, was investigated initially, but was found to be not as stable
276   as proton disordered or antiferroelectric variants of ice $I_h$. The
277   proton ordered variant of ice $I_h$ used here is a simple
278 < antiferroelectric version that we divised, and it has an 8 molecule
278 > antiferroelectric version that we devised, and it has an 8 molecule
279   unit cell similar to other predicted antiferroelectric $I_h$
280   crystals.\cite{Davidson84} The crystals contained 648 or 1728
281   molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
# Line 311 | Line 310 | models studied. With the free energy at these state po
310  
311   The free energy values computed for the studied polymorphs indicate
312   that Ice-{\it i} is the most stable state for all of the common water
313 < models studied. With the free energy at these state points, the
314 < Gibbs-Helmholtz equation was used to project to other state points and
315 < to build phase diagrams.  Figures
313 > models studied. With the calculated free energy at these state points,
314 > the Gibbs-Helmholtz equation was used to project to other state points
315 > and to build phase diagrams.  Figures
316   \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built
317   from the free energy results. All other models have similar structure,
318 < although the crossing points between the phases exist at slightly
318 > although the crossing points between the phases move to slightly
319   different temperatures and pressures. It is interesting to note that
320   ice $I$ does not exist in either cubic or hexagonal form in any of the
321   phase diagrams for any of the models. For purposes of this study, ice
322   B is representative of the dense ice polymorphs. A recent study by
323   Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and
324 < TIP4P in the high pressure regime.\cite{Sanz04}
324 > TIP4P at higher pressures than those studied here.\cite{Sanz04}
325  
326   \begin{figure}
327   \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
# Line 385 | Line 384 | conditions. While troubling, this behavior turned out
384   not exhibit a melting point at 1 atm, but it shows a sublimation point
385   at 355 K. This is due to the significant stability of Ice-{\it i} over
386   all other polymorphs for this particular model under these
387 < conditions. While troubling, this behavior turned out to be
388 < advantageous in that it facilitated the spontaneous crystallization of
389 < Ice-{\it i}. These observations provide a warning that simulations of
387 > conditions. While troubling, this behavior resulted in spontaneous
388 > crystallization of Ice-{\it i} and led us to investigate this
389 > structure. These observations provide a warning that simulations of
390   SSD/E as a ``liquid'' near 300 K are actually metastable and run the
391 < risk of spontaneous crystallization. However, this risk changes when
391 > risk of spontaneous crystallization. However, this risk lessens when
392   applying a longer cutoff.
393  
394   \begin{figure}
# Line 423 | Line 422 | energy of identical crystals with and without PME usin
422   was estimated by applying the potential energy difference do to its
423   inclusion in systems in the presence and absence of the
424   correction. This was accomplished by calculation of the potential
425 < energy of identical crystals with and without PME using TINKER. The
426 < free energies for the investigated polymorphs using the TIP3P and
427 < SPC/E water models are shown in Table \ref{pmeShift}. The same trend
428 < pointed out through increase of cutoff radius is observed in these PME
425 > energy of identical crystals both with and without PME. The free
426 > energies for the investigated polymorphs using the TIP3P and SPC/E
427 > water models are shown in Table \ref{pmeShift}. The same trend pointed
428 > out through increase of cutoff radius is observed in these PME
429   results. Ice-{\it i} is the preferred polymorph at ambient conditions
430   for both the TIP3P and SPC/E water models; however, the narrowing of
431   the free energy differences between the various solid forms is
432   significant enough that it becomes less clear that it is the most
433 < stable polymorph in the SPC/E.  The free energies of Ice-{\it i} and
434 < ice B nearly overlap within error, with ice $I_c$ just outside as
435 < well, indicating that Ice-{\it i} might be metastable with respect to
436 < ice B and possibly ice $I_c$ in the SPC/E water model. However, these
437 < results do not significantly alter the finding that the Ice-{\it i}
438 < polymorph is a stable crystal structure that should be considered when
439 < studying the phase behavior of water models.
433 > stable polymorph with the SPC/E model.  The free energies of Ice-{\it
434 > i} and ice B nearly overlap within error, with ice $I_c$ just outside
435 > as well, indicating that Ice-{\it i} might be metastable with respect
436 > to ice B and possibly ice $I_c$ with SPC/E. However, these results do
437 > not significantly alter the finding that the Ice-{\it i} polymorph is
438 > a stable crystal structure that should be considered when studying the
439 > phase behavior of water models.
440  
441   \begin{table*}
442   \begin{minipage}{\linewidth}
# Line 461 | Line 460 | $I$, ice B, and recently discovered Ice-{\it i} were c
460   \section{Conclusions}
461  
462   The free energy for proton ordered variants of hexagonal and cubic ice
463 < $I$, ice B, and recently discovered Ice-{\it i} were calculated under
464 < standard conditions for several common water models via thermodynamic
465 < integration. All the water models studied show Ice-{\it i} to be the
466 < minimum free energy crystal structure in the with a 9 \AA\ switching
467 < function cutoff. Calculated melting and boiling points show
468 < surprisingly good agreement with the experimental values; however, the
469 < solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of
470 < interaction truncation was investigated through variation of the
471 < cutoff radius, use of a reaction field parameterized model, and
463 > $I$, ice B, and our recently discovered Ice-{\it i} structure were
464 > calculated under standard conditions for several common water models
465 > via thermodynamic integration. All the water models studied show
466 > Ice-{\it i} to be the minimum free energy crystal structure with a 9
467 > \AA\ switching function cutoff. Calculated melting and boiling points
468 > show surprisingly good agreement with the experimental values;
469 > however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The
470 > effect of interaction truncation was investigated through variation of
471 > the cutoff radius, use of a reaction field parameterized model, and
472   estimation of the results in the presence of the Ewald
473   summation. Interaction truncation has a significant effect on the
474   computed free energy values, and may significantly alter the free
# Line 477 | Line 476 | Due to this relative stability of Ice-{\it i} in all m
476   these effects, these results show Ice-{\it i} to be an important ice
477   polymorph that should be considered in simulation studies.
478  
479 < Due to this relative stability of Ice-{\it i} in all manner of
480 < investigated simulation examples, the question arises as to possible
479 > Due to this relative stability of Ice-{\it i} in all of the
480 > investigated simulation conditions, the question arises as to possible
481   experimental observation of this polymorph.  The rather extensive past
482   and current experimental investigation of water in the low pressure
483   regime makes us hesitant to ascribe any relevance of this work outside
# Line 490 | Line 489 | and the structure factor ($S(\vec{q})$ for ice $I_c$ a
489   deposition environments, and in clathrate structures involving small
490   non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
491   our predictions for both the pair distribution function ($g_{OO}(r)$)
492 < and the structure factor ($S(\vec{q})$ for ice $I_c$ and for ice-{\it
493 < i} at a temperature of 77K.  In a quick comparison of the predicted
494 < S(q) for Ice-{\it i} and experimental studies of amorphous solid
495 < water, it is possible that some of the ``spurious'' peaks that could
496 < not be assigned in HDA could correspond to peaks labeled in this
497 < S(q).\cite{Bizid87} It should be noted that there is typically poor
498 < agreement on crystal densities between simulation and experiment, so
499 < such peak comparisons should be made with caution.  We will leave it
501 < to our experimental colleagues to determine whether this ice polymorph
502 < is named appropriately or if it should be promoted to Ice-0.
492 > and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
493 > ice-{\it i} at a temperature of 77K.  In studies of the high and low
494 > density forms of amorphous ice, ``spurious'' diffraction peaks have
495 > been observed experimentally.\cite{Bizid87} It is possible that a
496 > variant of Ice-{\it i} could explain some of this behavior; however,
497 > we will leave it to our experimental colleagues to make the final
498 > determination on whether this ice polymorph is named appropriately
499 > (i.e. with an imaginary number) or if it can be promoted to Ice-0.
500  
501   \begin{figure}
502   \includegraphics[width=\linewidth]{iceGofr.eps}
503 < \caption{Radial distribution functions of Ice-{\it i} and ice $I_c$
504 < calculated from from simulations of the SSD/RF water model at 77 K.}
503 > \caption{Radial distribution functions of ice $I_h$, $I_c$,
504 > Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
505 > of the SSD/RF water model at 77 K.}
506   \label{fig:gofr}
507   \end{figure}
508  
509   \begin{figure}
510   \includegraphics[width=\linewidth]{sofq.eps}
511 < \caption{Predicted structure factors for Ice-{\it i} and ice $I_c$ at
512 < 77 K.  The raw structure factors have been convoluted with a gaussian
513 < instrument function (0.075 \AA$^{-1}$ width) to compensate for the
514 < trunction effects in our finite size simulations. The labeled peaks
515 < compared favorably with ``spurious'' peaks observed in experimental
518 < studies of amorphous solid water.\cite{Bizid87}}
511 > \caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
512 > and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
513 > been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
514 > width) to compensate for the trunction effects in our finite size
515 > simulations.}
516   \label{fig:sofq}
517   \end{figure}
518  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines