--- trunk/ssdePaper/nptSSD.tex 2003/09/10 22:42:57 759 +++ trunk/ssdePaper/nptSSD.tex 2003/11/06 23:00:00 856 @@ -36,12 +36,12 @@ for the SSD water model, both with and without the use \begin{abstract} NVE and NPT molecular dynamics simulations were performed in order to investigate the density maximum and temperature dependent transport -for the SSD water model, both with and without the use of reaction -field. The constant pressure simulations of the melting of both $I_h$ -and $I_c$ ice showed a 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 +for SSD and related water models, both with and without the use of +reaction field. The constant pressure simulations of the melting of +both $I_h$ and $I_c$ ice showed a 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 density behavior, SSD was reparameterized for use both with and without a long-range interaction @@ -146,9 +146,9 @@ water.\cite{Ichiye96} In the original molecular dynami 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 -water.\cite{Ichiye96} In the original molecular dynamics studies of -SSD, it was shown that it actually improves upon the prediction of -water's dynamical properties 3 and 4-point models.\cite{Ichiye99} This +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 @@ -205,14 +205,14 @@ SSD more compatible with a reaction field. 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. - + Simulations were performed in both the isobaric-isothermal and microcanonical ensembles. The constant pressure simulations were implemented using an integral thermostat and barostat as outlined by -Hoover.\cite{Hoover85,Hoover86} For the constant pressure -simulations, the \emph{Q} parameter for the was set to 5.0 amu -\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at -100 ps. +Hoover.\cite{Hoover85,Hoover86} All particles were treated as +non-linear rigid bodies. Vibrational constraints are not necessary in +simulations of SSD, because there are no explicit hydrogen atoms, and +thus no molecular vibrational modes need to be considered. Integration of the equations of motion was carried out using the symplectic splitting method proposed by Dullweber \emph{et @@ -220,11 +220,11 @@ requirement that is actually quite sensitive to errors 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 -requirement that is actually quite sensitive to errors in the -equations of motion. The original implementation of this code utilized -quaternions for rotational motion propagation; however, a detailed -investigation showed that they resulted in a steady drift in the total -energy, something that has been observed by others.\cite{Laird97} +requirement that is quite sensitive to errors in the equations of +motion. The original implementation of this code utilized quaternions +for rotational motion propagation; however, a detailed investigation +showed that they resulted in a steady drift in the total energy, +something that has been observed by others.\cite{Laird97} The key difference in the integration method proposed by Dullweber \emph{et al.} is that the entire rotation matrix is propagated from @@ -244,11 +244,12 @@ simpler arithmetic quaternion propagation. On average, 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. On average, a 1000 SSD -particle simulation shows a 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}. +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} @@ -307,31 +308,26 @@ constant pressure and temperature dynamics. This invol \section{Results and discussion} Melting studies were performed on the randomized ice crystals using -constant pressure and temperature dynamics. This involved an initial -randomization of velocities about the starting temperature of 25 K for -varying amounts of time. The systems were all 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). In order to further improve statistics, five separate -simulation progressions were performed, and the averaged results from -the $I_h$ melting simulations are shown in figure \ref{dense1}. - -\begin{figure} -\includegraphics[width=65mm, angle=-90]{1hdense.epsi} -\caption{Average density of SSD water at increasing temperatures -starting from ice $I_h$ lattice.} -\label{dense1} -\end{figure} +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 regimes. 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). \subsection{Density Behavior} In the initial average density versus temperature plot, the density -maximum clearly 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 +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$ @@ -347,8 +343,11 @@ TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD wi \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}. } + 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{dense2} \end{figure} @@ -576,37 +575,38 @@ soft sticky dipole reaction field (SSD/RF). \begin{table} \caption{Parameters for the original and adjusted models} -\begin{tabular}{ l c c c } +\begin{tabular}{ l c c c c } \hline \\[-3mm] -\ Parameters & \ \ \ SSD$^\dagger$\ \ \ \ & \ SSD/E\ \ & \ SSD/RF\ \ \\ +\ \ \ Parameters\ \ \ & \ \ \ SSD$^\dagger$ \ \ \ & \ SSD1$^\ddagger$\ \ & \ SSD/E\ \ & \ SSD/RF \\ \hline \\[-3mm] -\ \ \ $\sigma$ (\AA) & 3.051 & 3.035 & 3.019\\ -\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\ -\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\ -\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\ -\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\ -\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\ -\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\ -\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\ -\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\ +\ \ \ $\sigma$ (\AA) & 3.051 & 3.016 & 3.035 & 3.019\\ +\ \ \ $\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\\ +\ \ \ $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\\ \\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96} +\\$^\ddagger$ ref. \onlinecite{Ichiye03} \end{tabular} \label{params} \end{table} \begin{figure} -\includegraphics[width=85mm]{gofrCompare.epsi} +\includegraphics[width=85mm]{GofRCompare.epsi} \caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E -and SSD without reaction field (top), as well as SSD/RF and SSD with +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 - SSD (with and without reaction field).} +and SSD/RF, and dotted line - SSD1 (with and without reaction field).} \label{grcompare} \end{figure} \begin{figure} \includegraphics[width=85mm]{dualsticky.ps} -\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \& +\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.} \label{isosurface} @@ -683,12 +683,13 @@ stated earlier in this paper. stated earlier in this paper. \begin{figure} -\includegraphics[width=85mm]{ssdecompare.epsi} +\includegraphics[width=62mm, angle=-90]{ssdeDense.epsi} \caption{Comparison of densities calculated with SSD/E to SSD without a -reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b}, -SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot -includes error bars, and the calculated results from the other -references were removed for clarity.} +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{figure} @@ -714,12 +715,13 @@ justify the modifications taken. justify the modifications taken. \begin{figure} -\includegraphics[width=85mm]{ssdrfcompare.epsi} +\includegraphics[width=62mm, angle=-90]{ssdrfDense.epsi} \caption{Comparison of densities calculated with SSD/RF to SSD with a -reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b}, -SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot -includes error bars, and the calculated results from the other -references were removed for clarity.} +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{figure} @@ -757,25 +759,29 @@ lower densities with increasing temperature as rapidly lower densities with increasing temperature as rapidly. \begin{figure} -\includegraphics[width=85mm]{ssdediffuse.epsi} -\caption{Plots of the diffusion constants calculated from SSD/E and SSD, - both without a reaction field along with experimental results from - Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The - upper plot is at densities calculated from the NPT simulations at a - pressure of 1 atm, while the lower plot is at the experimentally - calculated densities.} -\label{ssdediffuse} +\includegraphics[width=65mm, angle=-90]{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.} +\label{ssdrfdiffuse} \end{figure} \begin{figure} -\includegraphics[width=85mm]{ssdrfdiffuse.epsi} -\caption{Plots of the diffusion constants calculated from SSD/RF and SSD, - both with an active reaction field along with experimental results +\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 - upper plot is at densities calculated from the NPT simulations at a - pressure of 1 atm, while the lower plot is at the experimentally - calculated densities.} -\label{ssdrfdiffuse} + 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 @@ -902,10 +908,10 @@ The authors would like to thank the National Science F simulations of biochemical systems. \section{Acknowledgments} -The authors would like to thank the National Science Foundation for -funding under grant CHE-0134881. Computation time was provided by the -Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR -00 79647. +Support for this project was provided by the National Science +Foundation under grant CHE-0134881. Computation time was provided by +the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant +DMR 00 79647. \bibliographystyle{jcp}