| 23 |
|
|
| 24 |
|
\begin{document} |
| 25 |
|
|
| 26 |
< |
\title{On the temperature dependent structural and transport properties of the soft sticky dipole (SSD) and related single point water models} |
| 26 |
> |
\title{On the temperature dependent properties of the soft sticky dipole (SSD) and related single point water models} |
| 27 |
|
|
| 28 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote} |
| 29 |
|
\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}} |
| 146 |
|
simulations using this model, Ichiye \emph{et al.} reported a |
| 147 |
|
calculation speed up of up to an order of magnitude over other |
| 148 |
|
comparable models while maintaining the structural behavior of |
| 149 |
< |
water.\cite{Ichiye96} In the original molecular dynamics studies of |
| 150 |
< |
SSD, it was shown that it actually improves upon the prediction of |
| 151 |
< |
water's dynamical properties 3 and 4-point models.\cite{Ichiye99} This |
| 149 |
> |
water.\cite{Ichiye96} In the original molecular dynamics studies, it |
| 150 |
> |
was shown that SSD improves on the prediction of many of water's |
| 151 |
> |
dynamical properties over TIP3P and SPC/E.\cite{Ichiye99} This |
| 152 |
|
attractive combination of speed and accurate depiction of solvent |
| 153 |
|
properties makes SSD a model of interest for the simulation of large |
| 154 |
|
scale biological systems, such as membrane phase behavior, a specific |
| 160 |
|
systems, the Ewald summation and even particle-mesh Ewald become |
| 161 |
|
computational burdens with their respective ideal $N^\frac{3}{2}$ and |
| 162 |
|
$N\log N$ calculation scaling orders for $N$ particles.\cite{Darden99} |
| 163 |
+ |
In applying this water model in these types of systems, it would be |
| 164 |
+ |
useful to know its properties and behavior with the more |
| 165 |
+ |
computationally efficient reaction field (RF) technique, and even with |
| 166 |
+ |
a cutoff that lacks any form of long range correction. This study |
| 167 |
+ |
addresses these issues by looking at the structural and transport |
| 168 |
+ |
behavior of SSD over a variety of temperatures, with the purpose of |
| 169 |
+ |
utilizing the RF correction technique. Towards the end, we suggest |
| 170 |
+ |
alterations to the parameters that result in more water-like |
| 171 |
+ |
behavior. It should be noted that in a recent publication, some the |
| 172 |
+ |
original investigators of the SSD water model have put forth |
| 173 |
+ |
adjustments to the original SSD water model to address abnormal |
| 174 |
+ |
density behavior (also observed here), calling the corrected model |
| 175 |
+ |
SSD1.\cite{Ichiye03} This study will consider this new model's |
| 176 |
+ |
behavior as well, and hopefully improve upon its depiction of water |
| 177 |
+ |
under conditions without the Ewald Sum. |
| 178 |
|
|
| 164 |
– |
Up to this point, a detailed look at the model's structure and ion |
| 165 |
– |
solvation abilities has been performed.\cite{Ichiye96} In addition, a |
| 166 |
– |
thorough investigation of the dynamic properties of SSD was performed |
| 167 |
– |
by Chandra and Ichiye focusing on translational and orientational |
| 168 |
– |
properties at 298 K.\cite{Ichiye99} This study focuses on determining |
| 169 |
– |
the density maximum for SSD utilizing both microcanonical and |
| 170 |
– |
isobaric-isothermal ensemble molecular dynamics, while using the |
| 171 |
– |
reaction field method for handling long-ranged dipolar interactions. A |
| 172 |
– |
reaction field method has been previously implemented in Monte Carlo |
| 173 |
– |
simulations by Liu and Ichiye in order to study the static dielectric |
| 174 |
– |
constant for the model.\cite{Ichiye96b} This paper will expand the |
| 175 |
– |
scope of these original simulations to look on how the reaction field |
| 176 |
– |
affects the physical and dynamic properties of SSD systems. |
| 177 |
– |
|
| 179 |
|
\section{Methods} |
| 180 |
|
|
| 181 |
|
As stated previously, in this study the long-range dipole-dipole |
| 205 |
|
to the use of reaction field, simulations were also performed without |
| 206 |
|
a surrounding dielectric and suggestions are proposed on how to make |
| 207 |
|
SSD more compatible with a reaction field. |
| 208 |
< |
|
| 208 |
> |
|
| 209 |
|
Simulations were performed in both the isobaric-isothermal and |
| 210 |
|
microcanonical ensembles. The constant pressure simulations were |
| 211 |
|
implemented using an integral thermostat and barostat as outlined by |
| 212 |
< |
Hoover.\cite{Hoover85,Hoover86} For the constant pressure |
| 213 |
< |
simulations, the \emph{Q} parameter for the was set to 5.0 amu |
| 214 |
< |
\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at |
| 215 |
< |
100 ps. |
| 212 |
> |
Hoover.\cite{Hoover85,Hoover86} All particles were treated as |
| 213 |
> |
non-linear rigid bodies. Vibrational constraints are not necessary in |
| 214 |
> |
simulations of SSD, because there are no explicit hydrogen atoms, and |
| 215 |
> |
thus no molecular vibrational modes need to be considered. |
| 216 |
|
|
| 217 |
|
Integration of the equations of motion was carried out using the |
| 218 |
|
symplectic splitting method proposed by Dullweber \emph{et |
| 220 |
|
deals with poor energy conservation of rigid body systems using |
| 221 |
|
quaternions. While quaternions work well for orientational motion in |
| 222 |
|
alternate ensembles, the microcanonical ensemble has a constant energy |
| 223 |
< |
requirement that is actually quite sensitive to errors in the |
| 224 |
< |
equations of motion. The original implementation of this code utilized |
| 225 |
< |
quaternions for rotational motion propagation; however, a detailed |
| 226 |
< |
investigation showed that they resulted in a steady drift in the total |
| 227 |
< |
energy, something that has been observed by others.\cite{Laird97} |
| 223 |
> |
requirement that is quite sensitive to errors in the equations of |
| 224 |
> |
motion. The original implementation of this code utilized quaternions |
| 225 |
> |
for rotational motion propagation; however, a detailed investigation |
| 226 |
> |
showed that they resulted in a steady drift in the total energy, |
| 227 |
> |
something that has been observed by others.\cite{Laird97} |
| 228 |
|
|
| 229 |
|
The key difference in the integration method proposed by Dullweber |
| 230 |
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 233 |
|
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 234 |
|
quaternions respectively. System memory has become much less of an |
| 235 |
|
issue in recent times, and this has resulted in substantial benefits |
| 236 |
< |
in energy conservation. There is still the issue of an additional 5 or |
| 237 |
< |
6 additional elements for describing the orientation of each particle, |
| 238 |
< |
which will increase dump files substantially. Simply translating the |
| 239 |
< |
rotation matrix into its component Euler angles or quaternions for |
| 240 |
< |
storage purposes relieves this burden. |
| 236 |
> |
in energy conservation. There is still the issue of 5 or 6 additional |
| 237 |
> |
elements for describing the orientation of each particle, which will |
| 238 |
> |
increase dump files substantially. Simply translating the rotation |
| 239 |
> |
matrix into its component Euler angles or quaternions for storage |
| 240 |
> |
purposes relieves this burden. |
| 241 |
|
|
| 242 |
|
The symplectic splitting method allows for Verlet style integration of |
| 243 |
|
both linear and angular motion of rigid bodies. In the integration |
| 244 |
|
method, the orientational propagation involves a sequence of matrix |
| 245 |
|
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 246 |
|
matrix rotations end up being more costly computationally than the |
| 247 |
< |
simpler arithmetic quaternion propagation. On average, a 1000 SSD |
| 248 |
< |
particle simulation shows a 7\% increase in simulation time using the |
| 249 |
< |
symplectic step method in place of quaternions. This cost is more than |
| 250 |
< |
justified when comparing the energy conservation of the two methods as |
| 251 |
< |
illustrated in figure \ref{timestep}. |
| 247 |
> |
simpler arithmetic quaternion propagation. With the same time step, a |
| 248 |
> |
1000 SSD particle simulation shows an average 7\% increase in |
| 249 |
> |
computation time using the symplectic step method in place of |
| 250 |
> |
quaternions. This cost is more than justified when comparing the |
| 251 |
> |
energy conservation of the two methods as illustrated in figure |
| 252 |
> |
\ref{timestep}. |
| 253 |
|
|
| 254 |
|
\begin{figure} |
| 255 |
|
\includegraphics[width=61mm, angle=-90]{timeStep.epsi} |
| 271 |
|
with the quaternion method showing a slight energy drift over time in |
| 272 |
|
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 273 |
|
energy conservation benefits of the symplectic step method are clearly |
| 274 |
< |
demonstrated. |
| 274 |
> |
demonstrated. Thus, while maintaining the same degree of energy |
| 275 |
> |
conservation, one can take considerably longer time steps, leading to |
| 276 |
> |
an overall reduction in computation time. |
| 277 |
|
|
| 278 |
|
Energy drift in these SSD particle simulations was unnoticeable for |
| 279 |
|
time steps up to three femtoseconds. A slight energy drift on the |
| 317 |
|
increment was decreased from 25 K to 10 and then 5 K. The above |
| 318 |
|
equilibration and production times were sufficient in that the system |
| 319 |
|
volume fluctuations dampened out in all but the very cold simulations |
| 320 |
< |
(below 225 K). In order to further improve statistics, five separate |
| 321 |
< |
simulation progressions were performed, and the averaged results from |
| 322 |
< |
the $I_h$ melting simulations are shown in figure \ref{dense1}. |
| 319 |
< |
|
| 320 |
< |
\begin{figure} |
| 321 |
< |
\includegraphics[width=65mm, angle=-90]{1hdense.epsi} |
| 322 |
< |
\caption{Average density of SSD water at increasing temperatures |
| 323 |
< |
starting from ice $I_h$ lattice.} |
| 324 |
< |
\label{dense1} |
| 325 |
< |
\end{figure} |
| 320 |
> |
(below 225 K). In order to further improve statistics, an ensemble |
| 321 |
> |
average was accumulated from five separate simulation progressions, |
| 322 |
> |
each starting from a different ice crystal. |
| 323 |
|
|
| 324 |
|
\subsection{Density Behavior} |
| 325 |
|
In the initial average density versus temperature plot, the density |
| 896 |
|
simulations of biochemical systems. |
| 897 |
|
|
| 898 |
|
\section{Acknowledgments} |
| 899 |
< |
The authors would like to thank the National Science Foundation for |
| 900 |
< |
funding under grant CHE-0134881. Computation time was provided by the |
| 901 |
< |
Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR |
| 902 |
< |
00 79647. |
| 899 |
> |
Support for this project was provided by the National Science |
| 900 |
> |
Foundation under grant CHE-0134881. Computation time was provided by |
| 901 |
> |
the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant |
| 902 |
> |
DMR 00 79647. |
| 903 |
|
|
| 904 |
|
\bibliographystyle{jcp} |
| 905 |
|
|
