--- trunk/ssdePaper/nptSSD.tex	2004/01/12 16:20:53	921
+++ trunk/ssdePaper/nptSSD.tex	2004/02/04 20:19:21	1019
@@ -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})\  ,
@@ -591,29 +594,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 +617,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 +795,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 +813,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 +851,9 @@ throughout the temperature range shown and with only a
 In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
 compared to SSD1 with an active reaction field. Note that SSD/RF
 tracks the experimental results quantitatively, identical within error
-throughout the temperature range shown and with only a slight
-increasing trend at higher temperatures. SSD1 tends to diffuse more
-slowly at low temperatures and deviates to diffuse too rapidly at
+throughout most of the temperature range shown and exhibiting only a
+slight increasing trend at higher temperatures. SSD1 tends to diffuse
+more slowly at low temperatures and deviates to diffuse too rapidly at
 temperatures greater than 330 K.  As stated above, this deviation away
 from the ideal trend is due to a rapid decrease in density at higher
 temperatures. SSD/RF does not suffer from this problem as much as SSD1
@@ -869,6 +861,71 @@ 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. Again,
+SSD/E and SSD/RF show improved behavior over SSD1 both with and
+without an active reaction field. Numbers published from the original
+SSD dynamics studies appear closer to the experimental values, and we
+attribute this discrepancy to the implimentation of an Ewald sum
+versus a reaction field.
+
 \subsection{Additional Observations}
 
 \begin{figure}
