--- trunk/ssdePaper/nptSSD.tex	2004/01/12 16:20:53	921
+++ trunk/ssdePaper/nptSSD.tex	2004/02/04 22:13:36	1022
@@ -139,6 +139,7 @@ + s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i
 \frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)
 + s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf
 \Omega}_j)]\ .
+\label{stickyfunction}
 \end{equation}
 Here, $\nu_0$ is a strength parameter for the sticky potential, and
 $s$ and $s^\prime$ are cubic switching functions which turn off the
@@ -151,7 +152,7 @@ w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
 while the $w^\prime$ function counters the normal aligned and
 anti-aligned structures favored by point dipoles:
 \begin{equation}
-w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
+w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^\circ,
 \end{equation}
 It should be noted that $w$ is proportional to the sum of the $Y_3^2$
 and $Y_3^{-2}$ spherical harmonics (a linear combination which
@@ -208,8 +209,10 @@ cubic switching function at a cutoff radius.  Under th
 
 Long-range dipole-dipole interactions were accounted for in this study
 by using either the reaction field method or by resorting to a simple
-cubic switching function at a cutoff radius.  Under the first method,
-the magnitude of the reaction field acting on dipole $i$ is
+cubic switching function at a cutoff radius.  The reaction field
+method was actually first used in Monte Carlo simulations of liquid
+water.\cite{Barker73} Under this method, the magnitude of the reaction
+field acting on dipole $i$ is
 \begin{equation}
 \mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
 \frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij})\  ,
@@ -232,8 +235,8 @@ at the cutoff radius) and as a result we have two repa
  
 We have also performed a companion set of simulations {\it without} a
 surrounding dielectric (i.e. using a simple cubic switching function
-at the cutoff radius) and as a result we have two reparamaterizations
-of SSD which could be used either with or without the Reaction Field
+at the cutoff radius), and as a result we have two reparamaterizations
+of SSD which could be used either with or without the reaction field
 turned on.
 
 Simulations to obtain the preferred density were performed in the
@@ -251,12 +254,11 @@ traditional quaternion integration.\cite{Evans77,Evans
 symplectic splitting method proposed by Dullweber {\it et
 al.}\cite{Dullweber1997} Our reason for selecting this integrator
 centers on poor energy conservation of rigid body dynamics using
-traditional quaternion integration.\cite{Evans77,Evans77b} While quaternions
-may work well for orientational motion under NVT or NPT integrators,
-our limits on energy drift in the microcanonical ensemble were quite
-strict, and the drift under quaternions was substantially greater than
-in the symplectic splitting method.  This steady drift in the total
-energy has also been observed by Kol {\it et al.}\cite{Laird97}
+traditional quaternion integration.\cite{Evans77,Evans77b} In typical
+microcanonical ensemble simulations, the energy drift when using
+quaternions was substantially greater than when using the symplectic
+splitting method (fig. \ref{timestep}).  This steady drift in the
+total energy has also been observed by Kol {\it et al.}\cite{Laird97}
 
 The key difference in the integration method proposed by Dullweber
 \emph{et al.} is that the entire rotation matrix is propagated from
@@ -446,7 +448,7 @@ results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} 
 mean-square displacement as a function of time. The averaged results
 from five sets of NVE simulations are displayed in figure
 \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
-results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01} 
+results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01} 
 
 \begin{figure}
 \begin{center}
