--- trunk/ssdePaper/nptSSD.tex	2003/11/06 23:00:00	856
+++ trunk/ssdePaper/nptSSD.tex	2003/11/12 13:37:15	861
@@ -1,5 +1,5 @@
-\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
-%\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
+%\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
+\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
 \usepackage{amsmath}
 \usepackage{graphicx}
 
@@ -42,13 +42,12 @@ experimental very well in both the normal and super-co
 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 density behavior, SSD was
+experimental water very well in both the normal and super-cooled
+liquid regimes. In order to correct the density behavior, SSD was
 reparameterized for use both with and without a long-range interaction
-correction, SSD/RF and SSD/E respectively. In addition to correcting
-the abnormally low densities, these new versions were show to maintain
-or improve upon the transport and structural features of the original
-water model.
+correction, SSD/RF and SSD/E respectively. Compared to the density
+corrected version of SSD (SSD1), these modified models were shown to
+maintain or improve upon the structural and transport properties.
 \end{abstract}
 
 \maketitle
@@ -145,14 +144,13 @@ comparable models while maintaining the structural beh
 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
+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
 attractive combination of speed and accurate depiction of solvent
 properties makes SSD a model of interest for the simulation of large
-scale biological systems, such as membrane phase behavior, a specific
-interest within our group.
+scale biological systems, such as membrane phase behavior.
 
 One of the key limitations of this water model, however, is that it
 has been parameterized for use with the Ewald Sum technique for the
@@ -166,15 +164,15 @@ utilizing the RF correction technique. Towards the end
 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. Towards the end, we suggest
+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 original SSD water model to address abnormal
-density behavior (also observed here), calling the corrected model
-SSD1.\cite{Ichiye03} This study will consider this new model's
-behavior as well, and hopefully improve upon its depiction of water
-under conditions without the Ewald Sum.
+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.
 
 \section{Methods}
 
@@ -204,7 +202,7 @@ SSD more compatible with a 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
-SSD more compatible with a reaction field.
+SSD more accurate both with and without a reaction field.
 
 Simulations were performed in both the isobaric-isothermal and
 microcanonical ensembles. The constant pressure simulations were
@@ -310,35 +308,36 @@ the heat capacity, in addition to determining the dens
 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,
+the heat capacity, in addition to determining the density maximum -
 provided that the density maximum occurs in the liquid and not the
-supercooled regimes. An ensemble average from five separate melting
+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, ranging 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 then 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).
+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}
-In the initial average density versus temperature plot, the density
-maximum appears between 255 and 265 K. The calculated densities within
-this range were nearly indistinguishable, as can be seen in the zoom
-of this region of interest, shown in figure
-\ref{dense1}. The greater certainty 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 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 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.
+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, where the calculated densities
+within this range were nearly indistinguishable. The greater certainty
+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
+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
+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.
 
 \begin{figure}
 \includegraphics[width=65mm,angle=-90]{dense2.eps}
@@ -348,43 +347,19 @@ common liquid state density maximum with the $I_h$ sim
  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{dense2}
+\label{dense1}
 \end{figure}
 
 The density maximum for SSD actually compares quite favorably to other
-simple water models. Figure \ref{dense2} shows a plot of these
-findings with the density progression of several other models and
-experiment obtained from other
+simple water models. Figure \ref{dense1} also shows calculated
+densities of several other models and experiment obtained from other
 sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
 models, SSD has results closest to the experimentally observed water
 density maximum. Of the listed water models, TIP4P has a density
-maximum behavior most like that seen in SSD. Though not shown, it is
-useful to note that TIP5P has a water density maximum nearly identical
-to experiment.
+maximum behavior most like that seen in SSD. Though not included in
+this particular plot, it is useful to note that TIP5P has a water
+density maximum nearly identical to experiment.
 
-Possibly of more importance is the density scaling of SSD relative to
-other common models at any given temperature (Fig. \ref{dense2}). 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
-simulations is the presence of the reaction field.\cite{Berendsen98}
-In order to test the effect of the reaction field on the density of
-the systems, the simulations were repeated for the temperature region
-of interest without a reaction field present. The results of these
-simulations are also displayed in figure \ref{dense2}. Without
-reaction field, these 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 slight
-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 simulations.\cite{Clancy94,Jorgensen98b}
-
 It has been observed that densities are dependent on the cutoff radius
 used for a variety of water models in simulations both with and
 without the use of reaction field.\cite{Berendsen98} In order to
@@ -396,9 +371,37 @@ use of a longer cutoff radius results in a near doubli
 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 a near doubling of the time
-required to compute a single trajectory.
+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
+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
+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
+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
+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
+was what lead Ichiye \emph{et al.} to recently reparameterize SSD and
+make SSD1.\cite{Ichiye03} In discussing potential adjustments later in
+this paper, all comparisons were performed with this new model.
+
 \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
