| 1 | chrisfen | 1453 | %\documentclass[prb,aps,twocolumn,tabularx]{revtex4} | 
| 2 | gezelter | 1463 | \documentclass[11pt]{article} | 
| 3 | chrisfen | 1453 | %\documentclass[11pt]{article} | 
| 4 | gezelter | 1463 | \usepackage{endfloat} | 
| 5 | chrisfen | 1453 | \usepackage{amsmath} | 
| 6 |  |  | \usepackage{epsf} | 
| 7 |  |  | \usepackage{berkeley} | 
| 8 | gezelter | 1463 | \usepackage{setspace} | 
| 9 |  |  | \usepackage{tabularx} | 
| 10 | chrisfen | 1453 | \usepackage{graphicx} | 
| 11 | gezelter | 1463 | \usepackage[ref]{overcite} | 
| 12 |  |  | \pagestyle{plain} | 
| 13 |  |  | \pagenumbering{arabic} | 
| 14 |  |  | \oddsidemargin 0.0cm \evensidemargin 0.0cm | 
| 15 |  |  | \topmargin -21pt \headsep 10pt | 
| 16 |  |  | \textheight 9.0in \textwidth 6.5in | 
| 17 |  |  | \brokenpenalty=10000 | 
| 18 |  |  | \renewcommand{\baselinestretch}{1.2} | 
| 19 |  |  | \renewcommand\citemid{\ } % no comma in optional reference note | 
| 20 | chrisfen | 1453 |  | 
| 21 |  |  | \begin{document} | 
| 22 |  |  |  | 
| 23 | gezelter | 1465 | \title{Ice-{\it i}: a simulation-predicted ice polymorph which is more | 
| 24 |  |  | stable than Ice $I_h$ for point-charge and point-dipole water models} | 
| 25 | chrisfen | 1453 |  | 
| 26 | gezelter | 1463 | \author{Christopher J. Fennell and J. Daniel Gezelter \\ | 
| 27 |  |  | Department of Chemistry and Biochemistry\\ University of Notre Dame\\ | 
| 28 | chrisfen | 1453 | Notre Dame, Indiana 46556} | 
| 29 |  |  |  | 
| 30 |  |  | \date{\today} | 
| 31 |  |  |  | 
| 32 | gezelter | 1463 | \maketitle | 
| 33 | chrisfen | 1453 | %\doublespacing | 
| 34 |  |  |  | 
| 35 |  |  | \begin{abstract} | 
| 36 | chrisfen | 1459 | The free energies of several ice polymorphs in the low pressure regime | 
| 37 | gezelter | 1463 | were calculated using thermodynamic integration.  These integrations | 
| 38 |  |  | were done for most of the common water models. Ice-{\it i}, a | 
| 39 |  |  | structure we recently observed to be stable in one of the single-point | 
| 40 |  |  | water models, was determined to be the stable crystalline state (at 1 | 
| 41 |  |  | atm) for {\it all} the water models investigated.  Phase diagrams were | 
| 42 |  |  | generated, and phase coexistence lines were determined for all of the | 
| 43 |  |  | known low-pressure ice structures under all of the common water | 
| 44 |  |  | models.  Additionally, potential truncation was shown to have an | 
| 45 |  |  | effect on the calculated free energies, and can result in altered free | 
| 46 |  |  | energy landscapes. | 
| 47 | chrisfen | 1453 | \end{abstract} | 
| 48 |  |  |  | 
| 49 |  |  | %\narrowtext | 
| 50 |  |  |  | 
| 51 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 52 |  |  | %                              BODY OF TEXT | 
| 53 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 54 |  |  |  | 
| 55 |  |  | \section{Introduction} | 
| 56 |  |  |  | 
| 57 | chrisfen | 1471 | Computer simulations are a valuable tool for studying the phase | 
| 58 |  |  | behavior of systems ranging from small or simple molecules to complex | 
| 59 | chrisfen | 1472 | biological species.\cite{Matsumoto02,Li96,Marrink01} Useful techniques | 
| 60 | chrisfen | 1471 | have been developed to investigate the thermodynamic properites of | 
| 61 |  |  | model substances, providing both qualitative and quantitative | 
| 62 |  |  | comparisons between simulations and | 
| 63 |  |  | experiment.\cite{Widom63,Frenkel84} Investigation of these properties | 
| 64 |  |  | leads to the development of new and more accurate models, leading to | 
| 65 |  |  | better understanding and depiction of physical processes and intricate | 
| 66 |  |  | molecular systems. | 
| 67 | chrisfen | 1459 |  | 
| 68 |  |  | Water has proven to be a challenging substance to depict in | 
| 69 | gezelter | 1463 | simulations, and a variety of models have been developed to describe | 
| 70 |  |  | its behavior under varying simulation | 
| 71 | chrisfen | 1471 | conditions.\cite{Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Liu96,Berendsen98,Mahoney00,Fennell04} | 
| 72 | gezelter | 1463 | These models have been used to investigate important physical | 
| 73 | chrisfen | 1471 | phenomena like phase transitions, molecule transport, and the | 
| 74 |  |  | hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the | 
| 75 |  |  | choice of models available, it is only natural to compare the models | 
| 76 |  |  | under interesting thermodynamic conditions in an attempt to clarify | 
| 77 |  |  | the limitations of each of the | 
| 78 |  |  | models.\cite{Jorgensen83,Jorgensen98b,Clancy94,Mahoney01} Two | 
| 79 |  |  | important property to quantify are the Gibbs and Helmholtz free | 
| 80 |  |  | energies, particularly for the solid forms of water.  Difficulty in | 
| 81 |  |  | these types of studies typically arises from the assortment of | 
| 82 |  |  | possible crystalline polymorphs that water adopts over a wide range of | 
| 83 |  |  | pressures and temperatures.  There are currently 13 recognized forms | 
| 84 |  |  | of ice, and it is a challenging task to investigate the entire free | 
| 85 |  |  | energy landscape.\cite{Sanz04} Ideally, research is focused on the | 
| 86 |  |  | phases having the lowest free energy at a given state point, because | 
| 87 |  |  | these phases will dictate the true transition temperatures and | 
| 88 |  |  | pressures for the model. | 
| 89 | chrisfen | 1459 |  | 
| 90 | gezelter | 1465 | In this paper, standard reference state methods were applied to known | 
| 91 | chrisfen | 1471 | crystalline water polymorphs in the low pressure regime.  This work is | 
| 92 | gezelter | 1465 | unique in the fact that one of the crystal lattices was arrived at | 
| 93 |  |  | through crystallization of a computationally efficient water model | 
| 94 |  |  | under constant pressure and temperature conditions. Crystallization | 
| 95 |  |  | events are interesting in and of | 
| 96 |  |  | themselves;\cite{Matsumoto02,Yamada02} however, the crystal structure | 
| 97 |  |  | obtained in this case is different from any previously observed ice | 
| 98 |  |  | polymorphs in experiment or simulation.\cite{Fennell04} We have named | 
| 99 |  |  | this structure Ice-{\it i} to indicate its origin in computational | 
| 100 | chrisfen | 1459 | simulation. The unit cell (Fig. \ref{iceiCell}A) consists of eight | 
| 101 |  |  | water molecules that stack in rows of interlocking water | 
| 102 |  |  | tetramers. Proton ordering can be accomplished by orienting two of the | 
| 103 | gezelter | 1465 | molecules so that both of their donated hydrogen bonds are internal to | 
| 104 | chrisfen | 1459 | their tetramer (Fig. \ref{protOrder}). As expected in an ice crystal | 
| 105 |  |  | constructed of water tetramers, the hydrogen bonds are not as linear | 
| 106 |  |  | as those observed in ice $I_h$, however the interlocking of these | 
| 107 |  |  | subunits appears to provide significant stabilization to the overall | 
| 108 |  |  | crystal. The arrangement of these tetramers results in surrounding | 
| 109 |  |  | open octagonal cavities that are typically greater than 6.3 \AA\ in | 
| 110 |  |  | diameter. This relatively open overall structure leads to crystals | 
| 111 |  |  | that are 0.07 g/cm$^3$ less dense on average than ice $I_h$. | 
| 112 | gezelter | 1463 |  | 
| 113 | chrisfen | 1460 | \begin{figure} | 
| 114 | gezelter | 1463 | \includegraphics[width=\linewidth]{unitCell.eps} | 
| 115 | gezelter | 1465 | \caption{Unit cells for (A) Ice-{\it i} and (B) Ice-$i^\prime$, the | 
| 116 | chrisfen | 1471 | elongated variant of Ice-{\it i}.  The spheres represent the | 
| 117 |  |  | center-of-mass locations of the water molecules.  The $a$ to $c$ | 
| 118 |  |  | ratios for Ice-{\it i} and Ice-{\it i}$^\prime$ are given by | 
| 119 |  |  | $a:2.1214c$ and $a:1.7850c$ respectively.} | 
| 120 | chrisfen | 1460 | \label{iceiCell} | 
| 121 |  |  | \end{figure} | 
| 122 | gezelter | 1463 |  | 
| 123 | chrisfen | 1460 | \begin{figure} | 
| 124 | gezelter | 1463 | \includegraphics[width=\linewidth]{orderedIcei.eps} | 
| 125 | chrisfen | 1460 | \caption{Image of a proton ordered crystal of Ice-{\it i} looking | 
| 126 |  |  | down the (001) crystal face. The rows of water tetramers surrounded by | 
| 127 |  |  | octagonal pores leads to a crystal structure that is significantly | 
| 128 |  |  | less dense than ice $I_h$.} | 
| 129 |  |  | \label{protOrder} | 
| 130 |  |  | \end{figure} | 
| 131 | chrisfen | 1459 |  | 
| 132 | gezelter | 1465 | Results from our previous study indicated that Ice-{\it i} is the | 
| 133 |  |  | minimum energy crystal structure for the single point water models we | 
| 134 |  |  | investigated (for discussions on these single point dipole models, see | 
| 135 | chrisfen | 1471 | our previous work and related | 
| 136 |  |  | articles).\cite{Fennell04,Liu96,Bratko85} Those results only | 
| 137 | gezelter | 1465 | considered energetic stabilization and neglected entropic | 
| 138 |  |  | contributions to the overall free energy. To address this issue, the | 
| 139 |  |  | absolute free energy of this crystal was calculated using | 
| 140 |  |  | thermodynamic integration and compared to the free energies of cubic | 
| 141 |  |  | and hexagonal ice $I$ (the experimental low density ice polymorphs) | 
| 142 |  |  | and ice B (a higher density, but very stable crystal structure | 
| 143 |  |  | observed by B\`{a}ez and Clancy in free energy studies of | 
| 144 |  |  | SPC/E).\cite{Baez95b} This work includes results for the water model | 
| 145 |  |  | from which Ice-{\it i} was crystallized (SSD/E) in addition to several | 
| 146 |  |  | common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction | 
| 147 |  |  | field parametrized single point dipole water model (SSD/RF). It should | 
| 148 |  |  | be noted that a second version of Ice-{\it i} (Ice-$i^\prime$) was used | 
| 149 |  |  | in calculations involving SPC/E, TIP4P, and TIP5P. The unit cell of | 
| 150 |  |  | this crystal (Fig. \ref{iceiCell}B) is similar to the Ice-{\it i} unit | 
| 151 |  |  | it is extended in the direction of the (001) face and compressed along | 
| 152 |  |  | the other two faces. | 
| 153 | chrisfen | 1459 |  | 
| 154 | chrisfen | 1453 | \section{Methods} | 
| 155 |  |  |  | 
| 156 | chrisfen | 1454 | Canonical ensemble (NVT) molecular dynamics calculations were | 
| 157 | gezelter | 1465 | performed using the OOPSE molecular mechanics package.\cite{Meineke05} | 
| 158 |  |  | All molecules were treated as rigid bodies, with orientational motion | 
| 159 |  |  | propagated using the symplectic DLM integration method. Details about | 
| 160 | chrisfen | 1471 | the implementation of this technique can be found in a recent | 
| 161 | gezelter | 1468 | publication.\cite{Dullweber1997} | 
| 162 | chrisfen | 1454 |  | 
| 163 | chrisfen | 1471 | Thermodynamic integration is an established technique for | 
| 164 |  |  | determination of free energies of condensed phases of | 
| 165 |  |  | materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This | 
| 166 |  |  | method, implemented in the same manner illustrated by B\`{a}ez and | 
| 167 |  |  | Clancy, was utilized to calculate the free energy of several ice | 
| 168 |  |  | crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and | 
| 169 |  |  | SSD/E water models.\cite{Baez95a} Liquid state free energies at 300 | 
| 170 |  |  | and 400 K for all of these water models were also determined using | 
| 171 |  |  | this same technique in order to determine melting points and generate | 
| 172 |  |  | phase diagrams. All simulations were carried out at densities | 
| 173 |  |  | resulting in a pressure of approximately 1 atm at their respective | 
| 174 |  |  | temperatures. | 
| 175 | chrisfen | 1454 |  | 
| 176 | chrisfen | 1458 | A single thermodynamic integration involves a sequence of simulations | 
| 177 |  |  | over which the system of interest is converted into a reference system | 
| 178 | gezelter | 1465 | for which the free energy is known analytically. This transformation | 
| 179 |  |  | path is then integrated in order to determine the free energy | 
| 180 |  |  | difference between the two states: | 
| 181 | chrisfen | 1458 | \begin{equation} | 
| 182 |  |  | \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda | 
| 183 |  |  | )}{\partial\lambda}\right\rangle_\lambda d\lambda, | 
| 184 |  |  | \end{equation} | 
| 185 |  |  | where $V$ is the interaction potential and $\lambda$ is the | 
| 186 | chrisfen | 1459 | transformation parameter that scales the overall | 
| 187 | chrisfen | 1471 | potential. Simulations are distributed strategically along this path | 
| 188 |  |  | in order to sufficiently sample the regions of greatest change in the | 
| 189 | chrisfen | 1459 | potential. Typical integrations in this study consisted of $\sim$25 | 
| 190 |  |  | simulations ranging from 300 ps (for the unaltered system) to 75 ps | 
| 191 |  |  | (near the reference state) in length. | 
| 192 | chrisfen | 1458 |  | 
| 193 | chrisfen | 1454 | For the thermodynamic integration of molecular crystals, the Einstein | 
| 194 | chrisfen | 1471 | crystal was chosen as the reference system. In an Einstein crystal, | 
| 195 |  |  | the molecules are restrained at their ideal lattice locations and | 
| 196 |  |  | orientations. Using harmonic restraints, as applied by B\`{a}ez and | 
| 197 |  |  | Clancy, the total potential for this reference crystal | 
| 198 |  |  | ($V_\mathrm{EC}$) is the sum of all the harmonic restraints, | 
| 199 |  |  | \begin{equation} | 
| 200 |  |  | V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} + | 
| 201 |  |  | \frac{K_\omega\omega^2}{2}, | 
| 202 |  |  | \end{equation} | 
| 203 |  |  | where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are | 
| 204 |  |  | the spring constants restraining translational motion and deflection | 
| 205 |  |  | of and rotation around the principle axis of the molecule | 
| 206 |  |  | respectively.  It is clear from Fig. \ref{waterSpring} that the values | 
| 207 |  |  | of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges from | 
| 208 |  |  | $-\pi$ to $\pi$.  The partition function for a molecular crystal | 
| 209 | gezelter | 1465 | restrained in this fashion can be evaluated analytically, and the | 
| 210 |  |  | Helmholtz Free Energy ({\it A}) is given by | 
| 211 | chrisfen | 1454 | \begin{eqnarray} | 
| 212 |  |  | A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left | 
| 213 |  |  | [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right | 
| 214 |  |  | )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right | 
| 215 |  |  | )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right | 
| 216 |  |  | )^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi | 
| 217 |  |  | K_\omega K_\theta)^{\frac{1}{2}}}\exp\left | 
| 218 |  |  | (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right | 
| 219 |  |  | )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ], | 
| 220 |  |  | \label{ecFreeEnergy} | 
| 221 |  |  | \end{eqnarray} | 
| 222 | chrisfen | 1471 | where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum | 
| 223 |  |  | potential energy of the ideal crystal.\cite{Baez95a} | 
| 224 | gezelter | 1463 |  | 
| 225 | chrisfen | 1456 | \begin{figure} | 
| 226 | gezelter | 1463 | \includegraphics[width=\linewidth]{rotSpring.eps} | 
| 227 | chrisfen | 1456 | \caption{Possible orientational motions for a restrained molecule. | 
| 228 |  |  | $\theta$ angles correspond to displacement from the body-frame {\it | 
| 229 |  |  | z}-axis, while $\omega$ angles correspond to rotation about the | 
| 230 |  |  | body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring | 
| 231 |  |  | constants for the harmonic springs restraining motion in the $\theta$ | 
| 232 |  |  | and $\omega$ directions.} | 
| 233 |  |  | \label{waterSpring} | 
| 234 |  |  | \end{figure} | 
| 235 | chrisfen | 1454 |  | 
| 236 | chrisfen | 1471 | In the case of molecular liquids, the ideal vapor is chosen as the | 
| 237 |  |  | target reference state.  There are several examples of liquid state | 
| 238 |  |  | free energy calculations of water models present in the | 
| 239 |  |  | literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods | 
| 240 |  |  | typically differ in regard to the path taken for switching off the | 
| 241 |  |  | interaction potential to convert the system to an ideal gas of water | 
| 242 |  |  | molecules.  In this study, we apply of one of the most convenient | 
| 243 |  |  | methods and integrate over the $\lambda^4$ path, where all interaction | 
| 244 |  |  | parameters are scaled equally by this transformation parameter.  This | 
| 245 |  |  | method has been shown to be reversible and provide results in | 
| 246 |  |  | excellent agreement with other established methods.\cite{Baez95b} | 
| 247 |  |  |  | 
| 248 | chrisfen | 1456 | Charge, dipole, and Lennard-Jones interactions were modified by a | 
| 249 | chrisfen | 1462 | cubic switching between 100\% and 85\% of the cutoff value (9 \AA | 
| 250 |  |  | ). By applying this function, these interactions are smoothly | 
| 251 | gezelter | 1465 | truncated, thereby avoiding the poor energy conservation which results | 
| 252 | chrisfen | 1462 | from harsher truncation schemes. The effect of a long-range correction | 
| 253 |  |  | was also investigated on select model systems in a variety of | 
| 254 |  |  | manners. For the SSD/RF model, a reaction field with a fixed | 
| 255 |  |  | dielectric constant of 80 was applied in all | 
| 256 |  |  | simulations.\cite{Onsager36} For a series of the least computationally | 
| 257 |  |  | expensive models (SSD/E, SSD/RF, and TIP3P), simulations were | 
| 258 |  |  | performed with longer cutoffs of 12 and 15 \AA\ to compare with the 9 | 
| 259 |  |  | \AA\ cutoff results. Finally, results from the use of an Ewald | 
| 260 |  |  | summation were estimated for TIP3P and SPC/E by performing | 
| 261 | chrisfen | 1456 | calculations with Particle-Mesh Ewald (PME) in the TINKER molecular | 
| 262 | gezelter | 1465 | mechanics software package.\cite{Tinker} The calculated energy | 
| 263 | chrisfen | 1462 | difference in the presence and absence of PME was applied to the | 
| 264 | gezelter | 1465 | previous results in order to predict changes to the free energy | 
| 265 | chrisfen | 1462 | landscape. | 
| 266 | chrisfen | 1454 |  | 
| 267 | chrisfen | 1456 | \section{Results and discussion} | 
| 268 | chrisfen | 1454 |  | 
| 269 | chrisfen | 1456 | The free energy of proton ordered Ice-{\it i} was calculated and | 
| 270 |  |  | compared with the free energies of proton ordered variants of the | 
| 271 |  |  | experimentally observed low-density ice polymorphs, $I_h$ and $I_c$, | 
| 272 |  |  | as well as the higher density ice B, observed by B\`{a}ez and Clancy | 
| 273 |  |  | and thought to be the minimum free energy structure for the SPC/E | 
| 274 |  |  | model at ambient conditions (Table \ref{freeEnergy}).\cite{Baez95b} | 
| 275 | gezelter | 1465 | Ice XI, the experimentally-observed proton-ordered variant of ice | 
| 276 |  |  | $I_h$, was investigated initially, but was found to be not as stable | 
| 277 |  |  | as proton disordered or antiferroelectric variants of ice $I_h$. The | 
| 278 |  |  | proton ordered variant of ice $I_h$ used here is a simple | 
| 279 |  |  | antiferroelectric version that has an 8 molecule unit | 
| 280 |  |  | cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules | 
| 281 |  |  | for ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for | 
| 282 |  |  | ice $I_c$, or 1024 molecules for Ice-{\it i}. The larger crystal sizes | 
| 283 | chrisfen | 1456 | were necessary for simulations involving larger cutoff values. | 
| 284 | chrisfen | 1454 |  | 
| 285 | chrisfen | 1456 | \begin{table*} | 
| 286 |  |  | \begin{minipage}{\linewidth} | 
| 287 |  |  | \renewcommand{\thefootnote}{\thempfootnote} | 
| 288 |  |  | \begin{center} | 
| 289 |  |  | \caption{Calculated free energies for several ice polymorphs with a | 
| 290 |  |  | variety of common water models. All calculations used a cutoff radius | 
| 291 |  |  | of 9 \AA\ and were performed at 200 K and $\sim$1 atm. Units are | 
| 292 | chrisfen | 1466 | kcal/mol. Calculated error of the final digits is in parentheses. *Ice | 
| 293 |  |  | $I_c$ rapidly converts to a liquid at 200 K with the SSD/RF model.} | 
| 294 | chrisfen | 1456 | \begin{tabular}{ l  c  c  c  c } | 
| 295 | gezelter | 1463 | \hline | 
| 296 | chrisfen | 1466 | Water Model & $I_h$ & $I_c$ & B & Ice-{\it i}\\ | 
| 297 | gezelter | 1463 | \hline | 
| 298 | chrisfen | 1466 | TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\ | 
| 299 |  |  | TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\ | 
| 300 |  |  | TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\ | 
| 301 |  |  | SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\ | 
| 302 |  |  | SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\ | 
| 303 |  |  | SSD/RF & -11.51(4) & NA* & -12.08(5) & -12.29(4)\\ | 
| 304 | chrisfen | 1456 | \end{tabular} | 
| 305 |  |  | \label{freeEnergy} | 
| 306 |  |  | \end{center} | 
| 307 |  |  | \end{minipage} | 
| 308 |  |  | \end{table*} | 
| 309 | chrisfen | 1453 |  | 
| 310 | chrisfen | 1456 | The free energy values computed for the studied polymorphs indicate | 
| 311 |  |  | that Ice-{\it i} is the most stable state for all of the common water | 
| 312 |  |  | models studied. With the free energy at these state points, the | 
| 313 | gezelter | 1465 | Gibbs-Helmholtz equation was used to project to other state points and | 
| 314 |  |  | to build phase diagrams.  Figures | 
| 315 | chrisfen | 1456 | \ref{tp3phasedia} and \ref{ssdrfphasedia} are example diagrams built | 
| 316 |  |  | from the free energy results. All other models have similar structure, | 
| 317 | gezelter | 1465 | although the crossing points between the phases exist at slightly | 
| 318 |  |  | different temperatures and pressures. It is interesting to note that | 
| 319 |  |  | ice $I$ does not exist in either cubic or hexagonal form in any of the | 
| 320 |  |  | phase diagrams for any of the models. For purposes of this study, ice | 
| 321 |  |  | B is representative of the dense ice polymorphs. A recent study by | 
| 322 |  |  | Sanz {\it et al.} goes into detail on the phase diagrams for SPC/E and | 
| 323 |  |  | TIP4P in the high pressure regime.\cite{Sanz04} | 
| 324 | gezelter | 1463 |  | 
| 325 | chrisfen | 1456 | \begin{figure} | 
| 326 |  |  | \includegraphics[width=\linewidth]{tp3PhaseDia.eps} | 
| 327 |  |  | \caption{Phase diagram for the TIP3P water model in the low pressure | 
| 328 |  |  | regime. The displayed $T_m$ and $T_b$ values are good predictions of | 
| 329 |  |  | the experimental values; however, the solid phases shown are not the | 
| 330 |  |  | experimentally observed forms. Both cubic and hexagonal ice $I$ are | 
| 331 |  |  | higher in energy and don't appear in the phase diagram.} | 
| 332 |  |  | \label{tp3phasedia} | 
| 333 |  |  | \end{figure} | 
| 334 | gezelter | 1463 |  | 
| 335 | chrisfen | 1456 | \begin{figure} | 
| 336 |  |  | \includegraphics[width=\linewidth]{ssdrfPhaseDia.eps} | 
| 337 |  |  | \caption{Phase diagram for the SSD/RF water model in the low pressure | 
| 338 |  |  | regime. Calculations producing these results were done under an | 
| 339 |  |  | applied reaction field. It is interesting to note that this | 
| 340 |  |  | computationally efficient model (over 3 times more efficient than | 
| 341 |  |  | TIP3P) exhibits phase behavior similar to the less computationally | 
| 342 |  |  | conservative charge based models.} | 
| 343 |  |  | \label{ssdrfphasedia} | 
| 344 |  |  | \end{figure} | 
| 345 |  |  |  | 
| 346 |  |  | \begin{table*} | 
| 347 |  |  | \begin{minipage}{\linewidth} | 
| 348 |  |  | \renewcommand{\thefootnote}{\thempfootnote} | 
| 349 |  |  | \begin{center} | 
| 350 |  |  | \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) | 
| 351 | chrisfen | 1466 | temperatures at 1 atm for several common water models compared with | 
| 352 |  |  | experiment. The $T_m$ and $T_s$ values from simulation correspond to a | 
| 353 |  |  | transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the | 
| 354 |  |  | liquid or gas state.} | 
| 355 | chrisfen | 1456 | \begin{tabular}{ l  c  c  c  c  c  c  c } | 
| 356 | gezelter | 1463 | \hline | 
| 357 | chrisfen | 1466 | Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\ | 
| 358 | gezelter | 1463 | \hline | 
| 359 | chrisfen | 1466 | $T_m$ (K)  & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\ | 
| 360 |  |  | $T_b$ (K)  & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\ | 
| 361 |  |  | $T_s$ (K)  & - & - & - & - & 355(3) & - & -\\ | 
| 362 | chrisfen | 1456 | \end{tabular} | 
| 363 |  |  | \label{meltandboil} | 
| 364 |  |  | \end{center} | 
| 365 |  |  | \end{minipage} | 
| 366 |  |  | \end{table*} | 
| 367 |  |  |  | 
| 368 |  |  | Table \ref{meltandboil} lists the melting and boiling temperatures | 
| 369 |  |  | calculated from this work. Surprisingly, most of these models have | 
| 370 |  |  | melting points that compare quite favorably with experiment. The | 
| 371 |  |  | unfortunate aspect of this result is that this phase change occurs | 
| 372 |  |  | between Ice-{\it i} and the liquid state rather than ice $I_h$ and the | 
| 373 |  |  | liquid state. These results are actually not contrary to previous | 
| 374 |  |  | studies in the literature. Earlier free energy studies of ice $I$ | 
| 375 |  |  | using TIP4P predict a $T_m$ ranging from 214 to 238 K (differences | 
| 376 |  |  | being attributed to choice of interaction truncation and different | 
| 377 | chrisfen | 1466 | ordered and disordered molecular | 
| 378 |  |  | arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and | 
| 379 |  |  | Ice-{\it i} were omitted, a $T_m$ value around 210 K would be | 
| 380 | chrisfen | 1456 | predicted from this work. However, the $T_m$ from Ice-{\it i} is | 
| 381 |  |  | calculated at 265 K, significantly higher in temperature than the | 
| 382 |  |  | previous studies. Also of interest in these results is that SSD/E does | 
| 383 |  |  | not exhibit a melting point at 1 atm, but it shows a sublimation point | 
| 384 |  |  | at 355 K. This is due to the significant stability of Ice-{\it i} over | 
| 385 |  |  | all other polymorphs for this particular model under these | 
| 386 |  |  | conditions. While troubling, this behavior turned out to be | 
| 387 | chrisfen | 1459 | advantageous in that it facilitated the spontaneous crystallization of | 
| 388 | chrisfen | 1456 | Ice-{\it i}. These observations provide a warning that simulations of | 
| 389 |  |  | SSD/E as a ``liquid'' near 300 K are actually metastable and run the | 
| 390 |  |  | risk of spontaneous crystallization. However, this risk changes when | 
| 391 |  |  | applying a longer cutoff. | 
| 392 |  |  |  | 
| 393 | chrisfen | 1458 | \begin{figure} | 
| 394 |  |  | \includegraphics[width=\linewidth]{cutoffChange.eps} | 
| 395 |  |  | \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B) | 
| 396 |  |  | TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12 | 
| 397 |  |  | \AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9 | 
| 398 | chrisfen | 1466 | \AA . These crystals are unstable at 200 K and rapidly convert into | 
| 399 |  |  | liquids. The connecting lines are qualitative visual aid.} | 
| 400 | chrisfen | 1458 | \label{incCutoff} | 
| 401 |  |  | \end{figure} | 
| 402 |  |  |  | 
| 403 | chrisfen | 1457 | Increasing the cutoff radius in simulations of the more | 
| 404 |  |  | computationally efficient water models was done in order to evaluate | 
| 405 |  |  | the trend in free energy values when moving to systems that do not | 
| 406 |  |  | involve potential truncation. As seen in Fig. \ref{incCutoff}, the | 
| 407 |  |  | free energy of all the ice polymorphs show a substantial dependence on | 
| 408 |  |  | cutoff radius. In general, there is a narrowing of the free energy | 
| 409 | chrisfen | 1459 | differences while moving to greater cutoff radius. Interestingly, by | 
| 410 |  |  | increasing the cutoff radius, the free energy gap was narrowed enough | 
| 411 |  |  | in the SSD/E model that the liquid state is preferred under standard | 
| 412 |  |  | simulation conditions (298 K and 1 atm). Thus, it is recommended that | 
| 413 |  |  | simulations using this model choose interaction truncation radii | 
| 414 | gezelter | 1469 | greater than 9 \AA\ . This narrowing trend is much more subtle in the | 
| 415 | chrisfen | 1459 | case of SSD/RF, indicating that the free energies calculated with a | 
| 416 |  |  | reaction field present provide a more accurate picture of the free | 
| 417 |  |  | energy landscape in the absence of potential truncation. | 
| 418 | chrisfen | 1456 |  | 
| 419 | chrisfen | 1457 | To further study the changes resulting to the inclusion of a | 
| 420 |  |  | long-range interaction correction, the effect of an Ewald summation | 
| 421 |  |  | was estimated by applying the potential energy difference do to its | 
| 422 |  |  | inclusion in systems in the presence and absence of the | 
| 423 |  |  | correction. This was accomplished by calculation of the potential | 
| 424 |  |  | energy of identical crystals with and without PME using TINKER. The | 
| 425 |  |  | free energies for the investigated polymorphs using the TIP3P and | 
| 426 | chrisfen | 1471 | SPC/E water models are shown in Table \ref{pmeShift}. The same trend | 
| 427 |  |  | pointed out through increase of cutoff radius is observed in these PME | 
| 428 |  |  | results. Ice-{\it i} is the preferred polymorph at ambient conditions | 
| 429 |  |  | for both the TIP3P and SPC/E water models; however, the narrowing of | 
| 430 |  |  | the free energy differences between the various solid forms is | 
| 431 |  |  | significant enough that it becomes less clear that it is the most | 
| 432 |  |  | stable polymorph.  The free energies of Ice-{\it i} and ice B overlap | 
| 433 |  |  | within error, with ice $I_c$ just outside, indicating that Ice-{\it i} | 
| 434 |  |  | might be metastable with respect to ice B and possibly ice $I_c$ in | 
| 435 |  |  | the SPC/E water model. However, these results do not significantly | 
| 436 |  |  | alter the finding that the Ice-{\it i} polymorph is a stable crystal | 
| 437 |  |  | structure that should be considered when studying the phase behavior | 
| 438 |  |  | of water models. | 
| 439 | chrisfen | 1456 |  | 
| 440 | chrisfen | 1457 | \begin{table*} | 
| 441 |  |  | \begin{minipage}{\linewidth} | 
| 442 |  |  | \renewcommand{\thefootnote}{\thempfootnote} | 
| 443 |  |  | \begin{center} | 
| 444 | chrisfen | 1458 | \caption{The free energy of the studied ice polymorphs after applying | 
| 445 |  |  | the energy difference attributed to the inclusion of the PME | 
| 446 |  |  | long-range interaction correction. Units are kcal/mol.} | 
| 447 | chrisfen | 1457 | \begin{tabular}{ l  c  c  c  c } | 
| 448 | gezelter | 1463 | \hline | 
| 449 | chrisfen | 1457 | \ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\ | 
| 450 | gezelter | 1463 | \hline | 
| 451 | chrisfen | 1466 | TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\ | 
| 452 |  |  | SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\ | 
| 453 | chrisfen | 1457 | \end{tabular} | 
| 454 |  |  | \label{pmeShift} | 
| 455 |  |  | \end{center} | 
| 456 |  |  | \end{minipage} | 
| 457 |  |  | \end{table*} | 
| 458 |  |  |  | 
| 459 | chrisfen | 1453 | \section{Conclusions} | 
| 460 |  |  |  | 
| 461 | chrisfen | 1458 | The free energy for proton ordered variants of hexagonal and cubic ice | 
| 462 | gezelter | 1465 | $I$, ice B, and recently discovered Ice-{\it i} were calculated under | 
| 463 | chrisfen | 1458 | standard conditions for several common water models via thermodynamic | 
| 464 |  |  | integration. All the water models studied show Ice-{\it i} to be the | 
| 465 |  |  | minimum free energy crystal structure in the with a 9 \AA\ switching | 
| 466 |  |  | function cutoff. Calculated melting and boiling points show | 
| 467 |  |  | surprisingly good agreement with the experimental values; however, the | 
| 468 |  |  | solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of | 
| 469 |  |  | interaction truncation was investigated through variation of the | 
| 470 |  |  | cutoff radius, use of a reaction field parameterized model, and | 
| 471 | gezelter | 1465 | estimation of the results in the presence of the Ewald | 
| 472 |  |  | summation. Interaction truncation has a significant effect on the | 
| 473 | chrisfen | 1459 | computed free energy values, and may significantly alter the free | 
| 474 |  |  | energy landscape for the more complex multipoint water models. Despite | 
| 475 |  |  | these effects, these results show Ice-{\it i} to be an important ice | 
| 476 |  |  | polymorph that should be considered in simulation studies. | 
| 477 | chrisfen | 1458 |  | 
| 478 | chrisfen | 1459 | Due to this relative stability of Ice-{\it i} in all manner of | 
| 479 |  |  | investigated simulation examples, the question arises as to possible | 
| 480 | gezelter | 1465 | experimental observation of this polymorph.  The rather extensive past | 
| 481 | chrisfen | 1459 | and current experimental investigation of water in the low pressure | 
| 482 | gezelter | 1465 | regime makes us hesitant to ascribe any relevance of this work outside | 
| 483 |  |  | of the simulation community.  It is for this reason that we chose a | 
| 484 |  |  | name for this polymorph which involves an imaginary quantity.  That | 
| 485 |  |  | said, there are certain experimental conditions that would provide the | 
| 486 |  |  | most ideal situation for possible observation. These include the | 
| 487 |  |  | negative pressure or stretched solid regime, small clusters in vacuum | 
| 488 |  |  | deposition environments, and in clathrate structures involving small | 
| 489 | gezelter | 1469 | non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain | 
| 490 |  |  | our predictions for both the pair distribution function ($g_{OO}(r)$) | 
| 491 |  |  | and the structure factor ($S(\vec{q})$ for ice $I_c$ and for ice-{\it | 
| 492 | chrisfen | 1470 | i} at a temperature of 77K.  In a quick comparison of the predicted | 
| 493 |  |  | S(q) for Ice-{\it i} and experimental studies of amorphous solid | 
| 494 |  |  | water, it is possible that some of the ``spurious'' peaks that could | 
| 495 |  |  | not be assigned in HDA could correspond to peaks labeled in this | 
| 496 |  |  | S(q).\cite{Bizid87} It should be noted that there is typically poor | 
| 497 |  |  | agreement on crystal densities between simulation and experiment, so | 
| 498 |  |  | such peak comparisons should be made with caution.  We will leave it | 
| 499 |  |  | to our experimental colleagues to determine whether this ice polymorph | 
| 500 |  |  | is named appropriately or if it should be promoted to Ice-0. | 
| 501 | chrisfen | 1459 |  | 
| 502 | chrisfen | 1467 | \begin{figure} | 
| 503 |  |  | \includegraphics[width=\linewidth]{iceGofr.eps} | 
| 504 | chrisfen | 1470 | \caption{Radial distribution functions of Ice-{\it i} and ice $I_c$ | 
| 505 |  |  | calculated from from simulations of the SSD/RF water model at 77 K.} | 
| 506 | chrisfen | 1467 | \label{fig:gofr} | 
| 507 |  |  | \end{figure} | 
| 508 |  |  |  | 
| 509 | gezelter | 1469 | \begin{figure} | 
| 510 |  |  | \includegraphics[width=\linewidth]{sofq.eps} | 
| 511 |  |  | \caption{Predicted structure factors for Ice-{\it i} and ice $I_c$ at | 
| 512 |  |  | 77 K.  The raw structure factors have been convoluted with a gaussian | 
| 513 | chrisfen | 1470 | instrument function (0.075 \AA$^{-1}$ width) to compensate for the | 
| 514 |  |  | trunction effects in our finite size simulations. The labeled peaks | 
| 515 |  |  | compared favorably with ``spurious'' peaks observed in experimental | 
| 516 |  |  | studies of amorphous solid water.\cite{Bizid87}} | 
| 517 | gezelter | 1469 | \label{fig:sofq} | 
| 518 |  |  | \end{figure} | 
| 519 |  |  |  | 
| 520 | chrisfen | 1453 | \section{Acknowledgments} | 
| 521 |  |  | Support for this project was provided by the National Science | 
| 522 |  |  | Foundation under grant CHE-0134881. Computation time was provided by | 
| 523 | chrisfen | 1458 | the Notre Dame High Performance Computing Cluster and the Notre Dame | 
| 524 |  |  | Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647). | 
| 525 | chrisfen | 1453 |  | 
| 526 |  |  | \newpage | 
| 527 |  |  |  | 
| 528 |  |  | \bibliographystyle{jcp} | 
| 529 |  |  | \bibliography{iceiPaper} | 
| 530 |  |  |  | 
| 531 |  |  |  | 
| 532 |  |  | \end{document} |