--- trunk/ssdePaper/nptSSD.tex 2003/09/03 21:36:09 743 +++ trunk/ssdePaper/nptSSD.tex 2003/09/19 19:29:24 777 @@ -23,7 +23,7 @@ \begin{document} -\title{On the temperature dependent structural and transport properties of the soft sticky dipole (SSD) and related single point water models} +\title{On the temperature dependent properties of the soft sticky dipole (SSD) and related single point water models} \author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote} \footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}} @@ -146,27 +146,35 @@ 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 interest within our group. -Up to this point, a detailed look at the model's structure and ion -solvation abilities has been performed.\cite{Ichiye96} In addition, a -thorough investigation of the dynamic properties of SSD was performed -by Chandra and Ichiye focusing on translational and orientational -properties at 298 K.\cite{Ichiye99} This study focuses on determining -the density maximum for SSD utilizing both microcanonical and -isobaric-isothermal ensemble molecular dynamics, while using the -reaction field method for handling long-ranged dipolar interactions. A -reaction field method has been previously implemented in Monte Carlo -simulations by Liu and Ichiye in order to study the static dielectric -constant for the model.\cite{Ichiye96b} This paper will expand the -scope of these original simulations to look on how the reaction field -affects the physical and dynamic properties of SSD systems. +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 +handling of long-ranged interactions. When studying very large +systems, the Ewald summation and even particle-mesh Ewald become +computational burdens with their respective ideal $N^\frac{3}{2}$ and +$N\log N$ calculation scaling orders for $N$ particles.\cite{Darden99} +In applying this water model in these types of systems, it would be +useful to know its properties and behavior with the more +computationally efficient reaction field (RF) technique, and even with +a cutoff that lacks any form of long range correction. This study +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 +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. \section{Methods} @@ -197,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 @@ -212,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 @@ -225,22 +233,23 @@ in energy conservation. There is still the issue of an nine elements long as opposed to 3 or 4 elements for Euler angles and quaternions respectively. System memory has become much less of an issue in recent times, and this has resulted in substantial benefits -in energy conservation. There is still the issue of an additional 5 or -6 additional elements for describing the orientation of each particle, -which will increase dump files substantially. Simply translating the -rotation matrix into its component Euler angles or quaternions for -storage purposes relieves this burden. +in energy conservation. There is still the issue of 5 or 6 additional +elements for describing the orientation of each particle, which will +increase dump files substantially. Simply translating the rotation +matrix into its component Euler angles or quaternions for storage +purposes relieves this burden. The symplectic splitting method allows for Verlet style integration of both linear and angular motion of rigid bodies. In the integration 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 simulation 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} @@ -262,7 +271,9 @@ demonstrated. with the quaternion method showing a slight energy drift over time in the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the energy conservation benefits of the symplectic step method are clearly -demonstrated. +demonstrated. Thus, while maintaining the same degree of energy +conservation, one can take considerably longer time steps, leading to +an overall reduction in computation time. Energy drift in these SSD particle simulations was unnoticeable for time steps up to three femtoseconds. A slight energy drift on the @@ -306,16 +317,9 @@ volume fluctuations dampened out in all but the very c 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} +(below 225 K). In order to further improve statistics, an ensemble +average was accumulated from five separate simulation progressions, +each starting from a different ice crystal. \subsection{Density Behavior} In the initial average density versus temperature plot, the density @@ -892,10 +896,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}