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 |
155 |
|
interest within our group. |
156 |
|
|
157 |
< |
Up to this point, a detailed look at the model's structure and ion |
158 |
< |
solvation abilities has been performed.\cite{Ichiye96} In addition, a |
159 |
< |
thorough investigation of the dynamic properties of SSD was performed |
160 |
< |
by Chandra and Ichiye focusing on translational and orientational |
161 |
< |
properties at 298 K.\cite{Ichiye99} This study focuses on determining |
162 |
< |
the density maximum for SSD utilizing both microcanonical and |
163 |
< |
isobaric-isothermal ensemble molecular dynamics, while using the |
164 |
< |
reaction field method for handling long-ranged dipolar interactions. A |
165 |
< |
reaction field method has been previously implemented in Monte Carlo |
166 |
< |
simulations by Liu and Ichiye in order to study the static dielectric |
167 |
< |
constant for the model.\cite{Ichiye96b} This paper will expand the |
168 |
< |
scope of these original simulations to look on how the reaction field |
169 |
< |
affects the physical and dynamic properties of SSD systems. |
157 |
> |
One of the key limitations of this water model, however, is that it |
158 |
> |
has been parameterized for use with the Ewald Sum technique for the |
159 |
> |
handling of long-ranged interactions. When studying very large |
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 |
|
|
179 |
|
\section{Methods} |
180 |
|
|
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}. |
312 |
< |
|
313 |
< |
\begin{figure} |
314 |
< |
\includegraphics[width=65mm, angle=-90]{1hdense.epsi} |
315 |
< |
\caption{Average density of SSD water at increasing temperatures |
316 |
< |
starting from ice $I_h$ lattice.} |
317 |
< |
\label{dense1} |
318 |
< |
\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 |
|
|