ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/scienceIcei/iceiPaper.tex
Revision: 1888
Committed: Wed Dec 15 18:47:55 2004 UTC (19 years, 6 months ago) by chrisfen
Content type: application/x-tex
File size: 17208 byte(s)
Log Message:
fixed grammer and typos

File Contents

# User Rev Content
1 chrisfen 1845 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2     \documentclass[11pt]{article}
3     \usepackage{endfloat}
4     \usepackage{amsmath}
5     \usepackage{epsf}
6     \usepackage{berkeley}
7     \usepackage{setspace}
8     \usepackage{tabularx}
9     \usepackage{graphicx}
10     \usepackage[ref]{overcite}
11     \pagestyle{plain}
12     \pagenumbering{arabic}
13     \oddsidemargin 0.0cm \evensidemargin 0.0cm
14     \topmargin -21pt \headsep 10pt
15     \textheight 9.0in \textwidth 6.5in
16     \brokenpenalty=10000
17     \renewcommand{\baselinestretch}{1.2}
18     \renewcommand\citemid{\ } % no comma in optional reference note
19    
20     \begin{document}
21    
22     \title{Free Energy Analysis of Simulated Ice Polymorphs}
23    
24     \author{Christopher J. Fennell and J. Daniel Gezelter \\
25     Department of Chemistry and Biochemistry\\ University of Notre Dame\\
26     Notre Dame, Indiana 46556}
27    
28     \date{\today}
29    
30     \maketitle
31     %\doublespacing
32    
33     \begin{abstract}
34     The absolute free energies of several ice polymorphs which are stable
35     at low pressures were calculated using thermodynamic integration with
36     a variety of common water models. A recently discovered ice polymorph
37 chrisfen 1888 that has as yet only been observed in computer simulations (Ice-{\it
38     i}), was determined to be the stable crystalline state for {\it all}
39     the water models investigated. Phase diagrams were generated, and
40     phase coexistence lines were determined for all of the known
41     low-pressure ice structures. Additionally, potential truncation was
42     shown to play a role in the resulting shape of the free energy
43     landscape.
44 chrisfen 1845 \end{abstract}
45    
46     %\narrowtext
47    
48     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49     % BODY OF TEXT
50     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51    
52     \section{Introduction}
53    
54     Water has proven to be a challenging substance to depict in
55     simulations, and a variety of models have been developed to describe
56     its behavior under varying simulation
57     conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
58     These models have been used to investigate important physical
59     phenomena like phase transitions, transport properties, and the
60     hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
61     choice of models available, it is only natural to compare the models
62     under interesting thermodynamic conditions in an attempt to clarify
63 chrisfen 1888 the limitations of
64     each.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two important
65     properties to quantify are the Gibbs and Helmholtz free energies,
66     particularly for the solid forms of water. Difficulties in studies
67     addressing these thermodynamic quantities typically arise from the
68     assortment of possible crystalline polymorphs that water adopts over a
69     wide range of pressures and temperatures. It is a challenging task to
70     investigate the entire free energy landscape\cite{Sanz04}; and
71     ideally, research is focused on the phases having the lowest free
72 chrisfen 1845 energy at a given state point, because these phases will dictate the
73     relevant transition temperatures and pressures for the model.
74    
75     In this paper, standard reference state methods were applied to known
76 chrisfen 1860 crystalline water polymorphs to evaluate their free energy in the low
77     pressure regime. This work is unique in that one of the crystal
78     lattices was arrived at through crystallization of a computationally
79     efficient water model under constant pressure and temperature
80     conditions. Crystallization events are interesting in and of
81     themselves\cite{Matsumoto02,Yamada02}; however, the crystal structure
82     obtained in this case is different from any previously observed ice
83     polymorphs in experiment or simulation.\cite{Fennell04} We have named
84     this structure Ice-{\it i} to indicate its origin in computational
85     simulation. The unit cell of Ice-{\it i} and an extruded variant named
86     Ice-{\it i}$^\prime$ both consist of eight water molecules that stack
87     in rows of interlocking water tetramers as illustrated in figures
88     \ref{iCrystal}A and
89 chrisfen 1845 \ref{iCrystal}B. These tetramers make the crystal structure similar
90     in appearance to a recent two-dimensional ice tessellation simulated
91     on a silica surface.\cite{Yang04} As expected in an ice crystal
92     constructed of water tetramers, the hydrogen bonds are not as linear
93     as those observed in ice $I_h$, however the interlocking of these
94     subunits appears to provide significant stabilization to the overall
95     crystal. The arrangement of these tetramers results in surrounding
96     open octagonal cavities that are typically greater than 6.3 \AA\ in
97     diameter (Fig. \ref{iCrystal}C). This open structure leads to
98     crystals that are typically 0.07 g/cm$^3$ less dense than ice $I_h$.
99    
100     \begin{figure}
101     \includegraphics[width=4in]{iCrystal.eps}
102     \caption{(A) Unit cell for Ice-{\it i}, (B) Ice-{\it i}$^\prime$,
103     and (C) a rendering of a proton ordered crystal of Ice-{\it i} looking
104     down the (001) crystal face. In the unit cells, the spheres represent
105     the center-of-mass locations of the water molecules. The $a$ to $c$
106     ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
107     $a:2.1214c$ and $a:1.785c$ respectively. The presence of large
108     octagonal pores in both crystal forms lead to a polymorph that is less
109     dense than ice $I_h$.}
110     \label{iCrystal}
111     \end{figure}
112    
113     Results from our previous study indicated that Ice-{\it i} is the
114     minimum energy crystal structure for the single point water models
115     investigated (for discussions on these single point dipole models, see
116     our previous work and related
117     articles).\cite{Fennell04,Liu96,Bratko85} Those results only
118     considered energetic stabilization and neglected entropic
119     contributions to the overall free energy. To address this issue, we
120     have calculated the absolute free energy of this crystal using
121 chrisfen 1888 thermodynamic integration and compared it to the free energies of ice
122     $I_c$ and ice $I_h$ (the experimental low density ice polymorphs) and
123     ice B (a higher density, but very stable crystal structure observed by
124     B\`{a}ez and Clancy in free energy studies of SPC/E).\cite{Baez95b}
125     This work includes results for the water model from which Ice-{\it i}
126     was crystallized (SSD/E) in addition to several common water models
127     (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction field parametrized
128     single point dipole water model (SSD/RF). The extruded variant,
129     Ice-{\it i}$^\prime$, was used in calculations involving SPC/E, TIP4P,
130     and TIP5P. These models exhibit enhanced stability with Ice-{\it
131     i}$^\prime$ because of their more tetrahedrally arranged internal
132     charge distributions. Additionally, there is often a small distortion
133     of proton ordered Ice-{\it i}$^\prime$ that converts the normally
134     square tetramer into a rhombus with alternating approximately 85 and
135     95 degree angles. The degree of this distortion is model dependent
136     and significant enough to split the tetramer diagonal location peak in
137     the radial distribution function.
138 chrisfen 1845
139     Thermodynamic integration was utilized to calculate the free energies
140     of the listed water models at various state points using a modified
141     form of the OOPSE molecular dynamics package.\cite{Meineke05}
142     This calculation method involves a sequence of simulations during
143     which the system of interest is converted into a reference system for
144     which the free energy is known analytically. This transformation path
145 chrisfen 1872 is then integrated, in order to determine the free energy difference
146 chrisfen 1845 between the two states:
147     \begin{equation}
148     \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
149     )}{\partial\lambda}\right\rangle_\lambda d\lambda,
150     \end{equation}
151     where $V$ is the interaction potential and $\lambda$ is the
152     transformation parameter that scales the overall potential. For
153     liquid and solid phases, the ideal gas and harmonically restrained
154     crystal are chosen as the reference states respectively. Thermodynamic
155     integration is an established technique that has been used extensively
156     in the calculation of free energies for condensed phases of
157     materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}.
158    
159 chrisfen 1888 The calculated free energies of proton-ordered variants of three low
160 chrisfen 1845 density polymorphs ($I_h$, $I_c$, and Ice-{\it i} or Ice-{\it
161     i}$^\prime$) and the stable higher density ice B are listed in Table
162 chrisfen 1888 \ref{freeEnergy}. Ice B was included because it has been
163 chrisfen 1845 shown to be a minimum free energy structure for SPC/E at ambient
164     conditions.\cite{Baez95b} In addition to the free energies, the
165 chrisfen 1888 relevant transition temperatures at standard pressure are also
166 chrisfen 1845 displayed in Table \ref{freeEnergy}. These free energy values
167     indicate that Ice-{\it i} is the most stable state for all of the
168     investigated water models. With the free energy at these state
169     points, the Gibbs-Helmholtz equation was used to project to other
170 chrisfen 1888 state points and to build phase diagrams. Figure \ref{tp3PhaseDia} is
171     an example diagram built from the results for the TIP3P water model.
172     All other models have similar structure, although the crossing points
173     between the phases move to different temperatures and pressures as
174     indicated from the transition temperatures in Table \ref{freeEnergy}.
175     It is interesting to note that ice $I_h$ (and ice $I_c$ for that
176     matter) do not appear in any of the phase diagrams for any of the
177     models. For purposes of this study, ice B is representative of the
178     dense ice polymorphs. A recent study by Sanz {\it et al.} provides
179     details on the phase diagrams for SPC/E and TIP4P at higher pressures
180     than those studied here.\cite{Sanz04}
181 chrisfen 1845
182     \begin{table*}
183     \begin{minipage}{\linewidth}
184     \begin{center}
185     \caption{Calculated free energies for several ice polymorphs along
186     with the calculated melting (or sublimation) and boiling points for
187     the investigated water models. All free energy calculations used a
188     cutoff radius of 9.0 \AA\ and were performed at 200 K and $\sim$1 atm.
189     Units of free energy are kcal/mol, while transition temperature are in
190     Kelvin. Calculated error of the final digits is in parentheses.}
191     \begin{tabular}{lccccccc}
192     \hline
193     Water Model & $I_h$ & $I_c$ & B & Ice-{\it i} & Ice-{\it i}$^\prime$ & $T_m$ (*$T_s$) & $T_b$\\
194     \hline
195     TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(4) & 357(2)\\
196     TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 266(5) & 354(2)\\
197     TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 271(4) & 337(2)\\
198     SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 296(3) & 396(2)\\
199     SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(2) & -\\
200     SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2) & - & 278(4) & 349(2)\\
201     \end{tabular}
202     \label{freeEnergy}
203     \end{center}
204     \end{minipage}
205     \end{table*}
206    
207     \begin{figure}
208     \includegraphics[width=\linewidth]{tp3PhaseDia.eps}
209     \caption{Phase diagram for the TIP3P water model in the low pressure
210     regime. The displayed $T_m$ and $T_b$ values are good predictions of
211     the experimental values; however, the solid phases shown are not the
212     experimentally observed forms. Both cubic and hexagonal ice $I$ are
213     higher in energy and don't appear in the phase diagram.}
214     \label{tp3PhaseDia}
215     \end{figure}
216    
217     Most of the water models have melting points that compare quite
218     favorably with the experimental value of 273 K. The unfortunate
219     aspect of this result is that this phase change occurs between
220     Ice-{\it i} and the liquid state rather than ice $I_h$ and the liquid
221     state. Surprisingly, these results are not contrary to other studies.
222     Studies of ice $I_h$ using TIP4P predict a $T_m$ ranging from 214 to
223     238 K (differences being attributed to choice of interaction
224     truncation and different ordered and disordered molecular
225     arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
226     Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
227     predicted from this work. However, the $T_m$ from Ice-{\it i} is
228     calculated to be 265 K, indicating that these simulation based
229     structures ought to be included in studies probing phase transitions
230     with this model. Also of interest in these results is that SSD/E does
231     not exhibit a melting point at 1 atm, but it rather shows a
232     sublimation point at 355 K. This is due to the significant stability
233     of Ice-{\it i} over all other polymorphs for this particular model
234     under these conditions. While troubling, this behavior resulted in
235     spontaneous crystallization of Ice-{\it i} and led us to investigate
236     this structure. These observations provide a warning that simulations
237     of SSD/E as a ``liquid'' near 300 K are actually metastable and run
238     the risk of spontaneous crystallization. However, when applying a
239     longer cutoff, the liquid state is preferred under standard
240     conditions.
241    
242     \begin{figure}
243     \includegraphics[width=\linewidth]{cutoffChange.eps}
244     \caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
245     SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
246     with an added Ewald correction term. Error for the larger cutoff
247     points is equivalent to that observed at 9.0\AA\ (see Table
248     \ref{freeEnergy}). Data for ice I$_c$ with TIP3P using both 12 and
249     13.5 \AA\ cutoffs were omitted because the crystal was prone to
250     distortion and melting at 200 K. Ice-{\it i}$^\prime$ is the form of
251     Ice-{\it i} used in the SPC/E simulations.}
252     \label{incCutoff}
253     \end{figure}
254    
255     Increasing the cutoff radius in simulations of the more
256     computationally efficient water models was done in order to evaluate
257     the trend in free energy values when moving to systems that do not
258     involve potential truncation. As seen in Fig. \ref{incCutoff}, the
259     free energy of the ice polymorphs with water models lacking a
260     long-range correction show a significant cutoff radius dependence. In
261     general, there is a narrowing of the free energy differences while
262     moving to greater cutoff radii. As the free energies for the
263     polymorphs converge, the stability advantage that Ice-{\it i} exhibits
264     is reduced. Adjacent to each of these model plots is a system with an
265     applied or estimated long-range correction. SSD/RF was parametrized
266     for use with a reaction field, and the benefit provided by this
267     computationally inexpensive correction is apparent. Due to the
268     relative independence of the resultant free energies, calculations
269     performed with a small cutoff radius provide resultant properties
270     similar to what one would expect for the bulk material. In the cases
271     of TIP3P and SPC/E, the effect of an Ewald summation was estimated by
272     applying the potential energy difference do to its inclusion in
273     systems in the presence and absence of the correction. This was
274     accomplished by calculation of the potential energy of identical
275     crystals both with and without particle mesh Ewald (PME). Similar
276     behavior to that observed with reaction field is seen for both of
277     these models. The free energies show less dependence on cutoff radius
278     and span a more narrowed range for the various polymorphs. Like the
279     dipolar water models, TIP3P displays a relatively constant preference
280     for the Ice-{\it i} polymorph. Crystal preference is much more
281     difficult to determine for SPC/E. Without a long-range correction,
282     each of the polymorphs studied assumes the role of the preferred
283     polymorph under different cutoff conditions. The inclusion of the
284     Ewald correction flattens and narrows the sequences of free energies
285 chrisfen 1872 such that they often overlap within error, indicating that other
286     conditions, such as the density in fixed volume simulations, can
287 chrisfen 1860 influence the chosen polymorph upon crystallization.
288 chrisfen 1845
289 chrisfen 1860 So what is the preferred solid polymorph for simulated water? The
290 chrisfen 1888 answer appears to be dependent both on the conditions and the model
291 chrisfen 1860 used. In the case of short cutoffs without a long-range interaction
292     correction, Ice-{\it i} and Ice-{\it i}$^\prime$ have the lowest free
293     energy of the studied polymorphs with all the models. Ideally,
294     crystallization of each model under constant pressure conditions, as
295     was done with SSD/E, would aid in the identification of their
296     respective preferred structures. This work, however, helps illustrate
297     how studies involving one specific model can lead to insight about
298     important behavior of others. In general, the above results support
299     the finding that the Ice-{\it i} polymorph is a stable crystal
300     structure that should be considered when studying the phase behavior
301     of water models.
302    
303     Finally, due to the stability of Ice-{\it i} in the investigated
304     simulation conditions, the question arises as to possible experimental
305     observation of this polymorph. The rather extensive past and current
306     experimental investigation of water in the low pressure regime makes
307 chrisfen 1888 us hesitant to ascribe any relevance to this work outside of the
308 chrisfen 1860 simulation community. It is for this reason that we chose a name for
309     this polymorph which involves an imaginary quantity. That said, there
310     are certain experimental conditions that would provide the most ideal
311     situation for possible observation. These include the negative
312     pressure or stretched solid regime, small clusters in vacuum
313 chrisfen 1845 deposition environments, and in clathrate structures involving small
314 chrisfen 1860 non-polar molecules.
315 chrisfen 1845
316     \section{Acknowledgments}
317     Support for this project was provided by the National Science
318     Foundation under grant CHE-0134881. Computation time was provided by
319     the Notre Dame High Performance Computing Cluster and the Notre Dame
320     Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
321    
322     \newpage
323    
324     \bibliographystyle{jcp}
325     \bibliography{iceiPaper}
326    
327    
328     \end{document}