@@ -434,9 +437,7 @@ experiment at these temperatures, albeit not at standa
 SSD are lower than experimental water at temperatures higher than room
 temperature. 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. Results under these conditions can be found later in this
-paper.
+experiment at these temperatures, albeit not at standard pressure.
 
 \subsection{Structural Changes and Characterization}
 By starting the simulations from the crystalline state, the melting
@@ -449,8 +450,7 @@ surprising in that SSD is a simple rigid body model wi
 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 in that SSD is a simple rigid body model with a fixed
-dipole.
+surprising when considering the simplicity of the SSD model.
 
 \begin{figure} 
 \includegraphics[width=85mm]{fullContours.eps} 
@@ -462,10 +462,6 @@ Additional analyses for understanding the melting phas
 \label{contour} 
 \end{figure} 
 
-Additional analyses for understanding the melting phase-transition
-process were performed via two-dimensional structure and dipole angle
-correlations. Expressions for the correlations are as follows:
-
 \begin{figure}
 \includegraphics[width=45mm]{corrDiag.eps} 
 \caption{Two dimensional illustration of the angles involved in the 
@@ -473,6 +469,10 @@ 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\ ,
@@ -481,26 +481,27 @@ where $\theta$ and $\omega$ refer to the angles shown 
 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}
-where $\theta$ and $\omega$ refer to the angles shown in the above
-illustration. By binning over both distance and the cosine of the
+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
 dissected to determine the common dipole arrangements that constitute
 the peaks and troughs. Frames A and B of figure \ref{contour} show a
 relatively crystalline state of an ice $I_c$ simulation. The first
-peak of the g(\emph{r}) primarily consists of the preferred hydrogen
-bonding arrangements as dictated by the tetrahedral sticky potential,
+peak of the g(\emph{r}) consists primarily of the preferred hydrogen
+bonding arrangements as dictated by the tetrahedral sticky potential -
 one peak for the donating and the other for the accepting hydrogen
 bonds. Due to the high degree of crystallinity of the sample, the
 second and third solvation shells show a repeated peak arrangement
 which decays at distances around the fourth solvation shell, near the
 imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
-In the higher temperature simulation shown in frames C and D, the
-repeated peak features are significantly blurred. The first solvation
-shell still shows the strong effect of the sticky-potential, although
-it covers a larger area, extending to include a fraction of aligned
-dipole peaks within the first solvation shell. The latter peaks lose
-definition as thermal motion and the competing dipole force overcomes
-the sticky potential's tight tetrahedral structuring of the fluid.
+In the higher temperature simulation shown in frames C and D, these
+longer-ranged repeated peak features deteriorate rapidly. The first
+solvation shell still shows the strong effect of the sticky-potential,
+although it covers a larger area, extending to include a fraction of
+aligned dipole peaks within the first solvation shell. The latter
+peaks lose definition as thermal motion and the competing dipole force
+overcomes the sticky potential's tight tetrahedral structuring of the
+fluid.
 
 This complex interplay between dipole and sticky interactions was
 remarked upon as a possible reason for the split second peak in the
@@ -517,34 +518,35 @@ immediately outside of it's range (around 4 \AA), and 
 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 of it's range (around 4 \AA), and arrangements
-that seek to ideally satisfy both the sticky and dipole forces locate
-themselves just beyond this region (around 5 \AA).
+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).
 
 From these findings, the split second peak is primarily the product of
-the dipolar repulsion term of the sticky potential. In fact, the
-leading of the two peaks can be pushed out and merged with the outer
-split peak just by 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 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 previously. A better correction would be to include
-the quadrupole-quadrupole interactions for the water particles outside
-of the first solvation shell, but this reduces the simplicity and
-speed advantage of SSD, so it is not the most desirable path to take.
+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 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
+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
+previously. A better correction would be to include the
+quadrupole-quadrupole interactions for the water particles outside of
+the first solvation shell, but this reduces the simplicity and speed
+advantage of SSD.
 
-\subsection{Adjusted Potentials: SSD/E and SSD/RF}
+\subsection{Adjusted Potentials: SSD/RF and SSD/E}
 The propensity of SSD to adopt lower than expected densities under
 varying conditions is troubling, especially at higher temperatures. In
-order to correct this behavior, it's necessary to adjust the force
-field parameters for the primary intermolecular interactions. In
-undergoing a reparameterization, it is important not to focus on just
-one property and neglect the other important properties. In this case,
-it would be ideal to correct the densities while maintaining the
-accurate transport properties.
+order to correct this model for use with a reaction field, it is
+necessary to adjust the force field parameters for the primary
+intermolecular interactions. In undergoing a reparameterization, it is
+important not to focus on just one property and neglect the other
+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$
 Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
