ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/scienceIcei/iceiPaper.tex
Revision: 1872
Committed: Thu Dec 9 20:23:48 2004 UTC (19 years, 6 months ago) by chrisfen
Content type: application/x-tex
File size: 17236 byte(s)
Log Message:
adjustments and inclusion of supplemental

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