--- trunk/ssdePaper/nptSSD.tex 2004/01/12 16:20:53 921 +++ trunk/ssdePaper/nptSSD.tex 2004/02/05 18:42:59 1027 @@ -139,6 +139,7 @@ + s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i \frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) + s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)]\ . +\label{stickyfunction} \end{equation} Here, $\nu_0$ is a strength parameter for the sticky potential, and $s$ and $s^\prime$ are cubic switching functions which turn off the @@ -151,7 +152,7 @@ w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = while the $w^\prime$ function counters the normal aligned and anti-aligned structures favored by point dipoles: \begin{equation} -w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0, +w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^\circ, \end{equation} It should be noted that $w$ is proportional to the sum of the $Y_3^2$ and $Y_3^{-2}$ spherical harmonics (a linear combination which @@ -208,17 +209,19 @@ cubic switching function at a cutoff radius. Under th Long-range dipole-dipole interactions were accounted for in this study by using either the reaction field method or by resorting to a simple -cubic switching function at a cutoff radius. Under the first method, -the magnitude of the reaction field acting on dipole $i$ is +cubic switching function at a cutoff radius. The reaction field +method was actually first used in Monte Carlo simulations of liquid +water.\cite{Barker73} Under this method, the magnitude of the reaction +field acting on dipole $i$ is \begin{equation} \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1} -\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij})\ , +\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij}), \label{rfequation} \end{equation} where $\mathcal{R}$ is the cavity defined by the cutoff radius ($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole -moment vector of particle $j$ and $f(r_{ij})$ is a cubic switching +moment vector of particle $j$, and $f(r_{ij})$ is a cubic switching function.\cite{AllenTildesley} The reaction field contribution to the total energy by particle $i$ is given by $-\frac{1}{2}{\bf \mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf @@ -232,8 +235,8 @@ at the cutoff radius) and as a result we have two repa We have also performed a companion set of simulations {\it without} a surrounding dielectric (i.e. using a simple cubic switching function -at the cutoff radius) and as a result we have two reparamaterizations -of SSD which could be used either with or without the Reaction Field +at the cutoff radius), and as a result we have two reparamaterizations +of SSD which could be used either with or without the reaction field turned on. Simulations to obtain the preferred density were performed in the @@ -248,15 +251,14 @@ symplectic splitting method proposed by Dullweber {\it need to be considered. Integration of the equations of motion was carried out using the -symplectic splitting method proposed by Dullweber {\it et -al.}\cite{Dullweber1997} Our reason for selecting this integrator -centers on poor energy conservation of rigid body dynamics using -traditional quaternion integration.\cite{Evans77,Evans77b} While quaternions -may work well for orientational motion under NVT or NPT integrators, -our limits on energy drift in the microcanonical ensemble were quite -strict, and the drift under quaternions was substantially greater than -in the symplectic splitting method. This steady drift in the total -energy has also been observed by Kol {\it et al.}\cite{Laird97} +symplectic splitting method proposed by Dullweber, Leimkuhler, and +McLachlan (DLM).\cite{Dullweber1997} Our reason for selecting this +integrator centers on poor energy conservation of rigid body dynamics +using traditional quaternion integration.\cite{Evans77,Evans77b} In +typical microcanonical ensemble simulations, the energy drift when +using quaternions was substantially greater than when using the DLM +method (fig. \ref{timestep}). This steady drift in the total energy +has also been observed by Kol {\it et al.}\cite{Laird97} The key difference in the integration method proposed by Dullweber \emph{et al.} is that the entire rotation matrix is propagated from @@ -265,19 +267,19 @@ The symplectic splitting method allows for Verlet styl rotation matrix into quaternions for storage purposes makes trajectory data quite compact. -The symplectic splitting method allows for Verlet style integration of -both translational and orientational motion of rigid bodies. In this +The DML method allows for Verlet style integration of both +translational and orientational motion of rigid bodies. In this integration method, the orientational propagation involves a sequence of matrix evaluations to update the rotation matrix.\cite{Dullweber1997} These matrix rotations are more costly than the simpler arithmetic quaternion propagation. With the same time step, a 1000 SSD particle simulation shows an average 7\% increase in -computation time using the symplectic step method in place of -quaternions. The additional expense per step is justified when one -considers the ability to use time steps that are nearly twice as large -under symplectic splitting than would be usable under quaternion -dynamics. The energy conservation of the two methods using a number -of different time steps is illustrated in figure +computation time using the DML method in place of quaternions. The +additional expense per step is justified when one considers the +ability to use time steps that are nearly twice as large under DML +than would be usable under quaternion dynamics. The energy +conservation of the two methods using a number of different time steps +is illustrated in figure \ref{timestep}. \begin{figure} @@ -285,33 +287,33 @@ the symplectic step method proposed by Dullweber \emph \epsfxsize=6in \epsfbox{timeStep.epsi} \caption{Energy conservation using both quaternion based integration and -the symplectic step method proposed by Dullweber \emph{et al.} with -increasing time step. The larger time step plots are shifted from the -true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.} +the symplectic splitting method proposed by Dullweber \emph{et al.} +with increasing time step. The larger time step plots are shifted from +the true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.} \label{timestep} \end{center} \end{figure} In figure \ref{timestep}, the resulting energy drift at various time -steps for both the symplectic step and quaternion integration schemes -is compared. All of the 1000 SSD particle simulations started with -the same configuration, and the only difference was the method used to -handle orientational motion. At time steps of 0.1 and 0.5 fs, both -methods for propagating the orientational degrees of freedom conserve -energy fairly well, with the quaternion method showing a slight energy -drift over time in the 0.5 fs time step simulation. At time steps of 1 -and 2 fs, the energy conservation benefits of the symplectic step -method are clearly demonstrated. Thus, while maintaining the same -degree of energy conservation, one can take considerably longer time -steps, leading to an overall reduction in computation time. +steps for both the DML and quaternion integration schemes is compared. +All of the 1000 SSD particle simulations started with the same +configuration, and the only difference was the method used to handle +orientational motion. At time steps of 0.1 and 0.5 fs, both methods +for propagating the orientational degrees of freedom conserve energy +fairly well, with the quaternion method showing a slight energy drift +over time in the 0.5 fs time step simulation. At time steps of 1 and 2 +fs, the energy conservation benefits of the DML method are clearly +demonstrated. Thus, while maintaining the same degree of energy +conservation, one can take considerably longer time steps, leading to +an overall reduction in computation time. -Energy drift in the symplectic step simulations was unnoticeable for -time steps up to 3 fs. A slight energy drift on the -order of 0.012 kcal/mol per nanosecond was observed at a time step of -4 fs, and as expected, this drift increases dramatically -with increasing time step. To insure accuracy in our microcanonical -simulations, time steps were set at 2 fs and kept at this value for -constant pressure simulations as well. +Energy drift in the simulations using DML integration was unnoticeable +for time steps up to 3 fs. A slight energy drift on the order of 0.012 +kcal/mol per nanosecond was observed at a time step of 4 fs, and as +expected, this drift increases dramatically with increasing time +step. To insure accuracy in our microcanonical simulations, time steps +were set at 2 fs and kept at this value for constant pressure +simulations as well. Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices were generated as starting points for all simulations. The $I_h$ @@ -428,7 +430,7 @@ was what lead Ichiye {\it et al.} to recently reparame simulations.\cite{Clancy94,Jorgensen98b} However, even without the reaction field, the density around 300 K is still significantly lower than experiment and comparable water models. This anomalous behavior -was what lead Ichiye {\it et al.} to recently reparameterize +was what lead Tan {\it et al.} to recently reparameterize SSD.\cite{Ichiye03} Throughout the remainder of the paper our reparamaterizations of SSD will be compared with the newer SSD1 model. @@ -446,7 +448,7 @@ results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} mean-square displacement as a function of time. The averaged results from five sets of NVE simulations are displayed in figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P -results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} +results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} \begin{figure} \begin{center} @@ -454,7 +456,7 @@ and Experimental data [Refs. \citen{Gillen72} and \cit \epsfbox{betterDiffuse.epsi} \caption{Average self-diffusion constant as a function of temperature for SSD, SPC/E [Ref. \citen{Clancy94}], TIP5P [Ref. \citen{Jorgensen01}], -and Experimental data [Refs. \citen{Gillen72} and \citen{Mills73}]. Of +and Experimental data [Refs. \citen{Gillen72} and \citen{Holz00}]. Of the three water models shown, SSD has the least deviation from the experimental values. The rapidly increasing diffusion constants for TIP5P and SSD correspond to significant decrease in density at the @@ -591,29 +593,15 @@ The parameters available for tuning include the $\sigm important properties. In this case, it would be ideal to correct the densities while maintaining the accurate transport behavior. -The parameters available for tuning include the $\sigma$ and $\epsilon$ -Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky -attractive and dipole repulsive terms with their respective -cutoffs. To alter the attractive and repulsive terms of the sticky -potential independently, it is necessary to separate the terms as -follows: -\begin{equation} -u_{ij}^{sp} -({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = -\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)], -\end{equation} -where $\nu_0$ scales the strength of the tetrahedral attraction and -$\nu_0^\prime$ scales the dipole repulsion term independently. The -separation was performed for purposes of the reparameterization, but -the final parameters were adjusted so that it is not necessary to -separate the terms when implementing the adjusted water -potentials. The results of the reparameterizations are shown in table -\ref{params}. Note that the tetrahedral attractive and dipolar -repulsive terms do not share the same lower cutoff ($r_l$) in the -newly parameterized potentials. We are calling these -reparameterizations the Soft Sticky Dipole / Reaction Field +The parameters available for tuning include the $\sigma$ and +$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the +strength of the sticky potential ($\nu_0$), and the sticky attractive +and dipole repulsive cubic switching function cutoffs ($r_l$, $r_u$ +and $r_l^\prime$, $r_u^\prime$ respectively). The results of the +reparameterizations are shown in table \ref{params}. We are calling +these reparameterizations the Soft Sticky Dipole / Reaction Field (SSD/RF - for use with a reaction field) and Soft Sticky Dipole -Enhanced (SSD/E - an attempt to improve the liquid structure in +Extended (SSD/E - an attempt to improve the liquid structure in simulations without a long-range correction). \begin{table} @@ -628,9 +616,9 @@ simulations without a long-range correction). \ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ \ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\ \ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ +\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\ \ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\ \ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\ -\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\ \ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\ \ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\ \end{tabular} @@ -654,7 +642,7 @@ and broadened the first peak of SSD/E and SSD/RF.} \begin{figure} \begin{center} \epsfxsize=6in -\epsfbox{dualsticky.ps} +\epsfbox{dualsticky_bw.eps} \caption{Isosurfaces of the sticky potential for SSD1 (left) and SSD/E \& SSD/RF (right). Light areas correspond to the tetrahedral attractive component, and darker areas correspond to the dipolar repulsive @@ -806,13 +794,14 @@ K, shown by SSD and SSD1 respectively. \begin{center} \epsfxsize=6in \epsfbox{ssdeDiffuse.epsi} -\caption{Plots of the diffusion constants calculated from SSD/E and SSD1, -both without a reaction field, along with experimental results -[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations were -performed at the average densities observed in the 1 atm NPT -simulations for the respective models. SSD/E is slightly more fluid -than experiment at all of the temperatures, but it is closer than SSD1 -without a long-range correction.} +\caption{The diffusion constants calculated from SSD/E and SSD1, + both without a reaction field, along with experimental results + [Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations + were performed at the average densities observed in the 1 atm NPT + simulations for the respective models. SSD/E is slightly more mobile + than experiment at all of the temperatures, but it is closer to + experiment at biologically relevant temperatures than SSD1 without a + long-range correction.} \label{ssdediffuse} \end{center} \end{figure} @@ -823,35 +812,38 @@ without an active reaction field, both at the densitie the densities, it is important that the excellent diffusive behavior of SSD be maintained or improved. Figure \ref{ssdediffuse} compares the temperature dependence of the diffusion constant of SSD/E to SSD1 -without an active reaction field, both at the densities calculated at -1 atm and at the experimentally calculated densities for super-cooled -and liquid water. The diffusion constant for SSD/E is consistently -higher than experiment, while SSD1 remains lower than experiment until -relatively high temperatures (greater than 330 K). Both models follow -the shape of the experimental curve well below 300 K but tend to -diffuse too rapidly at higher temperatures, something that is -especially apparent with SSD1. This increasing diffusion relative to +without an active reaction field at the densities calculated from +their respective NPT simulations at 1 atm. The diffusion constant for +SSD/E is consistently higher than experiment, while SSD1 remains lower +than experiment until relatively high temperatures (around 360 +K). Both models follow the shape of the experimental curve well below +300 K but tend to diffuse too rapidly at higher temperatures, as seen +in SSD1's crossing above 360 K. This increasing diffusion relative to the experimental values is caused by the rapidly decreasing system -density with increasing temperature. The densities of SSD1 decay more -rapidly with temperature than do those of SSD/E, leading to more -visible deviation from the experimental diffusion trend. Thus, the +density with increasing temperature. Both SSD1 and SSD/E show this +deviation in particle mobility, but this trend has different +implications on the diffusive behavior of the models. While SSD1 +shows more experimentally accurate diffusive behavior in the high +temperature regimes, SSD/E shows more accurate behavior in the +supercooled and biologically relevant temperature ranges. Thus, the changes made to improve the liquid structure may have had an adverse affect on the density maximum, but they improve the transport behavior -of SSD/E relative to SSD1. +of SSD/E relative to SSD1 under the most commonly simulated +conditions. \begin{figure} \begin{center} \epsfxsize=6in \epsfbox{ssdrfDiffuse.epsi} -\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1, +\caption{The diffusion constants calculated from SSD/RF and SSD1, both with an active reaction field, along with experimental results - [Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations + [Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations were performed at the average densities observed in the 1 atm NPT simulations for both of the models. Note how accurately SSD/RF simulates the diffusion of water throughout this temperature range. The more rapidly increasing diffusion constants at high - temperatures for both models is attributed to the significantly lower - densities than observed in experiment.} + temperatures for both models is attributed to lower calculated + densities than those observed in experiment.} \label{ssdrfdiffuse} \end{center} \end{figure} @@ -859,9 +851,9 @@ throughout the temperature range shown and with only a In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are compared to SSD1 with an active reaction field. Note that SSD/RF tracks the experimental results quantitatively, identical within error -throughout the temperature range shown and with only a slight -increasing trend at higher temperatures. SSD1 tends to diffuse more -slowly at low temperatures and deviates to diffuse too rapidly at +throughout most of the temperature range shown and exhibiting only a +slight increasing trend at higher temperatures. SSD1 tends to diffuse +more slowly at low temperatures and deviates to diffuse too rapidly at temperatures greater than 330 K. As stated above, this deviation away from the ideal trend is due to a rapid decrease in density at higher temperatures. SSD/RF does not suffer from this problem as much as SSD1 @@ -869,12 +861,89 @@ reparameterization when using an altered long-range co values. These results again emphasize the importance of careful reparameterization when using an altered long-range correction. +\begin{table} +\begin{center} +\caption{Calculated and experimental properties of the single point waters and liquid water at 298 K and 1 atm. (a) Ref. [\citen{Mills73}]. (b) Calculated by integrating the data in ref. \citen{Head-Gordon00_1}. (c) Calculated by integrating the data in ref. \citen{Soper86}. (d) Calculated for 298 K from data in ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.} +\begin{tabular}{ l c c c c c } +\hline \\[-3mm] +\ \ \ \ \ \ & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ +\ & \ SSD/RF \ \ \ & \ Expt. \\ + \hline \\[-3mm] +\ \ \ $\rho$ (g/cm$^3$) & 0.999 $\pm$0.001 & 0.996 $\pm$0.001 & 0.972 $\pm$0.002 & 0.997 $\pm$0.001 & 0.997 \\ +\ \ \ $C_p$ (cal/mol K) & 28.80 $\pm$0.11 & 25.45 $\pm$0.09 & 28.28 $\pm$0.06 & 23.83 $\pm$0.16 & 17.98 \\ +\ \ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78 $\pm$0.07 & 2.51 $\pm$0.18 & 2.00 $\pm$0.17 & 2.32 $\pm$0.06 & 2.299$^\text{a}$ \\ +\ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\ +\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.5$^\text{c}$ \\ +\ \ \ $\tau_1$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 5.7$^\text{d}$ \\ +\ \ \ $\tau_2$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2 $\pm$0.2 & 2.3$^\text{e}$ \\ +\end{tabular} +\label{liquidproperties} +\end{center} +\end{table} + +Table \ref{liquidproperties} gives a synopsis of the liquid state +properties of the water models compared in this study along with the +experimental values for liquid water at ambient conditions. The +coordination number ($N_C$) and hydrogen bonds per particle ($N_H$) +were calculated by integrating the following relations: +\begin{equation} +N_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr, +\end{equation} +\begin{equation} +N_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr, +\end{equation} +where $\rho$ is the number density of the specified pair interactions, +$a$ and $b$ are the radial locations of the minima following the first +solvation shell peak in g$_\text{OO}(r)$ or g$_\text{OH}(r)$ +respectively. The number of hydrogen bonds stays relatively constant +across all of the models, but the coordination numbers of SSD/E and +SSD/RF show an improvement over SSD1. This improvement is primarily +due to the widening of the first solvation shell peak, allowing the +first minima to push outward. Comparing the coordination number with +the number of hydrogen bonds can lead to more insight into the +structural character of the liquid. Because of the near identical +values for SSD1, it appears to be a little too exclusive, in that all +molecules in the first solvation shell are involved in forming ideal +hydrogen bonds. The differing numbers for the newly parameterized +models indicate the allowance of more fluid configurations in addition +to the formation of an acceptable number of ideal hydrogen bonds. + +The time constants for the self orientational autocorrelation function +are also displayed in Table \ref{liquidproperties}. The dipolar +orientational time correlation function ($\Gamma_{l}$) is described +by: +\begin{equation} +\Gamma_{l}(t) = \langle P_l[\mathbf{u}_j(0)\cdot\mathbf{u}_j(t)]\rangle, +\end{equation} +where $P_l$ is a Legendre polynomial of order $l$ and $\mathbf{u}_j$ +is the unit vector of the particle dipole.\cite{Rahman71} From these +correlation functions, the orientational relaxation time of the dipole +vector can be calculated from an exponential fit in the long-time +regime ($t > \tau_l$).\cite{Rothschild84} Calculation of these +time constants were averaged from five detailed NVE simulations +performed at the STP density for each of the respective models. It +should be noted that the commonly cited value for $\tau_2$ of 1.9 ps +was determined from the NMR data in reference \citen{Krynicki66} at a +temperature near 34$^\circ$C.\cite{Rahman71} Because of the strong +temperature dependence of $\tau_2$, it is necessary to recalculate it +at 298 K to make proper comparisons. The value shown in Table +\ref{liquidproperties} was calculated from the same NMR data in the +fashion described in reference \citen{Krynicki66}. Similarly, $\tau_1$ +was recomputed for 298 K from the data in ref \citen{Eisenberg69}. +Again, SSD/E and SSD/RF show improved behavior over SSD1, both with +and without an active reaction field. Turning on the reaction field +leads to much improved time constants for SSD1; however, these results +also include a corresponding decrease in system density. Numbers +published from the original SSD dynamics studies are shorter than the +values observed here, and this difference can be attributed to the use +of the Ewald sum technique versus a reaction field.\cite{Ichiye99} + \subsection{Additional Observations} \begin{figure} \begin{center} \epsfxsize=6in -\epsfbox{povIce.ps} +\epsfbox{icei_bw.eps} \caption{A water lattice built from the crystal structure assumed by SSD/E when undergoing an extremely restricted temperature NPT simulation. This form of ice is referred to as ice-{\it i} to @@ -954,7 +1023,7 @@ super-cooled liquid regimes. significantly lower than the densities obtained from other water models (and experiment). Analysis of self-diffusion showed SSD to capture the transport properties of water well in both the liquid and -super-cooled liquid regimes. +supercooled liquid regimes. In order to correct the density behavior, the original SSD model was reparameterized for use both with and without a reaction field (SSD/RF @@ -970,8 +1039,8 @@ effectively simulations of super-cooled or metastable The existence of a novel low-density ice structure that is preferred by the SSD family of water models is somewhat troubling, since liquid simulations on this family of water models at room temperature are -effectively simulations of super-cooled or metastable liquids. One -way to de-stabilize this unphysical ice structure would be to make the +effectively simulations of supercooled or metastable liquids. One +way to destabilize this unphysical ice structure would be to make the range of angles preferred by the attractive part of the sticky potential much narrower. This would require extensive reparameterization to maintain the same level of agreement with the