--- trunk/iceiPaper/iceiPaper.tex 2004/09/14 21:55:24 1456 +++ trunk/iceiPaper/iceiPaper.tex 2004/09/15 21:44:27 1464 @@ -1,45 +1,50 @@ %\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{Ice-{\it i}: a novel ice polymorph predicted via computer simulation} -\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 free energies of several ice polymorphs in the low pressure regime +were calculated using thermodynamic integration. These integrations +were done for most of the common water models. Ice-{\it i}, a +structure we recently observed to be stable in 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 the common water +models. Additionally, potential truncation was shown to have an +effect on the calculated free energies, and can result in altered free +energy landscapes. \end{abstract} -\maketitle - -\newpage - %\narrowtext %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -48,12 +53,103 @@ Notre Dame, Indiana 46556} \section{Introduction} +Molecular dynamics is a valuable tool for studying the phase behavior +of systems ranging from small or simple +molecules\cite{Matsumoto02andOthers} to complex biological +species.\cite{bigStuff} Many techniques have been developed to +investigate the thermodynamic properites of model substances, +providing both qualitative and quantitative comparisons between +simulations and experiment.\cite{thermMethods} Investigation of these +properties leads to the development of new and more accurate models, +leading to better understanding and depiction of physical processes +and intricate molecular systems. + +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{Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Mahoney00,Fennell04} +These models have been used to investigate important physical +phenomena like phase transitions and the hydrophobic +effect.\cite{Yamada02} 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{modelProps} Two important property 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 true transition +temperatures and pressures for their respective model. + +In this paper, standard reference state methods were applied to the +study of crystalline water polymorphs in the low pressure regime. This +work is unique in the fact 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 was different from any previously observed ice +polymorphs, in experiment or simulation.\cite{Fennell04} This crystal +was termed Ice-{\it i} in homage to 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. Proton ordering can be accomplished by orienting two of the +waters 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-2{\it i}, the elongated variant of Ice-{\it i}. For Ice-{\it i}, the $a$ to $c$ relation is given by $a = 2.1214c$, while for Ice-2{\it i}, $a = 1.7850c$.} +\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 in the previous study indicated that Ice-{\it i} is the +minimum energy crystal structure for the single point water models +being studied (for discussions on these single point dipole models, +see the previous work and related +articles\cite{Fennell04,Ichiye96,Bratko85}). Those results only +consider energetic stabilization and neglect entropic contributions to +the overall free energy. To address this issue, the absolute free +energy of this crystal was calculated 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 +(soft sticky dipole extended, 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 (soft sticky dipole +reaction field, SSD/RF). In should be noted that a second version of +Ice-{\it i} (Ice-2{\it i}) 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. + \section{Methods} Canonical ensemble (NVT) molecular dynamics calculations were performed using the OOPSE (Object-Oriented Parallel Simulation Engine) molecular mechanics package. All molecules were treated as rigid -bodies, with orientational motion propogated using the symplectic DLM +bodies, with orientational motion propagated using the symplectic DLM integration method. Details about the implementation of these techniques can be found in a recent publication.\cite{Meineke05} @@ -66,6 +162,23 @@ For the thermodynamic integration of molecular crystal resulting in a pressure of approximately 1 atm at their respective temperatures. +A single thermodynamic integration involves a sequence of simulations +over which the system of interest is converted into a reference system +for which the free energy is known. 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 unevenly 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 is chosen as the reference state that the system is converted to over the course of the simulation. In an Einstein Crystal, the @@ -92,8 +205,9 @@ state. minimum potential energy of the ideal crystal. In the case of molecular liquids, the ideal vapor is chosen as the target reference state. + \begin{figure} -\includegraphics[scale=1.0]{rotSpring.eps} +\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 @@ -104,23 +218,25 @@ cubic switching between 100\% and 85\% of the cutoff v \end{figure} 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 poor energy conserving dynamics resulting 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 15 \AA\ to -compare with the 9 \AA\ cutoff results. Finally, results from the use -of an Ewald summation were estimated for TIP3P and SPC/E by performing +cubic switching between 100\% and 85\% of the cutoff value (9 \AA +). By applying this function, these interactions are smoothly +truncated, thereby avoiding poor energy conserving dynamics resulting +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 15 \AA\ to compare with the 9 +\AA\ cutoff results. Finally, results from the use of an Ewald +summation were estimated for TIP3P and SPC/E by performing calculations with Particle-Mesh Ewald (PME) in the TINKER molecular -mechanics software package. TINKER was chosen because it can also -propogate the motion of rigid-bodies, and provides the most direct -comparison to the results from OOPSE. The calculated energy difference -in the presence and absence of PME was applied to the previous results -in order to predict changes in the free energy landscape. +mechanics software package.\cite{Tinker} TINKER was chosen because it +can also propagate the motion of rigid-bodies, and provides the most +direct comparison to the results from OOPSE. The calculated energy +difference in the presence and absence of PME was applied to the +previous results in order to predict changes in the free energy +landscape. \section{Results and discussion} @@ -149,9 +265,9 @@ kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are kcal/mol. *Ice $I_c$ is unstable at 200 K using SSD/RF.} \begin{tabular}{ l c c c c } -\hline \\[-7mm] +\hline \ \quad \ Water Model\ \ & \ \quad \ \ \ \ $I_h$ \ \ & \ \quad \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \quad \ \ \ Ice-{\it i} \ \quad \ \\ -\hline \\[-3mm] +\hline \ \quad \ TIP3P & \ \quad \ -11.41 & \ \quad \ -11.23 & \ \quad \ -11.82 & \quad -12.30\\ \ \quad \ TIP4P & \ \quad \ -11.84 & \ \quad \ -12.04 & \ \quad \ -12.08 & \quad -12.33\\ \ \quad \ TIP5P & \ \quad \ -11.85 & \ \quad \ -11.86 & \ \quad \ -11.96 & \quad -12.29\\ @@ -178,6 +294,7 @@ TIP4P in the high pressure regime.\cite{Sanz04} 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 in the high pressure regime.\cite{Sanz04} + \begin{figure} \includegraphics[width=\linewidth]{tp3PhaseDia.eps} \caption{Phase diagram for the TIP3P water model in the low pressure @@ -187,6 +304,7 @@ higher in energy and don't appear in the phase diagram 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 @@ -205,9 +323,9 @@ temperatures of several common water models compared w \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) temperatures of several common water models compared with experiment.} \begin{tabular}{ l c c c c c c c } -\hline \\[-7mm] +\hline \ \ Equilibria Point\ \ & \ \ \ \ \ TIP3P \ \ & \ \ \ \ \ TIP4P \ \ & \ \quad \ \ \ \ TIP5P \ \ & \ \ \ \ \ SPC/E \ \ & \ \ \ \ \ SSD/E \ \ & \ \ \ \ \ SSD/RF \ \ & \ \ \ \ \ Exp. \ \ \\ -\hline \\[-3mm] +\hline \ \ $T_m$ (K) & \ \ 269 & \ \ 265 & \ \ 271 & 297 & \ \ - & \ \ 278 & \ \ 273\\ \ \ $T_b$ (K) & \ \ 357 & \ \ 354 & \ \ 337 & 396 & \ \ - & \ \ 349 & \ \ 373\\ \ \ $T_s$ (K) & \ \ - & \ \ - & \ \ - & - & \ \ 355 & \ \ - & \ \ -\\ @@ -235,21 +353,114 @@ advantagious in that it facilitated the spontaneous cr 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 turned out to be -advantagious in that it facilitated the spontaneous crystallization of +advantageous in that it facilitated the spontaneous crystallization of Ice-{\it i}. 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 changes when applying a longer cutoff. +\begin{figure} +\includegraphics[width=\linewidth]{cutoffChange.eps} +\caption{Free energy as a function of cutoff radius for (A) SSD/E, (B) +TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12 +\AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9 +\AA\. These crystals are unstable at 200 K and rapidly convert into a +liquid. The connecting lines are qualitative visual aid.} +\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 all the ice polymorphs show a substantial dependence on +cutoff radius. In general, there is a narrowing of the free energy +differences while moving to greater cutoff radius. Interestingly, by +increasing the cutoff radius, the free energy gap was narrowed enough +in the SSD/E model 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\. This narrowing trend is much more subtle in the +case of SSD/RF, indicating that the free energies calculated with a +reaction field present provide a more accurate picture of the free +energy landscape in the absence of potential truncation. +To further study the changes resulting to the inclusion of a +long-range interaction correction, 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 with and without PME using TINKER. The +free energies for the investigated polymorphs using the TIP3P and +SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P +are not fully supported in TINKER, so the results for these models +could not be estimated. The same trend pointed out through increase of +cutoff radius is observed in these PME results. Ice-{\it i} is the +preferred polymorph at ambient conditions for both the TIP3P and SPC/E +water models; however, there is a narrowing of the free energy +differences between the various solid forms. In the case of SPC/E this +narrowing is significant enough that it becomes less clear cut that +Ice-{\it i} is the most stable polymorph, and is possibly metastable +with respect to ice B and possibly ice $I_c$. However, these results +do not significantly alter 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. + +\begin{table*} +\begin{minipage}{\linewidth} +\renewcommand{\thefootnote}{\thempfootnote} +\begin{center} +\caption{The free energy of the studied ice polymorphs after applying +the energy difference attributed to the inclusion of the PME +long-range interaction correction. Units are kcal/mol.} +\begin{tabular}{ l c c c c } +\hline +\ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\ +\hline +\ \ TIP3P & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\ +\ \ SPC/E & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\ +\end{tabular} +\label{pmeShift} +\end{center} +\end{minipage} +\end{table*} + \section{Conclusions} +The free energy for proton ordered variants of hexagonal and cubic ice +$I$, ice B, and recently discovered Ice-{\it i} where 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 in the 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 +correction. 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 manner of +investigated simulation examples, 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 leads the authors to be hesitant in ascribing relevance outside +of computational models, hence the descriptive name presented. That +being 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. + \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