--- trunk/iceiPaper/iceiPaper.tex 2004/09/15 06:34:49 1458 +++ 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -47,13 +52,104 @@ 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} @@ -72,17 +168,16 @@ the two states: integrated in order to determine the free energy difference between the two states: \begin{equation} -\begin{center} \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda )}{\partial\lambda}\right\rangle_\lambda d\lambda, -\end{center} \end{equation} where $V$ is the interaction potential and $\lambda$ is the -transformation parameter. 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. +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 @@ -110,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 @@ -122,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} @@ -167,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\\ @@ -196,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 @@ -205,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 @@ -223,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 & \ \ - & \ \ -\\ @@ -253,7 +353,7 @@ 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 @@ -275,11 +375,15 @@ differences while moving to greater cutoff radius. Thi 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. This 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. +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 @@ -291,7 +395,7 @@ cutoff radius is observed in these results. Ice-{\it i 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 results. Ice-{\it i} is the +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 @@ -310,9 +414,9 @@ long-range interaction correction. Units are kcal/mol. 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 \\[-7mm] +\hline \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\ -\hline \\[-3mm] +\hline \ \ TIP3P & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\ \ \ SPC/E & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\ \end{tabular} @@ -335,9 +439,23 @@ computed free energy values, but Ice-{\it i} is still 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, but Ice-{\it i} is still observed to be a -relavent ice polymorph in simulation studies. +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