--- trunk/fennellDissertation/iceChapter.tex 2006/08/27 19:34:53 2978 +++ trunk/fennellDissertation/iceChapter.tex 2006/08/30 22:14:37 2986 @@ -1,18 +1,19 @@ \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS} -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 +As discussed in the previous chapter, water has proven to be a +challenging substance to depict in simulations, and a variety of +models have been developed to describe its behavior under varying +simulation conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04} These models have been used to investigate important physical phenomena like phase transitions and the hydrophobic effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models -available, it is only natural to compare the models under interesting +available, it is only natural to compare them under interesting thermodynamic conditions in an attempt to clarify the limitations of -each of the models.\cite{Jorgensen83,Jorgensen98,Baez94,Mahoney01} Two -important property to quantify are the Gibbs and Helmholtz free -energies, particularly for the solid forms of water, as these predict -the thermodynamic stability of the various phases. Water has a +each.\cite{Jorgensen83,Jorgensen98b,Baez94,Mahoney01} Two important +property to quantify are the Gibbs and Helmholtz free energies, +particularly for the solid forms of water, as these predict the +thermodynamic stability of the various phases. Water has a particularly rich phase diagram and takes on a number of different and stable crystalline structures as the temperature and pressure are varied. This complexity makes it a challenging task to investigate the @@ -23,31 +24,36 @@ computatuionally. In this chapter, standard reference The high-pressure phases of water (ice II-ice X as well as ice XII) have been studied extensively both experimentally and -computatuionally. In this chapter, standard reference state methods +computationally. In this chapter, standard reference state methods were applied in the {\it low} pressure regime to evaluate the free energies for a few known crystalline water polymorphs that might be stable at these pressures. 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 is different from any previously -observed ice polymorphs in experiment or simulation.\cite{Fennell04} -We have named this structure Ice-{\it i} to indicate its origin in -computational simulation. The unit cell of Ice-$i$ and an axially -elongated variant named Ice-$i^\prime$ both consist of eight water -molecules that stack in rows of interlocking water tetramers as -illustrated in figure \ref{fig:unitCell}A,B. These tetramers form a -crystal structure similar in appearance to a recent two-dimensional -surface tessellation simulated on silica.\cite{Yang04} As expected in -an ice crystal constructed of water tetramers, the hydrogen bonds are -not as linear as those observed in ice I$_\textrm{h}$; however, the -interlocking of these subunits appears to provide significant -stabilization to the overall crystal. The arrangement of these -tetramers results in open octagonal cavities that are typically -greater than 6.3\AA\ in diameter (see figure +temperature conditions. + +While performing a series of melting simulations on an early iteration +of SSD/E, we observed several recrystallization events at a constant +pressure of 1 atm. After melting from ice I$_\textrm{h}$ at 235~K, two +of five systems recrystallized near 245~K. Crystallization events are +interesting in and of themselves;\cite{Matsumoto02,Yamada02} however, +the crystal structure extracted from these systems is different from +any previously observed ice polymorphs in experiment or +simulation.\cite{Fennell04} We have named this structure Ice-{\it i} +to indicate its origin in computational simulation. The unit cell of +Ice-$i$ and an axially elongated variant named Ice-$i^\prime$ both +consist of eight water molecules that stack in rows of interlocking +water tetramers as illustrated in figure \ref{fig:iceiCell}A,B. These +tetramers form a crystal structure similar in appearance to a recent +two-dimensional surface tessellation simulated on silica.\cite{Yang04} +As expected in an ice crystal constructed of water tetramers, the +hydrogen bonds are not as linear as those observed in ice +I$_\textrm{h}$; however, the interlocking of these subunits appears to +provide significant stabilization to the overall crystal. The +arrangement of these tetramers results in open octagonal cavities that +are typically greater than 6.3~\AA\ in diameter (see figure \ref{fig:protOrder}). This open structure leads to crystals that are -typically 0.07 g/cm$^3$ less dense than ice I$_\textrm{h}$. +typically 0.07~g/cm$^3$ less dense than ice I$_\textrm{h}$. \begin{figure} \includegraphics[width=\linewidth]{./figures/unitCell.pdf} @@ -68,18 +74,18 @@ Results from our previous study indicated that Ice-{\i \label{fig:protOrder} \end{figure} -Results from our previous study indicated that Ice-{\it i} is the +Results from our initial studies indicated that Ice-{\it i} is the minimum energy crystal structure for the single point water models investigated (for discussions on these single point dipole models, see the previous work and related -articles\cite{Fennell04,Liu96,Bratko85}). Our earlier results only +articles\cite{Fennell04,Liu96,Bratko85}). These earlier results only considered energetic stabilization and neglected entropic contributions to the overall free energy. To address this issue, we have calculated the absolute free energy of this crystal using thermodynamic integration and compared to the free energies of ice I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common 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 +structure observed by B\'{a}ez and Clancy in free energy studies of SPC/E).\cite{Baez95b} This work includes results for the water model from which Ice-{\it i} was crystallized (SSD/E) in addition to several common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction @@ -93,68 +99,103 @@ to diagonal sites in the tetramers. to split the peak in the radial distribution function which corresponds to diagonal sites in the tetramers. -\section{Methods} +\section{Methods and Thermodynamic Integration} Canonical ensemble ({\it NVT}) molecular dynamics calculations were performed using the OOPSE molecular mechanics package.\cite{Meineke05} The densities chosen for the simulations were taken from isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and -200K. Each model (and each crystal structure) was allowed to relax for -300ps in the {\it NPT} ensemble before averaging the density to obtain +200~K. Each model (and each crystal structure) was allowed to relax for +300~ps in the {\it NPT} ensemble before averaging the density to obtain the volumes for the {\it NVT} simulations.All molecules were treated as rigid bodies, with orientational motion propagated using the symplectic DLM integration method described in section -\ref{sec:IntroIntegration}. +\ref{sec:IntroIntegrate}. + We used thermodynamic integration to calculate the Helmholtz free energies ({\it A}) of the listed water models at various state points. Thermodynamic integration is an established technique that has been used extensively in the calculation of free energies for condensed phases of -materials.\cite{Frenkel84,Hermans88,Meijer90,Baez95,Vlot99} This -method uses a sequence of simulations during which the system of +materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This +method uses a sequence of simulations over which the system of interest is converted into a reference system for which the free -energy is known analytically ($A_0$). This transformation path is then -integrated in order to determine the free energy difference between -the two states: +energy is known analytically ($A_0$). The difference in potential +energy between the reference system and the system of interest +($\Delta V$) is then integrated in order to determine the free energy +difference between the two states: \begin{equation} -\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda -)}{\partial\lambda}\right\rangle_\lambda d\lambda, + A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda. \end{equation} -where $V$ is the interaction potential and $\lambda$ is the -transformation parameter that scales the overall -potential. Simulations are distributed unevenly along this path in -order to sufficiently sample the regions of greatest change in the -potential. Typical integrations in this study consisted of $\sim$25 -simulations ranging from 300ps (for the unaltered system) to 75ps -(near the reference state) in length. +Here, $\lambda$ is the parameter that governs the transformation +between the reference system and the system of interest. For +crystalline phases, an harmonically-restrained (Einstein) crystal is +chosen as the reference state, while for liquid phases, the ideal gas +is taken as the reference state. Figure \ref{fig:integrationPath} +shows an example integration path for converting a crystalline system +to the Einstein crystal reference state. +\begin{figure} +\includegraphics[width=\linewidth]{./figures/integrationPath.pdf} +\caption{An example integration path to convert an unrestrained +crystal ($\lambda = 1$) to the Einstein crystal reference state +($\lambda = 0$). Note the increase in samples at either end of the +path to improve the smoothness of the curve. For reversible processes, +conversion of the Einstein crystal back to the system of interest will +give an identical plot, thereby integrating to the same result.} +\label{fig:integrationPath} +\end{figure} -For the thermodynamic integration of molecular crystals, the Einstein -crystal was chosen as the reference state. In an Einstein crystal, the -molecules are harmonically restrained at their ideal lattice locations -and orientations. The partition function for a molecular crystal -restrained in this fashion can be evaluated analytically, and the -Helmholtz Free Energy ({\it A}) is given by -\begin{eqnarray} -A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left -[\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right -)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right -)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right -)^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi -K_\omega K_\theta)^{\frac{1}{2}}}\exp\left -(-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right -)^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ], +In an Einstein crystal, the molecules are restrained at their ideal +lattice locations and orientations. Using harmonic restraints, as +applied by B\'{a}ez and Clancy, the total potential for this reference +crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints, +\begin{equation} +V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} + +\frac{K_\omega\omega^2}{2}, +\end{equation} +where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are +the spring constants restraining translational motion and deflection +of and rotation around the principle axis of the molecule +respectively. These spring constants are typically calculated from +the mean-square displacements of water molecules in an unrestrained +ice crystal at 200~K. For these studies, $K_\mathrm{v} = +4.29$~kcal~mol$^{-1}$~\AA$^{-2}$, $K_\theta\ = +13.88$~kcal~mol$^{-1}$~rad$^{-2}$, and $K_\omega\ = +17.75$~kcal~mol$^{-1}$~rad$^{-2}$. It is clear from +Fig. \ref{fig:waterSpring} that the values of $\theta$ range from $0$ +to $\pi$, while $\omega$ ranges from $-\pi$ to $\pi$. The partition +function for a molecular crystal restrained in this fashion can be +evaluated analytically, and the Helmholtz Free Energy ({\it A}) is +given by +\begin{equation} +\begin{split} +A = E_m &- kT\ln\left(\frac{kT}{h\nu}\right)^3 \\ + &- kT\ln\left[\pi^\frac{1}{2}\left( + \frac{8\pi^2I_\mathrm{A}kT}{h^2}\right)^\frac{1}{2} + \left(\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right)^\frac{1}{2} + \left(\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right)^\frac{1}{2} + \right] \\ + &- kT\ln\left[\frac{kT}{2(\pi K_\omega K_\theta)^{\frac{1}{2}}} + \exp\left(-\frac{kT}{2K_\theta}\right) + \int_0^{\left(\frac{kT}{2K_\theta}\right)^\frac{1}{2}} + \exp(t^2)\mathrm{d}t\right], +\end{split} \label{eq:ecFreeEnergy} -\end{eqnarray} -where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation -\ref{eq:ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and -$K_\mathrm{\omega}$ are the spring constants restraining translational -motion and deflection of and rotation around the principle axis of the -molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the -minimum potential energy of the ideal crystal. In the case of -molecular liquids, the ideal vapor is chosen as the target reference -state. - +\end{equation} +where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum +potential energy of the ideal crystal.\cite{Baez95a} The choice of an +Einstein crystal reference state is somewhat arbitrary. Any ideal +system for which the partition function is known exactly could be used +as a reference point as long as the system does not undergo a phase +transition during the integration path between the real and ideal +systems. Nada and van der Eerden have shown that the use of different +force constants in the Einstein crystal does not affect the total +free energy, and Gao {\it et al.} have shown that free energies +computed with the Debye crystal reference state differ from the +Einstein crystal by only a few tenths of a +kJ~mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can +lead to some uncertainty in the computed melting point of the solids. \begin{figure} \centering \includegraphics[width=3.5in]{./figures/rotSpring.pdf} @@ -164,87 +205,90 @@ and $\omega$ directions.} body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring constants for the harmonic springs restraining motion in the $\theta$ and $\omega$ directions.} -\label{waterSpring} +\label{fig:waterSpring} \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 the poor energy conservation which results from -harsher truncation schemes. The effect of a long-range correction was -also investigated on select model systems in a variety of manners. For -the SSD/RF model, a reaction field with a fixed dielectric constant of -80 was applied in all simulations.\cite{Onsager36} For a series of the -least computationally expensive models (SSD/E, SSD/RF, and TIP3P), -simulations were performed with longer cutoffs of 12 and 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.\cite{Tinker} The calculated energy -difference in the presence and absence of PME was applied to the -previous results in order to predict changes to the free energy -landscape. +In the case of molecular liquids, the ideal vapor is chosen as the +target reference state. There are several examples of liquid state +free energy calculations of water models present in the +literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods +typically differ in regard to the path taken for switching off the +interaction potential to convert the system to an ideal gas of water +molecules. In this study, we applied one of the most convenient +methods and integrated over the $\lambda^4$ path, where all +interaction parameters are scaled equally by this transformation +parameter. This method has been shown to be reversible and provide +results in excellent agreement with other established +methods.\cite{Baez95b} -\section{Results and discussion} +Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and +Lennard-Jones interactions were gradually reduced by a cubic switching +function. By applying this function, these interactions are smoothly +truncated, thereby avoiding the poor energy conservation which results +from harsher truncation schemes. The effect of a long-range +correction was also investigated on select model systems in a variety +of manners. For the SSD/RF model, a reaction field with a fixed +dielectric constant of 80 was applied in all +simulations.\cite{Onsager36} For a series of the least computationally +expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were +performed with longer cutoffs of 10.5, 12, 13.5, and 15~\AA\ to +compare with the 9~\AA\ cutoff results. Finally, the effects of using +the Ewald summation were estimated for TIP3P and SPC/E by performing +single configuration Particle-Mesh Ewald (PME) calculations for each +of the ice polymorphs.\cite{Ponder87} The calculated energy difference +in the presence and absence of PME was applied to the previous results +in order to predict changes to the free energy landscape. -The free energy of proton ordered Ice-{\it i} was calculated and -compared with the free energies of proton ordered variants of the -experimentally observed low-density ice polymorphs, I$_\textrm{h}$ and -I$_\textrm{c}$, as well as the higher density ice B, observed by -B\`{a}ez and Clancy and thought to be the minimum free energy -structure for the SPC/E model at ambient conditions.\cite{Baez95b} Ice -XI, the experimentally-observed proton-ordered variant of ice -I$_\textrm{h}$, was investigated initially, but was found to be not as -stable as proton disordered or antiferroelectric variants of ice -I$_\textrm{h}$. The proton ordered variant of ice I$_\textrm{h}$ used -here is a simple antiferroelectric version that has an 8 molecule unit -cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules -for ice B, 1024 or 1280 molecules for ice I$_\textrm{h}$, 1000 -molecules for ice I$_\textrm{c}$, or 1024 molecules for Ice-{\it -i}. The larger crystal sizes were necessary for simulations involving -larger cutoff values. +\section{Initial Free Energy Results} +The calculated free energies of proton-ordered variants of three low +density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it i} or +Ice-$i^\prime$) and the stable higher density ice B are listed in +table \ref{tab:freeEnergy}. Ice B was included because it has been +shown to be a minimum free energy structure for SPC/E at ambient +conditions.\cite{Baez95b} In addition to the free energies, the +relevant transition temperatures at standard pressure are also +displayed in table \ref{tab:freeEnergy}. These free energy values +indicate that Ice-{\it i} is the most stable state for all of the +investigated water models. With the free energy at these state +points, the Gibbs-Helmholtz equation was used to project to other +state points and to build phase diagrams. Figures +\ref{fig:tp3PhaseDia} and \ref{fig:ssdrfPhaseDia} are example diagrams +built from the results for the TIP3P and SSD/RF water models. All +other models have similar structure, although the crossing points +between the phases move to different temperatures and pressures as +indicated from the transition temperatures in table +\ref{tab:freeEnergy}. It is interesting to note that ice +I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear +in any of the phase diagrams for any of the models. For purposes of +this study, ice B is representative of the dense ice polymorphs. A +recent study by Sanz {\it et al.} provides details on the phase +diagrams for SPC/E and TIP4P at higher pressures than those studied +here.\cite{Sanz04} \begin{table} \centering -\caption{HELMHOLTZ FREE ENERGIES FOR SEVERAL ICE POLYMORPHS WITH A -VARIETY OF COMMON WATER MODELS AT 200 KELVIN AND 1 ATMOSPHERE} -\begin{tabular}{ l c c c c } +\caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1 +ATMOSPHERE FOR SEVERAL WATER MODELS} + +\footnotesize +\begin{tabular}{lccccccc} \toprule \toprule -Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or $i^\prime$) \\ -& kcal/mol & kcal/mol & kcal/mol & kcal/mol \\ +Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\ +\cmidrule(lr){2-6} +\cmidrule(l){7-8} +& \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\ \midrule -TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\ -TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\ -TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\ -SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\ -SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\ -SSD/RF & -11.51(4) & NA & -12.08(5) & -12.29(4)\\ +TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\ +TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\ +TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\ +SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\ +SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\ +SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\ \bottomrule \end{tabular} \label{tab:freeEnergy} \end{table} - -Table \ref{tab:freeEnergy} shows the results of the free energy -calculations with a cutoff radius of 9\AA. It should be noted that the -ice I$_\textrm{c}$ crystal polymorph is not stable at 200K and 1 atm -with the SSD/RF water model, hense omitted results for that cell. The -free energy values displayed in this table, it is clear that Ice-{\it -i} (or Ice-$i^\prime$ for TIP4P, TIP5P, and SPC/E) is the most stable -state for all of the common water models studied. - -With the free energy at these state points, the Gibbs-Helmholtz -equation was used to project to other state points and to build phase -diagrams. Figures \ref{fig:tp3phasedia} and \ref{fig:ssdrfphasedia} -are example diagrams built from the free energy results. All other -models have similar structure, although the crossing points between -the phases exist at slightly different temperatures and pressures. It -is interesting to note that ice I does not exist in either cubic or -hexagonal form in any of the phase diagrams for any of the models. For -purposes of this study, ice B is representative of the dense ice -polymorphs. A recent study by Sanz {\it et al.} goes into detail on -the phase diagrams for SPC/E and TIP4P in the high pressure -regime.\cite{Sanz04} \begin{figure} \centering @@ -254,7 +298,7 @@ higher in energy and don't appear in the phase diagram the experimental values; however, the solid phases shown are not the experimentally observed forms. Both cubic and hexagonal ice $I$ are higher in energy and don't appear in the phase diagram.} -\label{fig:tp3phasedia} +\label{fig:tp3PhaseDia} \end{figure} \begin{figure} @@ -266,178 +310,200 @@ conservative charge based models.} computationally efficient model (over 3 times more efficient than TIP3P) exhibits phase behavior similar to the less computationally conservative charge based models.} -\label{fig:ssdrfphasedia} +\label{fig:ssdrfPhaseDia} \end{figure} -\begin{table} -\centering -\caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$) -temperatures at 1 atm for several common water models compared with -experiment. The $T_m$ and $T_s$ values from simulation correspond to a -transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the -liquid or gas state.} -\begin{tabular}{ l c c c c c c c } -\toprule -\toprule -Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\ -\midrule -$T_m$ (K) & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\ -$T_b$ (K) & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\ -$T_s$ (K) & - & - & - & - & 355(3) & - & -\\ -\bottomrule -\end{tabular} -\label{tab:meltandboil} -\end{table} +We note that all of the crystals investigated in this study are ideal +proton-ordered antiferroelectric structures. All of the structures +obey the Bernal-Fowler rules and should be able to form stable +proton-{\it disordered} crystals which have the traditional +$k_\textrm{B}$ln(3/2) residual entropy at 0~K.\cite{Bernal33,Pauling35} +Simulations of proton-disordered structures are relatively unstable +with all but the most expensive water models.\cite{Nada03} Our +simulations have therefore been performed with the ordered +antiferroelectric structures which do not require the residual entropy +term to be accounted for in the free energies. This may result in some +discrepancies when comparing our melting temperatures to the melting +temperatures that have been calculated via thermodynamic integrations +of the disordered structures.\cite{Sanz04} -Table \ref{tab:meltandboil} lists the melting and boiling temperatures -calculated from this work. Surprisingly, most of these models have -melting points that compare quite favorably with experiment. The -unfortunate aspect of this result is that this phase change occurs -between Ice-{\it i} and the liquid state rather than ice -I$_\textrm{h}$ and the liquid state. These results are actually not -contrary to previous studies in the literature. Earlier free energy -studies of ice I using TIP4P predict a $T_m$ ranging from 214 to 238K +Most of the water models have melting points that compare quite +favorably with the experimental value of 273~K. The unfortunate +aspect of this result is that this phase change occurs between +Ice-{\it i} and the liquid state rather than ice I$_h$ and the liquid +state. These results do not contradict other studies. Studies of ice +I$_h$ using TIP4P predict a $T_m$ ranging from 191 to 238~K (differences being attributed to choice of interaction truncation and different ordered and disordered molecular -arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and -Ice-{\it i} were omitted, a $T_m$ value around 210K would be predicted -from this work. However, the $T_m$ from Ice-{\it i} is calculated at -265K, significantly higher in temperature than the previous -studies. Also of interest in these results is that SSD/E does not -exhibit a melting point at 1 atm, but it shows a sublimation point at -355K. 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 -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 300K are actually metastable and run the -risk of spontaneous crystallization. However, this risk changes when -applying a longer cutoff. +arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice +B and Ice-{\it i} were omitted, a $T_\textrm{m}$ value around 200~K +would be predicted from this work. However, the $T_\textrm{m}$ from +Ice-{\it i} is calculated to be 262~K, indicating that these +simulation based structures ought to be included in studies probing +phase transitions with this model. Also of interest in these results +is that SSD/E does not exhibit a melting point at 1 atm but does +sublime at 355~K. This is due to the significant stability of +Ice-{\it i} over all other polymorphs for this particular model under +these conditions. While troubling, this behavior resulted in the +spontaneous crystallization of Ice-{\it i} which led us to investigate +this structure. These observations provide a warning that simulations +of SSD/E as a ``liquid'' near 300~K are actually metastable and run +the risk of spontaneous crystallization. However, when a longer +cutoff radius is used, SSD/E prefers the liquid state under standard +temperature and pressure. +\section{Effects of Potential Truncation} + \begin{figure} \includegraphics[width=\linewidth]{./figures/cutoffChange.pdf} -\caption{Free energy as a function of cutoff radius for (A) SSD/E, (B) -TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: -I$_\textrm{c}$ 12\AA\, TIP3P: I$_\textrm{c}$ 12\AA\ and B 12\AA\, -and SSD/RF: I$_\textrm{c}$ 9\AA . These crystals are unstable at 200 K -and rapidly convert into liquids. The connecting lines are qualitative -visual aid.} +\caption{Free energy as a function of cutoff radius for SSD/E, TIP3P, +SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models +with an added Ewald correction term. Error for the larger cutoff +points is equivalent to that observed at 9.0~\AA\ (see Table +\ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using +both 12 and 13.5~\AA\ cutoffs were omitted because the crystal was +prone to distortion and melting at 200~K. Ice-$i^\prime$ is the +form of Ice-{\it i} used in the SPC/E simulations.} \label{fig:incCutoff} \end{figure} -Increasing the cutoff radius in simulations of the more -computationally efficient water models was done in order to evaluate -the trend in free energy values when moving to systems that do not -involve potential truncation. As seen in figure \ref{fig:incCutoff}, -the free energy of all the ice polymorphs show a substantial -dependence on cutoff radius. In general, there is a narrowing of the -free energy differences while moving to greater cutoff -radius. Interestingly, by increasing the cutoff radius, the free -energy gap was narrowed enough in the SSD/E model that the liquid -state is preferred under standard simulation conditions (298K 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. +For the more computationally efficient water models, we have also +investigated the effect of potential truncation on the computed free +energies as a function of the cutoff radius. As seen in +Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with +water models lacking a long-range correction show significant cutoff +dependence. In general, there is a narrowing of the free energy +differences while moving to greater cutoff radii. As the free +energies for the polymorphs converge, the stability advantage that +Ice-{\it i} exhibits is reduced. Adjacent to each of these plots are +results for systems with applied or estimated long-range corrections. +SSD/RF was parametrized for use with a reaction field, and the benefit +provided by this computationally inexpensive correction is apparent. +The free energies are largely independent of the size of the reaction +field cavity in this model, so small cutoff radii mimic bulk +calculations quite well under SSD/RF. + +Although TIP3P was parametrized for use without the Ewald summation, +we have estimated the effect of this method for computing long-range +electrostatics for both TIP3P and SPC/E. This was accomplished by +calculating the potential energy of identical crystals both with and +without particle mesh Ewald (PME). Similar behavior to that observed +with reaction field is seen for both of these models. The free +energies show reduced dependence on cutoff radius and span a narrower +range for the various polymorphs. Like the dipolar water models, +TIP3P displays a relatively constant preference for the Ice-{\it i} +polymorph. Crystal preference is much more difficult to determine for +SPC/E. Without a long-range correction, each of the polymorphs +studied assumes the role of the preferred polymorph under different +cutoff radii. The inclusion of the Ewald correction flattens and +narrows the gap in free energies such that the polymorphs are +isoenergetic within statistical uncertainty. This suggests that other +conditions, such as the density in fixed-volume simulations, can +influence the polymorph expressed upon crystallization. -To further study the changes resulting to the inclusion of a -long-range interaction correction, the effect of an Ewald summation -was estimated by applying the potential energy difference do to its -inclusion in systems in the presence and absence of the -correction. This was accomplished by calculation of the potential -energy of identical crystals with and without PME using TINKER. The -free energies for the investigated polymorphs using the TIP3P and -SPC/E water models are shown in Table \ref{tab:pmeShift}. TIP4P and -TIP5P are not fully supported in TINKER, so the results for these -models could not be estimated. The same trend pointed out through -increase of cutoff radius is observed in these PME results. Ice-{\it -i} is the preferred polymorph at ambient conditions for both the TIP3P -and SPC/E water models; however, there is a narrowing of the free -energy differences between the various solid forms. In the case of -SPC/E this narrowing is significant enough that it becomes less clear -that Ice-{\it i} is the most stable polymorph, and is possibly -metastable with respect to ice B and possibly ice -I$_\textrm{c}$. However, these results do not significantly alter the -finding that the Ice-{\it i} polymorph is a stable crystal structure -that should be considered when studying the phase behavior of water -models. +\section{Expanded Results Using Damped Shifted Force Electrostatics} +In chapter \ref{chap:electrostatics}, we discussed in detail a +pairwise method for handling electrostatics (shifted force, {\sc sf}) +that can be used as a simple and efficient replacement for the Ewald +summation. Answering the question of the free energies of these ice +polymorphs with varying water models would be an interesting +application of this technique. To this end, we set up thermodynamic +integrations of all of the previously discussed ice polymorphs using +the {\sc sf} technique with a cutoff radius of 12~\AA\ and an $\alpha$ +of 0.2125~\AA . These calculations were performed on TIP5P-E and +TIP4P-Ew (variants of the root models optimized for the Ewald +summation) as well as SPC/E, SSD/RF, and TRED (see section +\ref{sec:tredWater}). + \begin{table} \centering -\caption{The free energy of the studied ice polymorphs after applying -the energy difference attributed to the inclusion of the PME -long-range interaction correction. Units are kcal/mol.} -\begin{tabular}{ l c c c c } +\caption{HELMHOLTZ FREE ENERGIES OF ICE POLYMORPHS USING THE DAMPED +SHIFTED FORCE CORRECTION} +\begin{tabular}{ lccccc } \toprule \toprule -Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or Ice-$i^\prime$) \\ -\midrule -TIP3P & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\ -SPC/E & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\ -\bottomrule +Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\ +\cmidrule(lr){2-6} +& \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\ +\midrule +TIP5P-E & -10.76(4) & -10.72(4) & & - & -10.68(4) \\ +TIP4P-Ew & & -11.77(3) & & - & -11.60(3) \\ +SPC/E & -12.98(3) & -11.60(3) & & - & -12.93(3) \\ +SSD/RF & -11.81(4) & -11.65(3) & & -12.41(4) & - \\ +TRED & -12.58(3) & -12.44(3) & & -13.09(4) & - \\ \end{tabular} -\label{tab:pmeShift} +\label{tab:dampedFreeEnergy} \end{table} + \section{Conclusions} -The free energy for proton ordered variants of hexagonal and cubic ice -$I$, ice B, and recently discovered Ice-{\it i} were calculated under -standard conditions for several common water models via thermodynamic -integration. All the water models studied show Ice-{\it i} to be the -minimum free energy crystal structure in the with a 9\AA\ switching -function cutoff. Calculated melting and boiling points show -surprisingly good agreement with the experimental values; however, the -solid phase at 1 atm is Ice-{\it i}, not ice I$_\textrm{h}$. The -effect of interaction truncation was investigated through variation of -the cutoff radius, use of a reaction field parameterized model, and -estimation of the results in the presence of the Ewald -summation. Interaction truncation has a significant effect on the -computed free energy values, and may significantly alter the free -energy landscape for the more complex multipoint water models. Despite -these effects, these results show Ice-{\it i} to be an important ice -polymorph that should be considered in simulation studies. +In this work, thermodynamic integration was used to determine the +absolute free energies of several ice polymorphs. The new polymorph, +Ice-{\it i} was observed to be the stable crystalline state for {\it +all} the water models when using a 9.0~\AA\ cutoff. However, the free +energy partially depends on simulation conditions (particularly on the +choice of long range correction method). Regardless, Ice-{\it i} was +still observed to be a stable polymorph for all of the studied water +models. -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 makes us hesitant to ascribe any relevance of this work outside -of the simulation community. It is for this reason that we chose a -name for this polymorph which involves an imaginary quantity. That -said, there are certain experimental conditions that would provide the -most ideal situation for possible observation. These include the -negative pressure or stretched solid regime, small clusters in vacuum +So what is the preferred solid polymorph for simulated water? As +indicated above, the answer appears to be dependent both on the +conditions and the model used. In the case of short cutoffs without a +long-range interaction correction, Ice-{\it i} and Ice-$i^\prime$ have +the lowest free energy of the studied polymorphs with all the models. +Ideally, crystallization of each model under constant pressure +conditions, as was done with SSD/E, would aid in the identification of +their respective preferred structures. This work, however, helps +illustrate how studies involving one specific model can lead to +insight about important behavior of others. + +We also note that none of the water models used in this study are +polarizable or flexible models. It is entirely possible that the +polarizability of real water makes Ice-{\it i} substantially less +stable than ice I$_h$. However, the calculations presented above seem +interesting enough to communicate before the role of polarizability +(or flexibility) has been thoroughly investigated. + +Finally, due to the stability of Ice-{\it i} in the investigated +simulation conditions, the question arises as to possible experimental +observation of this polymorph. The rather extensive past and current +experimental investigation of water in the low pressure regime makes +us hesitant to ascribe any relevance to this work outside of the +simulation community. It is for this reason that we chose a name for +this polymorph which involves an imaginary quantity. That said, there +are certain experimental conditions that would provide the most ideal +situation for possible observation. These include the negative +pressure or stretched solid regime, small clusters in vacuum deposition environments, and in clathrate structures involving small -non-polar molecules. Figs. \ref{fig:gofr} and \ref{fig:sofq} contain -our predictions for both the pair distribution function ($g_{OO}(r)$) -and the structure factor ($S(\vec{q})$ for ice I$_\textrm{c}$ and for -ice-{\it i} at a temperature of 77K. In a quick comparison of the -predicted S(q) for Ice-{\it i} and experimental studies of amorphous -solid water, it is possible that some of the ``spurious'' peaks that -could not be assigned in HDA could correspond to peaks labeled in this -S(q).\cite{Bizid87} It should be noted that there is typically poor -agreement on crystal densities between simulation and experiment, so -such peak comparisons should be made with caution. We will leave it -to our experimental colleagues to determine whether this ice polymorph -is named appropriately or if it should be promoted to Ice-0. +non-polar molecules. For the purpose of comparison with experimental +results, we have calculated the oxygen-oxygen pair correlation +function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$ +for the two Ice-{\it i} variants (along with example ice +I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77~K, and they are shown in +figures \ref{fig:gofr} and \ref{fig:sofq} respectively. It is +interesting to note that the structure factors for Ice-$i^\prime$ and +Ice-I$_c$ are quite similar. The primary differences are small peaks +at 1.125, 2.29, and 2.53~\AA$^{-1}$, so particular attention to these +regions would be needed to identify the new $i^\prime$ variant from +the I$_\textrm{c}$ polymorph. + \begin{figure} \includegraphics[width=\linewidth]{./figures/iceGofr.pdf} \caption{Radial distribution functions of Ice-{\it i} and ice I$_\textrm{c}$ calculated from from simulations of the SSD/RF water -model at 77 K.} +model at 77~K.} \label{fig:gofr} \end{figure} \begin{figure} \includegraphics[width=\linewidth]{./figures/sofq.pdf} \caption{Predicted structure factors for Ice-{\it i} and ice -I$_\textrm{c}$ at 77 K. The raw structure factors have been -convoluted with a gaussian instrument function (0.075 \AA$^{-1}$ -width) to compensate for the trunction effects in our finite size +I$_\textrm{c}$ at 77~K. The raw structure factors have been +convoluted with a gaussian instrument function (0.075~\AA$^{-1}$ +width) to compensate for the truncation effects in our finite size simulations. The labeled peaks compared favorably with ``spurious'' peaks observed in experimental studies of amorphous solid water.\cite{Bizid87}}