@@ -570,8 +572,10 @@ parameterized potentials - soft sticky dipole enhanced
 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 enhanced (SSD/E) and
-soft sticky dipole reaction field (SSD/RF).
+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}
 \caption{Parameters for the original and adjusted models}
@@ -622,69 +626,69 @@ g(\emph{r}) by comparing the original SSD (with and wi
 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 original SSD (with and without reaction
-field), 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 a
-reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
-of the sticky potential strength and cutoffs. 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 the unphysical
-arrangements. This compares with the original SSD's excluding dipolar
-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 allowing the persistence of
-full dipolar character below the previous 4.0 \AA\ cutoff.
+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
+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 to some degree, 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 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
+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
+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 improve the
+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 available for
-adjusting the dipole
+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.418 and 2.48 D for SSD/E and SSD/RF
-respectively is moderate in the range of the experimental values;
-however, it leads to significant changes in the density and transport
-of the water models.
+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.
 
-In order to demonstrate the benefits of this reparameterization, a
+In order to demonstrate the benefits of these reparameterizations, a
 series of NPT and NVE simulations were performed to probe the density
 and transport properties of the adapted models and compare the results
 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 both self-consistent and experimental
-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 previously. Like
-before, all of the NPT simulations were equilibrated for 100 ps before
-a 200 ps data collection run, and they used the previous temperature's
-final configuration as a starting point. All of the NVE simulations
-had the same thermalization, equilibration, and data collection times
-stated earlier in this paper.
+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
+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.
 
 \begin{figure} 
 \includegraphics[width=62mm, angle=-90]{ssdeDense.epsi} 
-\caption{Comparison of densities calculated with SSD/E to SSD without a 
+\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
 expansion around 300 K with error bars included to clarify this region
@@ -693,30 +697,29 @@ Figure \ref{ssdedense} shows the density profile for t
 \label{ssdedense} 
 \end{figure} 
 
-Figure \ref{ssdedense} shows the density profile for the SSD/E water
-model in comparison to the original SSD without a reaction field,
-experiment, and the other common water models considered
-previously. The calculated densities have increased significantly over
-the original SSD model and match the experimental value just below 298
-K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
-compares well with the experimental value of 0.997 g/cm$^3$ and is
-considerably better than the SSD value of 0.967$\pm$0.003
-g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
-out the curve at higher temperatures, only the improvement is marginal
-at best. This steep drop in densities is due to the dipolar rather
-than charge based interactions which decay more rapidly at longer
-distances.
- 
-By monitoring C$\text{p}$ throughout these simulations, the melting
-transition for SSD/E was observed at 230 K, about 5 degrees lower than
-SSD. The resulting density maximum is located at 240 K, again about 5
-degrees lower than the SSD value of 245 K. Though there is a decrease
-in both of these values, the corrected densities near room temperature
-justify the modifications taken.
+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.
 
 \begin{figure} 
 \includegraphics[width=62mm, angle=-90]{ssdrfDense.epsi} 
-\caption{Comparison of densities calculated with SSD/RF to SSD with a 
+\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
 necessity of reparameterization when utilizing a reaction field
@@ -725,38 +728,58 @@ Figure \ref{ssdrfdense} shows a density comparison bet
 \label{ssdrfdense} 
 \end{figure} 
 
-Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
-SSD with an active reaction field. Like in the simulations of SSD/E,
-the densities show a dramatic increase over normal SSD. 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$. The melting point is observed at 240 K,
-which is 5 degrees lower than SSD with a reaction field, and the
-density maximum at 255 K, again 5 degrees lower than SSD. The density
-at higher temperature still drops off more rapidly than the charge
-based models but is in better agreement than SSD/E.
+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. 
 
+\begin{figure} 
+\includegraphics[width=65mm, angle=-90]{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.}
+\label{ssdediffuse} 
+\end{figure} 
+
 The reparameterization of the SSD water model, both for use with and
 without an applied long-range correction, brought the densities up to
 what is expected for simulating liquid water. In addition to improving
 the densities, it is important that particle transport be maintained
 or improved. Figure \ref{ssdediffuse} compares the temperature
-dependence of the diffusion constant of SSD/E to SSD without an active
-reaction field, both at the densities calculated at 1 atm and at the
-experimentally calculated densities for super-cooled and liquid
+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 SSD starts off
-slower than experiment and crosses to merge with SSD/E at high
-temperatures. Both models follow the experimental trend well, but
-diffuse too rapidly at higher temperatures. This abnormally fast
-diffusion is caused by the decreased system density. Since the
-densities of SSD/E don't deviate as much from experiment as those of
-SSD, it follows the experimental trend more closely. This observation
-is backed up by looking at the lower plot. The diffusion constants for
-SSD/E track with the experimental values while SSD deviates on the low
-side of the trend with increasing temperature. This is again a product
-of SSD/E having densities closer to experiment, and not deviating to
-lower densities with increasing temperature as rapidly.
+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.
 
 \begin{figure} 
 \includegraphics[width=65mm, angle=-90]{ssdrfDiffuse.epsi} 
@@ -772,101 +795,71 @@ lower densities with increasing temperature as rapidly
 \label{ssdrfdiffuse} 
 \end{figure} 
 
-\begin{figure} 
-\includegraphics[width=65mm, angle=-90]{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.}
-\label{ssdediffuse} 
-\end{figure} 
-
 In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
-compared with SSD with an active reaction field. In the upper plot,
-SSD/RF tracks with the experimental results incredibly well, identical
-within error throughout the temperature range and only showing a
-slight increasing trend at higher temperatures. SSD also tracks
-experiment well, only it tends to diffuse a little more slowly at low
-temperatures and deviates to diffuse too rapidly at high
-temperatures. 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 SSD, because the calculated densities are more true to
-experiment. This is again emphasized in the lower plot, where SSD/RF
-tracks the experimental diffusion exactly while SSD's diffusion
-constants are slightly too low due to its need for a lower density at
-the specified temperature.
+compared with 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
+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
+long-range correction.
 
 \subsection{Additional Observations}
 
-While performing the melting sequences of SSD/E, some interesting
-observations were made. After melting at 230 K, two of the 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. These simulations were excluded from the data
-set shown in figure \ref{ssdedense} and replaced with two additional
-melting sequences that did not undergo this anomalous phase
-transition, while this crystallization event was investigated
-separately.
-
 \begin{figure}
 \includegraphics[width=85mm]{povIce.ps} 
-\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
+\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.}
 \label{weirdice} 
 \end{figure}
 
-The final configurations of these two melting sequences shows 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 than a near ideal crystal structure could be
-obtained. Figure \ref{weirdice} shows the repeating crystal structure
-of a typical crystal at 5 K. The unit cell contains eight molecules,
-and figure \ref{unitcell} shows a unit cell built from the water
-particle center of masses that can be used to construct a repeating
-lattice, similar to figure \ref{weirdice}. Each molecule is hydrogen
-bonded to four other water molecules; 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 greater than room
-temperature.
+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$.
 
-\begin{figure}
-\includegraphics[width=65mm]{ice0cell.eps} 
-\caption{Simple unit cell for constructing ice 0. In this cell, $c$ is 
-equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
-\label{unitcell} 
-\end{figure}
-
-The initial simulations indicated that ice 0 is the preferred ice
+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 SSD. Near ideal versions of
-the three types of crystals were cooled to ~1 K, and the potential
+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 SSD, 6.5\% lower with SSD/E, and
-7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
-was observed to be ~2\% less stable than ice $I_h$. In addition to
-having the lowest potential energy, ice 0 was the most expanded of the
-three ice crystals, ~5\% less dense than ice $I_h$ with all of the
-water models. In all three water models, ice $I_c$ was observed to be
-~2\% more dense than ice $I_h$. 
+energy: 5\% lower than $I_h$ with SSD1, 6.5\% lower with SSD/E, and
+7.5\% lower with SSD/RF. 
 
-In addition to the low temperature comparisons, melting sequences were
-performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
-and SSD both with and without a reaction field. The melting
-transitions for both SSD/E and SSD without a reaction field occurred
-at temperature in excess of 375 K. SSD/RF and SSD with a reaction
+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
@@ -874,8 +867,8 @@ models, it is interesting to speculate on the favorabi
 
 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 favorability of this
-crystal structure with the different charge based models. As a quick
+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
@@ -884,7 +877,7 @@ continuation on work studing ice 0 with multipoint wat
 $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 studing ice 0 with multipoint water models will
+continuation on work studying ice 0 with multi-point water models will
 be published in a coming article.
 
 \section{Conclusions}
@@ -897,15 +890,16 @@ density behavior, SSD was reparameterized for use both
 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
-density behavior, SSD was reparameterized for use both with and
-without a long-range interaction correction, SSD/RF and SSD/E
-respectively. In addition to correcting the abnormally low densities,
-these new versions were show to maintain or improve upon the transport
-and structural features of the original water model, all while
-maintaining the fast performance of the SSD water model. This work
-shows these simple water models, and in particular SSD/E and SSD/RF,
-to be excellent choices to represent explicit water in future
-simulations of biochemical systems.
+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.
 
 \section{Acknowledgments}
 Support for this project was provided by the National Science