ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/ssdePaper/nptSSD.tex

Comparing trunk/ssdePaper/nptSSD.tex (file contents):
Revision 863 by chrisfen, Thu Nov 13 15:55:20 2003 UTC vs.
Revision 1019 by chrisfen, Wed Feb 4 20:19:21 2004 UTC

#	Line 21 | Line 21
21
22		\begin{document}
23
24	<	\title{On the temperature dependent properties of the soft sticky dipole (SSD) and related single point water models}
24	>	\title{On the structural and transport properties of the soft sticky
25	>	dipole (SSD) and related single point water models}
26
27		\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
28		Department of Chemistry and Biochemistry\\ University of Notre Dame\\
#	Line 32 | Line 33 | NVE and NPT molecular dynamics simulations were perfor
33		\maketitle
34
35		\begin{abstract}
36	<	NVE and NPT molecular dynamics simulations were performed in order to
37	<	investigate the density maximum and temperature dependent transport
38	<	for SSD and related water models, both with and without the use of
39	<	reaction field. The constant pressure simulations of the melting of
40	<	both $I_h$ and $I_c$ ice showed a density maximum near 260 K. In most
41	<	cases, the calculated densities were significantly lower than the
42	<	densities calculated in simulations of other water models. Analysis of
43	<	particle diffusion showed SSD to capture the transport properties of
36	>	The density maximum and temperature dependence of the self-diffusion
37	>	constant were investigated for the soft sticky dipole (SSD) water
38	>	model and two related re-parameterizations of this single-point model.
39	>	A combination of microcanonical and isobaric-isothermal molecular
40	>	dynamics simulations were used to calculate these properties, both
41	>	with and without the use of reaction field to handle long-range
42	>	electrostatics.  The isobaric-isothermal (NPT) simulations of the
43	>	melting of both ice-$I_h$ and ice-$I_c$ showed a density maximum near
44	>	260 K.  In most cases, the use of the reaction field resulted in
45	>	calculated densities which were were significantly lower than
46	>	experimental densities.  Analysis of self-diffusion constants shows
47	>	that the original SSD model captures the transport properties of
48		experimental water very well in both the normal and super-cooled
49	<	liquid regimes. In order to correct the density behavior, SSD was
50	<	reparameterized for use both with and without a long-range interaction
51	<	correction, SSD/RF and SSD/E respectively. Compared to the density
52	<	corrected version of SSD (SSD1), these modified models were shown to
53	<	maintain or improve upon the structural and transport properties.
49	>	liquid regimes.  We also present our re-parameterized versions of SSD
50	>	for use both with the reaction field or without any long-range
51	>	electrostatic corrections.  These are called the SSD/RF and SSD/E
52	>	models respectively.  These modified models were shown to maintain or
53	>	improve upon the experimental agreement with the structural and
54	>	transport properties that can be obtained with either the original SSD
55	>	or the density corrected version of the original model (SSD1).
56	>	Additionally, a novel low-density ice structure is presented
57	>	which appears to be the most stable ice structure for the entire SSD
58	>	family.
59		\end{abstract}
60
61		\newpage
#	Line 60 | Line 70 | systems is the proper depiction of water and water sol
70		\section{Introduction}
71
72		One of the most important tasks in the simulation of biochemical
73	<	systems is the proper depiction of water and water solvation. In fact,
74	<	the bulk of the calculations performed in solvated simulations are of
75	<	interactions with or between solvent molecules. Thus, the outcomes of
76	<	these types of simulations are highly dependent on the physical
77	<	properties of water, both as individual molecules and in clusters or
78	<	bulk. Due to the fact that explicit solvent accounts for a massive
79	<	portion of the calculations, it necessary to simplify the solvent to
70	<	some extent in order to complete simulations in a reasonable amount of
71	<	time. In the case of simulating water in biomolecular studies, the
72	<	balance between accurate properties and computational efficiency is
73	<	especially delicate, and it has resulted in a variety of different
74	<	water models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these
75	<	models predict specific properties more accurately than their
76	<	predecessors, but often at the cost of other properties or of computer
77	<	time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P improves
78	<	upon the structural and transport properties of water relative to the
79	<	previous TIP models, yet this comes at a greater than 50\% increase in
80	<	computational cost.\cite{Jorgensen01,Jorgensen00} One recently
81	<	developed model that succeeds in both retaining the accuracy of system
82	<	properties and simplifying calculations to increase computational
83	<	efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96}
73	>	systems is the proper depiction of the aqueous environment of the
74	>	molecules of interest.  In some cases (such as in the simulation of
75	>	phospholipid bilayers), the majority of the calculations that are
76	>	performed involve interactions with or between solvent molecules.
77	>	Thus, the properties one may observe in biochemical simulations are
78	>	going to be highly dependent on the physical properties of the water
79	>	model that is chosen.
80
81	<	The Soft Sticky Dipole (SSD)\ water model was developed by Ichiye
82	<	\emph{et al.} as a modified form of the hard-sphere water model
83	<	proposed by Bratko, Blum, and Luzar.\cite{Bratko85,Bratko95} SSD
84	<	consists of a single point dipole with a Lennard-Jones core and a
85	<	sticky potential that directs the particles to assume the proper
86	<	hydrogen bond orientation in the first solvation shell. Thus, the
87	<	interaction between two SSD water molecules \emph{i} and \emph{j} is
88	<	given by the potential
81	>	There is an especially delicate balance between computational
82	>	efficiency and the ability of the water model to accurately predict
83	>	the properties of bulk
84	>	water.\cite{Jorgensen83,Berendsen87,Jorgensen00} For example, the
85	>	TIP5P model improves on the structural and transport properties of
86	>	water relative to the previous TIP models, yet this comes at a greater
87	>	than 50\% increase in computational
88	>	cost.\cite{Jorgensen01,Jorgensen00}
89	>
90	>	One recently developed model that largely succeeds in retaining the
91	>	accuracy of bulk properties while greatly reducing the computational
92	>	cost is the Soft Sticky Dipole (SSD) water
93	>	model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The SSD model was
94	>	developed by Ichiye \emph{et al.} as a modified form of the
95	>	hard-sphere water model proposed by Bratko, Blum, and
96	>	Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model which
97	>	has an interaction site that is both a point dipole along with a
98	>	Lennard-Jones core.  However, since the normal aligned and
99	>	anti-aligned geometries favored by point dipoles are poor mimics of
100	>	local structure in liquid water, a short ranged ``sticky'' potential
101	>	is also added.  The sticky potential directs the molecules to assume
102	>	the proper hydrogen bond orientation in the first solvation
103	>	shell.
104	>
105	>	The interaction between two SSD water molecules \emph{i} and \emph{j}
106	>	is given by the potential
107		\begin{equation}
108		u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
109	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
109	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)\ +
110		u_{ij}^{sp}
111	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
111	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j),
112		\end{equation}
113	<	where the $\mathbf{r}_{ij}$ is the position vector between molecules
114	<	\emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and
115	<	$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
116	<	orientations of the respective molecules. The Lennard-Jones, dipole,
117	<	and sticky parts of the potential are giving by the following
104	<	equations:
113	>	where the ${\bf r}_{ij}$ is the position vector between molecules
114	>	\emph{i} and \emph{j} with magnitude $r_{ij}$, and
115	>	${\bf \Omega}_i$ and ${\bf \Omega}_j$ represent the orientations of
116	>	the two molecules. The Lennard-Jones and dipole interactions are given
117	>	by the following familiar forms:
118		\begin{equation}
119	<	u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
119	>	u_{ij}^{LJ}(r_{ij}) = 4\epsilon
120	>	\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]
121	>	\ ,
122		\end{equation}
123	+	and
124		\begin{equation}
125	<	u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ ,
125	>	u_{ij}^{dp} = \frac{|\mu_i||\mu_j|}{4 \pi \epsilon_0 r_{ij}^3} \left(
126	>	\hat{\bf u}_i \cdot \hat{\bf u}_j - 3(\hat{\bf u}_i\cdot\hat{\bf
127	>	r}_{ij})(\hat{\bf u}_j\cdot\hat{\bf r}_{ij}) \right)\ ,
128		\end{equation}
129	+	where $\hat{\bf u}_i$ and $\hat{\bf u}_j$ are the unit vectors along
130	+	the dipoles of molecules $i$ and $j$ respectively. $|\mu_i|$ and
131	+	$|\mu_j|$ are the strengths of the dipole moments, and $\hat{\bf
132	+	r}_{ij}$ is the unit vector pointing from molecule $j$ to molecule
133	+	$i$.
134	+
135	+	The sticky potential is somewhat less familiar:
136		\begin{equation}
137		u_{ij}^{sp}
138	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) =
139	<	\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) + s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
138	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
139	>	\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)
140	>	+ s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf
141	>	\Omega}_j)]\ .
142	>	\label{stickyfunction}
143		\end{equation}
144	<	where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
145	<	unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
146	<	$\nu_0$ scales the strength of the overall sticky potential, and $s$
147	<	and $s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
148	<	functions take the following forms:
144	>	Here, $\nu_0$ is a strength parameter for the sticky potential, and
145	>	$s$ and $s^\prime$ are cubic switching functions which turn off the
146	>	sticky interaction beyond the first solvation shell. The $w$ function
147	>	can be thought of as an attractive potential with tetrahedral
148	>	geometry:
149		\begin{equation}
150	<	w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
150	>	w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
151		\end{equation}
152	+	while the $w^\prime$ function counters the normal aligned and
153	+	anti-aligned structures favored by point dipoles:
154		\begin{equation}
155	<	w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
155	>	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,
156		\end{equation}
157	<	where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
158	<	term that promotes hydrogen bonding orientations within the first
159	<	solvation shell, and $w^\prime$ is a dipolar repulsion term that
160	<	repels unrealistic dipolar arrangements within the first solvation
161	<	shell. A more detailed description of the functional parts and
162	<	variables in this potential can be found in other
163	<	articles.\cite{Ichiye96,Ichiye99}
157	>	It should be noted that $w$ is proportional to the sum of the $Y_3^2$
158	>	and $Y_3^{-2}$ spherical harmonics (a linear combination which
159	>	enhances the tetrahedral geometry for hydrogen bonded structures),
160	>	while $w^\prime$ is a purely empirical function.  A more detailed
161	>	description of the functional parts and variables in this potential
162	>	can be found in the original SSD
163	>	articles.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03}
164
165	<	Being that this is a one-site point dipole model, the actual force
166	<	calculations are simplified significantly. In the original Monte Carlo
167	<	simulations using this model, Ichiye \emph{et al.} reported an
168	<	increase in calculation efficiency of up to an order of magnitude over
169	<	other comparable models, while maintaining the structural behavior of
170	<	water.\cite{Ichiye96} In the original molecular dynamics studies, it
171	<	was shown that SSD improves on the prediction of many of water's
172	<	dynamical properties over TIP3P and SPC/E.\cite{Ichiye99} This
173	<	attractive combination of speed and accurate depiction of solvent
174	<	properties makes SSD a model of interest for the simulation of large
175	<	scale biological systems, such as membrane phase behavior.
165	>	Since SSD is a single-point {\it dipolar} model, the force
166	>	calculations are simplified significantly relative to the standard
167	>	{\it charged} multi-point models. In the original Monte Carlo
168	>	simulations using this model, Ichiye {\it et al.} reported that using
169	>	SSD decreased computer time by a factor of 6-7 compared to other
170	>	models.\cite{Ichiye96} What is most impressive is that this savings
171	>	did not come at the expense of accurate depiction of the liquid state
172	>	properties.  Indeed, SSD maintains reasonable agreement with the Soper
173	>	data for the structural features of liquid
174	>	water.\cite{Soper86,Ichiye96} Additionally, the dynamical properties
175	>	exhibited by SSD agree with experiment better than those of more
176	>	computationally expensive models (like TIP3P and
177	>	SPC/E).\cite{Ichiye99} The combination of speed and accurate depiction
178	>	of solvent properties makes SSD a very attractive model for the
179	>	simulation of large scale biochemical simulations.
180
181	<	One of the key limitations of this water model, however, is that it
182	<	has been parameterized for use with the Ewald Sum technique for the
183	<	handling of long-ranged interactions.  When studying very large
184	<	systems, the Ewald summation and even particle-mesh Ewald become
185	<	computational burdens, with their respective ideal $N^\frac{3}{2}$ and
186	<	$N\log N$ calculation scaling orders for $N$ particles.\cite{Darden99}
187	<	In applying this water model in these types of systems, it would be
188	<	useful to know its properties and behavior with the more
189	<	computationally efficient reaction field (RF) technique, and even with
190	<	a cutoff that lacks any form of long-range correction. This study
191	<	addresses these issues by looking at the structural and transport
192	<	behavior of SSD over a variety of temperatures with the purpose of
193	<	utilizing the RF correction technique. We then suggest alterations to
194	<	the parameters that result in more water-like behavior. It should be
195	<	noted that in a recent publication, some of the original investigators of
196	<	the SSD water model have put forth adjustments to the SSD water model
197	<	to address abnormal density behavior (also observed here), calling the
198	<	corrected model SSD1.\cite{Ichiye03} This study will make comparisons
199	<	with SSD1's behavior with the goal of improving upon the
200	<	depiction of water under conditions without the Ewald Sum.
181	>	One feature of the SSD model is that it was parameterized for use with
182	>	the Ewald sum to handle long-range interactions.  This would normally
183	>	be the best way of handling long-range interactions in systems that
184	>	contain other point charges.  However, our group has recently become
185	>	interested in systems with point dipoles as mimics for neutral, but
186	>	polarized regions on molecules (e.g. the zwitterionic head group
187	>	regions of phospholipids).  If the system of interest does not contain
188	>	point charges, the Ewald sum and even particle-mesh Ewald become
189	>	computational bottlenecks.  Their respective ideal $N^\frac{3}{2}$ and
190	>	$N\log N$ calculation scaling orders for $N$ particles can become
191	>	prohibitive when $N$ becomes large.\cite{Darden99} In applying this
192	>	water model in these types of systems, it would be useful to know its
193	>	properties and behavior under the more computationally efficient
194	>	reaction field (RF) technique, or even with a simple cutoff. This
195	>	study addresses these issues by looking at the structural and
196	>	transport behavior of SSD over a variety of temperatures with the
197	>	purpose of utilizing the RF correction technique.  We then suggest
198	>	modifications to the parameters that result in more realistic bulk
199	>	phase behavior.  It should be noted that in a recent publication, some
200	>	of the original investigators of the SSD water model have suggested
201	>	adjustments to the SSD water model to address abnormal density
202	>	behavior (also observed here), calling the corrected model
203	>	SSD1.\cite{Ichiye03} In what follows, we compare our
204	>	reparamaterization of SSD with both the original SSD and SSD1 models
205	>	with the goal of improving the bulk phase behavior of an SSD-derived
206	>	model in simulations utilizing the Reaction Field.
207
208		\section{Methods}
209
210	<	As stated previously, the long-range dipole-dipole interactions were
211	<	accounted for in this study by using the reaction field method. The
212	<	magnitude of the reaction field acting on dipole \emph{i} is given by
210	>	Long-range dipole-dipole interactions were accounted for in this study
211	>	by using either the reaction field method or by resorting to a simple
212	>	cubic switching function at a cutoff radius.  The reaction field
213	>	method was actually first used in Monte Carlo simulations of liquid
214	>	water.\cite{Barker73} Under this method, the magnitude of the reaction
215	>	field acting on dipole $i$ is
216		\begin{equation}
217		\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
218	<	\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} \boldsymbol{\mu}_{j} f(r_{ij})\  ,
218	>	\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} f(r_{ij})\  ,
219		\label{rfequation}
220		\end{equation}
221		where $\mathcal{R}$ is the cavity defined by the cutoff radius
222		($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
223	<	system (80 in this case), $\boldsymbol{\mu}_{j}$ is the dipole moment
224	<	vector of particle \emph{j}, and $f(r_{ij})$ is a cubic switching
223	>	system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole
224	>	moment vector of particle $j$ and $f(r_{ij})$ is a cubic switching
225		function.\cite{AllenTildesley} The reaction field contribution to the
226	<	total energy by particle \emph{i} is given by
227	<	$-\frac{1}{2}\boldsymbol{\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque
228	<	on dipole \emph{i} by
229	<	$\boldsymbol{\mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley} Use
230	<	of reaction field is known to alter the orientational dynamic
231	<	properties, such as the dielectric relaxation time, based on changes
232	<	in the length of the cutoff radius.\cite{Berendsen98} This variable
233	<	behavior makes reaction field a less attractive method than other
234	<	methods, like the Ewald summation; however, for the simulation of
235	<	large-scale systems, the computational cost benefit of reaction field
236	<	is dramatic. To address some of the dynamical property alterations due
237	<	to the use of reaction field, simulations were also performed without
238	<	a surrounding dielectric and suggestions are presented on how to make
239	<	SSD more accurate both with and without a reaction field.
226	>	total energy by particle $i$ is given by $-\frac{1}{2}{\bf
227	>	\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf
228	>	\mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley}  Use of the reaction
229	>	field is known to alter the bulk orientational properties, such as the
230	>	dielectric relaxation time.  There is particular sensitivity of this
231	>	property on changes in the length of the cutoff
232	>	radius.\cite{Berendsen98} This variable behavior makes reaction field
233	>	a less attractive method than the Ewald sum.  However, for very large
234	>	systems, the computational benefit of reaction field is dramatic.
235	>
236	>	We have also performed a companion set of simulations {\it without} a
237	>	surrounding dielectric (i.e. using a simple cubic switching function
238	>	at the cutoff radius) and as a result we have two reparamaterizations
239	>	of SSD which could be used either with or without the Reaction Field
240	>	turned on.
241
242	<	Simulations were performed in both the isobaric-isothermal and
243	<	microcanonical ensembles. The constant pressure simulations were
244	<	implemented using an integral thermostat and barostat as outlined by
245	<	Hoover.\cite{Hoover85,Hoover86} All particles were treated as
246	<	non-linear rigid bodies. Vibrational constraints are not necessary in
247	<	simulations of SSD, because there are no explicit hydrogen atoms, and
248	<	thus no molecular vibrational modes need to be considered.
242	>	Simulations to obtain the preferred density were performed in the
243	>	isobaric-isothermal (NPT) ensemble, while all dynamical properties
244	>	were obtained from microcanonical (NVE) simulations done at densities
245	>	matching the NPT density for a particular target temperature.  The
246	>	constant pressure simulations were implemented using an integral
247	>	thermostat and barostat as outlined by Hoover.\cite{Hoover85,Hoover86}
248	>	All molecules were treated as non-linear rigid bodies. Vibrational
249	>	constraints are not necessary in simulations of SSD, because there are
250	>	no explicit hydrogen atoms, and thus no molecular vibrational modes
251	>	need to be considered.
252
253		Integration of the equations of motion was carried out using the
254	<	symplectic splitting method proposed by Dullweber \emph{et
255	<	al.}\cite{Dullweber1997} The reason for this integrator selection
256	<	deals with poor energy conservation of rigid body systems using
257	<	quaternions. While quaternions work well for orientational motion in
258	<	alternate ensembles, the microcanonical ensemble has a constant energy
259	<	requirement that is quite sensitive to errors in the equations of
260	<	motion. The original implementation of this code utilized quaternions
261	<	for rotational motion propagation; however, a detailed investigation
262	<	showed that they resulted in a steady drift in the total energy,
216	<	something that has been observed by others.\cite{Laird97}
254	>	symplectic splitting method proposed by Dullweber {\it et
255	>	al.}\cite{Dullweber1997} Our reason for selecting this integrator
256	>	centers on poor energy conservation of rigid body dynamics using
257	>	traditional quaternion integration.\cite{Evans77,Evans77b} While quaternions
258	>	may work well for orientational motion under NVT or NPT integrators,
259	>	our limits on energy drift in the microcanonical ensemble were quite
260	>	strict, and the drift under quaternions was substantially greater than
261	>	in the symplectic splitting method.  This steady drift in the total
262	>	energy has also been observed by Kol {\it et al.}\cite{Laird97}
263
264		The key difference in the integration method proposed by Dullweber
265		\emph{et al.} is that the entire rotation matrix is propagated from
266	<	one time step to the next. In the past, this would not have been as
267	<	feasible an option, being that the rotation matrix for a single body is
268	<	nine elements long as opposed to 3 or 4 elements for Euler angles and
269	<	quaternions respectively. System memory has become much less of an
224	<	issue in recent times, and this has resulted in substantial benefits
225	<	in energy conservation. There is still the issue of 5 or 6 additional
226	<	elements for describing the orientation of each particle, which will
227	<	increase dump files substantially. Simply translating the rotation
228	<	matrix into its component Euler angles or quaternions for storage
229	<	purposes relieves this burden.
266	>	one time step to the next.  The additional memory required by the
267	>	algorithm is inconsequential on modern computers, and translating the
268	>	rotation matrix into quaternions for storage purposes makes trajectory
269	>	data quite compact.
270
271		The symplectic splitting method allows for Verlet style integration of
272	<	both linear and angular motion of rigid bodies. In this integration
273	<	method, the orientational propagation involves a sequence of matrix
274	<	evaluations to update the rotation matrix.\cite{Dullweber1997} These
275	<	matrix rotations are more costly computationally than the simpler
276	<	arithmetic quaternion propagation. With the same time step, a 1000 SSD
277	<	particle simulation shows an average 7\% increase in computation time
278	<	using the symplectic step method in place of quaternions. This cost is
279	<	more than justified when comparing the energy conservation of the two
280	<	methods as illustrated in figure \ref{timestep}.
272	>	both translational and orientational motion of rigid bodies. In this
273	>	integration method, the orientational propagation involves a sequence
274	>	of matrix evaluations to update the rotation
275	>	matrix.\cite{Dullweber1997} These matrix rotations are more costly
276	>	than the simpler arithmetic quaternion propagation. With the same time
277	>	step, a 1000 SSD particle simulation shows an average 7\% increase in
278	>	computation time using the symplectic step method in place of
279	>	quaternions. The additional expense per step is justified when one
280	>	considers the ability to use time steps that are nearly twice as large
281	>	under symplectic splitting than would be usable under quaternion
282	>	dynamics.  The energy conservation of the two methods using a number
283	>	of different time steps is illustrated in figure
284	>	\ref{timestep}.
285
286		\begin{figure}
287		\begin{center}
288		\epsfxsize=6in
289		\epsfbox{timeStep.epsi}
290	<	\caption{Energy conservation using quaternion based integration versus
290	>	\caption{Energy conservation using both quaternion based integration and
291		the symplectic step method proposed by Dullweber \emph{et al.} with
292	<	increasing time step. The larger time step plots are shifted up from
293	<	the true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.}
292	>	increasing time step. The larger time step plots are shifted from the
293	>	true energy baseline (that of $\Delta t$ = 0.1 fs) for clarity.}
294		\label{timestep}
295		\end{center}
296		\end{figure}
297
298		In figure \ref{timestep}, the resulting energy drift at various time
299		steps for both the symplectic step and quaternion integration schemes
300	<	is compared. All of the 1000 SSD particle simulations started with the
301	<	same configuration, and the only difference was the method used to
302	<	handle rotational motion. At time steps of 0.1 and 0.5 fs, both
303	<	methods for propagating particle rotation conserve energy fairly well,
304	<	with the quaternion method showing a slight energy drift over time in
305	<	the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
306	<	energy conservation benefits of the symplectic step method are clearly
307	<	demonstrated. Thus, while maintaining the same degree of energy
308	<	conservation, one can take considerably longer time steps, leading to
309	<	an overall reduction in computation time.
300	>	is compared.  All of the 1000 SSD particle simulations started with
301	>	the same configuration, and the only difference was the method used to
302	>	handle orientational motion. At time steps of 0.1 and 0.5 fs, both
303	>	methods for propagating the orientational degrees of freedom conserve
304	>	energy fairly well, with the quaternion method showing a slight energy
305	>	drift over time in the 0.5 fs time step simulation. At time steps of 1
306	>	and 2 fs, the energy conservation benefits of the symplectic step
307	>	method are clearly demonstrated. Thus, while maintaining the same
308	>	degree of energy conservation, one can take considerably longer time
309	>	steps, leading to an overall reduction in computation time.
310
311		Energy drift in the symplectic step simulations was unnoticeable for
312	<	time steps up to three femtoseconds. A slight energy drift on the
312	>	time steps up to 3 fs. A slight energy drift on the
313		order of 0.012 kcal/mol per nanosecond was observed at a time step of
314	<	four femtoseconds, and as expected, this drift increases dramatically
315	<	with increasing time step. To insure accuracy in the constant energy
314	>	4 fs, and as expected, this drift increases dramatically
315	>	with increasing time step. To insure accuracy in our microcanonical
316		simulations, time steps were set at 2 fs and kept at this value for
317		constant pressure simulations as well.
318
319	<	Ice crystals in both the $I_h$ and $I_c$ lattices were generated as
320	<	starting points for all simulations. The $I_h$ crystals were formed by
321	<	first arranging the centers of mass of the SSD particles into a
322	<	``hexagonal'' ice lattice of 1024 particles. Because of the crystal
323	<	structure of $I_h$ ice, the simulation box assumed a rectangular shape
324	<	with an edge length ratio of approximately
319	>	Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
320	>	were generated as starting points for all simulations. The $I_h$
321	>	crystals were formed by first arranging the centers of mass of the SSD
322	>	particles into a ``hexagonal'' ice lattice of 1024 particles. Because
323	>	of the crystal structure of $I_h$ ice, the simulation box assumed an
324	>	orthorhombic shape with an edge length ratio of approximately
325		1.00$\times$1.06$\times$1.23. The particles were then allowed to
326		orient freely about fixed positions with angular momenta randomized at
327		400 K for varying times. The rotational temperature was then scaled
#	Line 297 | Line 341 | constant pressure and temperature dynamics. During mel
341		\section{Results and discussion}
342
343		Melting studies were performed on the randomized ice crystals using
344	<	constant pressure and temperature dynamics. During melting
345	<	simulations, the melting transition and the density maximum can both
346	<	be observed, provided that the density maximum occurs in the liquid
347	<	and not the supercooled regime. An ensemble average from five separate
348	<	melting simulations was acquired, each starting from different ice
349	<	crystals generated as described previously. All simulations were
350	<	equilibrated for 100 ps prior to a 200 ps data collection run at each
351	<	temperature setting. The temperature range of study spanned from 25 to
352	<	400 K, with a maximum degree increment of 25 K. For regions of
353	<	interest along this stepwise progression, the temperature increment
354	<	was decreased from 25 K to 10 and 5 K. The above equilibration and
355	<	production times were sufficient in that the system volume
356	<	fluctuations dampened out in all but the very cold simulations (below
313	<	225 K).
344	>	isobaric-isothermal (NPT) dynamics. During melting simulations, the
345	>	melting transition and the density maximum can both be observed,
346	>	provided that the density maximum occurs in the liquid and not the
347	>	supercooled regime. An ensemble average from five separate melting
348	>	simulations was acquired, each starting from different ice crystals
349	>	generated as described previously. All simulations were equilibrated
350	>	for 100 ps prior to a 200 ps data collection run at each temperature
351	>	setting. The temperature range of study spanned from 25 to 400 K, with
352	>	a maximum degree increment of 25 K. For regions of interest along this
353	>	stepwise progression, the temperature increment was decreased from 25
354	>	K to 10 and 5 K.  The above equilibration and production times were
355	>	sufficient in that fluctuations in the volume autocorrelation function
356	>	were damped out in all simulations in under 20 ps.
357
358		\subsection{Density Behavior}
316	–	Initial simulations focused on the original SSD water model, and an
317	–	average density versus temperature plot is shown in figure
318	–	\ref{dense1}. Note that the density maximum when using a reaction
319	–	field appears between 255 and 265 K, where the calculated densities
320	–	within this range were nearly indistinguishable. The greater certainty
321	–	of the average value at 260 K makes a good argument for the actual
322	–	density maximum residing at this midpoint value. Figure \ref{dense1}
323	–	was constructed using ice $I_h$ crystals for the initial
324	–	configuration; though not pictured, the simulations starting from ice
325	–	$I_c$ crystal configurations showed similar results, with a
326	–	liquid-phase density maximum in this same region (between 255 and 260
327	–	K). In addition, the $I_c$ crystals are more fragile than the $I_h$
328	–	crystals, leading to deformation into a dense glassy state at lower
329	–	temperatures. This resulted in an overall low temperature density
330	–	maximum at 200 K, while still retaining a liquid state density maximum
331	–	in common with the $I_h$ simulations.
359
360	+	Our initial simulations focused on the original SSD water model, and
361	+	an average density versus temperature plot is shown in figure
362	+	\ref{dense1}. Note that the density maximum when using a reaction
363	+	field appears between 255 and 265 K.  There were smaller fluctuations
364	+	in the density at 260 K than at either 255 or 265, so we report this
365	+	value as the location of the density maximum. Figure \ref{dense1} was
366	+	constructed using ice $I_h$ crystals for the initial configuration;
367	+	though not pictured, the simulations starting from ice $I_c$ crystal
368	+	configurations showed similar results, with a liquid-phase density
369	+	maximum in this same region (between 255 and 260 K).
370	+
371		\begin{figure}
372		\begin{center}
373		\epsfxsize=6in
374		\epsfbox{denseSSD.eps}
375	<	\caption{Density versus temperature for TIP4P,\cite{Jorgensen98b}
376	<	TIP3P,\cite{Jorgensen98b} SPC/E,\cite{Clancy94} SSD without Reaction
377	<	Field, SSD, and experiment.\cite{CRC80} The arrows indicate the
378	<	change in densities observed when turning off the reaction field. The
379	<	the lower than expected densities for the SSD model were what
380	<	prompted the original reparameterization.\cite{Ichiye03}}
375	>	\caption{Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}],
376	>	TIP3P [Ref. \citen{Jorgensen98b}], SPC/E [Ref. \citen{Clancy94}], SSD
377	>	without Reaction Field, SSD, and experiment [Ref. \citen{CRC80}]. The
378	>	arrows indicate the change in densities observed when turning off the
379	>	reaction field. The the lower than expected densities for the SSD
380	>	model were what prompted the original reparameterization of SSD1
381	>	[Ref. \citen{Ichiye03}].}
382		\label{dense1}
383		\end{center}
384		\end{figure}
385
386	<	The density maximum for SSD actually compares quite favorably to other
387	<	simple water models. Figure \ref{dense1} also shows calculated
388	<	densities of several other models and experiment obtained from other
386	>	The density maximum for SSD compares quite favorably to other simple
387	>	water models. Figure \ref{dense1} also shows calculated densities of
388	>	several other models and experiment obtained from other
389		sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
390	<	models, SSD has results closest to the experimentally observed water
391	<	density maximum. Of the listed water models, TIP4P has a density
392	<	maximum behavior most like that seen in SSD. Though not included in
393	<	this particular plot, it is useful to note that TIP5P has a water
394	<	density maximum nearly identical to experiment.
390	>	models, SSD has a temperature closest to the experimentally observed
391	>	density maximum. Of the {\it charge-based} models in
392	>	Fig. \ref{dense1}, TIP4P has a density maximum behavior most like that
393	>	seen in SSD. Though not included in this plot, it is useful
394	>	to note that TIP5P has a density maximum nearly identical to the
395	>	experimentally measured temperature.
396
397	<	It has been observed that densities are dependent on the cutoff radius
398	<	used for a variety of water models in simulations both with and
399	<	without the use of reaction field.\cite{Berendsen98} In order to
400	<	address the possible affect of cutoff radius, simulations were
401	<	performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the
402	<	previous SSD simulations, all performed with a cutoff of 9.0 \AA. All
403	<	of the resulting densities overlapped within error and showed no
404	<	significant trend toward lower or higher densities as a function of
405	<	cutoff radius, for simulations both with and without reaction
406	<	field. These results indicate that there is no major benefit in
407	<	choosing a longer cutoff radius in simulations using SSD. This is
408	<	advantageous in that the use of a longer cutoff radius results in
409	<	significant increases in the time required to obtain a single
370	<	trajectory.
397	>	It has been observed that liquid state densities in water are
398	>	dependent on the cutoff radius used both with and without the use of
399	>	reaction field.\cite{Berendsen98} In order to address the possible
400	>	effect of cutoff radius, simulations were performed with a dipolar
401	>	cutoff radius of 12.0 \AA\ to complement the previous SSD simulations,
402	>	all performed with a cutoff of 9.0 \AA. All of the resulting densities
403	>	overlapped within error and showed no significant trend toward lower
404	>	or higher densities as a function of cutoff radius, for simulations
405	>	both with and without reaction field. These results indicate that
406	>	there is no major benefit in choosing a longer cutoff radius in
407	>	simulations using SSD. This is advantageous in that the use of a
408	>	longer cutoff radius results in a significant increase in the time
409	>	required to obtain a single trajectory.
410
411		The key feature to recognize in figure \ref{dense1} is the density
412		scaling of SSD relative to other common models at any given
413	<	temperature. Note that the SSD model assumes a lower density than any
414	<	of the other listed models at the same pressure, behavior which is
415	<	especially apparent at temperatures greater than 300 K. Lower than
416	<	expected densities have been observed for other systems using a
417	<	reaction field for long-range electrostatic interactions, so the most
418	<	likely reason for the significantly lower densities seen in these
419	<	simulations is the presence of the reaction
420	<	field.\cite{Berendsen98,Nezbeda02} In order to test the effect of the
421	<	reaction field on the density of the systems, the simulations were
422	<	repeated without a reaction field present. The results of these
423	<	simulations are also displayed in figure \ref{dense1}. Without
424	<	reaction field, the densities increase considerably to more
425	<	experimentally reasonable values, especially around the freezing point
426	<	of liquid water. The shape of the curve is similar to the curve
427	<	produced from SSD simulations using reaction field, specifically the
428	<	rapidly decreasing densities at higher temperatures; however, a shift
429	<	in the density maximum location, down to 245 K, is observed. This is a
430	<	more accurate comparison to the other listed water models, in that no
431	<	long range corrections were applied in those
393	<	simulations.\cite{Clancy94,Jorgensen98b} However, even without a
413	>	temperature. SSD assumes a lower density than any of the other listed
414	>	models at the same pressure, behavior which is especially apparent at
415	>	temperatures greater than 300 K. Lower than expected densities have
416	>	been observed for other systems using a reaction field for long-range
417	>	electrostatic interactions, so the most likely reason for the
418	>	significantly lower densities seen in these simulations is the
419	>	presence of the reaction field.\cite{Berendsen98,Nezbeda02} In order
420	>	to test the effect of the reaction field on the density of the
421	>	systems, the simulations were repeated without a reaction field
422	>	present. The results of these simulations are also displayed in figure
423	>	\ref{dense1}. Without the reaction field, the densities increase
424	>	to more experimentally reasonable values, especially around the
425	>	freezing point of liquid water. The shape of the curve is similar to
426	>	the curve produced from SSD simulations using reaction field,
427	>	specifically the rapidly decreasing densities at higher temperatures;
428	>	however, a shift in the density maximum location, down to 245 K, is
429	>	observed. This is a more accurate comparison to the other listed water
430	>	models, in that no long range corrections were applied in those
431	>	simulations.\cite{Clancy94,Jorgensen98b} However, even without the
432		reaction field, the density around 300 K is still significantly lower
433		than experiment and comparable water models. This anomalous behavior
434	<	was what lead Ichiye \emph{et al.} to recently reparameterize SSD and
435	<	make SSD1.\cite{Ichiye03} In discussing potential adjustments later in
436	<	this paper, all comparisons were performed with this new model.
434	>	was what lead Ichiye {\it et al.} to recently reparameterize
435	>	SSD.\cite{Ichiye03} Throughout the remainder of the paper our
436	>	reparamaterizations of SSD will be compared with the newer SSD1 model.
437
438		\subsection{Transport Behavior}
401	–	Of importance in these types of studies are the transport properties
402	–	of the particles and their change in responce to altering
403	–	environmental conditions. In order to probe transport, constant energy
404	–	simulations were performed about the average density uncovered by the
405	–	constant pressure simulations. Simulations started with randomized
406	–	velocities and underwent 50 ps of temperature scaling and 50 ps of
407	–	constant energy equilibration before obtaining a 200 ps
408	–	trajectory. Diffusion constants were calculated via root-mean square
409	–	deviation analysis. The averaged results from five sets of NVE
410	–	simulations are displayed in figure \ref{diffuse}, alongside
411	–	experimental, SPC/E, and TIP5P
412	–	results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01}
439
440	+	Accurate dynamical properties of a water model are particularly
441	+	important when using the model to study permeation or transport across
442	+	biological membranes.  In order to probe transport in bulk water,
443	+	constant energy (NVE) simulations were performed at the average
444	+	density obtained by the NPT simulations at an identical target
445	+	temperature. Simulations started with randomized velocities and
446	+	underwent 50 ps of temperature scaling and 50 ps of constant energy
447	+	equilibration before a 200 ps data collection run. Diffusion constants
448	+	were calculated via linear fits to the long-time behavior of the
449	+	mean-square displacement as a function of time. The averaged results
450	+	from five sets of NVE simulations are displayed in figure
451	+	\ref{diffuse}, alongside experimental, SPC/E, and TIP5P
452	+	results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01}
453	+
454		\begin{figure}
455		\begin{center}
456		\epsfxsize=6in
457		\epsfbox{betterDiffuse.epsi}
458	<	\caption{Average diffusion coefficient over increasing temperature for
459	<	SSD, SPC/E,\cite{Clancy94} TIP5P,\cite{Jorgensen01} and Experimental
460	<	data.\cite{Gillen72,Mills73} Of the three water models shown, SSD has
461	<	the least deviation from the experimental values. The rapidly
462	<	increasing diffusion constants for TIP5P and SSD correspond to
463	<	significant decrease in density at the higher temperatures.}
458	>	\caption{Average self-diffusion constant as a function of temperature for
459	>	SSD, SPC/E [Ref. \citen{Clancy94}], TIP5P [Ref. \citen{Jorgensen01}],
460	>	and Experimental data [Refs. \citen{Gillen72} and \citen{Mills73}]. Of
461	>	the three water models shown, SSD has the least deviation from the
462	>	experimental values. The rapidly increasing diffusion constants for
463	>	TIP5P and SSD correspond to significant decrease in density at the
464	>	higher temperatures.}
465		\label{diffuse}
466		\end{center}
467		\end{figure}
468
469		The observed values for the diffusion constant point out one of the
470	<	strengths of the SSD model. Of the three experimental models shown,
471	<	the SSD model has the most accurate depiction of the diffusion trend
472	<	seen in experiment in both the supercooled and liquid temperature
473	<	regimes. SPC/E does a respectable job by producing values similar to
474	<	SSD and experiment around 290 K; however, it deviates at both higher
475	<	and lower temperatures, failing to predict the experimental
476	<	trend. TIP5P and SSD both start off low at colder temperatures and
477	<	tend to diffuse too rapidly at higher temperatures. This trend at
478	<	higher temperatures is not surprising in that the densities of both
479	<	TIP5P and SSD are lower than experimental water at these higher
480	<	temperatures. When calculating the diffusion coefficients for SSD at
481	<	experimental densities, the resulting values fall more in line with
482	<	experiment at these temperatures, albeit not at standard pressure.
470	>	strengths of the SSD model. Of the three models shown, the SSD model
471	>	has the most accurate depiction of self-diffusion in both the
472	>	supercooled and liquid regimes.  SPC/E does a respectable job by
473	>	reproducing values similar to experiment around 290 K; however, it
474	>	deviates at both higher and lower temperatures, failing to predict the
475	>	correct thermal trend. TIP5P and SSD both start off low at colder
476	>	temperatures and tend to diffuse too rapidly at higher temperatures.
477	>	This behavior at higher temperatures is not particularly surprising
478	>	since the densities of both TIP5P and SSD are lower than experimental
479	>	water densities at higher temperatures.  When calculating the
480	>	diffusion coefficients for SSD at experimental densities (instead of
481	>	the densities from the NPT simulations), the resulting values fall
482	>	more in line with experiment at these temperatures.
483
484		\subsection{Structural Changes and Characterization}
485	+
486		By starting the simulations from the crystalline state, the melting
487	<	transition and the ice structure can be studied along with the liquid
487	>	transition and the ice structure can be obtained along with the liquid
488		phase behavior beyond the melting point. The constant pressure heat
489		capacity (C$_\text{p}$) was monitored to locate the melting transition
490		in each of the simulations. In the melting simulations of the 1024
#	Line 450 | Line 492 | considerably lower than the experimental value, but th
492		at 245 K, indicating a first order phase transition for the melting of
493		these ice crystals. When the reaction field is turned off, the melting
494		transition occurs at 235 K.  These melting transitions are
495	<	considerably lower than the experimental value, but this is not a
454	<	surprise considering the simplicity of the SSD model.
495	>	considerably lower than the experimental value.
496
497		\begin{figure}
498		\begin{center}
499		\epsfxsize=6in
500		\epsfbox{corrDiag.eps}
501		\caption{Two dimensional illustration of angles involved in the
502	<	correlations observed in figure \ref{contour}.}
502	>	correlations observed in Fig. \ref{contour}.}
503		\label{corrAngle}
504		\end{center}
505		\end{figure}
#	Line 470 | Line 511 | depressions. White areas have g(\emph{r}) values below
511		\caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at
512		100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for
513		clarity: dark areas signify peaks while light areas signify
514	<	depressions. White areas have g(\emph{r}) values below 0.5 and black
514	>	depressions. White areas have $g(r)$ values below 0.5 and black
515		areas have values above 1.5.}
516		\label{contour}
517		\end{center}
518		\end{figure}
519
520	<	Additional analysis of the melting phase-transition process was
521	<	performed by using two-dimensional structure and dipole angle
522	<	correlations. Expressions for these correlations are as follows:
520	>	Additional analysis of the melting process was performed using
521	>	two-dimensional structure and dipole angle correlations. Expressions
522	>	for these correlations are as follows:
523
524		\begin{equation}
525	<	g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
525	>	g_{\text{AB}}(r,\cos\theta) = \frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\theta-\cos\theta_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
526		\end{equation}
527		\begin{equation}
528		g_{\text{AB}}(r,\cos\omega) =
529	<	\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|\mathbf{r}_{ij}\right|)\rangle\ ,
529	>	\frac{V}{N_\text{A}N_\text{B}}\langle\sum_{i\in\text{A}}\sum_{j\in\text{B}}\delta(\cos\omega-\cos\omega_{ij})\delta(r-\left|{\bf r}_{ij}\right|)\rangle\ ,
530		\end{equation}
531		where $\theta$ and $\omega$ refer to the angles shown in figure
532		\ref{corrAngle}. By binning over both distance and the cosine of the
533	<	desired angle between the two dipoles, the g(\emph{r}) can be
534	<	dissected to determine the common dipole arrangements that constitute
535	<	the peaks and troughs. Frames A and B of figure \ref{contour} show a
536	<	relatively crystalline state of an ice $I_c$ simulation. The first
537	<	peak of the g(\emph{r}) consists primarily of the preferred hydrogen
533	>	desired angle between the two dipoles, the $g(r)$ can be analyzed to
534	>	determine the common dipole arrangements that constitute the peaks and
535	>	troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
536	>	of figure \ref{contour} show results from an ice $I_c$ simulation. The
537	>	first peak in the $g(r)$ consists primarily of the preferred hydrogen
538		bonding arrangements as dictated by the tetrahedral sticky potential -
539	<	one peak for the donating and the other for the accepting hydrogen
540	<	bonds. Due to the high degree of crystallinity of the sample, the
541	<	second and third solvation shells show a repeated peak arrangement
539	>	one peak for the hydrogen bond donor and the other for the hydrogen
540	>	bond acceptor.  Due to the high degree of crystallinity of the sample,
541	>	the second and third solvation shells show a repeated peak arrangement
542		which decays at distances around the fourth solvation shell, near the
543		imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
544		In the higher temperature simulation shown in frames C and D, these
545	<	longer-ranged repeated peak features deteriorate rapidly. The first
546	<	solvation shell still shows the strong effect of the sticky-potential,
547	<	although it covers a larger area, extending to include a fraction of
548	<	aligned dipole peaks within the first solvation shell. The latter
549	<	peaks lose definition as thermal motion and the competing dipole force
550	<	overcomes the sticky potential's tight tetrahedral structuring of the
510	<	fluid.
545	>	long-range features deteriorate rapidly. The first solvation shell
546	>	still shows the strong effect of the sticky-potential, although it
547	>	covers a larger area, extending to include a fraction of aligned
548	>	dipole peaks within the first solvation shell. The latter peaks lose
549	>	due to thermal motion and as the competing dipole force overcomes the
550	>	sticky potential's tight tetrahedral structuring of the crystal.
551
552		This complex interplay between dipole and sticky interactions was
553		remarked upon as a possible reason for the split second peak in the
554	<	oxygen-oxygen g(\emph{r}).\cite{Ichiye96} At low temperatures, the
555	<	second solvation shell peak appears to have two distinct components
556	<	that blend together to form one observable peak. At higher
554	>	oxygen-oxygen $g_\mathrm{OO}(r)$.\cite{Ichiye96} At low temperatures,
555	>	the second solvation shell peak appears to have two distinct
556	>	components that blend together to form one observable peak. At higher
557		temperatures, this split character alters to show the leading 4 \AA\
558		peak dominated by equatorial anti-parallel dipole orientations. There
559		is also a tightly bunched group of axially arranged dipoles that most
#	Line 522 | Line 562 | dipolar repulsion term. Primary energetically favorabl
562		dipoles that assume hydrogen bond arrangements similar to those seen
563		in the first solvation shell. This evidence indicates that the dipole
564		pair interaction begins to dominate outside of the range of the
565	<	dipolar repulsion term. Primary energetically favorable dipole
565	>	dipolar repulsion term.  The energetically favorable dipole
566		arrangements populate the region immediately outside this repulsion
567	<	region (around 4 \AA), while arrangements that seek to ideally satisfy
568	<	both the sticky and dipole forces locate themselves just beyond this
567	>	region (around 4 \AA), while arrangements that seek to satisfy both
568	>	the sticky and dipole forces locate themselves just beyond this
569		initial buildup (around 5 \AA).
570
571		From these findings, the split second peak is primarily the product of
572		the dipolar repulsion term of the sticky potential. In fact, the inner
573		peak can be pushed out and merged with the outer split peak just by
574	<	extending the switching function cutoff ($s^\prime(r_{ij})$) from its
575	<	normal 4.0 \AA\ to values of 4.5 or even 5 \AA. This type of
574	>	extending the switching function ($s^\prime(r_{ij})$) from its normal
575	>	4.0 \AA\ cutoff to values of 4.5 or even 5 \AA. This type of
576		correction is not recommended for improving the liquid structure,
577		since the second solvation shell would still be shifted too far
578		out. In addition, this would have an even more detrimental effect on
579		the system densities, leading to a liquid with a more open structure
580	<	and a density considerably lower than the normal SSD behavior shown
581	<	previously. A better correction would be to include the
582	<	quadrupole-quadrupole interactions for the water particles outside of
583	<	the first solvation shell, but this reduces the simplicity and speed
584	<	advantage of SSD.
580	>	and a density considerably lower than the already low SSD density.  A
581	>	better correction would be to include the quadrupole-quadrupole
582	>	interactions for the water particles outside of the first solvation
583	>	shell, but this would remove the simplicity and speed advantage of
584	>	SSD.
585
586		\subsection{Adjusted Potentials: SSD/RF and SSD/E}
587	+
588		The propensity of SSD to adopt lower than expected densities under
589		varying conditions is troubling, especially at higher temperatures. In
590		order to correct this model for use with a reaction field, it is
#	Line 551 | Line 592 | densities while maintaining the accurate transport pro
592		intermolecular interactions. In undergoing a reparameterization, it is
593		important not to focus on just one property and neglect the other
594		important properties. In this case, it would be ideal to correct the
595	<	densities while maintaining the accurate transport properties.
595	>	densities while maintaining the accurate transport behavior.
596
597	<	The parameters available for tuning include the $\sigma$ and $\epsilon$
598	<	Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
599	<	attractive and dipole repulsive terms with their respective
600	<	cutoffs. To alter the attractive and repulsive terms of the sticky
601	<	potential independently, it is necessary to separate the terms as
602	<	follows:
603	<	\begin{equation}
604	<	u_{ij}^{sp}
605	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) =
565	<	\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)] + \frac{\nu_0^\prime}{2} [s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)],
566	<	\end{equation}
567	<
568	<	where $\nu_0$ scales the strength of the tetrahedral attraction and
569	<	$\nu_0^\prime$ acts in an identical fashion on the dipole repulsion
570	<	term. The separation was performed for purposes of the
571	<	reparameterization, but the final parameters were adjusted so that it
572	<	is unnecessary to separate the terms when implementing the adjusted
573	<	water potentials. The results of the reparameterizations are shown in
574	<	table \ref{params}. Note that the tetrahedral attractive and dipolar
575	<	repulsive terms do not share the same lower cutoff ($r_l$) in the
576	<	newly parameterized potentials - soft sticky dipole reaction field
577	<	(SSD/RF - for use with a reaction field) and soft sticky dipole
578	<	enhanced (SSD/E - an attempt to improve the liquid structure in
597	>	The parameters available for tuning include the $\sigma$ and
598	>	$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the
599	>	strength of the sticky potential ($\nu_0$), and the sticky attractive
600	>	and dipole repulsive cubic switching function cutoffs ($r_l$, $r_u$
601	>	and $r_l^\prime$, $r_u^\prime$ respectively). The results of the
602	>	reparameterizations are shown in table \ref{params}. We are calling
603	>	these reparameterizations the Soft Sticky Dipole / Reaction Field
604	>	(SSD/RF - for use with a reaction field) and Soft Sticky Dipole
605	>	Extended (SSD/E - an attempt to improve the liquid structure in
606		simulations without a long-range correction).
607
608		\begin{table}
#	Line 583 | Line 610 | simulations without a long-range correction).
610		\caption{Parameters for the original and adjusted models}
611		\begin{tabular}{ l  c  c  c  c }
612		\hline \\[-3mm]
613	<	\ \ \ Parameters\ \ \  & \ \ \ SSD\cite{Ichiye96} \ \ \ & \ SSD1\cite{Ichiye03}\ \  & \ SSD/E\ \  & \ SSD/RF \\
613	>	\ \ \ Parameters\ \ \  & \ \ \ SSD [Ref. \citen{Ichiye96}] \ \ \
614	>	& \ SSD1 [Ref. \citen{Ichiye03}]\ \  & \ SSD/E\ \  & \ SSD/RF \\
615		\hline \\[-3mm]
616		\ \ \ $\sigma$ (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
617		\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
618		\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\
619		\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
620	+	\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\
621		\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
622		\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
594	–	\ \ \ $\nu_0^\prime$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
623		\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\
624		\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
625		\end{tabular}
#	Line 603 | Line 631 | simulations without a long-range correction).
631		\begin{center}
632		\epsfxsize=5in
633		\epsfbox{GofRCompare.epsi}
634	<	\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E
634	>	\caption{Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with SSD/E
635		and SSD1 without reaction field (top), as well as SSD/RF and SSD1 with
636		reaction field turned on (bottom). The insets show the respective
637		first peaks in detail. Note how the changes in parameters have lowered
#	Line 624 | Line 652 | In the paper detailing the development of SSD, Liu and
652		\end{center}
653		\end{figure}
654
655	<	In the paper detailing the development of SSD, Liu and Ichiye placed
656	<	particular emphasis on an accurate description of the first solvation
657	<	shell. This resulted in a somewhat tall and narrow first peak in the
658	<	g(\emph{r}) that integrated to give similar coordination numbers to
655	>	In the original paper detailing the development of SSD, Liu and Ichiye
656	>	placed particular emphasis on an accurate description of the first
657	>	solvation shell. This resulted in a somewhat tall and narrow first
658	>	peak in $g(r)$ that integrated to give similar coordination numbers to
659		the experimental data obtained by Soper and
660		Phillips.\cite{Ichiye96,Soper86} New experimental x-ray scattering
661		data from the Head-Gordon lab indicates a slightly lower and shifted
662	<	first peak in the g$_\mathrm{OO}(r)$, so adjustments to SSD were made
663	<	while taking into consideration the new experimental
662	>	first peak in the g$_\mathrm{OO}(r)$, so our adjustments to SSD were
663	>	made while taking into consideration the new experimental
664		findings.\cite{Head-Gordon00_1} Figure \ref{grcompare} shows the
665	<	relocation of the first peak of the oxygen-oxygen g(\emph{r}) by
666	<	comparing the revised SSD model (SSD1), SSD-E, and SSD-RF to the new
665	>	relocation of the first peak of the oxygen-oxygen $g(r)$ by comparing
666	>	the revised SSD model (SSD1), SSD/E, and SSD/RF to the new
667		experimental results. Both modified water models have shorter peaks
668	<	that are brought in more closely to the experimental peak (as seen in
669	<	the insets of figure \ref{grcompare}).  This structural alteration was
668	>	that match more closely to the experimental peak (as seen in the
669	>	insets of figure \ref{grcompare}).  This structural alteration was
670		accomplished by the combined reduction in the Lennard-Jones $\sigma$
671	<	variable and adjustment of the sticky potential strength and
672	<	cutoffs. As can be seen in table \ref{params}, the cutoffs for the
673	<	tetrahedral attractive and dipolar repulsive terms were nearly swapped
674	<	with each other. Isosurfaces of the original and modified sticky
675	<	potentials are shown in figure \ref{isosurface}. In these isosurfaces,
676	<	it is easy to see how altering the cutoffs changes the repulsive and
677	<	attractive character of the particles. With a reduced repulsive
678	<	surface (darker region), the particles can move closer to one another,
679	<	increasing the density for the overall system. This change in
680	<	interaction cutoff also results in a more gradual orientational motion
681	<	by allowing the particles to maintain preferred dipolar arrangements
682	<	before they begin to feel the pull of the tetrahedral
683	<	restructuring. As the particles move closer together, the dipolar
684	<	repulsion term becomes active and excludes unphysical nearest-neighbor
685	<	arrangements. This compares with how SSD and SSD1 exclude preferred
686	<	dipole alignments before the particles feel the pull of the ``hydrogen
687	<	bonds''. Aside from improving the shape of the first peak in the
688	<	g(\emph{r}), this modification improves the densities considerably by
689	<	allowing the persistence of full dipolar character below the previous
690	<	4.0 \AA\ cutoff.
671	>	variable and adjustment of the sticky potential strength and cutoffs.
672	>	As can be seen in table \ref{params}, the cutoffs for the tetrahedral
673	>	attractive and dipolar repulsive terms were nearly swapped with each
674	>	other.  Isosurfaces of the original and modified sticky potentials are
675	>	shown in figure \ref{isosurface}. In these isosurfaces, it is easy to
676	>	see how altering the cutoffs changes the repulsive and attractive
677	>	character of the particles. With a reduced repulsive surface (darker
678	>	region), the particles can move closer to one another, increasing the
679	>	density for the overall system.  This change in interaction cutoff also
680	>	results in a more gradual orientational motion by allowing the
681	>	particles to maintain preferred dipolar arrangements before they begin
682	>	to feel the pull of the tetrahedral restructuring. As the particles
683	>	move closer together, the dipolar repulsion term becomes active and
684	>	excludes unphysical nearest-neighbor arrangements. This compares with
685	>	how SSD and SSD1 exclude preferred dipole alignments before the
686	>	particles feel the pull of the ``hydrogen bonds''. Aside from
687	>	improving the shape of the first peak in the g(\emph{r}), this
688	>	modification improves the densities considerably by allowing the
689	>	persistence of full dipolar character below the previous 4.0 \AA\
690	>	cutoff.
691
692	<	While adjusting the location and shape of the first peak of
693	<	g(\emph{r}) improves the densities, these changes alone are
694	<	insufficient to bring the system densities up to the values observed
695	<	experimentally. To further increase the densities, the dipole moments
696	<	were increased in both of the adjusted models. Since SSD is a dipole
697	<	based model, the structure and transport are very sensitive to changes
698	<	in the dipole moment. The original SSD simply used the dipole moment
699	<	calculated from the TIP3P water model, which at 2.35 D is
700	<	significantly greater than the experimental gas phase value of 1.84
701	<	D. The larger dipole moment is a more realistic value and improves the
702	<	dielectric properties of the fluid. Both theoretical and experimental
703	<	measurements indicate a liquid phase dipole moment ranging from 2.4 D
704	<	to values as high as 3.11 D, providing a substantial range of
705	<	reasonable values for a dipole
678	<	moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
692	>	While adjusting the location and shape of the first peak of $g(r)$
693	>	improves the densities, these changes alone are insufficient to bring
694	>	the system densities up to the values observed experimentally.  To
695	>	further increase the densities, the dipole moments were increased in
696	>	both of our adjusted models. Since SSD is a dipole based model, the
697	>	structure and transport are very sensitive to changes in the dipole
698	>	moment. The original SSD simply used the dipole moment calculated from
699	>	the TIP3P water model, which at 2.35 D is significantly greater than
700	>	the experimental gas phase value of 1.84 D. The larger dipole moment
701	>	is a more realistic value and improves the dielectric properties of
702	>	the fluid. Both theoretical and experimental measurements indicate a
703	>	liquid phase dipole moment ranging from 2.4 D to values as high as
704	>	3.11 D, providing a substantial range of reasonable values for a
705	>	dipole moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
706		increasing the dipole moments to 2.42 and 2.48 D for SSD/E and SSD/RF,
707		respectively, leads to significant changes in the density and
708		transport of the water models.
#	Line 693 | Line 720 | collection times as stated earlier in this paper.
720		run at each temperature step, and the final configuration from the
721		previous temperature simulation was used as a starting point. All NVE
722		simulations had the same thermalization, equilibration, and data
723	<	collection times as stated earlier in this paper.
723	>	collection times as stated previously.
724
725		\begin{figure}
726		\begin{center}
727		\epsfxsize=6in
728		\epsfbox{ssdeDense.epsi}
729		\caption{Comparison of densities calculated with SSD/E to SSD1 without a
730	<	reaction field, TIP3P,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00}
731	<	SPC/E,\cite{Clancy94} and experiment.\cite{CRC80} The window shows a
732	<	expansion around 300 K with error bars included to clarify this region
733	<	of interest. Note that both SSD1 and SSD/E show good agreement with
730	>	reaction field, TIP3P [Ref. \citen{Jorgensen98b}], TIP5P
731	>	[Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}] and
732	>	experiment [Ref. \citen{CRC80}]. The window shows a expansion around
733	>	300 K with error bars included to clarify this region of
734	>	interest. Note that both SSD1 and SSD/E show good agreement with
735		experiment when the long-range correction is neglected.}
736		\label{ssdedense}
737		\end{center}
738		\end{figure}
739
740	<	Figure \ref{ssdedense} shows the density profile for the SSD/E model
740	>	Fig. \ref{ssdedense} shows the density profile for the SSD/E model
741		in comparison to SSD1 without a reaction field, other common water
742		models, and experimental results. The calculated densities for both
743		SSD/E and SSD1 have increased significantly over the original SSD
744	<	model (see figure \ref{dense1}) and are in better agreement with the
744	>	model (see fig. \ref{dense1}) and are in better agreement with the
745		experimental values. At 298 K, the densities of SSD/E and SSD1 without
746		a long-range correction are 0.996$\pm$0.001 g/cm$^3$ and
747		0.999$\pm$0.001 g/cm$^3$ respectively.  These both compare well with
#	Line 725 | Line 753 | little effect on the melting transition. By monitoring
753		comes about via an increase in the liquid disorder through the
754		weakening of the sticky potential and strengthening of the dipolar
755		character. However, this increasing disorder in the SSD/E model has
756	<	little effect on the melting transition. By monitoring C$\text{p}$
756	>	little effect on the melting transition. By monitoring $C_p$
757		throughout these simulations, the melting transition for SSD/E was
758	<	shown to occur at 235 K, the same transition temperature observed with
759	<	SSD and SSD1.
758	>	shown to occur at 235 K.  The same transition temperature observed
759	>	with SSD and SSD1.
760
761		\begin{figure}
762		\begin{center}
763		\epsfxsize=6in
764		\epsfbox{ssdrfDense.epsi}
765		\caption{Comparison of densities calculated with SSD/RF to SSD1 with a
766	<	reaction field, TIP3P,\cite{Jorgensen98b} TIP5P,\cite{Jorgensen00}
767	<	SPC/E,\cite{Clancy94} and experiment.\cite{CRC80} The inset shows the
768	<	necessity of reparameterization when utilizing a reaction field
769	<	long-ranged correction - SSD/RF provides significantly more accurate
770	<	densities than SSD1 when performing room temperature simulations.}
766	>	reaction field, TIP3P [Ref. \citen{Jorgensen98b}], TIP5P
767	>	[Ref. \citen{Jorgensen00}], SPC/E [Ref. \citen{Clancy94}], and
768	>	experiment [Ref. \citen{CRC80}]. The inset shows the necessity of
769	>	reparameterization when utilizing a reaction field long-ranged
770	>	correction - SSD/RF provides significantly more accurate densities
771	>	than SSD1 when performing room temperature simulations.}
772		\label{ssdrfdense}
773		\end{center}
774		\end{figure}
775
776		Including the reaction field long-range correction in the simulations
777	<	results in a more interesting comparison. A density profile including
777	>	results in a more interesting comparison.  A density profile including
778		SSD/RF and SSD1 with an active reaction field is shown in figure
779		\ref{ssdrfdense}.  As observed in the simulations without a reaction
780		field, the densities of SSD/RF and SSD1 show a dramatic increase over
781		normal SSD (see figure \ref{dense1}). At 298 K, SSD/RF has a density
782		of 0.997$\pm$0.001 g/cm$^3$, directly in line with experiment and
783	<	considerably better than the SSD value of 0.941$\pm$0.001 g/cm$^3$ and
784	<	the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results further
785	<	emphasize the importance of reparameterization in order to model the
786	<	density properly under different simulation conditions. Again, these
787	<	changes have only a minor effect on the melting point, which observed
788	<	at 245 K for SSD/RF, is identical to SSD and only 5 K lower than SSD1
789	<	with a reaction field. Additionally, the difference in density maxima
790	<	is not as extreme, with SSD/RF showing a density maximum at 255 K,
791	<	fairly close to the density maxima of 260 K and 265 K, shown by SSD
792	<	and SSD1 respectively.
783	>	considerably better than the original SSD value of 0.941$\pm$0.001
784	>	g/cm$^3$ and the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results
785	>	further emphasize the importance of reparameterization in order to
786	>	model the density properly under different simulation conditions.
787	>	Again, these changes have only a minor effect on the melting point,
788	>	which observed at 245 K for SSD/RF, is identical to SSD and only 5 K
789	>	lower than SSD1 with a reaction field. Additionally, the difference in
790	>	density maxima is not as extreme, with SSD/RF showing a density
791	>	maximum at 255 K, fairly close to the density maxima of 260 K and 265
792	>	K, shown by SSD and SSD1 respectively.
793
794		\begin{figure}
795		\begin{center}
796		\epsfxsize=6in
797		\epsfbox{ssdeDiffuse.epsi}
798	<	\caption{Plots of the diffusion constants calculated from SSD/E and SSD1,
799	<	both without a reaction field, along with experimental
800	<	results.\cite{Gillen72,Mills73} The NVE calculations were performed
801	<	at the average densities observed in the 1 atm NPT simulations for
802	<	the respective models. SSD/E is slightly more fluid than experiment
803	<	at all of the temperatures, but it is closer than SSD1 without a
798	>	\caption{The diffusion constants calculated from SSD/E and SSD1,
799	>	both without a reaction field, along with experimental results
800	>	[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations
801	>	were performed at the average densities observed in the 1 atm NPT
802	>	simulations for the respective models. SSD/E is slightly more mobile
803	>	than experiment at all of the temperatures, but it is closer to
804	>	experiment at biologically relavent temperatures than SSD1 without a
805		long-range correction.}
806		\label{ssdediffuse}
807		\end{center}
#	Line 780 | Line 810 | the densities, it is important that particle transport
810		The reparameterization of the SSD water model, both for use with and
811		without an applied long-range correction, brought the densities up to
812		what is expected for simulating liquid water. In addition to improving
813	<	the densities, it is important that particle transport be maintained
814	<	or improved. Figure \ref{ssdediffuse} compares the temperature
815	<	dependence of the diffusion constant of SSD/E to SSD1 without an
816	<	active reaction field, both at the densities calculated at 1 atm and
817	<	at the experimentally calculated densities for super-cooled and liquid
818	<	water. The diffusion constant for SSD/E is consistently a little
819	<	higher than experiment, while SSD1 remains lower than experiment until
820	<	relatively high temperatures (greater than 330 K). Both models follow
821	<	the shape of the experimental curve well below 300 K but tend to
822	<	diffuse too rapidly at higher temperatures, something that is
823	<	especially apparent with SSD1. This accelerated increasing of
824	<	diffusion is caused by the rapidly decreasing system density with
825	<	increasing temperature. Though it is difficult to see in figure
826	<	\ref{ssdedense}, the densities of SSD1 decay more rapidly with
827	<	temperature than do those of SSD/E, leading to more visible deviation
828	<	from the experimental diffusion trend. Thus, the changes made to
813	>	the densities, it is important that the excellent diffusive behavior
814	>	of SSD be maintained or improved. Figure \ref{ssdediffuse} compares
815	>	the temperature dependence of the diffusion constant of SSD/E to SSD1
816	>	without an active reaction field at the densities calculated from the
817	>	NPT simulations at 1 atm. The diffusion constant for SSD/E is
818	>	consistently higher than experiment, while SSD1 remains lower than
819	>	experiment until relatively high temperatures (around 360 K). Both
820	>	models follow the shape of the experimental curve well below 300 K but
821	>	tend to diffuse too rapidly at higher temperatures, as seen in SSD1's
822	>	crossing above 360 K.  This increasing diffusion relative to the
823	>	experimental values is caused by the rapidly decreasing system density
824	>	with increasing temperature.  Both SSD1 and SSD/E show this deviation
825	>	in diffusive behavior, but this trend has different implications on
826	>	the diffusive behavior of the models.  While SSD1 shows more
827	>	experimentally accurate diffusive behavior in the high temperature
828	>	regimes, SSD/E shows more accurate behavior in the supercooled and
829	>	biologically relavent temperature ranges.  Thus, the changes made to
830		improve the liquid structure may have had an adverse affect on the
831		density maximum, but they improve the transport behavior of SSD/E
832	<	relative to SSD1.
832	>	relative to SSD1 under the most commonly simulated conditions.
833
834		\begin{figure}
835		\begin{center}
836		\epsfxsize=6in
837		\epsfbox{ssdrfDiffuse.epsi}
838	<	\caption{Plots of the diffusion constants calculated from SSD/RF and SSD1,
839	<	both with an active reaction field, along with experimental
840	<	results.\cite{Gillen72,Mills73} The NVE calculations were performed
841	<	at the average densities observed in the 1 atm NPT simulations for
842	<	both of the models. Note how accurately SSD/RF simulates the
843	<	diffusion of water throughout this temperature range. The more
844	<	rapidly increasing diffusion constants at high temperatures for both
845	<	models is attributed to the significantly lower densities than
846	<	observed in experiment.}
838	>	\caption{The diffusion constants calculated from SSD/RF and SSD1,
839	>	both with an active reaction field, along with experimental results
840	>	[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations
841	>	were performed at the average densities observed in the 1 atm NPT
842	>	simulations for both of the models. Note how accurately SSD/RF
843	>	simulates the diffusion of water throughout this temperature
844	>	range. The more rapidly increasing diffusion constants at high
845	>	temperatures for both models is attributed to lower calculated
846	>	densities than those observed in experiment.}
847		\label{ssdrfdiffuse}
848		\end{center}
849		\end{figure}
850
851		In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
852		compared to SSD1 with an active reaction field. Note that SSD/RF
853	<	tracks the experimental results incredibly well, identical within
854	<	error throughout the temperature range shown and with only a slight
855	<	increasing trend at higher temperatures. SSD1 tends to diffuse more
856	<	slowly at low temperatures and deviates to diffuse too rapidly at
857	<	temperatures greater than 330 K. As stated in relation to SSD/E, this
858	<	deviation away from the ideal trend is due to a rapid decrease in
859	<	density at higher temperatures. SSD/RF does not suffer from this
860	<	problem as much as SSD1, because the calculated densities are closer
861	<	to the experimental value. These results again emphasize the
862	<	importance of careful reparameterization when using an altered
832	<	long-range correction.
853	>	tracks the experimental results quantitatively, identical within error
854	>	throughout most of the temperature range shown and exhibiting only a
855	>	slight increasing trend at higher temperatures. SSD1 tends to diffuse
856	>	more slowly at low temperatures and deviates to diffuse too rapidly at
857	>	temperatures greater than 330 K.  As stated above, this deviation away
858	>	from the ideal trend is due to a rapid decrease in density at higher
859	>	temperatures. SSD/RF does not suffer from this problem as much as SSD1
860	>	because the calculated densities are closer to the experimental
861	>	values. These results again emphasize the importance of careful
862	>	reparameterization when using an altered long-range correction.
863
864	+	\begin{table}
865	+	\begin{center}
866	+	\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}.}
867	+	\begin{tabular}{ l  c  c  c  c  c }
868	+	\hline \\[-3mm]
869	+	\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ SSD/E \ \ \ & \ SSD1 (RF) \ \
870	+	\ & \ SSD/RF \ \ \ & \ Expt. \\
871	+	\hline \\[-3mm]
872	+	\ \ \ $\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 \\
873	+	\ \ \ $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 \\
874	+	\ \ \ $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}$ \\
875	+	\ \ \ Coordination Number & 3.9 & 4.3 & 3.8 & 4.4 & 4.7$^\text{b}$ \\
876	+	\ \ \ H-bonds per particle & 3.7 & 3.6 & 3.7 & 3.7 & 3.4$^\text{c}$ \\
877	+	\ \ \ $\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}$ \\
878	+	\ \ \ $\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}$ \\
879	+	\end{tabular}
880	+	\label{liquidproperties}
881	+	\end{center}
882	+	\end{table}
883	+
884	+	Table \ref{liquidproperties} gives a synopsis of the liquid state
885	+	properties of the water models compared in this study along with the
886	+	experimental values for liquid water at ambient conditions. The
887	+	coordination number and hydrogen bonds per particle were calculated by
888	+	integrating the following relation:
889	+	\begin{equation}
890	+	4\pi\rho\int_{0}^{a}r^2\text{g}(r)dr,
891	+	\end{equation}
892	+	where $\rho$ is the number density of pair interactions, $a$ is the
893	+	radial location of the minima following the first solvation shell
894	+	peak, and g$(r)$ is either g$_\text{OO}(r)$ or g$_\text{OH}(r)$ for
895	+	calculation of the coordination number or hydrogen bonds per particle
896	+	respectively. The number of hydrogen bonds stays relatively constant
897	+	across all of the models, but the coordination numbers of SSD/E and
898	+	SSD/RF show an improvement over SSD1. This improvement is primarily
899	+	due to the widening of the first solvation shell peak, allowing the
900	+	first minima to push outward. Comparing the coordination number with
901	+	the number of hydrogen bonds can lead to more insight into the
902	+	structural character of the liquid.  Because of the near identical
903	+	values for SSD1, it appears to be a little too exclusive, in that all
904	+	molecules in the first solvation shell are involved in forming ideal
905	+	hydrogen bonds.  The differing numbers for the newly parameterized
906	+	models indicate the allowance of more fluid configurations in addition
907	+	to the formation of an acceptable number of ideal hydrogen bonds.
908	+
909	+	The time constants for the self orientational autocorrelation function
910	+	are also displayed in Table \ref{liquidproperties}. The dipolar
911	+	orientational time correlation function ($\Gamma_{l}$) is described
912	+	by:
913	+	\begin{equation}
914	+	\Gamma_{l}(t) = \langle P_l[\mathbf{u}_j(0)\cdot\mathbf{u}_j(t)]\rangle,
915	+	\end{equation}
916	+	where $P_l$ is a Legendre polynomial of order $l$ and $\mathbf{u}_j$
917	+	is the unit vector of the particle dipole.\cite{Rahman71} From these
918	+	correlation functions, the orientational relaxation time of the dipole
919	+	vector can be calculated from an exponential fit in the long-time
920	+	regime ($t > \tau_l^\mu$).\cite{Rothschild84} Calculation of these
921	+	time constants were averaged from five detailed NVE simulations
922	+	performed at the STP density for each of the respective models. Again,
923	+	SSD/E and SSD/RF show improved behavior over SSD1 both with and
924	+	without an active reaction field. Numbers published from the original
925	+	SSD dynamics studies appear closer to the experimental values, and we
926	+	attribute this discrepancy to the implimentation of an Ewald sum
927	+	versus a reaction field.
928	+
929		\subsection{Additional Observations}
930
931		\begin{figure}
#	Line 838 | Line 933 | long-range correction.
933		\epsfxsize=6in
934		\epsfbox{povIce.ps}
935		\caption{A water lattice built from the crystal structure assumed by
936	<	SSD/E when undergoing an extremely restricted temperature NPT
937	<	simulation. This form of ice is referred to as ice \emph{i} to
938	<	emphasize its simulation origins. This image was taken of the (001)
939	<	face of the crystal.}
936	>	SSD/E when undergoing an extremely restricted temperature NPT
937	>	simulation. This form of ice is referred to as ice-{\it i} to
938	>	emphasize its simulation origins. This image was taken of the (001)
939	>	face of the crystal.}
940		\label{weirdice}
941		\end{center}
942		\end{figure}
943
944	<	While performing restricted temperature melting sequences of SSD/E not
945	<	previously discussed, some interesting observations were made. After
946	<	melting at 235 K, two of five systems underwent crystallization events
947	<	near 245 K. As the heating process continued, the two systems remained
948	<	crystalline until finally melting between 320 and 330 K. The final
949	<	configurations of these two melting sequences show an expanded
950	<	zeolite-like crystal structure that does not correspond to any known
951	<	form of ice. For convenience, and to help distinguish it from the
952	<	experimentally observed forms of ice, this crystal structure will
953	<	henceforth be referred to as ice $\sqrt{\smash[b]{-\text{I}}}$ (ice
954	<	\emph{i}). The crystallinity was extensive enough that a near ideal
955	<	crystal structure of ice \emph{i} could be obtained. Figure
956	<	\ref{weirdice} shows the repeating crystal structure of a typical
957	<	crystal at 5 K. Each water molecule is hydrogen bonded to four others;
958	<	however, the hydrogen bonds are flexed rather than perfectly
959	<	straight. This results in a skewed tetrahedral geometry about the
960	<	central molecule. Referring to figure \ref{isosurface}, these flexed
961	<	hydrogen bonds are allowed due to the conical shape of the attractive
962	<	regions, with the greatest attraction along the direct hydrogen bond
944	>	While performing a series of melting simulations on an early iteration
945	>	of SSD/E not discussed in this paper, we observed recrystallization
946	>	into a novel structure not previously known for water.  After melting
947	>	at 235 K, two of five systems underwent crystallization events near
948	>	245 K.  The two systems remained crystalline up to 320 and 330 K,
949	>	respectively.  The crystal exhibits an expanded zeolite-like structure
950	>	that does not correspond to any known form of ice.  This appears to be
951	>	an artifact of the point dipolar models, so to distinguish it from the
952	>	experimentally observed forms of ice, we have denoted the structure
953	>	Ice-$\sqrt{\smash[b]{-\text{I}}}$ (ice-{\it i}).  A large enough
954	>	portion of the sample crystallized that we have been able to obtain a
955	>	near ideal crystal structure of ice-{\it i}. Figure \ref{weirdice}
956	>	shows the repeating crystal structure of a typical crystal at 5
957	>	K. Each water molecule is hydrogen bonded to four others; however, the
958	>	hydrogen bonds are bent rather than perfectly straight. This results
959	>	in a skewed tetrahedral geometry about the central molecule.  In
960	>	figure \ref{isosurface}, it is apparent that these flexed hydrogen
961	>	bonds are allowed due to the conical shape of the attractive regions,
962	>	with the greatest attraction along the direct hydrogen bond
963		configuration. Though not ideal, these flexed hydrogen bonds are
964	<	favorable enough to stabilize an entire crystal generated around
870	<	them. In fact, the imperfect ice \emph{i} crystals were so stable that
871	<	they melted at temperatures nearly 100 K greater than both ice I$_c$
872	<	and I$_h$.
964	>	favorable enough to stabilize an entire crystal generated around them.
965
966	<	These initial simulations indicated that ice \emph{i} is the preferred
967	<	ice structure for at least the SSD/E model. To verify this, a
968	<	comparison was made between near ideal crystals of ice $I_h$, ice
969	<	$I_c$, and ice 0 at constant pressure with SSD/E, SSD/RF, and
970	<	SSD1. Near ideal versions of the three types of crystals were cooled
971	<	to 1 K, and the potential energies of each were compared using all
972	<	three water models. With every water model, ice \emph{i} turned out to
973	<	have the lowest potential energy: 5\% lower than $I_h$ with SSD1,
974	<	6.5\% lower with SSD/E, and 7.5\% lower with SSD/RF.
966	>	Initial simulations indicated that ice-{\it i} is the preferred ice
967	>	structure for at least the SSD/E model. To verify this, a comparison
968	>	was made between near ideal crystals of ice~$I_h$, ice~$I_c$, and
969	>	ice-{\it i} at constant pressure with SSD/E, SSD/RF, and
970	>	SSD1. Near-ideal versions of the three types of crystals were cooled
971	>	to 1 K, and the enthalpies of each were compared using all three water
972	>	models. With every model in the SSD family, ice-{\it i} had the lowest
973	>	calculated enthalpy: 5\% lower than $I_h$ with SSD1, 6.5\% lower with
974	>	SSD/E, and 7.5\% lower with SSD/RF.  The enthalpy data is summarized
975	>	in Table \ref{iceenthalpy}.
976
977	<	In addition to these low temperature comparisons, melting sequences
978	<	were performed with ice \emph{i} as the initial configuration using
979	<	SSD/E, SSD/RF, and SSD1 both with and without a reaction field. The
980	<	melting transitions for both SSD/E and SSD1 without a reaction field
981	<	occurred at temperature in excess of 375 K. SSD/RF and SSD1 with a
982	<	reaction field showed more reasonable melting transitions near 325
983	<	K. These melting point observations emphasize the preference for this
984	<	crystal structure over the most common types of ice when using these
985	<	single point water models.
977	>	\begin{table}
978	>	\begin{center}
979	>	\caption{Enthalpies (in kcal / mol) of the three crystal structures (at 1
980	>	K) exhibited by the SSD family of water models}
981	>	\begin{tabular}{ l  c  c  c  }
982	>	\hline \\[-3mm]
983	>	\ \ \ Water Model \ \ \  & \ \ \ Ice-$I_h$ \ \ \ & \ Ice-$I_c$\ \  & \
984	>	Ice-{\it i} \\
985	>	\hline \\[-3mm]
986	>	\ \ \ SSD/E & -12.286 & -12.292 & -13.590 \\
987	>	\ \ \ SSD/RF & -12.935 & -12.917 & -14.022 \\
988	>	\ \ \ SSD1 & -12.496 & -12.411 & -13.417 \\
989	>	\ \ \ SSD1 (RF) & -12.504 & -12.411 & -13.134 \\
990	>	\end{tabular}
991	>	\label{iceenthalpy}
992	>	\end{center}
993	>	\end{table}
994
995	<	Recognizing that the above tests show ice \emph{i} to be both the most
996	<	stable and lowest density crystal structure for these single point
997	<	water models, it is interesting to speculate on the relative stability
998	<	of this crystal structure with charge based water models. As a quick
999	<	test, these 3 crystal types were converted from SSD type particles to
1000	<	TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
1001	<	minimizations were performed on the crystals to compare the system
1002	<	energies. Again, ice \emph{i} was observed to have the lowest total
902	<	system energy. The total energy of ice \emph{i} was ~2\% lower than
903	<	ice $I_h$, which was in turn ~3\% lower than ice $I_c$. Based on these
904	<	initial studies, it would not be surprising if results from the other
905	<	common water models show ice \emph{i} to be the lowest energy crystal
906	<	structure. A continuation of this work studying ice \emph{i} with
907	<	multi-point water models will be published in a coming article.
995	>	In addition to these energetic comparisons, melting simulations were
996	>	performed with ice-{\it i} as the initial configuration using SSD/E,
997	>	SSD/RF, and SSD1 both with and without a reaction field. The melting
998	>	transitions for both SSD/E and SSD1 without reaction field occurred at
999	>	temperature in excess of 375~K.  SSD/RF and SSD1 with a reaction field
1000	>	showed more reasonable melting transitions near 325~K.  These melting
1001	>	point observations clearly show that all of the SSD-derived models
1002	>	prefer the ice-{\it i} structure.
1003
1004		\section{Conclusions}
910	–	The density maximum and temperature dependent transport for the SSD
911	–	water model, both with and without the use of reaction field, were
912	–	studied via a series of NPT and NVE simulations. The constant pressure
913	–	simulations of the melting of both $I_h$ and $I_c$ ice showed a
914	–	density maximum near 260 K. In most cases, the calculated densities
915	–	were significantly lower than the densities calculated in simulations
916	–	of other water models. Analysis of particle diffusion showed SSD to
917	–	capture the transport properties of experimental water well in both
918	–	the liquid and super-cooled liquid regimes. In order to correct the
919	–	density behavior, the original SSD model was reparameterized for use
920	–	both with and without a reaction field (SSD/RF and SSD/E), and
921	–	comparison simulations were performed with SSD1, the density corrected
922	–	version of SSD. Both models improve the liquid structure, density
923	–	values, and diffusive properties under their respective conditions,
924	–	indicating the necessity of reparameterization when altering the
925	–	long-range correction specifics. When taking into account the
926	–	appropriate considerations, these simple water models are excellent
927	–	choices for representing explicit water in large scale simulations of
928	–	biochemical systems.
1005
1006	+	The density maximum and temperature dependence of the self-diffusion
1007	+	constant were studied for the SSD water model, both with and without
1008	+	the use of reaction field, via a series of NPT and NVE
1009	+	simulations. The constant pressure simulations showed a density
1010	+	maximum near 260 K. In most cases, the calculated densities were
1011	+	significantly lower than the densities obtained from other water
1012	+	models (and experiment). Analysis of self-diffusion showed SSD to
1013	+	capture the transport properties of water well in both the liquid and
1014	+	super-cooled liquid regimes.
1015	+
1016	+	In order to correct the density behavior, the original SSD model was
1017	+	reparameterized for use both with and without a reaction field (SSD/RF
1018	+	and SSD/E), and comparisons were made with SSD1, Ichiye's density
1019	+	corrected version of SSD. Both models improve the liquid structure,
1020	+	densities, and diffusive properties under their respective simulation
1021	+	conditions, indicating the necessity of reparameterization when
1022	+	changing the method of calculating long-range electrostatic
1023	+	interactions.  In general, however, these simple water models are
1024	+	excellent choices for representing explicit water in large scale
1025	+	simulations of biochemical systems.
1026	+
1027	+	The existence of a novel low-density ice structure that is preferred
1028	+	by the SSD family of water models is somewhat troubling, since liquid
1029	+	simulations on this family of water models at room temperature are
1030	+	effectively simulations of super-cooled or metastable liquids.  One
1031	+	way to de-stabilize this unphysical ice structure would be to make the
1032	+	range of angles preferred by the attractive part of the sticky
1033	+	potential much narrower.  This would require extensive
1034	+	reparameterization to maintain the same level of agreement with the
1035	+	experiments.
1036	+
1037	+	Additionally, our initial calculations show that the ice-{\it i}
1038	+	structure may also be a preferred crystal structure for at least one
1039	+	other popular multi-point water model (TIP3P), and that much of the
1040	+	simulation work being done using this popular model could also be at
1041	+	risk for crystallization into this unphysical structure.  A future
1042	+	publication will detail the relative stability of the known ice
1043	+	structures for a wide range of popular water models.
1044	+
1045		\section{Acknowledgments}
1046		Support for this project was provided by the National Science
1047		Foundation under grant CHE-0134881. Computation time was provided by
1048		the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
1049	<	DMR 00 79647.
1049	>	DMR-0079647.
1050
936	–
1051		\newpage
1052
1053		\bibliographystyle{jcp}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines