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