@@ -454,7 +456,7 @@ and Experimental data [Refs. \citen{Gillen72} and \cit
 \epsfbox{betterDiffuse.epsi} 
 \caption{Average self-diffusion constant as a function of temperature for 
 SSD, SPC/E [Ref. \citen{Clancy94}], TIP5P [Ref. \citen{Jorgensen01}],
-and Experimental data [Refs. \citen{Gillen72} and \citen{Mills73}]. Of
+and Experimental data [Refs. \citen{Gillen72} and \citen{Holz00}]. Of
 the three water models shown, SSD has the least deviation from the
 experimental values. The rapidly increasing diffusion constants for
 TIP5P and SSD correspond to significant decrease in density at the
@@ -591,29 +593,15 @@ The parameters available for tuning include the $\sigm
 important properties. In this case, it would be ideal to correct the
 densities while maintaining the accurate transport behavior.
 
-The parameters available for tuning include the $\sigma$ and $\epsilon$
-Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
-attractive and dipole repulsive terms with their respective
-cutoffs. To alter the attractive and repulsive terms of the sticky
-potential independently, it is necessary to separate the terms as
-follows:
-\begin{equation}
-u_{ij}^{sp}
-({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
-\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)],
-\end{equation}
-where $\nu_0$ scales the strength of the tetrahedral attraction and
-$\nu_0^\prime$ scales the dipole repulsion term independently. The
-separation was performed for purposes of the reparameterization, but
-the final parameters were adjusted so that it is not necessary to
-separate the terms when implementing the adjusted water
-potentials. The results of the reparameterizations are shown in table
-\ref{params}. Note that the tetrahedral attractive and dipolar
-repulsive terms do not share the same lower cutoff ($r_l$) in the
-newly parameterized potentials.  We are calling these
-reparameterizations the Soft Sticky Dipole / Reaction Field
+The parameters available for tuning include the $\sigma$ and
+$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the
+strength of the sticky potential ($\nu_0$), and the sticky attractive
+and dipole repulsive cubic switching function cutoffs ($r_l$, $r_u$
+and $r_l^\prime$, $r_u^\prime$ respectively). The results of the
+reparameterizations are shown in table \ref{params}. We are calling
+these reparameterizations the Soft Sticky Dipole / Reaction Field
 (SSD/RF - for use with a reaction field) and Soft Sticky Dipole
-Enhanced (SSD/E - an attempt to improve the liquid structure in
+Extended (SSD/E - an attempt to improve the liquid structure in
 simulations without a long-range correction).
 
 \begin{table}
@@ -628,9 +616,9 @@ simulations without a long-range correction).
 \ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
 \ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\
 \ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
+\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\
 \ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
 \ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
-\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
 \ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\
 \ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
 \end{tabular}
@@ -806,13 +794,14 @@ K, shown by SSD and SSD1 respectively.
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdeDiffuse.epsi} 
-\caption{Plots of the diffusion constants calculated from SSD/E and SSD1,
-both without a reaction field, along with experimental results
-[Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations were
-performed at the average densities observed in the 1 atm NPT
-simulations for the respective models. SSD/E is slightly more fluid
-than experiment at all of the temperatures, but it is closer than SSD1
-without a long-range correction.}
+\caption{The diffusion constants calculated from SSD/E and SSD1,
+ both without a reaction field, along with experimental results
+ [Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations
+ were performed at the average densities observed in the 1 atm NPT
+ simulations for the respective models. SSD/E is slightly more mobile
+ than experiment at all of the temperatures, but it is closer to
+ experiment at biologically relavent temperatures than SSD1 without a
+ long-range correction.}
 \label{ssdediffuse} 
 \end{center}
 \end{figure} 
@@ -823,35 +812,37 @@ without an active reaction field, both at the densitie
 the densities, it is important that the excellent diffusive behavior
 of SSD be maintained or improved. Figure \ref{ssdediffuse} compares
 the temperature dependence of the diffusion constant of SSD/E to SSD1
-without an active reaction field, both at the densities calculated at
-1 atm and at the experimentally calculated densities for super-cooled
-and liquid water. The diffusion constant for SSD/E is consistently
-higher than experiment, while SSD1 remains lower than experiment until
-relatively high temperatures (greater than 330 K). Both models follow
-the shape of the experimental curve well below 300 K but tend to
-diffuse too rapidly at higher temperatures, something that is
-especially apparent with SSD1.  This increasing diffusion relative to
-the experimental values is caused by the rapidly decreasing system
-density with increasing temperature.  The densities of SSD1 decay more
-rapidly with temperature than do those of SSD/E, leading to more
-visible deviation from the experimental diffusion trend.  Thus, the
-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.
+without an active reaction field at the densities calculated from the
+NPT simulations at 1 atm. The diffusion constant for SSD/E is
+consistently higher than experiment, while SSD1 remains lower than
+experiment until relatively high temperatures (around 360 K). Both
+models follow the shape of the experimental curve well below 300 K but
+tend to diffuse too rapidly at higher temperatures, as seen in SSD1's
+crossing above 360 K.  This increasing diffusion relative to the
+experimental values is caused by the rapidly decreasing system density
+with increasing temperature.  Both SSD1 and SSD/E show this deviation
+in diffusive behavior, but this trend has different implications on
+the diffusive behavior of the models.  While SSD1 shows more
+experimentally accurate diffusive behavior in the high temperature
+regimes, SSD/E shows more accurate behavior in the supercooled and
+biologically relavent temperature ranges.  Thus, the changes made to
+improve the liquid structure may have had an adverse affect on the
+density maximum, but they improve the transport behavior of SSD/E
+relative to SSD1 under the most commonly simulated conditions.
 
 \begin{figure}
 \begin{center} 
 \epsfxsize=6in
 \epsfbox{ssdrfDiffuse.epsi} 
-\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1,
+\caption{The diffusion constants calculated from SSD/RF and SSD1,
  both with an active reaction field, along with experimental results
- [Refs. \citen{Gillen72} and \citen{Mills73}]. The NVE calculations
+ [Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations
  were performed at the average densities observed in the 1 atm NPT
  simulations for both of the models. Note how accurately SSD/RF
  simulates the diffusion of water throughout this temperature
  range. The more rapidly increasing diffusion constants at high
- temperatures for both models is attributed to the significantly lower
- densities than observed in experiment.}
+ temperatures for both models is attributed to lower calculated
+ densities than those observed in experiment.}
 \label{ssdrfdiffuse} 
 \end{center}
 \end{figure} 
@@ -859,9 +850,9 @@ throughout the temperature range shown and with only a
 In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
 compared to SSD1 with an active reaction field. Note that SSD/RF
 tracks the experimental results quantitatively, identical within error
-throughout the temperature range shown and with only a slight
-increasing trend at higher temperatures. SSD1 tends to diffuse more
-slowly at low temperatures and deviates to diffuse too rapidly at
+throughout most of the temperature range shown and exhibiting only a
+slight increasing trend at higher temperatures. SSD1 tends to diffuse
+more slowly at low temperatures and deviates to diffuse too rapidly at
 temperatures greater than 330 K.  As stated above, this deviation away
 from the ideal trend is due to a rapid decrease in density at higher
 temperatures. SSD/RF does not suffer from this problem as much as SSD1
@@ -869,6 +860,80 @@ reparameterization when using an altered long-range co
 values. These results again emphasize the importance of careful
 reparameterization when using an altered long-range correction.
 
+\begin{table}
+\begin{center}
+\caption{Calculated and experimental properties of the single point waters and liquid water at 298 K and 1 atm. (a) Ref. [\citen{Mills73}]. (b) Calculated by integrating the data in ref. \citen{Head-Gordon00_1}. (c) Calculated by integrating the data in ref. \citen{Soper86}. (d) Ref. [\citen{Eisenberg69}]. (e) Calculated for 298 K from data in ref. \citen{Krynicki66}.}
+\begin{tabular}{ l  c  c  c  c  c } 
+\hline \\[-3mm]
+\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \ 
+\ & \ SSD/RF \ \ \ & \ Expt. \\
+ \hline \\[-3mm]
+\ \ \ $\rho$ (g/cm$^3$) & 0.999 $\pm$0.001 & 0.996 $\pm$0.001 & 0.972 $\pm$0.002 & 0.997 $\pm$0.001 & 0.997 \\
+\ \ \ $C_p$ (cal/mol K) & 28.80 $\pm$0.11 & 25.45 $\pm$0.09 & 28.28 $\pm$0.06 & 23.83 $\pm$0.16 & 17.98 \\
+\ \ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78 $\pm$0.07 & 2.51 $\pm$0.18 & 2.00 $\pm$0.17 & 2.32 $\pm$0.06 & 2.299$^\text{a}$ \\
+\ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\
+\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.4$^\text{c}$ \\
+\ \ \ $\tau_1^\mu$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 & 7.2 $\pm$0.4 & 4.76$^\text{d}$ \\
+\ \ \ $\tau_2^\mu$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2 $\pm$0.2 & 2.3$^\text{e}$ \\
+\end{tabular}
+\label{liquidproperties}
+\end{center}
+\end{table}
+
+Table \ref{liquidproperties} gives a synopsis of the liquid state
+properties of the water models compared in this study along with the
+experimental values for liquid water at ambient conditions. The
+coordination number and hydrogen bonds per particle were calculated by
+integrating the following relation:
+\begin{equation}
+4\pi\rho\int_{0}^{a}r^2\text{g}(r)dr,
+\end{equation}
+where $\rho$ is the number density of pair interactions, $a$ is the
+radial location of the minima following the first solvation shell
+peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for
+calculation of the coordination number or hydrogen bonds per particle
+respectively. The number of hydrogen bonds stays relatively constant
+across all of the models, but the coordination numbers of SSD/E and
+SSD/RF show an improvement over SSD1. This improvement is primarily
+due to the widening of the first solvation shell peak, allowing the
+first minima to push outward. Comparing the coordination number with
+the number of hydrogen bonds can lead to more insight into the
+structural character of the liquid.  Because of the near identical
+values for SSD1, it appears to be a little too exclusive, in that all
+molecules in the first solvation shell are involved in forming ideal
+hydrogen bonds.  The differing numbers for the newly parameterized
+models indicate the allowance of more fluid configurations in addition
+to the formation of an acceptable number of ideal hydrogen bonds.
+
+The time constants for the self orientational autocorrelation function
+are also displayed in Table \ref{liquidproperties}. The dipolar
+orientational time correlation function ($\Gamma_{l}$) is described
+by:
+\begin{equation}
+\Gamma_{l}(t) = \langle P_l[\mathbf{u}_j(0)\cdot\mathbf{u}_j(t)]\rangle,
+\end{equation}
+where $P_l$ is a Legendre polynomial of order $l$ and $\mathbf{u}_j$
+is the unit vector of the particle dipole.\cite{Rahman71} From these
+correlation functions, the orientational relaxation time of the dipole
+vector can be calculated from an exponential fit in the long-time
+regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculation of these
+time constants were averaged from five detailed NVE simulations
+performed at the STP density for each of the respective models. It
+should be noted that the commonly cited value for $\tau_2$ of 1.9 ps
+was determined from the NMR data in reference \citen{Krynicki66} at a
+temperature near 34$^\circ$C.\cite{Rahman73} Because of the strong
+temperature dependence of $\tau_2$, it is necessary to recalculate it
+at 298 K to make proper comparisons. The value shown in Table
+\ref{liquidproperties} was calculated from the same NMR data in the 
+fashion described in reference \citen{Krynicki66}. Again, SSD/E and
+SSD/RF show improved behavior over SSD1, both with and without an
+active reaction field. Turning on the reaction field leads to much
+improved time constants for SSD1; however, these results also include
+a corresponding decrease in system density. Numbers published from the
+original SSD dynamics studies appear closer to the experimental
+values, and this difference can be attributed to the use of the Ewald
+sum technique versus a reaction field.\cite{Ichiye99}
+
 \subsection{Additional Observations}
 
 \begin{figure}
