--- trunk/ssdePaper/nptSSD.tex 2004/02/05 22:54:57 1030 +++ trunk/ssdePaper/nptSSD.tex 2004/02/06 21:43:00 1036 @@ -1,4 +1,5 @@ %\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4} +%\documentclass[preprint,aps,endfloat]{revtex4} \documentclass[11pt]{article} \usepackage{endfloat} \usepackage{amsmath} @@ -8,45 +9,45 @@ \usepackage{tabularx} \usepackage{graphicx} \usepackage[ref]{overcite} -%\usepackage{berkeley} -%\usepackage{curves} \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight 9.0in \textwidth 6.5in \brokenpenalty=10000 -\renewcommand{\baselinestretch}{1.2} -\renewcommand\citemid{\ } % no comma in optional reference note +%\renewcommand\citemid{\ } % no comma in optional reference note + \begin{document} \title{On the structural and transport properties of the soft sticky dipole (SSD) and related single point water models} -\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ -Department of Chemistry and Biochemistry\\ University of Notre Dame\\ +\author{Christopher J. Fennell and J. Daniel +Gezelter\footnote{Corresponding author. \ Electronic mail: +gezelter@nd.edu} \\ Department of Chemistry and Biochemistry\\ University of Notre Dame\\ Notre Dame, Indiana 46556} \date{\today} \maketitle +\doublespacing \begin{abstract} The density maximum and temperature dependence of the self-diffusion constant were investigated for the soft sticky dipole (SSD) water -model and two related re-parameterizations of this single-point model. +model and two related reparameterizations of this single-point model. A combination of microcanonical and isobaric-isothermal molecular dynamics simulations were used to calculate these properties, both with and without the use of reaction field to handle long-range electrostatics. The isobaric-isothermal (NPT) simulations of the melting of both ice-$I_h$ and ice-$I_c$ showed a density maximum near -260 K. In most cases, the use of the reaction field resulted in +260~K. In most cases, the use of the reaction field resulted in calculated densities which were were significantly lower than experimental densities. Analysis of self-diffusion constants shows that the original SSD model captures the transport properties of experimental water very well in both the normal and super-cooled -liquid regimes. We also present our re-parameterized versions of SSD +liquid regimes. We also present our reparameterized versions of SSD for use both with the reaction field or without any long-range electrostatic corrections. These are called the SSD/RF and SSD/E models respectively. These modified models were shown to maintain or @@ -62,7 +63,6 @@ family. %\narrowtext - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % BODY OF TEXT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -94,7 +94,7 @@ which has an interaction site that is both a point dip was developed by Ichiye \emph{et al.} as a modified form of the hard-sphere water model proposed by Bratko, Blum, and Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model -which has an interaction site that is both a point dipole along with a +which has an interaction site that is both a point dipole and a Lennard-Jones core. However, since the normal aligned and anti-aligned geometries favored by point dipoles are poor mimics of local structure in liquid water, a short ranged ``sticky'' potential @@ -164,12 +164,12 @@ simulations using this model, Ichiye {\it et al.} repo Since SSD is a single-point {\it dipolar} model, the force calculations are simplified significantly relative to the standard {\it charged} multi-point models. In the original Monte Carlo -simulations using this model, Ichiye {\it et al.} reported that using -SSD decreased computer time by a factor of 6-7 compared to other +simulations using this model, Liu and Ichiye reported that using SSD +decreased computer time by a factor of 6-7 compared to other models.\cite{Ichiye96} What is most impressive is that this savings did not come at the expense of accurate depiction of the liquid state -properties. Indeed, SSD maintains reasonable agreement with the -Soper data for the structural features of liquid +properties. Indeed, SSD maintains reasonable agreement with the Soper +data for the structural features of liquid water.\cite{Soper86,Ichiye96} Additionally, the dynamical properties exhibited by SSD agree with experiment better than those of more computationally expensive models (like TIP3P and @@ -203,16 +203,16 @@ utilizing the Reaction Field. follows, we compare our reparamaterization of SSD with both the original SSD and SSD1 models with the goal of improving the bulk phase behavior of an SSD-derived model in simulations -utilizing the Reaction Field. +utilizing the reaction field. \section{Methods} 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. 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 +by using either the reaction field technique or by resorting to a +simple cubic switching function at a cutoff radius. One of the early +applications of a reaction field was actually 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} s(r_{ij}), @@ -226,12 +226,12 @@ field is known to alter the bulk orientational propert 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 \mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley} Use of the reaction -field is known to alter the bulk orientational properties, such as the -dielectric relaxation time. There is particular sensitivity of this -property on changes in the length of the cutoff -radius.\cite{Berendsen98} This variable behavior makes reaction field -a less attractive method than the Ewald sum. However, for very large -systems, the computational benefit of reaction field is dramatic. +field is known to alter the bulk orientational properties of simulated +water, and there is particular sensitivity of these properties on +changes in the length of the cutoff radius.\cite{Berendsen98} This +variable behavior makes reaction field a less attractive method than +the Ewald sum. However, for very large systems, the computational +benefit of reaction field is dramatic. We have also performed a companion set of simulations {\it without} a surrounding dielectric (i.e. using a simple cubic switching function @@ -287,10 +287,10 @@ time steps is illustrated in figure \begin{center} \epsfxsize=6in \epsfbox{timeStep.epsi} -\caption{Energy conservation using both quaternion-based integration and -the {\sc dlm} method 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.} +\caption{Energy conservation using both quaternion-based integration and the +{\sc dlm} method 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} @@ -299,38 +299,38 @@ handle orientational motion. At time steps of 0.1 and steps for both the {\sc dlm} 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 +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 {\sc dlm} method are +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 {\sc dlm} 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 simulations using {\sc dlm} 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 +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 +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$ -crystals were formed by first arranging the centers of mass of the -SSD particles into a ``hexagonal'' ice lattice of 1024 -particles. Because of the crystal structure of $I_h$ ice, the -simulation box assumed an orthorhombic shape with an edge length ratio -of approximately 1.00$\times$1.06$\times$1.23. The particles were then -allowed to orient freely about fixed positions with angular momenta -randomized at 400 K for varying times. The rotational temperature was -then scaled down in stages to slowly cool the crystals to 25 K. The -particles were then allowed to translate with fixed orientations at a -constant pressure of 1 atm for 50 ps at 25 K. Finally, all constraints -were removed and the ice crystals were allowed to equilibrate for 50 -ps at 25 K and a constant pressure of 1 atm. This procedure resulted -in structurally stable $I_h$ ice crystals that obey the Bernal-Fowler +crystals were formed by first arranging the centers of mass of the SSD +particles into a ``hexagonal'' ice lattice of 1024 particles. Because +of the crystal structure of $I_h$ ice, the simulation box assumed an +orthorhombic shape with an edge length ratio of approximately +1.00$\times$1.06$\times$1.23. The particles were then allowed to +orient freely about fixed positions with angular momenta randomized at +400~K for varying times. The rotational temperature was then scaled +down in stages to slowly cool the crystals to 25~K. The particles were +then allowed to translate with fixed orientations at a constant +pressure of 1 atm for 50~ps at 25~K. Finally, all constraints were +removed and the ice crystals were allowed to equilibrate for 50~ps at +25~K and a constant pressure of 1~atm. This procedure resulted in +structurally stable $I_h$ ice crystals that obey the Bernal-Fowler rules.\cite{Bernal33,Rahman72} This method was also utilized in the making of diamond lattice $I_c$ ice crystals, with each cubic simulation box consisting of either 512 or 1000 particles. Only @@ -347,32 +347,32 @@ for 100 ps prior to a 200 ps data collection run at ea supercooled regime. An ensemble average from five separate melting simulations was acquired, each starting from different ice crystals generated as described previously. All simulations were equilibrated -for 100 ps prior to a 200 ps data collection run at each temperature -setting. The temperature range of study spanned from 25 to 400 K, with -a maximum degree increment of 25 K. For regions of interest along this -stepwise progression, the temperature increment was decreased from 25 -K to 10 and 5 K. The above equilibration and production times were +for 100~ps prior to a 200~ps data collection run at each temperature +setting. The temperature range of study spanned from 25 to 400~K, with +a maximum degree increment of 25~K. For regions of interest along this +stepwise progression, the temperature increment was decreased from +25~K to 10 and 5~K. The above equilibration and production times were sufficient in that fluctuations in the volume autocorrelation function -were damped out in all simulations in under 20 ps. +were damped out in all simulations in under 20~ps. \subsection{Density Behavior} Our initial simulations focused on the original SSD water model, and an average density versus temperature plot is shown in figure \ref{dense1}. Note that the density maximum when using a reaction -field appears between 255 and 265 K. There were smaller fluctuations -in the density at 260 K than at either 255 or 265, so we report this +field appears between 255 and 265~K. There were smaller fluctuations +in the density at 260~K than at either 255 or 265~K, so we report this value as the location of the density maximum. Figure \ref{dense1} was constructed using ice $I_h$ crystals for the initial configuration; though not pictured, the simulations starting from ice $I_c$ crystal configurations showed similar results, with a liquid-phase density -maximum in this same region (between 255 and 260 K). +maximum in this same region (between 255 and 260~K). \begin{figure} \begin{center} \epsfxsize=6in \epsfbox{denseSSDnew.eps} -\caption{Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}], +\caption{ Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}], TIP3P [Ref. \citen{Jorgensen98b}], SPC/E [Ref. \citen{Clancy94}], SSD without Reaction Field, SSD, and experiment [Ref. \citen{CRC80}]. The arrows indicate the change in densities observed when turning off the @@ -398,8 +398,8 @@ cutoff radius of 12.0 \AA\ to complement the previous dependent on the cutoff radius used both with and without the use of reaction field.\cite{Berendsen98} In order to address the possible effect of cutoff radius, simulations were performed with a dipolar -cutoff radius of 12.0 \AA\ to complement the previous SSD -simulations, all performed with a cutoff of 9.0 \AA. All of the +cutoff radius of 12.0~\AA\ to complement the previous SSD +simulations, all performed with a cutoff of 9.0~\AA. All of the resulting densities overlapped within error and showed no significant trend toward lower or higher densities as a function of cutoff radius, for simulations both with and without reaction field. These results @@ -412,7 +412,7 @@ apparent at temperatures greater than 300 K. Lower tha scaling of SSD relative to other common models at any given temperature. SSD assumes a lower density than any of the other listed models at the same pressure, behavior which is especially -apparent at temperatures greater than 300 K. Lower than expected +apparent at temperatures greater than 300~K. Lower than expected densities have been observed for other systems using a reaction field for long-range electrostatic interactions, so the most likely reason for the significantly lower densities seen in these simulations is the @@ -425,11 +425,11 @@ however, a shift in the density maximum location, down freezing point of liquid water. The shape of the curve is similar to the curve produced from SSD simulations using reaction field, specifically the rapidly decreasing densities at higher temperatures; -however, a shift in the density maximum location, down to 245 K, is +however, a shift in the density maximum location, down to 245~K, is observed. This is a more accurate comparison to the other listed water models, in that no long range corrections were applied in those simulations.\cite{Clancy94,Jorgensen98b} However, even without the -reaction field, the density around 300 K is still significantly lower +reaction field, the density around 300~K is still significantly lower than experiment and comparable water models. This anomalous behavior was what lead Tan {\it et al.} to recently reparameterize SSD.\cite{Ichiye03} Throughout the remainder of the paper our @@ -444,8 +444,8 @@ underwent 50 ps of temperature scaling and 50 ps of co constant energy (NVE) simulations were performed at the average density obtained by the NPT simulations at an identical target temperature. Simulations started with randomized velocities and -underwent 50 ps of temperature scaling and 50 ps of constant energy -equilibration before a 200 ps data collection run. Diffusion constants +underwent 50~ps of temperature scaling and 50~ps of constant energy +equilibration before a 200~ps data collection run. Diffusion constants were calculated via linear fits to the long-time behavior of the mean-square displacement as a function of time. The averaged results from five sets of NVE simulations are displayed in figure @@ -456,7 +456,7 @@ results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} \begin{center} \epsfxsize=6in \epsfbox{betterDiffuse.epsi} -\caption{Average self-diffusion constant as a function of temperature for +\caption{ Average self-diffusion constant as a function of temperature for SSD, SPC/E [Ref. \citen{Clancy94}], and TIP5P [Ref. \citen{Jorgensen01}] compared with experimental data [Refs. \citen{Gillen72} and \citen{Holz00}]. Of the three water models @@ -471,7 +471,7 @@ reproducing values similar to experiment around 290 K; strengths of the SSD model. Of the three models shown, the SSD model has the most accurate depiction of self-diffusion in both the supercooled and liquid regimes. SPC/E does a respectable job by -reproducing values similar to experiment around 290 K; however, it +reproducing values similar to experiment around 290~K; however, it deviates at both higher and lower temperatures, failing to predict the correct thermal trend. TIP5P and SSD both start off low at colder temperatures and tend to diffuse too rapidly at higher temperatures. @@ -490,26 +490,17 @@ at 245 K, indicating a first order phase transition fo capacity (C$_\text{p}$) was monitored to locate the melting transition in each of the simulations. In the melting simulations of the 1024 particle ice $I_h$ simulations, a large spike in C$_\text{p}$ occurs -at 245 K, indicating a first order phase transition for the melting of +at 245~K, indicating a first order phase transition for the melting of these ice crystals. When the reaction field is turned off, the melting -transition occurs at 235 K. These melting transitions are +transition occurs at 235~K. These melting transitions are considerably lower than the experimental value. - -\begin{figure} -\begin{center} -\epsfxsize=6in -\epsfbox{corrDiag.eps} -\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.} -\label{corrAngle} -\end{center} -\end{figure} \begin{figure} \begin{center} \epsfxsize=6in \epsfbox{fullContours.eps} -\caption{Contour plots of 2D angular pair correlation functions for -512 SSD molecules at 100 K (A \& B) and 300 K (C \& D). Dark areas +\caption{ Contour plots of 2D angular pair correlation functions for +512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas signify regions of enhanced density while light areas signify depletion relative to the bulk density. White areas have pair correlation values below 0.5 and black areas have values above 1.5.} @@ -517,6 +508,15 @@ Additional analysis of the melting process was perform \end{center} \end{figure} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{corrDiag.eps} +\caption{ An illustration of angles involved in the correlations observed in Fig. \ref{contour}.} +\label{corrAngle} +\end{center} +\end{figure} + Additional analysis of the melting process was performed using two-dimensional structure and dipole angle correlations. Expressions for these correlations are as follows: @@ -555,25 +555,25 @@ this split character alters to show the leading 4 \AA\ $g_\mathrm{OO}(r)$.\cite{Ichiye96} At low temperatures, the second solvation shell peak appears to have two distinct components that blend together to form one observable peak. At higher temperatures, -this split character alters to show the leading 4 \AA\ peak dominated +this split character alters to show the leading 4~\AA\ peak dominated by equatorial anti-parallel dipole orientations. There is also a tightly bunched group of axially arranged dipoles that most likely consist of the smaller fraction of aligned dipole pairs. The trailing -component of the split peak at 5 \AA\ is dominated by aligned dipoles +component of the split peak at 5~\AA\ is dominated by aligned dipoles that assume hydrogen bond arrangements similar to those seen in the first solvation shell. This evidence indicates that the dipole pair interaction begins to dominate outside of the range of the dipolar repulsion term. The energetically favorable dipole arrangements populate the region immediately outside this repulsion region (around -4 \AA), while arrangements that seek to satisfy both the sticky and +4~\AA), while arrangements that seek to satisfy both the sticky and dipole forces locate themselves just beyond this initial buildup -(around 5 \AA). +(around 5~\AA). From these findings, the split second peak is primarily the product of the dipolar repulsion term of the sticky potential. In fact, the inner peak can be pushed out and merged with the outer split peak just by extending the switching function ($s^\prime(r_{ij})$) from its normal -4.0 \AA\ cutoff to values of 4.5 or even 5 \AA. This type of +4.0~\AA\ cutoff to values of 4.5 or even 5~\AA. This type of correction is not recommended for improving the liquid structure, since the second solvation shell would still be shifted too far out. In addition, this would have an even more detrimental effect on @@ -608,11 +608,11 @@ the liquid structure in simulations without a long-ran \begin{table} \begin{center} -\caption{Parameters for the original and adjusted models} +\caption{ Parameters for the original and adjusted models} \begin{tabular}{ l c c c c } \hline \\[-3mm] \ \ \ Parameters\ \ \ & \ \ \ SSD [Ref. \citen{Ichiye96}] \ \ \ -& \ SSD1 [Ref. \citen{Ichiye03}]\ \ & \ SSD/E\ \ & \ SSD/RF \\ +& \ SSD1 [Ref. \citen{Ichiye03}]\ \ & \ SSD/E\ \ & \ \ SSD/RF \\ \hline \\[-3mm] \ \ \ $\sigma$ (\AA) & 3.051 & 3.016 & 3.035 & 3.019\\ \ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\ @@ -632,7 +632,7 @@ the liquid structure in simulations without a long-ran \begin{center} \epsfxsize=5in \epsfbox{GofRCompare.epsi} -\caption{Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with +\caption{ Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with SSD/E and SSD1 without reaction field (top), as well as SSD/RF and SSD1 with reaction field turned on (bottom). The insets show the respective first peaks in detail. Note @@ -646,7 +646,7 @@ peak of SSD/E and SSD/RF.} \begin{center} \epsfxsize=6in \epsfbox{dualsticky_bw.eps} -\caption{Positive and negative isosurfaces of the sticky potential for +\caption{ Positive and negative 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 component.} @@ -688,27 +688,26 @@ persistence of full dipolar character below the previo particles feel the pull of the ``hydrogen bonds''. Aside from improving the shape of the first peak in the g(\emph{r}), this modification improves the densities considerably by allowing the -persistence of full dipolar character below the previous 4.0 \AA\ +persistence of full dipolar character below the previous 4.0~\AA\ cutoff. While adjusting the location and shape of the first peak of $g(r)$ improves the densities, these changes alone are insufficient to bring the system densities up to the values observed experimentally. To further increase the densities, the dipole moments were increased in -both of our adjusted models. Since SSD is a dipole based model, -the structure and transport are very sensitive to changes in the -dipole moment. The original SSD simply used the dipole moment -calculated from the TIP3P water model, which at 2.35 D is -significantly greater than the experimental gas phase value of 1.84 -D. The larger dipole moment is a more realistic value and improves the -dielectric properties of the fluid. Both theoretical and experimental -measurements indicate a liquid phase dipole moment ranging from 2.4 D -to values as high as 3.11 D, providing a substantial range of -reasonable values for a dipole -moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately -increasing the dipole moments to 2.42 and 2.48 D for SSD/E and -SSD/RF, respectively, leads to significant changes in the -density and transport of the water models. +both of our adjusted models. Since SSD is a dipole based model, the +structure and transport are very sensitive to changes in the dipole +moment. The original SSD simply used the dipole moment calculated from +the TIP3P water model, which at 2.35~D is significantly greater than +the experimental gas phase value of 1.84~D. The larger dipole moment +is a more realistic value and improves the dielectric properties of +the fluid. Both theoretical and experimental measurements indicate a +liquid phase dipole moment ranging from 2.4~D to values as high as +3.11~D, providing a substantial range of reasonable values for a +dipole moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately +increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF, +respectively, leads to significant changes in the density and +transport of the water models. In order to demonstrate the benefits of these reparameterizations, a series of NPT and NVE simulations were performed to probe the density @@ -719,7 +718,7 @@ simulation was equilibrated for 100 ps before a 200 ps results are obtained from five separate simulations of 1024 particle systems, and the melting sequences were started from different ice $I_h$ crystals constructed as described previously. Each NPT -simulation was equilibrated for 100 ps before a 200 ps data collection +simulation was equilibrated for 100~ps before a 200~ps data collection run at each temperature step, and the final configuration from the previous temperature simulation was used as a starting point. All NVE simulations had the same thermalization, equilibration, and data @@ -729,7 +728,7 @@ collection times as stated previously. \begin{center} \epsfxsize=6in \epsfbox{ssdeDense.epsi} -\caption{Comparison of densities calculated with SSD/E to +\caption{ Comparison of densities calculated with SSD/E to SSD1 without a reaction field, TIP3P [Ref. \citen{Jorgensen98b}], TIP5P [Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}] and experiment [Ref. \citen{CRC80}]. The window shows a expansion around @@ -759,14 +758,14 @@ melting transition for SSD/E was shown to occur at 235 strengthening of the dipolar character. However, this increasing disorder in the SSD/E model has little effect on the melting transition. By monitoring $C_p$ throughout these simulations, the -melting transition for SSD/E was shown to occur at 235 K. The +melting transition for SSD/E was shown to occur at 235~K. The same transition temperature observed with SSD and SSD1. \begin{figure} \begin{center} \epsfxsize=6in \epsfbox{ssdrfDense.epsi} -\caption{Comparison of densities calculated with SSD/RF to +\caption{ Comparison of densities calculated with SSD/RF to SSD1 with a reaction field, TIP3P [Ref. \citen{Jorgensen98b}], TIP5P [Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}], and experiment [Ref. \citen{CRC80}]. The inset shows the necessity of @@ -790,17 +789,17 @@ which observed at 245 K for SSD/RF, is identical to SS further emphasize the importance of reparameterization in order to model the density properly under different simulation conditions. Again, these changes have only a minor effect on the melting point, -which observed at 245 K for SSD/RF, is identical to SSD and only 5 K +which observed at 245~K for SSD/RF, is identical to SSD and only 5~K lower than SSD1 with a reaction field. Additionally, the difference in density maxima is not as extreme, with SSD/RF showing a density -maximum at 255 K, fairly close to the density maxima of 260 K and 265 -K, shown by SSD and SSD1 respectively. +maximum at 255~K, fairly close to the density maxima of 260~K and +265~K, shown by SSD and SSD1 respectively. \begin{figure} \begin{center} \epsfxsize=6in \epsfbox{ssdeDiffuse.epsi} -\caption{The diffusion constants calculated from SSD/E and +\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 @@ -823,8 +822,8 @@ K). Both models follow the shape of the experimental c 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 +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. Both SSD1 and SSD/E show this deviation in particle mobility, but this trend has different @@ -841,7 +840,7 @@ conditions. \begin{center} \epsfxsize=6in \epsfbox{ssdrfDiffuse.epsi} -\caption{The diffusion constants calculated from SSD/RF and +\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{Holz00}]. The NVE calculations were performed at the average densities observed in @@ -860,7 +859,7 @@ temperatures greater than 330 K. As stated above, thi 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 +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 because the calculated densities are closer to the experimental @@ -871,27 +870,25 @@ reparameterization when using an altered long-range co \begin{minipage}{\linewidth} \renewcommand{\thefootnote}{\thempfootnote} \begin{center} -\caption{Properties of the single-point water models compared with -experimental data at ambient conditions} +\caption{ Properties of the single-point water models compared with +experimental data at ambient conditions. Deviations of the of the +averages are given in parentheses.} \begin{tabular}{ l c c c c c } \hline \\[-3mm] -\ \ \ \ \ \ & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ -\ & \ SSD/RF \ \ \ & \ Expt. \\ +\ \ \ \ \ \ & \ \ \ 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\cite{Mills73} \\ -\ \ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 & +\ \ $\rho$ (g/cm$^3$) & 0.999(0.001) & 0.996(0.001) & 0.972(0.002) & 0.997(0.001) & 0.997 \\ +\ \ $C_p$ (cal/mol K) & 28.80(0.11) & 25.45(0.09) & 28.28(0.06) & 23.83(0.16) & 17.98 \\ +\ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78(0.7) & 2.51(0.18) & 2.00(0.17) & 2.32(0.06) & 2.299\cite{Mills73} \\ +\ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 & 4.7\footnote{Calculated by integrating $g_{\text{OO}}(r)$ in Ref. \citen{Head-Gordon00_1}} \\ -\ \ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 & +\ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 & 3.5\footnote{Calculated by integrating $g_{\text{OH}}(r)$ in Ref. \citen{Soper86}} \\ -\ \ \ $\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\footnote{Calculated for 298 K from data in Ref. \citen{Eisenberg69}} \\ -\ \ \ $\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\footnote{Calculated for 298 K from data in +\ \ $\tau_1$ (ps) & 10.9(0.6) & 7.3(0.4) & 7.5(0.7) & 7.2(0.4) & 5.7\footnote{Calculated for 298 K from data in Ref. \citen{Eisenberg69}} \\ +\ \ $\tau_2$ (ps) & 4.7(0.4) & 3.1(0.2) & 3.5(0.3) & 3.2(0.2) & 2.3\footnote{Calculated for 298 K from data in Ref. \citen{Krynicki66}} \end{tabular} \label{liquidproperties} @@ -943,11 +940,11 @@ dependence of $\tau_2$, it is necessary to recalculate the commonly cited value of 1.9 ps for $\tau_2$ was determined from the NMR data in Ref. \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 +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 Ref. \citen{Krynicki66}. Similarly, $\tau_1$ was -recomputed for 298 K from the data in Ref. \citen{Eisenberg69}. +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 @@ -962,7 +959,7 @@ can be attributed to the use of the Ewald sum.\cite{Ic \begin{center} \epsfxsize=6in \epsfbox{icei_bw.eps} -\caption{The most stable crystal structure assumed by the SSD family +\caption{ The most stable crystal structure assumed by the SSD family of water models. We refer to this structure as Ice-{\it i} to indicate its origins in computer simulation. This image was taken of the (001) face of the crystal.} @@ -973,9 +970,9 @@ water. After melting at 235 K, two of five systems un While performing a series of melting simulations on an early iteration of SSD/E not discussed in this paper, we observed recrystallization into a novel structure not previously known for -water. After melting at 235 K, two of five systems underwent -crystallization events near 245 K. The two systems remained -crystalline up to 320 and 330 K, respectively. The crystal exhibits +water. After melting at 235~K, two of five systems underwent +crystallization events near 245~K. The two systems remained +crystalline up to 320 and 330~K, respectively. The crystal exhibits an expanded zeolite-like structure that does not correspond to any known form of ice. This appears to be an artifact of the point dipolar models, so to distinguish it from the experimentally observed @@ -1009,12 +1006,12 @@ family, Ice-{\it i} had the lowest calculated enthalpy \begin{table} \begin{center} -\caption{Enthalpies of Formation (in kcal / mol) of the three crystal +\caption{ Enthalpies of Formation (in kcal / mol) of the three crystal structures (at 1 K) exhibited by the SSD family of water models} \begin{tabular}{ l c c c } \hline \\[-3mm] -\ \ \ Water Model \ \ \ & \ \ \ Ice-$I_h$ \ \ \ & \ Ice-$I_c$\ \ & \ -Ice-{\it i} \\ +\ \ \ Water Model \ \ \ & \ \ \ Ice-$I_h$ \ \ \ & \ \ \ Ice-$I_c$ \ \ \ & +\ \ \ \ Ice-{\it i} \\ \hline \\[-3mm] \ \ \ SSD/E & -72.444 & -72.450 & -73.748 \\ \ \ \ SSD/RF & -73.093 & -73.075 & -74.180 \\ @@ -1026,7 +1023,7 @@ performed with ice-{\it i} as the initial configuratio \end{table} In addition to these energetic comparisons, melting simulations were -performed with ice-{\it i} as the initial configuration using SSD/E, +performed with Ice-{\it i} as the initial configuration using SSD/E, SSD/RF, and SSD1 both with and without a reaction field. The melting transitions for both SSD/E and SSD1 without reaction field occurred at temperature in excess of 375~K. SSD/RF and SSD1 with a reaction field @@ -1084,8 +1081,7 @@ DMR-0079647. \newpage \bibliographystyle{jcp} -\bibliography{nptSSD} +\bibliography{nptSSD} -%\pagebreak \end{document}