[ViewVC] Diff of: group/trunk/iceiPaper/iceiPaper.tex

Comparing trunk/iceiPaper/iceiPaper.tex (file contents):
Revision 1453 by chrisfen, Mon Sep 13 18:39:53 2004 UTC vs.
Revision 1834 by chrisfen, Thu Dec 2 18:58:25 2004 UTC

+%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
-<
+\documentclass[preprint,aps,endfloats]{revtex4}
->
+\documentclass[11pt]{article}
+%\documentclass[11pt]{article}
-<
+%\usepackage{endfloat}
->
+\usepackage{endfloat}
+\usepackage{amsmath}
+\usepackage{epsf}
+\usepackage{berkeley}
-<
+%\usepackage{setspace}
-<
+%\usepackage{tabularx}
->
+\usepackage{setspace}
->
+\usepackage{tabularx}
+\usepackage{graphicx}
-<
+%\usepackage[ref]{overcite}
-<
+%\pagestyle{plain}
-<
+%\pagenumbering{arabic}
-<
+%\oddsidemargin 0.0cm \evensidemargin 0.0cm
-<
+%\topmargin -21pt \headsep 10pt
-<
+%\textheight 9.0in \textwidth 6.5in
-<
+%\brokenpenalty=10000
->
+\usepackage[ref]{overcite}
->
+\pagestyle{plain}
->
+\pagenumbering{arabic}
->
+\oddsidemargin 0.0cm \evensidemargin 0.0cm
->
+\topmargin -21pt \headsep 10pt
->
+\textheight 9.0in \textwidth 6.5in
->
+\brokenpenalty=10000
->
+\renewcommand{\baselinestretch}{1.2}
->
+\renewcommand\citemid{\ } % no comma in optional reference note
-–
+%\renewcommand\citemid{\ } % no comma in optional reference note
-–
+\begin{document}
-<
+\title{A Free Energy Study of Low Temperature and Anomolous Ice}
->
+\title{Free Energy Analysis of Simulated Ice Polymorphs Using Simple
->
+Dipolar and Charge Based Water Models}
-<
+\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
-<
+\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
-<
-<
+\address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
->
+\author{Christopher J. Fennell and J. Daniel Gezelter \\
->
+Department of Chemistry and Biochemistry\\ University of Notre Dame\\
+Notre Dame, Indiana 46556}
+\date{\today}
-<
+%\maketitle
->
+\maketitle
+%\doublespacing
+\begin{abstract}
-+
+The absolute free energies of several ice polymorphs which are stable
-+
+at low pressures were calculated using thermodynamic integration to a
-+
+reference system (the Einstein crystal).  These integrations were
-+
+performed for most of the common water models (SPC/E, TIP3P, TIP4P,
-+
+TIP5P, SSD/E, SSD/RF). Ice-{\it i}, a structure we recently observed
-+
+crystallizing at room temperature for one of the single-point water
-+
+models, was determined to be the stable crystalline state (at 1 atm)
-+
+for {\it all} the water models investigated.  Phase diagrams were
-+
+generated, and phase coexistence lines were determined for all of the
-+
+known low-pressure ice structures under all of these water models.
-+
+Additionally, potential truncation was shown to have an effect on the
-+
+calculated free energies, and can result in altered free energy
-+
+landscapes.  Structure factor predictions for the new crystal were
-+
+generated and we await experimental confirmation of the existence of
-+
+this new polymorph.
+\end{abstract}
-–
+\maketitle
-–
-–
+\newpage
-–
+%\narrowtext
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Introduction}
-+
+Water has proven to be a challenging substance to depict in
-+
+simulations, and a variety of models have been developed to describe
-+
+its behavior under varying simulation
-+
+conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,Berendsen98,Dill00,Mahoney00,Fennell04}
-+
+These models have been used to investigate important physical
-+
+phenomena like phase transitions, transport properties, and the
-+
+hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the
-+
+choice of models available, it is only natural to compare the models
-+
+under interesting thermodynamic conditions in an attempt to clarify
-+
+the limitations of each of the
-+
+models.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two
-+
+important properties to quantify are the Gibbs and Helmholtz free
-+
+energies, particularly for the solid forms of water.  Difficulty in
-+
+these types of studies typically arises from the assortment of
-+
+possible crystalline polymorphs that water adopts over a wide range of
-+
+pressures and temperatures.  There are currently 13 recognized forms
-+
+of ice, and it is a challenging task to investigate the entire free
-+
+energy landscape.\cite{Sanz04} Ideally, research is focused on the
-+
+phases having the lowest free energy at a given state point, because
-+
+these phases will dictate the relevant transition temperatures and
-+
+pressures for the model.
-+
-+
+In this paper, standard reference state methods were applied to known
-+
+crystalline water polymorphs in the low pressure regime.  This work is
-+
+unique in that one of the crystal lattices was arrived at through
-+
+crystallization of a computationally efficient water model under
-+
+constant pressure and temperature conditions.  Crystallization events
-+
+are interesting in and of themselves;\cite{Matsumoto02,Yamada02}
-+
+however, the crystal structure obtained in this case is different from
-+
+any previously observed ice polymorphs in experiment or
-+
+simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
-+
+to indicate its origin in computational simulation. The unit cell
-+
+(Fig. \ref{iceiCell}A) consists of eight water molecules that stack in
-+
+rows of interlocking water tetramers.  This crystal structure has a
-+
+limited resemblence to a recent two-dimensional ice tessellation
-+
+simulated on a silica surface.\cite{Yang04} Proton ordering can be
-+
+accomplished by orienting two of the molecules so that both of their
-+
+donated hydrogen bonds are internal to their tetramer
-+
+(Fig. \ref{protOrder}).  As expected in an ice crystal constructed of
-+
+water tetramers, the hydrogen bonds are not as linear as those
-+
+observed in ice $I_h$, however the interlocking of these subunits
-+
+appears to provide significant stabilization to the overall crystal.
-+
+The arrangement of these tetramers results in surrounding open
-+
+octagonal cavities that are typically greater than 6.3 \AA\ in
-+
+diameter.  This relatively open overall structure leads to crystals
-+
+that are 0.07 g/cm$^3$ less dense on average than ice $I_h$.
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{unitCell.eps}
-+
+\caption{Unit cells for (A) Ice-{\it i} and (B) Ice-{\it i}$^\prime$,
-+
+the elongated variant of Ice-{\it i}.  The spheres represent the
-+
+center-of-mass locations of the water molecules.  The $a$ to $c$
-+
+ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by
-+
+$a:2.1214c$ and $a:1.7850c$ respectively.}
-+
+\label{iceiCell}
-+
+\end{figure}
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{orderedIcei.eps}
-+
+\caption{Image of a proton ordered crystal of Ice-{\it i} looking
-+
+down the (001) crystal face.  The rows of water tetramers surrounded
-+
+by octagonal pores leads to a crystal structure that is significantly
-+
+less dense than ice $I_h$.}
-+
+\label{protOrder}
-+
+\end{figure}
-+
-+
+Results from our previous study indicated that Ice-{\it i} is the
-+
+minimum energy crystal structure for the single point water models we
-+
+had investigated (for discussions on these single point dipole models,
-+
+see our previous work and related
-+
+articles).\cite{Fennell04,Liu96,Bratko85} Those results only
-+
+considered energetic stabilization and neglected entropic
-+
+contributions to the overall free energy.  To address this issue, we
-+
+have calculated the absolute free energy of this crystal using
-+
+thermodynamic integration and compared to the free energies of cubic
-+
+and hexagonal ice $I$ (the experimental low density ice polymorphs)
-+
+and ice B (a higher density, but very stable crystal structure
-+
+observed by B\`{a}ez and Clancy in free energy studies of
-+
+SPC/E).\cite{Baez95b} This work includes results for the water model
-+
+from which Ice-{\it i} was crystallized (SSD/E) in addition to several
-+
+common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
-+
+field parametrized single point dipole water model (SSD/RF). It should
-+
+be noted that a second version of Ice-{\it i} (Ice-{\it i}$^\prime$)
-+
+was used in calculations involving SPC/E, TIP4P, and TIP5P.  The unit
-+
+cell of this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it
-+
+i} unit it is extended in the direction of the (001) face and
-+
+compressed along the other two faces.  There is typically a small
-+
+distortion of proton ordered Ice-{\it i}$^\prime$ that converts the
-+
+normally square tetramer into a rhombus with alternating approximately
-+
+and 95 degree angles.  The degree of this distortion is model
-+
+dependent and significant enough to split the tetramer diagonal
-+
+location peak in the radial distribution function.
-+
+\section{Methods}
-+
-+
+Canonical ensemble (NVT) molecular dynamics calculations were
-+
+performed using the OOPSE molecular mechanics package.\cite{Meineke05}
-+
+All molecules were treated as rigid bodies, with orientational motion
-+
+propagated using the symplectic DLM integration method.  Details about
-+
+the implementation of this technique can be found in a recent
-+
+publication.\cite{Dullweber1997}
-+
-+
+Thermodynamic integration is an established technique for
-+
+determination of free energies of condensed phases of
-+
+materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
-+
+method, implemented in the same manner illustrated by B\`{a}ez and
-+
+Clancy, was utilized to calculate the free energy of several ice
-+
+crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and
-+
+SSD/E water models.\cite{Baez95a} Liquid state free energies at 300
-+
+and 400 K for all of these water models were also determined using
-+
+this same technique in order to determine melting points and to
-+
+generate phase diagrams.  All simulations were carried out at
-+
+densities which correspond to a pressure of approximately 1 atm at
-+
+their respective temperatures.
-+
-+
+Thermodynamic integration involves a sequence of simulations during
-+
+which the system of interest is converted into a reference system for
-+
+which the free energy is known analytically.  This transformation path
-+
+is then integrated in order to determine the free energy difference
-+
+between the two states:
-+
+\begin{equation}
-+
+\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
-+
+)}{\partial\lambda}\right\rangle_\lambda d\lambda,
-+
+\end{equation}
-+
+where $V$ is the interaction potential and $\lambda$ is the
-+
+transformation parameter that scales the overall potential.
-+
+Simulations are distributed strategically along this path in order to
-+
+sufficiently sample the regions of greatest change in the potential.
-+
+Typical integrations in this study consisted of $\sim$25 simulations
-+
+ranging from 300 ps (for the unaltered system) to 75 ps (near the
-+
+reference state) in length.
-+
-+
+For the thermodynamic integration of molecular crystals, the Einstein
-+
+crystal was chosen as the reference system.  In an Einstein crystal,
-+
+the molecules are restrained at their ideal lattice locations and
-+
+orientations. Using harmonic restraints, as applied by B\`{a}ez and
-+
+Clancy, the total potential for this reference crystal
-+
+($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
-+
+\begin{equation}
-+
+V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
-+
+\frac{K_\omega\omega^2}{2},
-+
+\end{equation}
-+
+where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
-+
+the spring constants restraining translational motion and deflection
-+
+of and rotation around the principle axis of the molecule
-+
+respectively.  These spring constants are typically calculated from
-+
+the mean-square displacements of water molecules in an unrestrained
-+
+ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
-+
+mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
-+
+.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
-+
+the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
-+
+from $-\pi$ to $\pi$.  The partition function for a molecular crystal
-+
+restrained in this fashion can be evaluated analytically, and the
-+
+Helmholtz Free Energy ({\it A}) is given by
-+
+\begin{eqnarray}
-+
+A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
-+
+[\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
-+
+)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
-+
+)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
-+
+)^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
-+
+K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
-+
+(-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
-+
+)^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
-+
+\label{ecFreeEnergy}
-+
+\end{eqnarray}
-+
+where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
-+
+potential energy of the ideal crystal.\cite{Baez95a}
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{rotSpring.eps}
-+
+\caption{Possible orientational motions for a restrained molecule.
-+
+$\theta$ angles correspond to displacement from the body-frame {\it
-+
+z}-axis, while $\omega$ angles correspond to rotation about the
-+
+body-frame {\it z}-axis.  $K_\theta$ and $K_\omega$ are spring
-+
+constants for the harmonic springs restraining motion in the $\theta$
-+
+and $\omega$ directions.}
-+
+\label{waterSpring}
-+
+\end{figure}
-+
-+
+In the case of molecular liquids, the ideal vapor is chosen as the
-+
+target reference state.  There are several examples of liquid state
-+
+free energy calculations of water models present in the
-+
+literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
-+
+typically differ in regard to the path taken for switching off the
-+
+interaction potential to convert the system to an ideal gas of water
-+
+molecules.  In this study, we applied of one of the most convenient
-+
+methods and integrated over the $\lambda^4$ path, where all
-+
+interaction parameters are scaled equally by this transformation
-+
+parameter.  This method has been shown to be reversible and provide
-+
+results in excellent agreement with other established
-+
+methods.\cite{Baez95b}
-+
-+
+Charge, dipole, and Lennard-Jones interactions were modified by a
-+
+cubic switching between 100\% and 85\% of the cutoff value (9 \AA).
-+
+By applying this function, these interactions are smoothly truncated,
-+
+thereby avoiding the poor energy conservation which results from
-+
+harsher truncation schemes.  The effect of a long-range correction was
-+
+also investigated on select model systems in a variety of manners.
-+
+For the SSD/RF model, a reaction field with a fixed dielectric
-+
+constant of 80 was applied in all simulations.\cite{Onsager36} For a
-+
+series of the least computationally expensive models (SSD/E, SSD/RF,
-+
+and TIP3P), simulations were performed with longer cutoffs of 12 and
-+
+\AA\ to compare with the 9 \AA\ cutoff results.  Finally, the
-+
+effects of utilizing an Ewald summation were estimated for TIP3P and
-+
+SPC/E by performing single configuration calculations with
-+
+Particle-Mesh Ewald (PME) in the TINKER molecular mechanics software
-+
+package.\cite{Tinker} The calculated energy difference in the presence
-+
+and absence of PME was applied to the previous results in order to
-+
+predict changes to the free energy landscape.
+\section{Results and discussion}
-+
+The free energy of proton-ordered Ice-{\it i} was calculated and
-+
+compared with the free energies of proton ordered variants of the
-+
+experimentally observed low-density ice polymorphs, $I_h$ and $I_c$,
-+
+as well as the higher density ice B, observed by B\`{a}ez and Clancy
-+
+and thought to be the minimum free energy structure for the SPC/E
-+
+model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b}
-+
+Ice XI, the experimentally-observed proton-ordered variant of ice
-+
+$I_h$, was investigated initially, but was found to be not as stable
-+
+as proton disordered or antiferroelectric variants of ice $I_h$.  The
-+
+proton ordered variant of ice $I_h$ used here is a simple
-+
+antiferroelectric version that we devised, and it has an 8 molecule
-+
+unit cell similar to other predicted antiferroelectric $I_h$
-+
+crystals.\cite{Davidson84} The crystals contained 648 or 1728
-+
+molecules for ice B, 1024 or 1280 molecules for ice $I_h$, 1000
-+
+molecules for ice $I_c$, or 1024 molecules for Ice-{\it i}.  The
-+
+larger crystal sizes were necessary for simulations involving larger
-+
+cutoff values.
-+
-+
+\begin{table*}
-+
+\begin{minipage}{\linewidth}
-+
+\begin{center}
-+
-+
+\caption{Calculated free energies for several ice polymorphs with a
-+
+variety of common water models.  All calculations used a cutoff radius
-+
+of 9.0 \AA\ and were performed at 200 K and $\sim$1 atm.  Units are
-+
+kcal/mol.  Calculated error of the final digits is in parentheses.}
-+
-+
+\begin{tabular}{lccccc}
-+
+\hline
-+
+Water Model & $I_h$ & $I_c$ & B & Ice-{\it i} & Ice-{\it i}$^\prime$\\
-+
+\hline
-+
+TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & -\\
-+
+TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3)\\
-+
+TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2)\\
-+
+SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2)\\
-+
+SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & -\\
-+
+SSD/RF & -11.51(2) & -11.47(2) & -12.08(3) & -12.29(2) & -\\
-+
+\end{tabular}
-+
+\label{freeEnergy}
-+
+\end{center}
-+
+\end{minipage}
-+
+\end{table*}
-+
-+
+The free energy values computed for the studied polymorphs indicate
-+
+that Ice-{\it i} is the most stable state for all of the common water
-+
+models studied.  With the calculated free energy at these state
-+
+points, the Gibbs-Helmholtz equation was used to project to other
-+
+state points and to build phase diagrams.  Figures \ref{tp3phasedia}
-+
+and \ref{ssdrfphasedia} are example diagrams built from the free
-+
+energy results.  All other models have similar structure, although the
-+
+crossing points between the phases move to slightly different
-+
+temperatures and pressures.  It is interesting to note that ice $I$
-+
+does not exist in either cubic or hexagonal form in any of the phase
-+
+diagrams for any of the models.  For purposes of this study, ice B is
-+
+representative of the dense ice polymorphs.  A recent study by Sanz
-+
+{\it et al.} goes into detail on the phase diagrams for SPC/E and
-+
+TIP4P at higher pressures than those studied here.\cite{Sanz04}
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{tp3PhaseDia.eps}
-+
+\caption{Phase diagram for the TIP3P water model in the low pressure
-+
+regime.  The displayed $T_m$ and $T_b$ values are good predictions of
-+
+the experimental values; however, the solid phases shown are not the
-+
+experimentally observed forms.  Both cubic and hexagonal ice $I$ are
-+
+higher in energy and don't appear in the phase diagram.}
-+
+\label{tp3phasedia}
-+
+\end{figure}
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{ssdrfPhaseDia.eps}
-+
+\caption{Phase diagram for the SSD/RF water model in the low pressure
-+
+regime.  Calculations producing these results were done under an
-+
+applied reaction field.  It is interesting to note that this
-+
+computationally efficient model (over 3 times more efficient than
-+
+TIP3P) exhibits phase behavior similar to the less computationally
-+
+conservative charge based models.}
-+
+\label{ssdrfphasedia}
-+
+\end{figure}
-+
-+
+\begin{table*}
-+
+\begin{minipage}{\linewidth}
-+
+\begin{center}
-+
-+
+\caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
-+
+temperatures at 1 atm for several common water models compared with
-+
+experiment.  The $T_m$ and $T_s$ values from simulation correspond to
-+
+a transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
-+
+liquid or gas state.}
-+
-+
+\begin{tabular}{lccccccc}
-+
+\hline
-+
+Equilibrium Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
-+
+\hline
-+
+$T_m$ (K)  & 269(4) & 266(5) & 271(4) & 296(3) & - & 278(4) & 273\\
-+
+$T_b$ (K)  & 357(2) & 354(2) & 337(2) & 396(2) & - & 348(2) & 373\\
-+
+$T_s$ (K)  & - & - & - & - & 355(2) & - & -\\
-+
+\end{tabular}
-+
+\label{meltandboil}
-+
+\end{center}
-+
+\end{minipage}
-+
+\end{table*}
-+
-+
+Table \ref{meltandboil} lists the melting and boiling temperatures
-+
+calculated from this work.  Surprisingly, most of these models have
-+
+melting points that compare quite favorably with experiment.  The
-+
+unfortunate aspect of this result is that this phase change occurs
-+
+between Ice-{\it i} and the liquid state rather than ice $I_h$ and the
-+
+liquid state.  These results are actually not contrary to other
-+
+studies.  Studies of ice $I_h$ using TIP4P predict a $T_m$ ranging
-+
+from 214 to 238 K (differences being attributed to choice of
-+
+interaction truncation and different ordered and disordered molecular
-+
+arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
-+
+Ice-{\it i} were omitted, a $T_m$ value around 210 K would be
-+
+predicted from this work.  However, the $T_m$ from Ice-{\it i} is
-+
+calculated to be 265 K, indicating that these simulation based
-+
+structures ought to be included in studies probing phase transitions
-+
+with this model.  Also of interest in these results is that SSD/E does
-+
+not exhibit a melting point at 1 atm, but it shows a sublimation point
-+
+at 355 K.  This is due to the significant stability of Ice-{\it i}
-+
+over all other polymorphs for this particular model under these
-+
+conditions.  While troubling, this behavior resulted in spontaneous
-+
+crystallization of Ice-{\it i} and led us to investigate this
-+
+structure.  These observations provide a warning that simulations of
-+
+SSD/E as a ``liquid'' near 300 K are actually metastable and run the
-+
+risk of spontaneous crystallization.  However, this risk lessens when
-+
+applying a longer cutoff.
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{cutoffChange.eps}
-+
+\caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
-+
+SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
-+
+with an added Ewald correction term.  Calculations performed without a
-+
+long-range correction show noticable free energy dependence on the
-+
+cutoff radius and show some degree of converge at large cutoff radii.
-+
+Inclusion of a long-range correction reduces the cutoff radius
-+
+dependence of the free energy for all the models.  Error for the
-+
+larger cutoff points is equivalent to that observed at 9.0 \AA\ (see
-+
+Table \ref{freeEnergy}).  Data for ice I$_c$ with TIP3P using both 12
-+
+and 13.5 \AA\ cutoffs were omitted because the crystal was prone to
-+
+distortion and melting at 200 K.  Ice-{\it i}$^\prime$ is the form of
-+
+Ice-{\it i} used in the SPC/E simulations.}
-+
+\label{incCutoff}
-+
+\end{figure}
-+
-+
+Increasing the cutoff radius in simulations of the more
-+
+computationally efficient water models was done in order to evaluate
-+
+the trend in free energy values when moving to systems that do not
-+
+involve potential truncation.  As seen in Fig. \ref{incCutoff}, the
-+
+free energy of the ice polymorphs with water models lacking a
-+
+long-range correction show a significant cutoff radius dependence.  In
-+
+general, there is a narrowing of the free energy differences while
-+
+moving to greater cutoff radii.  As the free energies for the
-+
+polymorphs converge, the stability advantage that Ice-{\it i} exhibits
-+
+is reduced.  Interestingly, increasing the cutoff radius a mere 1.5
-+
+\AA\ with the SSD/E model destabilizes the Ice-{\it i} polymorph
-+
+enough that the liquid state is preferred under standard simulation
-+
+conditions (298 K and 1 atm).  Thus, it is recommended that
-+
+simulations using this model choose interaction truncation radii
-+
+greater than 9 \AA.  Considering the stabilization of Ice-{\it i} with
-+
+smaller cutoffs, it is not surprising that crystallization was
-+
+observed with SSD/E.  The choice of a 9 \AA\ cutoff in the previous
-+
+simulations gives the Ice-{\it i} polymorph a greater than 1 kcal/mol
-+
+lower free energy than the ice $I_\textrm{h}$ starting configurations.
-+
+Additionally, it should be noted that ice $I_c$ is not stable with
-+
+cutoff radii of 12 and 13.5 \AA\ using the TIP3P water model.  These
-+
+simulations showed bulk distortions of the simulation cell that
-+
+rapidly deteriorated crystal integrity.
-+
-+
+Adjacent to each of these model plots is a system with an applied or
-+
+estimated long-range correction.  SSD/RF was parametrized for use with
-+
+a reaction field, and the benefit provided by this computationally
-+
+inexpensive correction is apparent.  Due to the relative independence
-+
+of the resultant free energies, calculations performed with a small
-+
+cutoff radius provide resultant properties similar to what one would
-+
+expect for the bulk material.  In the cases of TIP3P and SPC/E, the
-+
+effect of an Ewald summation was estimated by applying the potential
-+
+energy difference do to its inclusion in systems in the presence and
-+
+absence of the correction.  This was accomplished by calculation of
-+
+the potential energy of identical crystals both with and without
-+
+particle mesh Ewald (PME).  Similar behavior to that observed with
-+
+reaction field is seen for both of these models.  The free energies
-+
+show less dependence on cutoff radius and span a more narrowed range
-+
+for the various polymorphs.  Like the dipolar water models, TIP3P
-+
+displays a relatively constant preference for the Ice-{\it i}
-+
+polymorph.  Crystal preference is much more difficult to determine for
-+
+SPC/E.  Without a long-range correction, each of the polymorphs
-+
+studied assumes the role of the preferred polymorph under different
-+
+cutoff conditions.  The inclusion of the Ewald correction flattens and
-+
+narrows the sequences of free energies so much that they often overlap
-+
+within error, indicating that other conditions, such as cell volume in
-+
+microcanonical simulations, can influence the chosen polymorph upon
-+
+crystallization.  All of these results support the finding that the
-+
+Ice-{\it i} polymorph is a stable crystal structure that should be
-+
+considered when studying the phase behavior of water models.
-+
+\section{Conclusions}
-+
+The free energy for proton ordered variants of hexagonal and cubic ice
-+
+$I$, ice B, and our recently discovered Ice-{\it i} structure were
-+
+calculated under standard conditions for several common water models
-+
+via thermodynamic integration.  All the water models studied show
-+
+Ice-{\it i} to be the minimum free energy crystal structure with a 9
-+
+\AA\ switching function cutoff.  Calculated melting and boiling points
-+
+show surprisingly good agreement with the experimental values;
-+
+however, the solid phase at 1 atm is Ice-{\it i}, not ice $I_h$.  The
-+
+effect of interaction truncation was investigated through variation of
-+
+the cutoff radius, use of a reaction field parameterized model, and
-+
+estimation of the results in the presence of the Ewald summation.
-+
+Interaction truncation has a significant effect on the computed free
-+
+energy values, and may significantly alter the free energy landscape
-+
+for the more complex multipoint water models.  Despite these effects,
-+
+these results show Ice-{\it i} to be an important ice polymorph that
-+
+should be considered in simulation studies.
-+
-+
+Due to this relative stability of Ice-{\it i} in all of the
-+
+investigated simulation conditions, the question arises as to possible
-+
+experimental observation of this polymorph.  The rather extensive past
-+
+and current experimental investigation of water in the low pressure
-+
+regime makes us hesitant to ascribe any relevance of this work outside
-+
+of the simulation community.  It is for this reason that we chose a
-+
+name for this polymorph which involves an imaginary quantity.  That
-+
+said, there are certain experimental conditions that would provide the
-+
+most ideal situation for possible observation. These include the
-+
+negative pressure or stretched solid regime, small clusters in vacuum
-+
+deposition environments, and in clathrate structures involving small
-+
+non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
-+
+our predictions for both the pair distribution function ($g_{OO}(r)$)
-+
+and the structure factor ($S(\vec{q})$ for ice $I_h$, $I_c$, and for
-+
+ice-{\it i} at a temperature of 77K.  In studies of the high and low
-+
+density forms of amorphous ice, ``spurious'' diffraction peaks have
-+
+been observed experimentally.\cite{Bizid87} It is possible that a
-+
+variant of Ice-{\it i} could explain some of this behavior; however,
-+
+we will leave it to our experimental colleagues to make the final
-+
+determination on whether this ice polymorph is named appropriately
-+
+(i.e. with an imaginary number) or if it can be promoted to Ice-0.
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{iceGofr.eps}
-+
+\caption{Radial distribution functions of ice $I_h$, $I_c$,
-+
+Ice-{\it i}, and Ice-{\it i}$^\prime$ calculated from from simulations
-+
+of the SSD/RF water model at 77 K.}
-+
+\label{fig:gofr}
-+
+\end{figure}
-+
-+
+\begin{figure}
-+
+\includegraphics[width=\linewidth]{sofq.eps}
-+
+\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
-+
+ and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
-+
+ been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
-+
+ width) to compensate for the trunction effects in our finite size
-+
+ simulations.}
-+
+\label{fig:sofq}
-+
+\end{figure}
-+
+\section{Acknowledgments}
+Support for this project was provided by the National Science
+Foundation under grant CHE-0134881. Computation time was provided by
-<
+the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
-<
+DMR-0079647.
->
+the Notre Dame High Performance Computing Cluster and the Notre Dame
->
+Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647).
+\newpage

Diff Legend

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

Comparing trunk/iceiPaper/iceiPaper.tex (file contents): Revision 1453 by chrisfen, Mon Sep 13 18:39:53 2004 UTC vs. Revision 1834 by chrisfen, Thu Dec 2 18:58:25 2004 UTC

Diff Legend

Comparing trunk/iceiPaper/iceiPaper.tex (file contents):
Revision 1453 by chrisfen, Mon Sep 13 18:39:53 2004 UTC vs.
Revision 1834 by chrisfen, Thu Dec 2 18:58:25 2004 UTC