--- trunk/ssdePaper/nptSSD.tex 2003/11/12 13:37:15 861 +++ trunk/ssdePaper/nptSSD.tex 2003/11/13 15:55:20 863 @@ -1,38 +1,36 @@ %\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4} -\documentclass[prb,aps,times,tabularx,preprint]{revtex4} +\documentclass[11pt]{article} +\usepackage{endfloat} \usepackage{amsmath} +\usepackage{epsf} +\usepackage{berkeley} +\usepackage{setspace} +\usepackage{tabularx} \usepackage{graphicx} - -%\usepackage{endfloat} +\usepackage[ref]{overcite} %\usepackage{berkeley} -%\usepackage{epsf} -%\usepackage[ref]{overcite} -%\usepackage{setspace} -%\usepackage{tabularx} -%\usepackage{graphicx} %\usepackage{curves} -%\usepackage{amsmath} -%\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 +\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 \begin{document} \title{On the temperature dependent properties of the soft sticky dipole (SSD) and related single point water models} -\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote} -\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}} - -\address{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 + \begin{abstract} NVE and NPT molecular dynamics simulations were performed in order to investigate the density maximum and temperature dependent transport @@ -50,7 +48,7 @@ maintain or improve upon the structural and transport maintain or improve upon the structural and transport properties. \end{abstract} -\maketitle +\newpage %\narrowtext @@ -61,27 +59,26 @@ One of the most important tasks in simulations of bioc \section{Introduction} -One of the most important tasks in simulations of biochemical systems -is the proper depiction of water and water solvation. In fact, the -bulk of the calculations performed in solvated simulations are of +One of the most important tasks in the simulation of biochemical +systems is the proper depiction of water and water solvation. In fact, +the bulk of the calculations performed in solvated simulations are of interactions with or between solvent molecules. Thus, the outcomes of these types of simulations are highly dependent on the physical -properties of water, both as individual molecules and in -groups/bulk. Due to the fact that explicit solvent accounts for a -massive portion of the calculations, it necessary to simplify the -solvent to some extent in order to complete simulations in a -reasonable amount of time. In the case of simulating water in -bio-molecular studies, the balance between accurate properties and -computational efficiency is especially delicate, and it has resulted -in a variety of different water -models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these models -get specific properties correct or better than their predecessors, but -this is often at a cost of some other properties or of computer -time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P succeeds -in improving the structural and transport properties over its -predecessors, yet this comes at a greater than 50\% increase in +properties of water, both as individual molecules and in clusters or +bulk. Due to the fact that explicit solvent accounts for a massive +portion of the calculations, it necessary to simplify the solvent to +some extent in order to complete simulations in a reasonable amount of +time. In the case of simulating water in biomolecular studies, the +balance between accurate properties and computational efficiency is +especially delicate, and it has resulted in a variety of different +water models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these +models predict specific properties more accurately than their +predecessors, but often at the cost of other properties or of computer +time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P improves +upon the structural and transport properties of water relative to the +previous TIP models, yet this comes at a greater than 50\% increase in computational cost.\cite{Jorgensen01,Jorgensen00} One recently -developed model that succeeds in both retaining accuracy of system +developed model that succeeds in both retaining the accuracy of system properties and simplifying calculations to increase computational efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96} @@ -100,11 +97,11 @@ where the $\mathbf{r}_{ij}$ is the position vector bet (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j), \end{equation} where the $\mathbf{r}_{ij}$ is the position vector between molecules -\emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and +\emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the orientations of the respective molecules. The Lennard-Jones, dipole, and sticky parts of the potential are giving by the following -equations, +equations: \begin{equation} u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right], \end{equation} @@ -112,20 +109,15 @@ u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ , \end{equation} \begin{equation} -\begin{split} u_{ij}^{sp} -(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) -&= -\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\ -& \quad \ + -s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ , -\end{split} +(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = +\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) + s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ , \end{equation} where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D, -$\nu_0$ scales the strength of the overall sticky potential, $s$ and -$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$ -functions take the following forms, +$\nu_0$ scales the strength of the overall sticky potential, and $s$ +and $s^\prime$ are cubic switching functions. The $w$ and $w^\prime$ +functions take the following forms: \begin{equation} w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij}, \end{equation} @@ -142,9 +134,9 @@ simulations using this model, Ichiye \emph{et al.} rep Being that this is a one-site point dipole model, the actual force calculations are simplified significantly. In the original Monte Carlo -simulations using this model, Ichiye \emph{et al.} reported a -calculation speed up of up to an order of magnitude over other -comparable models, while maintaining the structural behavior of +simulations using this model, Ichiye \emph{et al.} reported an +increase in calculation efficiency of up to an order of magnitude over +other comparable models, while maintaining the structural behavior of water.\cite{Ichiye96} In the original molecular dynamics studies, it was shown that SSD improves on the prediction of many of water's dynamical properties over TIP3P and SPC/E.\cite{Ichiye99} This @@ -156,28 +148,27 @@ computational burdens with their respective ideal $N^\ has been parameterized for use with the Ewald Sum technique for the handling of long-ranged interactions. When studying very large systems, the Ewald summation and even particle-mesh Ewald become -computational burdens with their respective ideal $N^\frac{3}{2}$ and +computational burdens, with their respective ideal $N^\frac{3}{2}$ and $N\log N$ calculation scaling orders for $N$ particles.\cite{Darden99} In applying this water model in these types of systems, it would be useful to know its properties and behavior with the more computationally efficient reaction field (RF) technique, and even with -a cutoff that lacks any form of long range correction. This study +a cutoff that lacks any form of long-range correction. This study addresses these issues by looking at the structural and transport -behavior of SSD over a variety of temperatures, with the purpose of -utilizing the RF correction technique. Toward the end, we suggest -alterations to the parameters that result in more water-like -behavior. It should be noted that in a recent publication, some the -original investigators of the SSD water model have put forth -adjustments to the SSD water model to address abnormal density -behavior (also observed here), calling the corrected model -SSD1.\cite{Ichiye03} This study will make comparisons with this new -model's behavior with the goal of improving upon the depiction of -water under conditions without the Ewald Sum. +behavior of SSD over a variety of temperatures with the purpose of +utilizing the RF correction technique. We then suggest alterations to +the parameters that result in more water-like behavior. It should be +noted that in a recent publication, some of the original investigators of +the SSD water model have put forth adjustments to the SSD water model +to address abnormal density behavior (also observed here), calling the +corrected model SSD1.\cite{Ichiye03} This study will make comparisons +with SSD1's behavior with the goal of improving upon the +depiction of water under conditions without the Ewald Sum. \section{Methods} -As stated previously, in this study the long-range dipole-dipole -interactions were accounted for using the reaction field method. The +As stated previously, the long-range dipole-dipole interactions were +accounted for in this study by using the reaction field method. The magnitude of the reaction field acting on dipole \emph{i} is given by \begin{equation} \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1} @@ -198,10 +189,10 @@ large-scale system, the computational cost benefit of in the length of the cutoff radius.\cite{Berendsen98} This variable behavior makes reaction field a less attractive method than other methods, like the Ewald summation; however, for the simulation of -large-scale system, the computational cost benefit of reaction field +large-scale systems, the computational cost benefit of reaction field is dramatic. To address some of the dynamical property alterations due to the use of reaction field, simulations were also performed without -a surrounding dielectric and suggestions are proposed on how to make +a surrounding dielectric and suggestions are presented on how to make SSD more accurate both with and without a reaction field. Simulations were performed in both the isobaric-isothermal and @@ -214,7 +205,7 @@ al.}.\cite{Dullweber1997} The reason for this integrat Integration of the equations of motion was carried out using the symplectic splitting method proposed by Dullweber \emph{et -al.}.\cite{Dullweber1997} The reason for this integrator selection +al.}\cite{Dullweber1997} The reason for this integrator selection deals with poor energy conservation of rigid body systems using quaternions. While quaternions work well for orientational motion in alternate ensembles, the microcanonical ensemble has a constant energy @@ -227,7 +218,7 @@ feasible a option, being that the rotation matrix for The key difference in the integration method proposed by Dullweber \emph{et al.} is that the entire rotation matrix is propagated from one time step to the next. In the past, this would not have been as -feasible a option, being that the rotation matrix for a single body is +feasible an option, being that the rotation matrix for a single body is nine elements long as opposed to 3 or 4 elements for Euler angles and quaternions respectively. System memory has become much less of an issue in recent times, and this has resulted in substantial benefits @@ -238,33 +229,33 @@ both linear and angular motion of rigid bodies. In the purposes relieves this burden. The symplectic splitting method allows for Verlet style integration of -both linear and angular motion of rigid bodies. In the integration +both linear and angular 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 end up being more costly computationally 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. This cost is more than justified when comparing the -energy conservation of the two methods as illustrated in figure -\ref{timestep}. +matrix rotations are more costly computationally 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. This cost is +more than justified when comparing the energy conservation of the two +methods as illustrated in figure \ref{timestep}. \begin{figure} -\includegraphics[width=61mm, angle=-90]{timeStep.epsi} +\begin{center} +\epsfxsize=6in +\epsfbox{timeStep.epsi} \caption{Energy conservation using quaternion based integration versus the symplectic step method proposed by Dullweber \emph{et al.} with -increasing time step. For each time step, the dotted line is total -energy using the symplectic step integrator, and the solid line comes -from the quaternion integrator. The larger time step plots are shifted -up from the true energy baseline for clarity.} +increasing time step. The larger time step plots are shifted up 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 for -handling rotational motion. At time steps of 0.1 and 0.5 fs, both +same configuration, and the only difference was the method used to +handle rotational motion. At time steps of 0.1 and 0.5 fs, both methods for propagating particle rotation 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 @@ -273,7 +264,7 @@ Energy drift in these SSD particle simulations was unn conservation, one can take considerably longer time steps, leading to an overall reduction in computation time. -Energy drift in these SSD particle simulations was unnoticeable for +Energy drift in the symplectic step simulations was unnoticeable for time steps up to three femtoseconds. A slight energy drift on the order of 0.012 kcal/mol per nanosecond was observed at a time step of four femtoseconds, and as expected, this drift increases dramatically @@ -282,21 +273,21 @@ starting points for all the simulations. The $I_h$ cry constant pressure simulations as well. Ice crystals in both the $I_h$ and $I_c$ lattices were generated as -starting points for all the simulations. The $I_h$ crystals were -formed by first arranging the center of masses 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 a -rectangular shape with a edge length ratio of approximately +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 a rectangular 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 down to 25 K. The particles -were then allowed translate with fixed orientations at a constant +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 +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 isotropic volume fluctuations were performed under constant pressure, @@ -306,20 +297,20 @@ constant pressure and temperature dynamics. By perform \section{Results and discussion} Melting studies were performed on the randomized ice crystals using -constant pressure and temperature dynamics. By performing melting -simulations, the melting transition can be determined by monitoring -the heat capacity, in addition to determining the density maximum - -provided that the density maximum occurs in the liquid and not the -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 -sufficient in that the system volume fluctuations dampened out in all -but the very cold simulations (below 225 K). +constant pressure and temperature dynamics. During melting +simulations, the melting transition and the density maximum can both +be observed, provided that the density maximum occurs in the liquid +and not the 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 sufficient in that the system volume +fluctuations dampened out in all but the very cold simulations (below +225 K). \subsection{Density Behavior} Initial simulations focused on the original SSD water model, and an @@ -330,24 +321,27 @@ configuration; and though not pictured, the simulation of the average value at 260 K makes a good argument for the actual density maximum residing at this midpoint value. Figure \ref{dense1} was constructed using ice $I_h$ crystals for the initial -configuration; and though not pictured, the simulations starting from -ice $I_c$ crystal configurations showed similar results, with a +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). In addition, the $I_c$ crystals are more fragile than the $I_h$ -crystals, leading them to deform into a dense glassy state at lower +crystals, leading to deformation into a dense glassy state at lower temperatures. This resulted in an overall low temperature density -maximum at 200 K, but they still retained a common liquid state -density maximum with the $I_h$ simulations. +maximum at 200 K, while still retaining a liquid state density maximum +in common with the $I_h$ simulations. \begin{figure} -\includegraphics[width=65mm,angle=-90]{dense2.eps} -\caption{Density versus temperature for TIP4P\cite{Jorgensen98b}, - TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction - Field, SSD, and Experiment\cite{CRC80}. The arrows indicate the +\begin{center} +\epsfxsize=6in +\epsfbox{denseSSD.eps} +\caption{Density versus temperature for TIP4P,\cite{Jorgensen98b} + TIP3P,\cite{Jorgensen98b} SPC/E,\cite{Clancy94} SSD without Reaction + Field, SSD, and experiment.\cite{CRC80} The arrows indicate the change in densities observed when turning off the reaction field. The the lower than expected densities for the SSD model were what prompted the original reparameterization.\cite{Ichiye03}} \label{dense1} +\end{center} \end{figure} The density maximum for SSD actually compares quite favorably to other @@ -366,35 +360,36 @@ the resulting densities overlapped within error and sh address the possible affect of cutoff radius, simulations were performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the previous SSD simulations, all performed with a cutoff of 9.0 \AA. All -the resulting densities overlapped within error and showed no -significant trend in lower or higher densities as a function of cutoff -radius, both for simulations with and without reaction field. These -results indicate that there is no major benefit in choosing a longer -cutoff radius in simulations using SSD. This is comforting in that the -use of a longer cutoff radius results in significant increases in the -time required to obtain a single trajectory. +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 indicate that there is no major benefit in +choosing a longer cutoff radius in simulations using SSD. This is +advantageous in that the use of a longer cutoff radius results in +significant increases in the time required to obtain a single +trajectory. -The most important thing to recognize in figure \ref{dense1} is the -density scaling of SSD relative to other common models at any given +The key feature to recognize in figure \ref{dense1} is the density +scaling of SSD relative to other common models at any given temperature. Note that the SSD model 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 densities have been observed for other systems with the use -of a reaction field for long-range electrostatic interactions, so the -most likely reason for these significantly lower densities in these +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 presence of the reaction field.\cite{Berendsen98,Nezbeda02} In order to test the effect of the reaction field on the density of the systems, the simulations were repeated without a reaction field present. The results of these simulations are also displayed in figure \ref{dense1}. Without -reaction field, these densities increase considerably to more +reaction field, the densities increase considerably to more experimentally reasonable values, especially around the 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 observed. This is -probably a more accurate comparison to the other listed water models, -in that no long range corrections were applied in those +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 a reaction field, the density around 300 K is still significantly lower than experiment and comparable water models. This anomalous behavior @@ -404,83 +399,94 @@ of the particles and how they change when altering the \subsection{Transport Behavior} Of importance in these types of studies are the transport properties -of the particles and how they change when altering the environmental -conditions. In order to probe transport, constant energy simulations -were performed about the average density uncovered by the constant -pressure simulations. Simulations started with randomized velocities -and underwent 50 ps of temperature scaling and 50 ps of constant -energy equilibration before obtaining a 200 ps trajectory. Diffusion -constants were calculated via root-mean square deviation analysis. The -averaged results from 5 sets of these NVE simulations is displayed in -figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P -results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} +of the particles and their change in responce to altering +environmental conditions. In order to probe transport, constant energy +simulations were performed about the average density uncovered by the +constant pressure simulations. Simulations started with randomized +velocities and underwent 50 ps of temperature scaling and 50 ps of +constant energy equilibration before obtaining a 200 ps +trajectory. Diffusion constants were calculated via root-mean square +deviation analysis. 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} \begin{figure} -\includegraphics[width=65mm, angle=-90]{betterDiffuse.epsi} +\begin{center} +\epsfxsize=6in +\epsfbox{betterDiffuse.epsi} \caption{Average diffusion coefficient over increasing temperature for -SSD, SPC/E\cite{Clancy94}, TIP5P\cite{Jorgensen01}, and Experimental -data from Gillen \emph{et al.}\cite{Gillen72}, and from -Mills\cite{Mills73}.} +SSD, SPC/E,\cite{Clancy94} TIP5P,\cite{Jorgensen01} and Experimental +data.\cite{Gillen72,Mills73} 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 higher temperatures.} \label{diffuse} +\end{center} \end{figure} The observed values for the diffusion constant point out one of the strengths of the SSD model. Of the three experimental models shown, the SSD model has the most accurate depiction of the diffusion trend -seen in experiment in both the supercooled and normal regimes. SPC/E -does a respectable job by getting similar values as SSD and experiment -around 290 K; however, it deviates at both higher and lower -temperatures, failing to predict the experimental trend. TIP5P and SSD -both start off low at the colder temperatures and tend to diffuse too -rapidly at the higher temperatures. This type of trend at the higher -temperatures is not surprising in that the densities of both TIP5P and -SSD are lower than experimental water at temperatures higher than room -temperature. When calculating the diffusion coefficients for SSD at +seen in experiment in both the supercooled and liquid temperature +regimes. SPC/E does a respectable job by producing values similar to +SSD and experiment around 290 K; however, it deviates at both higher +and lower temperatures, failing to predict the experimental +trend. TIP5P and SSD both start off low at colder temperatures and +tend to diffuse too rapidly at higher temperatures. This trend at +higher temperatures is not surprising in that the densities of both +TIP5P and SSD are lower than experimental water at these higher +temperatures. When calculating the diffusion coefficients for SSD at experimental densities, the resulting values fall more in line with experiment at these temperatures, albeit not at standard pressure. \subsection{Structural Changes and Characterization} By starting the simulations from the crystalline state, the melting transition and the ice structure can be studied along with the liquid -phase behavior beyond the melting point. To locate the melting -transition, the constant pressure heat capacity (C$_\text{p}$) was -monitored 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 these ice crystals. When the reaction field is turned off, -the melting transition occurs at 235 K. These melting transitions are -considerably lower than the experimental value, but this is not -surprising when considering the simplicity of the SSD model. +phase behavior beyond the melting point. The constant pressure heat +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 +these ice crystals. When the reaction field is turned off, the melting +transition occurs at 235 K. These melting transitions are +considerably lower than the experimental value, but this is not a +surprise considering the simplicity of the SSD model. -\begin{figure} -\includegraphics[width=85mm]{fullContours.eps} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{corrDiag.eps} +\caption{Two dimensional illustration of angles involved in the +correlations observed in figure \ref{contour}.} +\label{corrAngle} +\end{center} +\end{figure} + +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{fullContours.eps} \caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at 100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for clarity: dark areas signify peaks while light areas signify depressions. White areas have g(\emph{r}) values below 0.5 and black areas have values above 1.5.} \label{contour} +\end{center} \end{figure} -\begin{figure} -\includegraphics[width=45mm]{corrDiag.eps} -\caption{Two dimensional illustration of the angles involved in the -correlations observed in figure \ref{contour}.} -\label{corrAngle} -\end{figure} - Additional analysis of the melting phase-transition process was performed by using two-dimensional structure and dipole angle correlations. Expressions for these correlations are as follows: -\begin{multline} -g_{\text{AB}}(r,\cos\theta) = \\ -\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ , -\end{multline} -\begin{multline} -g_{\text{AB}}(r,\cos\omega) = \\ +\begin{equation} +g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ , +\end{equation} +\begin{equation} +g_{\text{AB}}(r,\cos\omega) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ , -\end{multline} +\end{equation} where $\theta$ and $\omega$ refer to the angles shown in figure \ref{corrAngle}. By binning over both distance and the cosine of the desired angle between the two dipoles, the g(\emph{r}) can be @@ -506,22 +512,21 @@ second solvation shell peak appears to have two distin This complex interplay between dipole and sticky interactions was remarked upon as a possible reason for the split second peak in the oxygen-oxygen g(\emph{r}).\cite{Ichiye96} At low temperatures, the -second solvation shell peak appears to have two distinct parts that -blend together to form one observable peak. At higher temperatures, -this split character alters to show the leading 4 \AA\ peak dominated -by equatorial anti-parallel dipole orientations, and there is tightly -bunched group of axially arranged dipoles that most likely consist of -the smaller fraction aligned dipole pairs. The trailing part of the -split peak at 5 \AA\ is dominated by aligned dipoles that range -primarily within the axial to the chief 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, with the primary -energetically favorable dipole arrangements populating the region -immediately outside this repulsion region (around 4 \AA), and -arrangements that seek to ideally satisfy both the sticky and dipole -forces locate themselves just beyond this initial buildup (around 5 -\AA). +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 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 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. Primary energetically favorable dipole +arrangements populate the region immediately outside this repulsion +region (around 4 \AA), while arrangements that seek to ideally satisfy +both the sticky and dipole forces locate themselves just beyond this +initial buildup (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 @@ -529,7 +534,7 @@ because the second solvation shell will still be shift extending the switching function cutoff ($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even 5 \AA. This type of correction is not recommended for improving the liquid structure, -because the second solvation shell will still be shifted too far +since the second solvation shell would still be shifted too far out. In addition, this would have an even more detrimental effect on the system densities, leading to a liquid with a more open structure and a density considerably lower than the normal SSD behavior shown @@ -548,40 +553,37 @@ The possible parameters for tuning include the $\sigma important properties. In this case, it would be ideal to correct the densities while maintaining the accurate transport properties. -The possible parameters for tuning include the $\sigma$ and $\epsilon$ +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} -\begin{split} u_{ij}^{sp} -(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) &= -\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\\ -& \quad \ + \frac{\nu_0^\prime}{2} -[s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)], -\end{split} +(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = +\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)], \end{equation} where $\nu_0$ scales the strength of the tetrahedral attraction and $\nu_0^\prime$ acts in an identical fashion on the dipole repulsion -term. For purposes of the reparameterization, the separation was -performed, but the final parameters were adjusted so that it is -unnecessary to separate the terms when implementing the adjusted water -potentials. The results of the reparameterizations are shown in table -\ref{params}. Note that both the tetrahedral attractive and dipolar -repulsive don't share the same lower cutoff ($r_l$) in the newly -parameterized potentials - 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 simulations without a -long-range correction). +term. The separation was performed for purposes of the +reparameterization, but the final parameters were adjusted so that it +is unnecessary 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 - 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 +simulations without a long-range correction). \begin{table} +\begin{center} \caption{Parameters for the original and adjusted models} \begin{tabular}{ l c c c c } \hline \\[-3mm] -\ \ \ Parameters\ \ \ & \ \ \ SSD$^\dagger$ \ \ \ & \ SSD1$^\ddagger$\ \ & \ SSD/E\ \ & \ SSD/RF \\ +\ \ \ Parameters\ \ \ & \ \ \ SSD\cite{Ichiye96} \ \ \ & \ SSD1\cite{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\\ @@ -592,71 +594,78 @@ long-range correction). \ \ \ $\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\\ -\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96} -\\$^\ddagger$ ref. \onlinecite{Ichiye03} \end{tabular} \label{params} +\end{center} \end{table} -\begin{figure} -\includegraphics[width=85mm]{GofRCompare.epsi} +\begin{figure} +\begin{center} +\epsfxsize=5in +\epsfbox{GofRCompare.epsi} \caption{Plots comparing experiment\cite{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. Solid Line - experiment, dashed line - SSD/E -and SSD/RF, and dotted line - SSD1 (with and without reaction field).} +first peaks in detail. Note how the changes in parameters have lowered +and broadened the first peak of SSD/E and SSD/RF.} \label{grcompare} +\end{center} \end{figure} -\begin{figure} -\includegraphics[width=85mm]{dualsticky.ps} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{dualsticky.ps} \caption{Isosurfaces of the sticky potential for SSD1 (left) and SSD/E \& SSD/RF (right). Light areas correspond to the tetrahedral attractive -part, and the darker areas correspond to the dipolar repulsive part.} +component, and darker areas correspond to the dipolar repulsive +component.} \label{isosurface} +\end{center} \end{figure} In the paper detailing the development of SSD, Liu and Ichiye placed particular emphasis on an accurate description of the first solvation -shell. This resulted in a somewhat tall and sharp first peak that -integrated to give similar coordination numbers to the experimental -data obtained by Soper and Phillips.\cite{Ichiye96,Soper86} New -experimental x-ray scattering data from the Head-Gordon lab indicates -a slightly lower and shifted first peak in the g$_\mathrm{OO}(r)$, so -adjustments to SSD were made while taking into consideration the new -experimental findings.\cite{Head-Gordon00_1} Figure \ref{grcompare} -shows the relocation of the first peak of the oxygen-oxygen -g(\emph{r}) by comparing the revised SSD model (SSD1), SSD-E, and -SSD-RF to the new experimental results. Both the modified water models -have shorter peaks that are brought in more closely to the -experimental peak (as seen in the insets of figure \ref{grcompare}). -This structural alteration was accomplished by the combined reduction -in the Lennard-Jones $\sigma$ variable and adjustment of the sticky -potential strength and cutoffs. As can be seen in table \ref{params}, -the cutoffs for the tetrahedral attractive and dipolar repulsive terms -were nearly swapped with each other. Isosurfaces of the original and -modified sticky potentials are shown in figure \cite{isosurface}. In -these isosurfaces, it is easy to see how altering the cutoffs changes -the repulsive and attractive character of the particles. With a -reduced repulsive surface (the darker region), the particles can move -closer to one another, increasing the density for the overall -system. This change in interaction cutoff also results in a more -gradual orientational motion by allowing the particles to maintain -preferred dipolar arrangements before they begin to feel the pull of -the tetrahedral restructuring. Upon moving closer together, the -dipolar repulsion term becomes active and excludes unphysical -nearest-neighbor arrangements. This compares with how SSD and SSD1 -exclude preferred dipole alignments before the particles feel the pull -of the ``hydrogen bonds''. Aside from improving the shape of the first -peak in the g(\emph{r}), this improves the densities considerably by +shell. This resulted in a somewhat tall and narrow first peak in the +g(\emph{r}) that integrated to give similar coordination numbers to +the experimental data obtained by Soper and +Phillips.\cite{Ichiye96,Soper86} New experimental x-ray scattering +data from the Head-Gordon lab indicates a slightly lower and shifted +first peak in the g$_\mathrm{OO}(r)$, so adjustments to SSD were made +while taking into consideration the new experimental +findings.\cite{Head-Gordon00_1} Figure \ref{grcompare} shows the +relocation of the first peak of the oxygen-oxygen g(\emph{r}) by +comparing the revised SSD model (SSD1), SSD-E, and SSD-RF to the new +experimental results. Both modified water models have shorter peaks +that are brought in more closely to the experimental peak (as seen in +the insets of figure \ref{grcompare}). This structural alteration was +accomplished by the combined reduction in the Lennard-Jones $\sigma$ +variable and adjustment of the sticky potential strength and +cutoffs. As can be seen in table \ref{params}, the cutoffs for the +tetrahedral attractive and dipolar repulsive terms were nearly swapped +with each other. Isosurfaces of the original and modified sticky +potentials are shown in figure \ref{isosurface}. In these isosurfaces, +it is easy to see how altering the cutoffs changes the repulsive and +attractive character of the particles. With a reduced repulsive +surface (darker region), the particles can move closer to one another, +increasing the density for the overall system. This change in +interaction cutoff also results in a more gradual orientational motion +by allowing the particles to maintain preferred dipolar arrangements +before they begin to feel the pull of the tetrahedral +restructuring. As the particles move closer together, the dipolar +repulsion term becomes active and excludes unphysical nearest-neighbor +arrangements. This compares with how SSD and SSD1 exclude preferred +dipole alignments before the 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\ cutoff. While adjusting the location and shape of the first peak of g(\emph{r}) improves the densities, these changes alone are insufficient to bring the system densities up to the values observed -experimentally. To finish bringing up the densities, the dipole -moments were increased in both the adjusted models. Being a dipole +experimentally. To further increase the densities, the dipole moments +were increased in both of the 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 @@ -664,12 +673,12 @@ to values as high as 3.11 D, so there is quite a range 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, so there is quite a range of available -values for a reasonable dipole -moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of -the dipole moments to 2.42 and 2.48 D for SSD/E and SSD/RF -respectively is moderate in this range; however, it leads to -significant changes in the density and transport of the water models. +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 @@ -677,85 +686,95 @@ results come from five separate simulations of 1024 pa to the original SSD model. This comparison involved full NPT melting sequences for both SSD/E and SSD/RF, as well as NVE transport calculations at the calculated self-consistent densities. Again, the -results come from five separate simulations of 1024 particle systems, -and the melting sequences were started from different ice $I_h$ -crystals constructed as stated earlier. Like before, each NPT +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 -run at each temperature step, and they used the final configuration -from the previous temperature simulation as a starting point. All of -the NVE simulations had the same thermalization, equilibration, and -data collection times stated earlier in this paper. +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 +collection times as stated earlier in this paper. -\begin{figure} -\includegraphics[width=62mm, angle=-90]{ssdeDense.epsi} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{ssdeDense.epsi} \caption{Comparison of densities calculated with SSD/E to SSD1 without a -reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00}, -SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The window shows a +reaction field, TIP3P,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00} +SPC/E,\cite{Clancy94} and experiment.\cite{CRC80} The window shows a expansion around 300 K with error bars included to clarify this region of interest. Note that both SSD1 and SSD/E show good agreement with experiment when the long-range correction is neglected.} \label{ssdedense} +\end{center} \end{figure} Figure \ref{ssdedense} shows the density profile for the SSD/E model -in comparison to SSD1 without a reaction field, experiment, and other -common water models. The calculated densities for both SSD/E and SSD1 -have increased significantly over the original SSD model (see figure -\ref{dense1} and are in significantly better agreement with the -experimental values. At 298 K, the density of SSD/E and SSD1 without a -long-range correction are 0.996$\pm$0.001 g/cm$^3$ and 0.999$\pm$0.001 -g/cm$^3$ respectively. These both compare well with the experimental -value of 0.997 g/cm$^3$, and they are considerably better than the SSD -value of 0.967$\pm$0.003 g/cm$^3$. The changes to the dipole moment -and sticky switching functions have improved the structuring of the -liquid (as seen in figure \ref{grcompare}, but they have shifted the -density maximum to much lower temperatures. This comes about via an -increase of the liquid disorder through the weakening of the sticky -potential and strengthening of the dipolar character. However, this -increasing disorder in the SSD/E model has little affect on the -melting transition. By monitoring C$\text{p}$ throughout these -simulations, the melting transition for SSD/E occurred at 235 K, the -same transition temperature observed with SSD and SSD1. +in comparison to SSD1 without a reaction field, other common water +models, and experimental results. The calculated densities for both +SSD/E and SSD1 have increased significantly over the original SSD +model (see figure \ref{dense1}) and are in better agreement with the +experimental values. At 298 K, the densities of SSD/E and SSD1 without +a long-range correction are 0.996$\pm$0.001 g/cm$^3$ and +0.999$\pm$0.001 g/cm$^3$ respectively. These both compare well with +the experimental value of 0.997 g/cm$^3$, and they are considerably +better than the SSD value of 0.967$\pm$0.003 g/cm$^3$. The changes to +the dipole moment and sticky switching functions have improved the +structuring of the liquid (as seen in figure \ref{grcompare}, but they +have shifted the density maximum to much lower temperatures. This +comes about via an increase in the liquid disorder through the +weakening of the sticky potential and strengthening of the dipolar +character. However, this increasing disorder in the SSD/E model has +little effect on the melting transition. By monitoring C$\text{p}$ +throughout these simulations, the melting transition for SSD/E was +shown to occur at 235 K, the same transition temperature observed with +SSD and SSD1. -\begin{figure} -\includegraphics[width=62mm, angle=-90]{ssdrfDense.epsi} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{ssdrfDense.epsi} \caption{Comparison of densities calculated with SSD/RF to SSD1 with a -reaction field, TIP3P\cite{Jorgensen98b}, TIP5P\cite{Jorgensen00}, -SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The inset shows the +reaction field, TIP3P,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00} +SPC/E,\cite{Clancy94} and experiment.\cite{CRC80} The inset shows the necessity of reparameterization when utilizing a reaction field long-ranged correction - SSD/RF provides significantly more accurate densities than SSD1 when performing room temperature simulations.} \label{ssdrfdense} +\end{center} \end{figure} -Including the reaction field long-range correction results in a more -interesting comparison. A density profile including SSD/RF and SSD1 -with an active reaction field is shown in figure \ref{ssdrfdense}. As -observed in the simulations without a reaction field, the densities of -SSD/RF and SSD1 show a dramatic increase over normal SSD (see figure -\ref{dense1}). At 298 K, SSD/RF has a density of 0.997$\pm$0.001 -g/cm$^3$, right in line with experiment and considerably better than -the SSD value of 0.941$\pm$0.001 g/cm$^3$ and the SSD1 value of -0.972$\pm$0.002 g/cm$^3$. These results further emphasize the -importance of reparameterization in order to model the density -properly under different simulation conditions. Again, these changes -don't have that profound an effect on the melting point which is -observed at 245 K for SSD/RF, identical to SSD and only 5 K lower than -SSD1 with a reaction field. However, the difference in density maxima -is not quite as extreme with SSD/RF showing a density maximum at 255 -K, fairly close to 260 and 265 K, the density maxima for SSD and SSD1 -respectively. +Including the reaction field long-range correction in the simulations +results in a more interesting comparison. A density profile including +SSD/RF and SSD1 with an active reaction field is shown in figure +\ref{ssdrfdense}. As observed in the simulations without a reaction +field, the densities of SSD/RF and SSD1 show a dramatic increase over +normal SSD (see figure \ref{dense1}). At 298 K, SSD/RF has a density +of 0.997$\pm$0.001 g/cm$^3$, directly in line with experiment and +considerably better than the SSD value of 0.941$\pm$0.001 g/cm$^3$ and +the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results 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 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. -\begin{figure} -\includegraphics[width=65mm, angle=-90]{ssdeDiffuse.epsi} +\begin{figure} +\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 are - from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{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.} + both without a reaction field, along with experimental + results.\cite{Gillen72,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.} \label{ssdediffuse} +\end{center} \end{figure} The reparameterization of the SSD water model, both for use with and @@ -766,119 +785,126 @@ water. In the upper plot, the diffusion constant for S 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. In the upper plot, the diffusion constant for SSD/E is -consistently a little faster than experiment, while SSD1 remains -slower than experiment until relatively high temperatures (greater -than 330 K). Both models follow the shape of the experimental trend -well below 300 K, but the trend leans toward diffusing too rapidly at -higher temperatures, something that is especially apparent with -SSD1. This accelerated increasing of diffusion is caused by the -rapidly decreasing system density with increasing temperature. Though -it is difficult to see in figure \ref{ssdedense}, 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 changes made to improve the liquid structure may have had an -adverse affect on the density maximum, but they improve the transport -behavior of the water model. +water. The diffusion constant for SSD/E is consistently a little +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 accelerated increasing of +diffusion is caused by the rapidly decreasing system density with +increasing temperature. Though it is difficult to see in figure +\ref{ssdedense}, 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 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. -\begin{figure} -\includegraphics[width=65mm, angle=-90]{ssdrfDiffuse.epsi} +\begin{figure} +\begin{center} +\epsfxsize=6in +\epsfbox{ssdrfDiffuse.epsi} \caption{Plots of the diffusion constants calculated from SSD/RF and SSD1, - both with an active reaction field, along with experimental results - from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. 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.} + both with an active reaction field, along with experimental + results.\cite{Gillen72,Mills73} 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.} \label{ssdrfdiffuse} +\end{center} \end{figure} In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are -compared with SSD1 with an active reaction field. Note that SSD/RF +compared to SSD1 with an active reaction field. Note that SSD/RF tracks the experimental results incredibly well, identical within -error throughout the temperature range shown and only showing a slight +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 -temperatures greater than 330 K. As was stated in the SSD/E -comparisons, this deviation away from the ideal trend is due to a -rapid decrease in density at higher temperatures. SSD/RF doesn't -suffer from this problem as much as SSD1, because the calculated -densities are more true to experiment. These results again emphasize -the importance of careful reparameterization when using an altered +temperatures greater than 330 K. As stated in relation to SSD/E, 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 value. These results again emphasize the +importance of careful reparameterization when using an altered long-range correction. \subsection{Additional Observations} \begin{figure} -\includegraphics[width=85mm]{povIce.ps} -\caption{A water lattice built from the crystal structure that SSD/E - assumed when undergoing an extremely restricted temperature NPT - simulation. This form of ice is referred to as ice 0 to emphasize its - simulation origins. This image was taken of the (001) face of the - crystal.} +\begin{center} +\epsfxsize=6in +\epsfbox{povIce.ps} +\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 \emph{i} to + emphasize its simulation origins. This image was taken of the (001) + face of the crystal.} \label{weirdice} +\end{center} \end{figure} While performing restricted temperature melting sequences of SSD/E not -discussed earlier in this paper, some interesting observations were -made. After melting at 235 K, two of five systems underwent -crystallization events near 245 K. As the heating process continued, -the two systems remained crystalline until finally melting between 320 -and 330 K. The final configurations of these two melting sequences -show an expanded zeolite-like crystal structure that does not -correspond to any known form of ice. For convenience and to help -distinguish it from the experimentally observed forms of ice, this -crystal structure will henceforth be referred to as ice-zero (ice -0). The crystallinity was extensive enough that a near ideal crystal -structure could be obtained. Figure \ref{weirdice} shows the repeating -crystal structure of a typical crystal at 5 K. Each water molecule is -hydrogen bonded to four others; however, the hydrogen bonds are flexed -rather than perfectly straight. This results in a skewed tetrahedral -geometry about the central molecule. Looking back at figure -\ref{isosurface}, it is easy to see how these flexed hydrogen bonds -are allowed in that the attractive regions are conical in shape, with -the greatest attraction in the central region. Though not ideal, these -flexed hydrogen bonds are favorable enough to stabilize an entire -crystal generated around them. In fact, the imperfect ice 0 crystals -were so stable that they melted at temperatures nearly 100 K greater -than both ice I$_c$ and I$_h$. +previously discussed, some interesting observations were made. After +melting at 235 K, two of five systems underwent crystallization events +near 245 K. As the heating process continued, the two systems remained +crystalline until finally melting between 320 and 330 K. The final +configurations of these two melting sequences show an expanded +zeolite-like crystal structure that does not correspond to any known +form of ice. For convenience, and to help distinguish it from the +experimentally observed forms of ice, this crystal structure will +henceforth be referred to as ice $\sqrt{\smash[b]{-\text{I}}}$ (ice +\emph{i}). The crystallinity was extensive enough that a near ideal +crystal structure of ice \emph{i} could be obtained. Figure +\ref{weirdice} shows the repeating crystal structure of a typical +crystal at 5 K. Each water molecule is hydrogen bonded to four others; +however, the hydrogen bonds are flexed rather than perfectly +straight. This results in a skewed tetrahedral geometry about the +central molecule. Referring to figure \ref{isosurface}, these flexed +hydrogen bonds are allowed due to the conical shape of the attractive +regions, with the greatest attraction along the direct hydrogen bond +configuration. Though not ideal, these flexed hydrogen bonds are +favorable enough to stabilize an entire crystal generated around +them. In fact, the imperfect ice \emph{i} crystals were so stable that +they melted at temperatures nearly 100 K greater than both ice I$_c$ +and I$_h$. -These initial simulations indicated that ice 0 is the preferred ice -structure for at least SSD/E. To verify this, a comparison was made -between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at -constant pressure with SSD/E, SSD/RF, and SSD1. Near ideal versions of -the three types of crystals were cooled to 1 K, and the potential -energies of each were compared using all three water models. With -every water model, ice 0 turned out to have the lowest potential -energy: 5\% lower than $I_h$ with SSD1, 6.5\% lower with SSD/E, and -7.5\% lower with SSD/RF. +These initial simulations indicated that ice \emph{i} is the preferred +ice structure for at least the SSD/E model. To verify this, a +comparison was made between near ideal crystals of ice $I_h$, ice +$I_c$, and ice 0 at constant pressure with SSD/E, SSD/RF, and +SSD1. Near ideal versions of the three types of crystals were cooled +to 1 K, and the potential energies of each were compared using all +three water models. With every water model, ice \emph{i} turned out to +have the lowest potential energy: 5\% lower than $I_h$ with SSD1, +6.5\% lower with SSD/E, and 7.5\% lower with SSD/RF. In addition to these low temperature comparisons, melting sequences -were performed with ice 0 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 a reaction field occurred -at temperature in excess of 375 K. SSD/RF and SSD1 with a reaction -field had more reasonable melting transitions, down near 325 K. These -melting point observations emphasize how preferred this crystal -structure is over the most common types of ice when using these single -point water models. +were performed with ice \emph{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 a reaction field +occurred at temperature in excess of 375 K. SSD/RF and SSD1 with a +reaction field showed more reasonable melting transitions near 325 +K. These melting point observations emphasize the preference for this +crystal structure over the most common types of ice when using these +single point water models. -Recognizing that the above tests show ice 0 to be both the most stable -and lowest density crystal structure for these single point water -models, it is interesting to speculate on the relative stability of -this crystal structure with charge based water models. As a quick +Recognizing that the above tests show ice \emph{i} to be both the most +stable and lowest density crystal structure for these single point +water models, it is interesting to speculate on the relative stability +of this crystal structure with charge based water models. As a quick test, these 3 crystal types were converted from SSD type particles to TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy -minimizations were performed on all of these crystals to compare the -system energies. Again, ice 0 was observed to have the lowest total -system energy. The total energy of ice 0 was ~2\% lower than ice -$I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial -results, we would not be surprised if results from the other common -water models show ice 0 to be the lowest energy crystal structure. A -continuation on work studying ice 0 with multi-point water models will -be published in a coming article. +minimizations were performed on the crystals to compare the system +energies. Again, ice \emph{i} was observed to have the lowest total +system energy. The total energy of ice \emph{i} was ~2\% lower than +ice $I_h$, which was in turn ~3\% lower than ice $I_c$. Based on these +initial studies, it would not be surprising if results from the other +common water models show ice \emph{i} to be the lowest energy crystal +structure. A continuation of this work studying ice \emph{i} with +multi-point water models will be published in a coming article. \section{Conclusions} The density maximum and temperature dependent transport for the SSD @@ -888,18 +914,18 @@ capture the transport properties of experimental very density maximum near 260 K. In most cases, the calculated densities were significantly lower than the densities calculated in simulations of other water models. Analysis of particle diffusion showed SSD to -capture the transport properties of experimental very well in both the -normal and super-cooled liquid regimes. In order to correct the +capture the transport properties of experimental water well in both +the liquid and super-cooled 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 and SSD/E), and comparison simulations were performed with SSD1, the density corrected version of SSD. Both models improve the liquid structure, density values, and diffusive properties under their respective conditions, indicating the necessity of reparameterization when altering the -long-range correction specifics. When taking the appropriate -considerations, these simple water models are excellent choices for -representing explicit water in large scale simulations of biochemical -systems. +long-range correction specifics. When taking into account the +appropriate considerations, these simple water models are excellent +choices for representing explicit water in large scale simulations of +biochemical systems. \section{Acknowledgments} Support for this project was provided by the National Science @@ -907,8 +933,10 @@ DMR 00 79647. the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR 00 79647. -\bibliographystyle{jcp} +\newpage + +\bibliographystyle{jcp} \bibliography{nptSSD} %\pagebreak