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 743 by chrisfen, Wed Sep 3 21:36:09 2003 UTC vs.
Revision 1029 by gezelter, Thu Feb 5 20:47:50 2004 UTC

#	Line 1 | Line 1
1	<	\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
2	<	%\documentclass[prb,aps,times,tabularx,preprint]{revtex4}
1	>	%\documentclass[prb,aps,times,twocolumn,tabularx]{revtex4}
2	>	\documentclass[11pt]{article}
3	>	\usepackage{endfloat}
4		\usepackage{amsmath}
5	+	\usepackage{epsf}
6	+	\usepackage{berkeley}
7	+	\usepackage{setspace}
8	+	\usepackage{tabularx}
9		\usepackage{graphicx}
10	<
6	<	%\usepackage{endfloat}
10	>	\usepackage[ref]{overcite}
11		%\usepackage{berkeley}
8	–	%\usepackage{epsf}
9	–	%\usepackage[ref]{overcite}
10	–	%\usepackage{setspace}
11	–	%\usepackage{tabularx}
12	–	%\usepackage{graphicx}
12		%\usepackage{curves}
13	<	%\usepackage{amsmath}
14	<	%\pagestyle{plain}
15	<	%\pagenumbering{arabic}
16	<	%\oddsidemargin 0.0cm \evensidemargin 0.0cm
17	<	%\topmargin -21pt \headsep 10pt
18	<	%\textheight 9.0in \textwidth 6.5in
19	<	%\brokenpenalty=10000
20	<	%\renewcommand{\baselinestretch}{1.2}
22	<	%\renewcommand\citemid{\ } % no comma in optional reference note
13	>	\pagestyle{plain}
14	>	\pagenumbering{arabic}
15	>	\oddsidemargin 0.0cm \evensidemargin 0.0cm
16	>	\topmargin -21pt \headsep 10pt
17	>	\textheight 9.0in \textwidth 6.5in
18	>	\brokenpenalty=10000
19	>	\renewcommand{\baselinestretch}{1.2}
20	>	\renewcommand\citemid{\ } % no comma in optional reference note
21
22		\begin{document}
23
24	<	\title{On the temperature dependent structural and transport 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 ({\sc ssd}) and related single point water models}
26
27	<	\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
28	<	\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
30	<
31	<	\address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
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\\
29		Notre Dame, Indiana 46556}
30
31		\date{\today}
32
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 the SSD water model, both with and without the use of reaction
39	<	field. The constant pressure simulations of the melting of both $I_h$
40	<	and $I_c$ ice showed a density maximum near 260 K. In most cases, the
41	<	calculated densities were significantly lower than the densities
42	<	calculated in simulations of other water models. Analysis of particle
43	<	diffusion showed SSD to capture the transport properties of
44	<	experimental very well in both the normal and super-cooled liquid
45	<	regimes. In order to correct the density behavior, SSD was
46	<	reparameterized for use both with and without a long-range interaction
47	<	correction, SSD/RF and SSD/E respectively. In addition to correcting
48	<	the abnormally low densities, these new versions were show to maintain
49	<	or improve upon the transport and structural features of the original
50	<	water model.
36	>	The density maximum and temperature dependence of the self-diffusion
37	>	constant were investigated for the soft sticky dipole ({\sc 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 {\sc ssd} model captures the transport properties of
48	>	experimental water very well in both the normal and super-cooled
49	>	liquid regimes.  We also present our re-parameterized versions of {\sc ssd}
50	>	for use both with the reaction field or without any long-range
51	>	electrostatic corrections.  These are called the {\sc ssd/rf} and {\sc 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 {\sc ssd}
55	>	or the density corrected version of the original model ({\sc ssd1}).
56	>	Additionally, a novel low-density ice structure is presented
57	>	which appears to be the most stable ice structure for the entire {\sc ssd}
58	>	family.
59		\end{abstract}
60
61	<	\maketitle
61	>	\newpage
62
63		%\narrowtext
64
#	Line 62 | Line 69 | One of the most important tasks in simulations of bioc
69
70		\section{Introduction}
71
72	<	One of the most important tasks in simulations of biochemical systems
73	<	is the proper depiction of water and water solvation. In fact, the
74	<	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
78	<	groups/bulk. Due to the fact that explicit solvent accounts for a
79	<	massive portion of the calculations, it necessary to simplify the
73	<	solvent to some extent in order to complete simulations in a
74	<	reasonable amount of time. In the case of simulating water in
75	<	bio-molecular studies, the balance between accurate properties and
76	<	computational efficiency is especially delicate, and it has resulted
77	<	in a variety of different water
78	<	models.\cite{Jorgensen83,Berendsen87,Jorgensen00} Many of these models
79	<	get specific properties correct or better than their predecessors, but
80	<	this is often at a cost of some other properties or of computer
81	<	time. As an example, compare TIP3P or TIP4P to TIP5P. TIP5P succeeds
82	<	in improving the structural and transport properties over its
83	<	predecessors, yet this comes at a greater than 50\% increase in
84	<	computational cost.\cite{Jorgensen01,Jorgensen00} One recently
85	<	developed model that succeeds in both retaining accuracy of system
86	<	properties and simplifying calculations to increase computational
87	<	efficiency is the Soft Sticky Dipole water model.\cite{Ichiye96}
72	>	One of the most important tasks in the simulation of biochemical
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 ({\sc ssd}) water
93	>	model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The {\sc 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} {\sc 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 {\sc 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
108	<	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}
116	–	\begin{split}
137		u_{ij}^{sp}
138	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)
139	<	&=
140	<	\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
141	<	& \quad \ +
142	<	s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
123	<	\end{split}
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, $s$ and
147	<	$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 {\sc 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 a
168	<	calculation speed up of up to an order of magnitude over other
169	<	comparable models while maintaining the structural behavior of
170	<	water.\cite{Ichiye96} In the original molecular dynamics studies of
171	<	SSD, it was shown that it actually improves upon the prediction of
172	<	water's dynamical properties 3 and 4-point models.\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, a specific
176	<	interest within our group.
165	>	Since {\sc 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	>	{\sc 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, {\sc 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 {\sc 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 {\sc ssd} a very attractive model for the
179	>	simulation of large scale biochemical simulations.
180
181	<	Up to this point, a detailed look at the model's structure and ion
182	<	solvation abilities has been performed.\cite{Ichiye96} In addition, a
183	<	thorough investigation of the dynamic properties of SSD was performed
184	<	by Chandra and Ichiye focusing on translational and orientational
185	<	properties at 298 K.\cite{Ichiye99} This study focuses on determining
186	<	the density maximum for SSD utilizing both microcanonical and
187	<	isobaric-isothermal ensemble molecular dynamics, while using the
188	<	reaction field method for handling long-ranged dipolar interactions. A
189	<	reaction field method has been previously implemented in Monte Carlo
190	<	simulations by Liu and Ichiye in order to study the static dielectric
191	<	constant for the model.\cite{Ichiye96b} This paper will expand the
192	<	scope of these original simulations to look on how the reaction field
193	<	affects the physical and dynamic properties of SSD systems.
181	>	One feature of the {\sc 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 {\sc 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 {\sc ssd} water model have suggested
201	>	adjustments to the {\sc ssd} water model to address abnormal density
202	>	behavior (also observed here), calling the corrected model
203	>	{\sc ssd1}.\cite{Ichiye03} In what follows, we compare our
204	>	reparamaterization of {\sc ssd} with both the original {\sc ssd} and {\sc ssd1} models
205	>	with the goal of improving the bulk phase behavior of an {\sc ssd}-derived
206	>	model in simulations utilizing the Reaction Field.
207
208		\section{Methods}
209
210	<	As stated previously, in this study the long-range dipole-dipole
211	<	interactions were accounted for 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} s(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 $s(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
195	<	large-scale system, the computational cost benefit of reaction field
196	<	is dramatic. To address some of the dynamical property alterations due
197	<	to the use of reaction field, simulations were also performed without
198	<	a surrounding dielectric and suggestions are proposed on how to make
199	<	SSD more compatible with 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	<	Simulations were performed in both the isobaric-isothermal and
237	<	microcanonical ensembles. The constant pressure simulations were
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 {\sc ssd} which could be used either with or without the reaction field
240	>	turned on.
241	>
242	>	Simulations to obtain the preferred densities of the models were
243	>	performed in the isobaric-isothermal (NPT) ensemble, while all
244	>	dynamical properties were obtained from microcanonical (NVE)
245	>	simulations done at densities matching the NPT density for a
246	>	particular target temperature.  The constant pressure simulations were
247		implemented using an integral thermostat and barostat as outlined by
248	<	Hoover.\cite{Hoover85,Hoover86} For the constant pressure
249	<	simulations, the \emph{Q} parameter for the was set to 5.0 amu
250	<	\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at
251	<	100 ps.
248	>	Hoover.\cite{Hoover85,Hoover86} All molecules were treated as
249	>	non-linear rigid bodies. Vibrational constraints are not necessary in
250	>	simulations of {\sc ssd}, because there are no explicit hydrogen atoms, and
251	>	thus no molecular vibrational modes 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 actually quite sensitive to errors in the
260	<	equations of motion. The original implementation of this code utilized
261	<	quaternions for rotational motion propagation; however, a detailed
218	<	investigation showed that they resulted in a steady drift in the total
219	<	energy, something that has been observed by others.\cite{Laird97}
254	>	symplectic splitting method proposed by Dullweber, Leimkuhler, and
255	>	McLachlan ({\sc dlm}).\cite{Dullweber1997} Our reason for selecting this
256	>	integrator centers on poor energy conservation of rigid body dynamics
257	>	using traditional quaternion integration.\cite{Evans77,Evans77b} In
258	>	typical microcanonical ensemble simulations, the energy drift when
259	>	using quaternions was substantially greater than when using the {\sc dlm}
260	>	method (fig. \ref{timestep}).  This steady drift in the total energy
261	>	has also been observed by Kol {\it et al.}\cite{Laird97}
262
263		The key difference in the integration method proposed by Dullweber
264		\emph{et al.} is that the entire rotation matrix is propagated from
265	<	one time step to the next. In the past, this would not have been as
266	<	feasible a option, being that the rotation matrix for a single body is
267	<	nine elements long as opposed to 3 or 4 elements for Euler angles and
268	<	quaternions respectively. System memory has become much less of an
227	<	issue in recent times, and this has resulted in substantial benefits
228	<	in energy conservation. There is still the issue of an additional 5 or
229	<	6 additional elements for describing the orientation of each particle,
230	<	which will increase dump files substantially. Simply translating the
231	<	rotation matrix into its component Euler angles or quaternions for
232	<	storage purposes relieves this burden.
265	>	one time step to the next.  The additional memory required by the
266	>	algorithm is inconsequential on modern computers, and translating the
267	>	rotation matrix into quaternions for storage purposes makes trajectory
268	>	data quite compact.
269
270	<	The symplectic splitting method allows for Verlet style integration of
271	<	both linear and angular motion of rigid bodies. In the integration
272	<	method, the orientational propagation involves a sequence of matrix
273	<	evaluations to update the rotation matrix.\cite{Dullweber1997} These
274	<	matrix rotations end up being more costly computationally than the
275	<	simpler arithmetic quaternion propagation. On average, a 1000 SSD
276	<	particle simulation shows a 7\% increase in simulation time using the
277	<	symplectic step method in place of quaternions. This cost is more than
278	<	justified when comparing the energy conservation of the two methods as
279	<	illustrated in figure \ref{timestep}.
270	>	The {\sc dlm} method allows for Verlet style integration of both
271	>	translational and orientational motion of rigid bodies. In this
272	>	integration method, the orientational propagation involves a sequence
273	>	of matrix evaluations to update the rotation
274	>	matrix.\cite{Dullweber1997} These matrix rotations are more costly
275	>	than the simpler arithmetic quaternion propagation. With the same time
276	>	step, a 1000 {\sc ssd} particle simulation shows an average 7\% increase in
277	>	computation time using the {\sc dlm} method in place of quaternions. The
278	>	additional expense per step is justified when one considers the
279	>	ability to use time steps that are nearly twice as large under {\sc dlm}
280	>	than would be usable under quaternion dynamics.  The energy
281	>	conservation of the two methods using a number of different time steps
282	>	is illustrated in figure
283	>	\ref{timestep}.
284
285		\begin{figure}
286	<	\includegraphics[width=61mm, angle=-90]{timeStep.epsi}
287	<	\caption{Energy conservation using quaternion based integration versus
288	<	the symplectic step method proposed by Dullweber \emph{et al.} with
289	<	increasing time step. For each time step, the dotted line is total
290	<	energy using the symplectic step integrator, and the solid line comes
291	<	from the quaternion integrator. The larger time step plots are shifted
292	<	up from the true energy baseline for clarity.}
286	>	\begin{center}
287	>	\epsfxsize=6in
288	>	\epsfbox{timeStep.epsi}
289	>	\caption{Energy conservation using both quaternion-based integration and
290	>	the {\sc dlm} method with increasing time step. The larger time step plots
291	>	are shifted from the true energy baseline (that of $\Delta t$ = 0.1
292	>	fs) for clarity.}
293		\label{timestep}
294	+	\end{center}
295		\end{figure}
296
297		In figure \ref{timestep}, the resulting energy drift at various time
298	<	steps for both the symplectic step and quaternion integration schemes
299	<	is compared. All of the 1000 SSD particle simulations started with the
300	<	same configuration, and the only difference was the method for
301	<	handling rotational motion. At time steps of 0.1 and 0.5 fs, both
302	<	methods for propagating particle rotation conserve energy fairly well,
303	<	with the quaternion method showing a slight energy drift over time in
304	<	the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
305	<	energy conservation benefits of the symplectic step method are clearly
306	<	demonstrated.
307	<
308	<	Energy drift in these SSD particle simulations was unnoticeable for
268	<	time steps up to three femtoseconds. A slight energy drift on the
269	<	order of 0.012 kcal/mol per nanosecond was observed at a time step of
270	<	four femtoseconds, and as expected, this drift increases dramatically
271	<	with increasing time step. To insure accuracy in the constant energy
272	<	simulations, time steps were set at 2 fs and kept at this value for
273	<	constant pressure simulations as well.
298	>	steps for both the {\sc dlm} and quaternion integration schemes is compared.
299	>	All of the 1000 {\sc ssd} particle simulations started with the same
300	>	configuration, and the only difference was the method used to handle
301	>	orientational motion. At time steps of 0.1 and 0.5 fs, both methods
302	>	for propagating the orientational degrees of freedom conserve energy
303	>	fairly well, with the quaternion method showing a slight energy drift
304	>	over time in the 0.5 fs time step simulation. At time steps of 1 and 2
305	>	fs, the energy conservation benefits of the {\sc dlm} method are clearly
306	>	demonstrated. Thus, while maintaining the same degree of energy
307	>	conservation, one can take considerably longer time steps, leading to
308	>	an overall reduction in computation time.
309
310	<	Ice crystals in both the $I_h$ and $I_c$ lattices were generated as
311	<	starting points for all the simulations. The $I_h$ crystals were
312	<	formed by first arranging the center of masses of the SSD particles
313	<	into a ``hexagonal'' ice lattice of 1024 particles. Because of the
314	<	crystal structure of $I_h$ ice, the simulation box assumed a
315	<	rectangular shape with a edge length ratio of approximately
310	>	Energy drift in the simulations using {\sc dlm} integration was unnoticeable
311	>	for time steps up to 3 fs. A slight energy drift on the order of 0.012
312	>	kcal/mol per nanosecond was observed at a time step of 4 fs, and as
313	>	expected, this drift increases dramatically with increasing time
314	>	step. To insure accuracy in our microcanonical simulations, time steps
315	>	were set at 2 fs and kept at this value for constant pressure
316	>	simulations as well.
317	>
318	>	Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
319	>	were generated as starting points for all simulations. The $I_h$
320	>	crystals were formed by first arranging the centers of mass of the {\sc ssd}
321	>	particles into a ``hexagonal'' ice lattice of 1024 particles. Because
322	>	of the crystal structure of $I_h$ ice, the simulation box assumed an
323	>	orthorhombic shape with an edge length ratio of approximately
324		1.00$\times$1.06$\times$1.23. The particles were then allowed to
325		orient freely about fixed positions with angular momenta randomized at
326		400 K for varying times. The rotational temperature was then scaled
327	<	down in stages to slowly cool the crystals down to 25 K. The particles
328	<	were then allowed translate with fixed orientations at a constant
327	>	down in stages to slowly cool the crystals to 25 K. The particles were
328	>	then allowed to translate with fixed orientations at a constant
329		pressure of 1 atm for 50 ps at 25 K. Finally, all constraints were
330		removed and the ice crystals were allowed to equilibrate for 50 ps at
331		25 K and a constant pressure of 1 atm.  This procedure resulted in
332		structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
333	<	rules\cite{Bernal33,Rahman72}.  This method was also utilized in the
333	>	rules.\cite{Bernal33,Rahman72} This method was also utilized in the
334		making of diamond lattice $I_c$ ice crystals, with each cubic
335		simulation box consisting of either 512 or 1000 particles. Only
336		isotropic volume fluctuations were performed under constant pressure,
#	Line 297 | Line 340 | constant pressure and temperature dynamics. This invol
340		\section{Results and discussion}
341
342		Melting studies were performed on the randomized ice crystals using
343	<	constant pressure and temperature dynamics. This involved an initial
344	<	randomization of velocities about the starting temperature of 25 K for
345	<	varying amounts of time. The systems were all equilibrated for 100 ps
346	<	prior to a 200 ps data collection run at each temperature setting,
347	<	ranging from 25 to 400 K, with a maximum degree increment of 25 K. For
348	<	regions of interest along this stepwise progression, the temperature
349	<	increment was decreased from 25 K to 10 and then 5 K. The above
350	<	equilibration and production times were sufficient in that the system
351	<	volume fluctuations dampened out in all but the very cold simulations
352	<	(below 225 K). In order to further improve statistics, five separate
353	<	simulation progressions were performed, and the averaged results from
354	<	the $I_h$ melting simulations are shown in figure \ref{dense1}.
343	>	isobaric-isothermal (NPT) dynamics. During melting simulations, the
344	>	melting transition and the density maximum can both be observed,
345	>	provided that the density maximum occurs in the liquid and not the
346	>	supercooled regime. An ensemble average from five separate melting
347	>	simulations was acquired, each starting from different ice crystals
348	>	generated as described previously. All simulations were equilibrated
349	>	for 100 ps prior to a 200 ps data collection run at each temperature
350	>	setting. The temperature range of study spanned from 25 to 400 K, with
351	>	a maximum degree increment of 25 K. For regions of interest along this
352	>	stepwise progression, the temperature increment was decreased from 25
353	>	K to 10 and 5 K.  The above equilibration and production times were
354	>	sufficient in that fluctuations in the volume autocorrelation function
355	>	were damped out in all simulations in under 20 ps.
356
313	–	\begin{figure}
314	–	\includegraphics[width=65mm, angle=-90]{1hdense.epsi}
315	–	\caption{Average density of SSD water at increasing temperatures
316	–	starting from ice $I_h$ lattice.}
317	–	\label{dense1}
318	–	\end{figure}
319	–
357		\subsection{Density Behavior}
321	–	In the initial average density versus temperature plot, the density
322	–	maximum clearly appears between 255 and 265 K. The calculated
323	–	densities within this range were nearly indistinguishable, as can be
324	–	seen in the zoom of this region of interest, shown in figure
325	–	\ref{dense1}. The greater certainty of the average value at 260 K makes
326	–	a good argument for the actual density maximum residing at this
327	–	midpoint value. Figure \ref{dense1} was constructed using ice $I_h$
328	–	crystals for the initial configuration; and though not pictured, the
329	–	simulations starting from ice $I_c$ crystal configurations showed
330	–	similar results, with a liquid-phase density maximum in this same
331	–	region (between 255 and 260 K). In addition, the $I_c$ crystals are
332	–	more fragile than the $I_h$ crystals, leading them to deform into a
333	–	dense glassy state at lower temperatures. This resulted in an overall
334	–	low temperature density maximum at 200 K, but they still retained a
335	–	common liquid state density maximum with the $I_h$ simulations.
358
359	+	Our initial simulations focused on the original {\sc ssd} water model, and
360	+	an average density versus temperature plot is shown in figure
361	+	\ref{dense1}. Note that the density maximum when using a reaction
362	+	field appears between 255 and 265 K.  There were smaller fluctuations
363	+	in the density at 260 K than at either 255 or 265, so we report this
364	+	value as the location of the density maximum. Figure \ref{dense1} was
365	+	constructed using ice $I_h$ crystals for the initial configuration;
366	+	though not pictured, the simulations starting from ice $I_c$ crystal
367	+	configurations showed similar results, with a liquid-phase density
368	+	maximum in this same region (between 255 and 260 K).
369	+
370		\begin{figure}
371	<	\includegraphics[width=65mm,angle=-90]{dense2.eps}
372	<	\caption{Density versus temperature for TIP4P\cite{Jorgensen98b},
373	<	TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction
374	<	Field, SSD, and Experiment\cite{CRC80}. }
375	<	\label{dense2}
371	>	\begin{center}
372	>	\epsfxsize=6in
373	>	\epsfbox{denseSSD.eps}
374	>	\caption{Density versus temperature for TIP4P [Ref. \citen{Jorgensen98b}],
375	>	TIP3P [Ref. \citen{Jorgensen98b}], SPC/E [Ref. \citen{Clancy94}], {\sc ssd}
376	>	without Reaction Field, {\sc ssd}, and experiment [Ref. \citen{CRC80}]. The
377	>	arrows indicate the change in densities observed when turning off the
378	>	reaction field. The the lower than expected densities for the {\sc ssd}
379	>	model were what prompted the original reparameterization of {\sc ssd1}
380	>	[Ref. \citen{Ichiye03}].}
381	>	\label{dense1}
382	>	\end{center}
383		\end{figure}
384
385	<	The density maximum for SSD actually compares quite favorably to other
386	<	simple water models. Figure \ref{dense2} shows a plot of these
387	<	findings with the density progression of several other models and
348	<	experiment obtained from other
385	>	The density maximum for {\sc ssd} compares quite favorably to other simple
386	>	water models. Figure \ref{dense1} also shows calculated densities of
387	>	several other models and experiment obtained from other
388		sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
389	<	models, SSD has results closest to the experimentally observed water
390	<	density maximum. Of the listed water models, TIP4P has a density
391	<	maximum behavior most like that seen in SSD. Though not shown, it is
392	<	useful to note that TIP5P has a water density maximum nearly identical
393	<	to experiment.
389	>	models, {\sc ssd} has a temperature closest to the experimentally observed
390	>	density maximum. Of the {\it charge-based} models in
391	>	Fig. \ref{dense1}, TIP4P has a density maximum behavior most like that
392	>	seen in {\sc ssd}. Though not included in this plot, it is useful
393	>	to note that TIP5P has a density maximum nearly identical to the
394	>	experimentally measured temperature.
395
396	<	Possibly of more importance is the density scaling of SSD relative to
397	<	other common models at any given temperature (Fig. \ref{dense2}). Note
398	<	that the SSD model assumes a lower density than any of the other
399	<	listed models at the same pressure, behavior which is especially
400	<	apparent at temperatures greater than 300 K. Lower than expected
401	<	densities have been observed for other systems with the use of a
402	<	reaction field for long-range electrostatic interactions, so the most
403	<	likely reason for these significantly lower densities in these
404	<	simulations is the presence of the reaction field.\cite{Berendsen98}
405	<	In order to test the effect of the reaction field on the density of
406	<	the systems, the simulations were repeated for the temperature region
407	<	of interest without a reaction field present. The results of these
408	<	simulations are also displayed in figure \ref{dense2}. Without
369	<	reaction field, these densities increase considerably to more
370	<	experimentally reasonable values, especially around the freezing point
371	<	of liquid water. The shape of the curve is similar to the curve
372	<	produced from SSD simulations using reaction field, specifically the
373	<	rapidly decreasing densities at higher temperatures; however, a slight
374	<	shift in the density maximum location, down to 245 K, is
375	<	observed. This is probably a more accurate comparison to the other
376	<	listed water models in that no long range corrections were applied in
377	<	those simulations.\cite{Clancy94,Jorgensen98b}
396	>	It has been observed that liquid state densities in water are
397	>	dependent on the cutoff radius used both with and without the use of
398	>	reaction field.\cite{Berendsen98} In order to address the possible
399	>	effect of cutoff radius, simulations were performed with a dipolar
400	>	cutoff radius of 12.0 \AA\ to complement the previous {\sc ssd} simulations,
401	>	all performed with a cutoff of 9.0 \AA. All of the resulting densities
402	>	overlapped within error and showed no significant trend toward lower
403	>	or higher densities as a function of cutoff radius, for simulations
404	>	both with and without reaction field. These results indicate that
405	>	there is no major benefit in choosing a longer cutoff radius in
406	>	simulations using {\sc ssd}. This is advantageous in that the use of a
407	>	longer cutoff radius results in a significant increase in the time
408	>	required to obtain a single trajectory.
409
410	<	It has been observed that densities are dependent on the cutoff radius
411	<	used for a variety of water models in simulations both with and
412	<	without the use of reaction field.\cite{Berendsen98} In order to
413	<	address the possible affect of cutoff radius, simulations were
414	<	performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the
415	<	previous SSD simulations, all performed with a cutoff of 9.0 \AA. All
416	<	the resulting densities overlapped within error and showed no
417	<	significant trend in lower or higher densities as a function of cutoff
418	<	radius, both for simulations with and without reaction field. These
419	<	results indicate that there is no major benefit in choosing a longer
420	<	cutoff radius in simulations using SSD. This is comforting in that the
421	<	use of a longer cutoff radius results in a near doubling of the time
422	<	required to compute a single trajectory.
410	>	The key feature to recognize in figure \ref{dense1} is the density
411	>	scaling of {\sc ssd} relative to other common models at any given
412	>	temperature. {\sc ssd} assumes a lower density than any of the other listed
413	>	models at the same pressure, behavior which is especially apparent at
414	>	temperatures greater than 300 K. Lower than expected densities have
415	>	been observed for other systems using a reaction field for long-range
416	>	electrostatic interactions, so the most likely reason for the
417	>	significantly lower densities seen in these simulations is the
418	>	presence of the reaction field.\cite{Berendsen98,Nezbeda02} In order
419	>	to test the effect of the reaction field on the density of the
420	>	systems, the simulations were repeated without a reaction field
421	>	present. The results of these simulations are also displayed in figure
422	>	\ref{dense1}. Without the reaction field, the densities increase
423	>	to more experimentally reasonable values, especially around the
424	>	freezing point of liquid water. The shape of the curve is similar to
425	>	the curve produced from {\sc ssd} simulations using reaction field,
426	>	specifically the rapidly decreasing densities at higher temperatures;
427	>	however, a shift in the density maximum location, down to 245 K, is
428	>	observed. This is a more accurate comparison to the other listed water
429	>	models, in that no long range corrections were applied in those
430	>	simulations.\cite{Clancy94,Jorgensen98b} However, even without the
431	>	reaction field, the density around 300 K is still significantly lower
432	>	than experiment and comparable water models. This anomalous behavior
433	>	was what lead Tan {\it et al.} to recently reparameterize
434	>	{\sc ssd}.\cite{Ichiye03} Throughout the remainder of the paper our
435	>	reparamaterizations of {\sc ssd} will be compared with their newer {\sc ssd1}
436	>	model.
437
438		\subsection{Transport Behavior}
439	<	Of importance in these types of studies are the transport properties
440	<	of the particles and how they change when altering the environmental
441	<	conditions. In order to probe transport, constant energy simulations
442	<	were performed about the average density uncovered by the constant
443	<	pressure simulations. Simulations started with randomized velocities
444	<	and underwent 50 ps of temperature scaling and 50 ps of constant
445	<	energy equilibration before obtaining a 200 ps trajectory. Diffusion
446	<	constants were calculated via root-mean square deviation analysis. The
447	<	averaged results from 5 sets of these NVE simulations is displayed in
448	<	figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
449	<	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,Holz00,Clancy94,Jorgensen01}
453
454		\begin{figure}
455	<	\includegraphics[width=65mm, angle=-90]{betterDiffuse.epsi}
456	<	\caption{Average diffusion coefficient over increasing temperature for
457	<	SSD, SPC/E\cite{Clancy94}, TIP5P\cite{Jorgensen01}, and Experimental
458	<	data from Gillen \emph{et al.}\cite{Gillen72}, and from
459	<	Mills\cite{Mills73}.}
455	>	\begin{center}
456	>	\epsfxsize=6in
457	>	\epsfbox{betterDiffuse.epsi}
458	>	\caption{Average self-diffusion constant as a function of temperature for
459	>	{\sc ssd}, SPC/E [Ref. \citen{Clancy94}], and TIP5P
460	>	[Ref. \citen{Jorgensen01}] compared with experimental data
461	>	[Refs. \citen{Gillen72} and \citen{Holz00}]. Of the three water models
462	>	shown, {\sc ssd} has the least deviation from the experimental values. The
463	>	rapidly increasing diffusion constants for TIP5P and {\sc ssd} correspond to
464	>	significant decreases in density at the 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 normal regimes. SPC/E
473	<	does a respectable job by getting similar values as SSD and experiment
474	<	around 290 K; however, it deviates at both higher and lower
475	<	temperatures, failing to predict the experimental trend. TIP5P and SSD
476	<	both start off low at the colder temperatures and tend to diffuse too
477	<	rapidly at the higher temperatures. This type of trend at the higher
478	<	temperatures is not surprising in that the densities of both TIP5P and
479	<	SSD are lower than experimental water at temperatures higher than room
480	<	temperature. 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
429	<	pressure. Results under these conditions can be found later in this
430	<	paper.
470	>	strengths of the {\sc ssd} model. Of the three models shown, the {\sc 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 {\sc 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 {\sc ssd} are lower than experimental
479	>	water densities at higher temperatures.  When calculating the
480	>	diffusion coefficients for {\sc 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
488	<	phase behavior beyond the melting point. To locate the melting
489	<	transition, the constant pressure heat capacity (C$_\text{p}$) was
490	<	monitored in each of the simulations. In the melting simulations of
491	<	the 1024 particle ice $I_h$ simulations, a large spike in C$_\text{p}$
492	<	occurs at 245 K, indicating a first order phase transition for the
493	<	melting of these ice crystals. When the reaction field is turned off,
494	<	the melting transition occurs at 235 K.  These melting transitions are
495	<	considerably lower than the experimental value, but this is not
443	<	surprising in that SSD is a simple rigid body model with a fixed
444	<	dipole.
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
491	>	particle ice $I_h$ simulations, a large spike in C$_\text{p}$ occurs
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.
496
446	–	\begin{figure}
447	–	\includegraphics[width=85mm]{fullContours.eps}
448	–	\caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at
449	–	100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for
450	–	clarity: dark areas signify peaks while light areas signify
451	–	depressions. White areas have g(\emph{r}) values below 0.5 and black
452	–	areas have values above 1.5.}
453	–	\label{contour}
454	–	\end{figure}
455	–
456	–	Additional analyses for understanding the melting phase-transition
457	–	process were performed via two-dimensional structure and dipole angle
458	–	correlations. Expressions for the correlations are as follows:
459	–
497		\begin{figure}
498	<	\includegraphics[width=45mm]{corrDiag.eps}
499	<	\caption{Two dimensional illustration of the angles involved in the
500	<	correlations observed in figure \ref{contour}.}
498	>	\begin{center}
499	>	\epsfxsize=6in
500	>	\epsfbox{corrDiag.eps}
501	>	\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
502		\label{corrAngle}
503	+	\end{center}
504		\end{figure}
505
506	<	\begin{multline}
507	<	g_{\text{AB}}(r,\cos\theta) = \\
508	<	\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\ ,
509	<	\end{multline}
510	<	\begin{multline}
511	<	g_{\text{AB}}(r,\cos\omega) = \\
512	<	\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\ ,
513	<	\end{multline}
514	<	where $\theta$ and $\omega$ refer to the angles shown in the above
515	<	illustration. By binning over both distance and the cosine of the
516	<	desired angle between the two dipoles, the g(\emph{r}) can be
517	<	dissected to determine the common dipole arrangements that constitute
518	<	the peaks and troughs. Frames A and B of figure \ref{contour} show a
519	<	relatively crystalline state of an ice $I_c$ simulation. The first
520	<	peak of the g(\emph{r}) primarily consists of the preferred hydrogen
521	<	bonding arrangements as dictated by the tetrahedral sticky potential,
522	<	one peak for the donating and the other for the accepting hydrogen
523	<	bonds. Due to the high degree of crystallinity of the sample, the
524	<	second and third solvation shells show a repeated peak arrangement
506	>	\begin{figure}
507	>	\begin{center}
508	>	\epsfxsize=6in
509	>	\epsfbox{fullContours.eps}
510	>	\caption{Contour plots of 2D angular pair correlation functions for
511	>	512 {\sc ssd} molecules at 100 K (A \& B) and 300 K (C \& D). Dark areas
512	>	signify regions of enhanced density while light areas signify
513	>	depletion relative to the bulk density. White areas have pair
514	>	correlation values below 0.5 and black areas have values above 1.5.}
515	>	\label{contour}
516	>	\end{center}
517	>	\end{figure}
518	>
519	>	Additional analysis of the melting process was performed using
520	>	two-dimensional structure and dipole angle correlations. Expressions
521	>	for these correlations are as follows:
522	>
523	>	\begin{equation}
524	>	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\ ,
525	>	\end{equation}
526	>	\begin{equation}
527	>	g_{\text{AB}}(r,\cos\omega) =
528	>	\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\ ,
529	>	\end{equation}
530	>	where $\theta$ and $\omega$ refer to the angles shown in figure
531	>	\ref{corrAngle}. By binning over both distance and the cosine of the
532	>	desired angle between the two dipoles, the $g(r)$ can be analyzed to
533	>	determine the common dipole arrangements that constitute the peaks and
534	>	troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
535	>	of figure \ref{contour} show results from an ice $I_c$ simulation. The
536	>	first peak in the $g(r)$ consists primarily of the preferred hydrogen
537	>	bonding arrangements as dictated by the tetrahedral sticky potential -
538	>	one peak for the hydrogen bond donor and the other for the hydrogen
539	>	bond acceptor.  Due to the high degree of crystallinity of the sample,
540	>	the second and third solvation shells show a repeated peak arrangement
541		which decays at distances around the fourth solvation shell, near the
542		imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
543	<	In the higher temperature simulation shown in frames C and D, the
544	<	repeated peak features are significantly blurred. The first solvation
545	<	shell still shows the strong effect of the sticky-potential, although
546	<	it covers a larger area, extending to include a fraction of aligned
543	>	In the higher temperature simulation shown in frames C and D, these
544	>	long-range features deteriorate rapidly. The first solvation shell
545	>	still shows the strong effect of the sticky-potential, although it
546	>	covers a larger area, extending to include a fraction of aligned
547		dipole peaks within the first solvation shell. The latter peaks lose
548	<	definition as thermal motion and the competing dipole force overcomes
549	<	the sticky potential's tight tetrahedral structuring of the fluid.
548	>	due to thermal motion and as the competing dipole force overcomes the
549	>	sticky potential's tight tetrahedral structuring of the crystal.
550
551		This complex interplay between dipole and sticky interactions was
552		remarked upon as a possible reason for the split second peak in the
553	<	oxygen-oxygen g(\emph{r}).\cite{Ichiye96} At low temperatures, the
554	<	second solvation shell peak appears to have two distinct parts that
553	>	oxygen-oxygen pair correlation function,
554	>	$g_\mathrm{OO}(r)$.\cite{Ichiye96} At low temperatures, the second
555	>	solvation shell peak appears to have two distinct components that
556		blend together to form one observable peak. At higher temperatures,
557		this split character alters to show the leading 4 \AA\ peak dominated
558	<	by equatorial anti-parallel dipole orientations, and there is tightly
559	<	bunched group of axially arranged dipoles that most likely consist of
560	<	the smaller fraction aligned dipole pairs. The trailing part of the
561	<	split peak at 5 \AA\ is dominated by aligned dipoles that range
562	<	primarily within the axial to the chief hydrogen bond arrangements
563	<	similar to those seen in the first solvation shell. This evidence
564	<	indicates that the dipole pair interaction begins to dominate outside
565	<	of the range of the dipolar repulsion term, with the primary
566	<	energetically favorable dipole arrangements populating the region
567	<	immediately outside of it's range (around 4 \AA), and arrangements
568	<	that seek to ideally satisfy both the sticky and dipole forces locate
569	<	themselves just beyond this region (around 5 \AA).
558	>	by equatorial anti-parallel dipole orientations. There is also a
559	>	tightly bunched group of axially arranged dipoles that most likely
560	>	consist of the smaller fraction of aligned dipole pairs. The trailing
561	>	component of the split peak at 5 \AA\ is dominated by aligned dipoles
562	>	that assume hydrogen bond arrangements similar to those seen in the
563	>	first solvation shell. This evidence indicates that the dipole pair
564	>	interaction begins to dominate outside of the range of the dipolar
565	>	repulsion term.  The energetically favorable dipole arrangements
566	>	populate the region immediately outside this repulsion region (around
567	>	4 \AA), while arrangements that seek to satisfy both the sticky and
568	>	dipole forces locate themselves just beyond this initial buildup
569	>	(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
573	<	leading of the two peaks can be pushed out and merged with the outer
574	<	split peak just by extending the switching function cutoff
575	<	($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even
576	<	5 \AA. This type of correction is not recommended for improving the
577	<	liquid structure, because the second solvation shell will still be
578	<	shifted too far out. In addition, this would have an even more
579	<	detrimental effect on the system densities, leading to a liquid with a
580	<	more open structure and a density considerably lower than the normal
581	<	SSD behavior shown previously. A better correction would be to include
582	<	the quadrupole-quadrupole interactions for the water particles outside
583	<	of the first solvation shell, but this reduces the simplicity and
584	<	speed advantage of SSD, so it is not the most desirable path to take.
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 ($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 already low {\sc 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	>	{\sc ssd}.
585
586	<	\subsection{Adjusted Potentials: SSD/E and SSD/RF}
587	<	The propensity of SSD to adopt lower than expected densities under
586	>	\subsection{Adjusted Potentials: {\sc ssd/rf} and {\sc ssd/e}}
587	>
588	>	The propensity of {\sc ssd} to adopt lower than expected densities under
589		varying conditions is troubling, especially at higher temperatures. In
590	<	order to correct this behavior, it's necessary to adjust the force
591	<	field parameters for the primary intermolecular interactions. In
592	<	undergoing a reparameterization, it is important not to focus on just
593	<	one property and neglect the other important properties. In this case,
594	<	it would be ideal to correct the densities while maintaining the
595	<	accurate transport properties.
590	>	order to correct this model for use with a reaction field, it is
591	>	necessary to adjust the force field parameters for the primary
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 behavior.
596
597	<	The possible parameters 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	<	\begin{split}
605	<	u_{ij}^{sp}
606	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) &=
550	<	\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\\
551	<	& \quad \ + \frac{\nu_0^\prime}{2}
552	<	[s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)],
553	<	\end{split}
554	<	\end{equation}
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 cutoff distances
600	>	for the sticky attractive and dipole repulsive cubic switching
601	>	function cutoffs ($r_l$, $r_u$ and $r_l^\prime$, $r_u^\prime$
602	>	respectively). The results of the reparameterizations are shown in
603	>	table \ref{params}. We are calling these reparameterizations the Soft
604	>	Sticky Dipole / Reaction Field ({\sc ssd/rf} - for use with a reaction
605	>	field) and Soft Sticky Dipole Extended ({\sc ssd/e} - an attempt to improve
606	>	the liquid structure in simulations without a long-range correction).
607
556	–	where $\nu_0$ scales the strength of the tetrahedral attraction and
557	–	$\nu_0^\prime$ acts in an identical fashion on the dipole repulsion
558	–	term. For purposes of the reparameterization, the separation was
559	–	performed, but the final parameters were adjusted so that it is
560	–	unnecessary to separate the terms when implementing the adjusted water
561	–	potentials. The results of the reparameterizations are shown in table
562	–	\ref{params}. Note that both the tetrahedral attractive and dipolar
563	–	repulsive don't share the same lower cutoff ($r_l$) in the newly
564	–	parameterized potentials - soft sticky dipole enhanced (SSD/E) and
565	–	soft sticky dipole reaction field (SSD/RF).
566	–
608		\begin{table}
609	+	\begin{center}
610		\caption{Parameters for the original and adjusted models}
611	<	\begin{tabular}{ l  c  c  c }
611	>	\begin{tabular}{ l  c  c  c  c }
612		\hline \\[-3mm]
613	<	\ Parameters & \ \ \  SSD$^\dagger$\ \ \ \  & \ SSD/E\ \  & \ SSD/RF\ \ \\
613	>	\ \ \ Parameters\ \ \  & \ \ \ {\sc ssd} [Ref. \citen{Ichiye96}] \ \ \
614	>	& \ {\sc ssd1} [Ref. \citen{Ichiye03}]\ \  & \ {\sc ssd/e}\ \  & \ {\sc ssd/rf} \\
615		\hline \\[-3mm]
616	<	\ \ \ $\sigma$ (\AA)  & 3.051 & 3.035 & 3.019\\
617	<	\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\
618	<	\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\
619	<	\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
620	<	\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\
621	<	\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\
622	<	\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
623	<	\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\
624	<	\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\
582	<	\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96}
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\\
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}
626		\label{params}
627	+	\end{center}
628		\end{table}
629
630	<	\begin{figure}
631	<	\includegraphics[width=85mm]{gofrCompare.epsi}
632	<	\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E
633	<	and SSD without reaction field (top), as well as SSD/RF and SSD with
630	>	\begin{figure}
631	>	\begin{center}
632	>	\epsfxsize=5in
633	>	\epsfbox{GofRCompare.epsi}
634	>	\caption{Plots comparing experiment [Ref. \citen{Head-Gordon00_1}] with {\sc ssd/e}
635	>	and {\sc ssd1} without reaction field (top), as well as {\sc ssd/rf} and {\sc ssd1} with
636		reaction field turned on (bottom). The insets show the respective
637	<	first peaks in detail. Solid Line - experiment, dashed line - SSD/E
638	<	and SSD/RF, and dotted line - SSD (with and without reaction field).}
637	>	first peaks in detail. Note how the changes in parameters have lowered
638	>	and broadened the first peak of {\sc ssd/e} and {\sc ssd/rf}.}
639		\label{grcompare}
640	+	\end{center}
641		\end{figure}
642
643	<	\begin{figure}
644	<	\includegraphics[width=85mm]{dualsticky.ps}
645	<	\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \&
646	<	SSD/RF (right). Light areas correspond to the tetrahedral attractive
647	<	part, and the darker areas correspond to the dipolar repulsive part.}
643	>	\begin{figure}
644	>	\begin{center}
645	>	\epsfxsize=6in
646	>	\epsfbox{dualsticky_bw.eps}
647	>	\caption{Positive and negative isosurfaces of the sticky potential for
648	>	{\sc ssd1} (left) and {\sc ssd/e} \& {\sc ssd/rf} (right). Light areas correspond to the
649	>	tetrahedral attractive component, and darker areas correspond to the
650	>	dipolar repulsive component.}
651		\label{isosurface}
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 sharp first peak that
658	<	integrated to give similar coordination numbers to the experimental
659	<	data obtained by Soper and Phillips.\cite{Ichiye96,Soper86} New
660	<	experimental x-ray scattering data from the Head-Gordon lab indicates
661	<	a slightly lower and shifted first peak in the g$_\mathrm{OO}(r)$, so
662	<	adjustments to SSD were made while taking into consideration the new
663	<	experimental findings.\cite{Head-Gordon00_1} Figure \ref{grcompare}
664	<	shows the relocation of the first peak of the oxygen-oxygen
665	<	g(\emph{r}) by comparing the original SSD (with and without reaction
666	<	field), SSD-E, and SSD-RF to the new experimental results. Both the
667	<	modified water models have shorter peaks that are brought in more
668	<	closely to the experimental peak (as seen in the insets of figure
669	<	\ref{grcompare}). This structural alteration was accomplished by a
670	<	reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
671	<	of the sticky potential strength and cutoffs. The cutoffs for the
672	<	tetrahedral attractive and dipolar repulsive terms were nearly swapped
673	<	with each other. Isosurfaces of the original and modified sticky
674	<	potentials are shown in figure \cite{isosurface}. In these
675	<	isosurfaces, it is easy to see how altering the cutoffs changes the
676	<	repulsive and attractive character of the particles. With a reduced
677	<	repulsive surface (the darker region), the particles can move closer
678	<	to one another, increasing the density for the overall system. This
679	<	change in interaction cutoff also results in a more gradual
680	<	orientational motion by allowing the particles to maintain preferred
681	<	dipolar arrangements before they begin to feel the pull of the
682	<	tetrahedral restructuring. Upon moving closer together, the dipolar
683	<	repulsion term becomes active and excludes the unphysical
684	<	arrangements. This compares with the original SSD's excluding dipolar
685	<	before the particles feel the pull of the ``hydrogen bonds''. Aside
686	<	from improving the shape of the first peak in the g(\emph{r}), this
687	<	improves the densities considerably by allowing the persistence of
688	<	full dipolar character below the previous 4.0 \AA\ cutoff.
655	>	In the original paper detailing the development of {\sc 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 our adjustments to {\sc ssd} were
663	>	made after 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(r)$ by comparing
666	>	the revised {\sc ssd} model ({\sc ssd1}), {\sc ssd/e}, and {\sc ssd/rf} to the new
667	>	experimental results. Both modified water models have shorter peaks
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 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 {\sc ssd} and {\sc 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 to some degree, these changes alone
694	<	are insufficient to bring the system densities up to the values
695	<	observed experimentally. To finish bringing up the densities, the
696	<	dipole moments were increased in both the adjusted models. Being a
697	<	dipole based model, the structure and transport are very sensitive to
698	<	changes in the dipole moment. The original SSD simply used the dipole
699	<	moment 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 improve 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, so there is quite a range available for
705	<	adjusting the dipole
706	<	moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of
707	<	the dipole moments to 2.418 and 2.48 D for SSD/E and SSD/RF
708	<	respectively is moderate in the range of the experimental values;
657	<	however, it leads to significant changes in the density and transport
658	<	of the water models.
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 {\sc ssd} is a dipole based model, the
697	>	structure and transport are very sensitive to changes in the dipole
698	>	moment. The original {\sc 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 {\sc ssd/e} and {\sc ssd/rf},
707	>	respectively, leads to significant changes in the density and
708	>	transport of the water models.
709
710	<	In order to demonstrate the benefits of this reparameterization, a
710	>	In order to demonstrate the benefits of these reparameterizations, a
711		series of NPT and NVE simulations were performed to probe the density
712		and transport properties of the adapted models and compare the results
713	<	to the original SSD model. This comparison involved full NPT melting
714	<	sequences for both SSD/E and SSD/RF, as well as NVE transport
715	<	calculations at both self-consistent and experimental
716	<	densities. Again, the results come from five separate simulations of
717	<	1024 particle systems, and the melting sequences were started from
718	<	different ice $I_h$ crystals constructed as stated previously. Like
719	<	before, all of the NPT simulations were equilibrated for 100 ps before
720	<	a 200 ps data collection run, and they used the previous temperature's
721	<	final configuration as a starting point. All of the NVE simulations
722	<	had the same thermalization, equilibration, and data collection times
723	<	stated earlier in this paper.
713	>	to the original {\sc ssd} model. This comparison involved full NPT melting
714	>	sequences for both {\sc ssd/e} and {\sc ssd/rf}, as well as NVE transport
715	>	calculations at the calculated self-consistent densities. Again, the
716	>	results are obtained from five separate simulations of 1024 particle
717	>	systems, and the melting sequences were started from different ice
718	>	$I_h$ crystals constructed as described previously. Each NPT
719	>	simulation was equilibrated for 100 ps before a 200 ps data collection
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 previously.
724
725	<	\begin{figure}
726	<	\includegraphics[width=85mm]{ssdecompare.epsi}
727	<	\caption{Comparison of densities calculated with SSD/E to SSD without a
728	<	reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
729	<	SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
730	<	includes error bars, and the calculated results from the other
731	<	references were removed for clarity.}
725	>	\begin{figure}
726	>	\begin{center}
727	>	\epsfxsize=6in
728	>	\epsfbox{ssdeDense.epsi}
729	>	\caption{Comparison of densities calculated with {\sc ssd/e} to {\sc ssd1} without a
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 {\sc ssd1} and {\sc 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 water
741	<	model in comparison to the original SSD without a reaction field,
742	<	experiment, and the other common water models considered
743	<	previously. The calculated densities have increased significantly over
744	<	the original SSD model and match the experimental value just below 298
745	<	K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
746	<	compares well with the experimental value of 0.997 g/cm$^3$ and is
747	<	considerably better than the SSD value of 0.967$\pm$0.003
748	<	g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
749	<	out the curve at higher temperatures, only the improvement is marginal
750	<	at best. This steep drop in densities is due to the dipolar rather
751	<	than charge based interactions which decay more rapidly at longer
752	<	distances.
753	<
754	<	By monitoring C$\text{p}$ throughout these simulations, the melting
755	<	transition for SSD/E was observed at 230 K, about 5 degrees lower than
756	<	SSD. The resulting density maximum is located at 240 K, again about 5
757	<	degrees lower than the SSD value of 245 K. Though there is a decrease
758	<	in both of these values, the corrected densities near room temperature
759	<	justify the modifications taken.
740	>	Fig. \ref{ssdedense} shows the density profile for the {\sc ssd/e} model
741	>	in comparison to {\sc ssd1} without a reaction field, other common water
742	>	models, and experimental results. The calculated densities for both
743	>	{\sc ssd/e} and {\sc ssd1} have increased significantly over the original {\sc ssd}
744	>	model (see fig. \ref{dense1}) and are in better agreement with the
745	>	experimental values. At 298 K, the densities of {\sc ssd/e} and {\sc 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
748	>	the experimental value of 0.997 g/cm$^3$, and they are considerably
749	>	better than the {\sc ssd} value of 0.967$\pm$0.003 g/cm$^3$. The changes to
750	>	the dipole moment and sticky switching functions have improved the
751	>	structuring of the liquid (as seen in figure \ref{grcompare}, but they
752	>	have shifted the density maximum to much lower temperatures. This
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 {\sc ssd/e} model has
756	>	little effect on the melting transition. By monitoring $C_p$
757	>	throughout these simulations, the melting transition for {\sc ssd/e} was
758	>	shown to occur at 235 K.  The same transition temperature observed
759	>	with {\sc ssd} and {\sc ssd1}.
760
761	<	\begin{figure}
762	<	\includegraphics[width=85mm]{ssdrfcompare.epsi}
763	<	\caption{Comparison of densities calculated with SSD/RF to SSD with a
764	<	reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
765	<	SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
766	<	includes error bars, and the calculated results from the other
767	<	references were removed for clarity.}
761	>	\begin{figure}
762	>	\begin{center}
763	>	\epsfxsize=6in
764	>	\epsfbox{ssdrfDense.epsi}
765	>	\caption{Comparison of densities calculated with {\sc ssd/rf} to {\sc ssd1} with a
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 - {\sc ssd/rf} provides significantly more accurate densities
771	>	than {\sc ssd1} when performing room temperature simulations.}
772		\label{ssdrfdense}
773	+	\end{center}
774		\end{figure}
775
776	<	Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
777	<	SSD with an active reaction field. Like in the simulations of SSD/E,
778	<	the densities show a dramatic increase over normal SSD. At 298 K,
779	<	SSD/RF has a density of 0.997$\pm$0.001 g/cm$^3$, right in line with
780	<	experiment and considerably better than the SSD value of
781	<	0.941$\pm$0.001 g/cm$^3$. The melting point is observed at 240 K,
782	<	which is 5 degrees lower than SSD with a reaction field, and the
783	<	density maximum at 255 K, again 5 degrees lower than SSD. The density
784	<	at higher temperature still drops off more rapidly than the charge
785	<	based models but is in better agreement than SSD/E.
776	>	Including the reaction field long-range correction in the simulations
777	>	results in a more interesting comparison.  A density profile including
778	>	{\sc ssd/rf} and {\sc 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 {\sc ssd/rf} and {\sc ssd1} show a dramatic increase over
781	>	normal {\sc ssd} (see figure \ref{dense1}). At 298 K, {\sc 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 original {\sc ssd} value of 0.941$\pm$0.001
784	>	g/cm$^3$ and the {\sc 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 {\sc ssd/rf}, is identical to {\sc ssd} and only 5 K
789	>	lower than {\sc ssd1} with a reaction field. Additionally, the difference in
790	>	density maxima is not as extreme, with {\sc 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 {\sc ssd} and {\sc ssd1} respectively.
793
794	<	The reparameterization of the SSD water model, both for use with and
795	<	without an applied long-range correction, brought the densities up to
796	<	what is expected for simulating liquid water. In addition to improving
797	<	the densities, it is important that particle transport be maintained
798	<	or improved. Figure \ref{ssdediffuse} compares the temperature
799	<	dependence of the diffusion constant of SSD/E to SSD without an active
800	<	reaction field, both at the densities calculated at 1 atm and at the
801	<	experimentally calculated densities for super-cooled and liquid
802	<	water. In the upper plot, the diffusion constant for SSD/E is
803	<	consistently a little faster than experiment, while SSD starts off
804	<	slower than experiment and crosses to merge with SSD/E at high
805	<	temperatures. Both models follow the experimental trend well, but
739	<	diffuse too rapidly at higher temperatures. This abnormally fast
740	<	diffusion is caused by the decreased system density. Since the
741	<	densities of SSD/E don't deviate as much from experiment as those of
742	<	SSD, it follows the experimental trend more closely. This observation
743	<	is backed up by looking at the lower plot. The diffusion constants for
744	<	SSD/E track with the experimental values while SSD deviates on the low
745	<	side of the trend with increasing temperature. This is again a product
746	<	of SSD/E having densities closer to experiment, and not deviating to
747	<	lower densities with increasing temperature as rapidly.
748	<
749	<	\begin{figure}
750	<	\includegraphics[width=85mm]{ssdediffuse.epsi}
751	<	\caption{Plots of the diffusion constants calculated from SSD/E and SSD,
752	<	both without a reaction field along with experimental results from
753	<	Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
754	<	upper plot is at densities calculated from the NPT simulations at a
755	<	pressure of 1 atm, while the lower plot is at the experimentally
756	<	calculated densities.}
794	>	\begin{figure}
795	>	\begin{center}
796	>	\epsfxsize=6in
797	>	\epsfbox{ssdeDiffuse.epsi}
798	>	\caption{The diffusion constants calculated from {\sc ssd/e} and {\sc ssd1} (both
799	>	without a reaction field) along with experimental results
800	>	[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations were
801	>	performed at the average densities observed in the 1 atm NPT
802	>	simulations for the respective models. {\sc ssd/e} is slightly more mobile
803	>	than experiment at all of the temperatures, but it is closer to
804	>	experiment at biologically relevant temperatures than {\sc ssd1} without a
805	>	long-range correction.}
806		\label{ssdediffuse}
807	+	\end{center}
808		\end{figure}
809
810	<	\begin{figure}
811	<	\includegraphics[width=85mm]{ssdrfdiffuse.epsi}
812	<	\caption{Plots of the diffusion constants calculated from SSD/RF and SSD,
813	<	both with an active reaction field along with experimental results
814	<	from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
815	<	upper plot is at densities calculated from the NPT simulations at a
816	<	pressure of 1 atm, while the lower plot is at the experimentally
817	<	calculated densities.}
810	>	The reparameterization of the {\sc 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 the diffusive behavior of {\sc ssd} be
814	>	maintained or improved. Figure \ref{ssdediffuse} compares the
815	>	temperature dependence of the diffusion constant of {\sc ssd/e} to {\sc ssd1}
816	>	without an active reaction field at the densities calculated from
817	>	their respective NPT simulations at 1 atm. The diffusion constant for
818	>	{\sc ssd/e} is consistently higher than experiment, while {\sc ssd1} remains lower
819	>	than experiment until relatively high temperatures (around 360
820	>	K). Both models follow the shape of the experimental curve well below
821	>	300 K but tend to diffuse too rapidly at higher temperatures, as seen
822	>	in {\sc ssd1}'s crossing above 360 K.  This increasing diffusion relative to
823	>	the experimental values is caused by the rapidly decreasing system
824	>	density with increasing temperature.  Both {\sc ssd1} and {\sc ssd/e} show this
825	>	deviation in particle mobility, but this trend has different
826	>	implications on the diffusive behavior of the models.  While {\sc ssd1}
827	>	shows more experimentally accurate diffusive behavior in the high
828	>	temperature regimes, {\sc ssd/e} shows more accurate behavior in the
829	>	supercooled and biologically relevant temperature ranges.  Thus, the
830	>	changes made to improve the liquid structure may have had an adverse
831	>	affect on the density maximum, but they improve the transport behavior
832	>	of {\sc ssd/e} relative to {\sc ssd1} under the most commonly simulated
833	>	conditions.
834	>
835	>	\begin{figure}
836	>	\begin{center}
837	>	\epsfxsize=6in
838	>	\epsfbox{ssdrfDiffuse.epsi}
839	>	\caption{The diffusion constants calculated from {\sc ssd/rf} and {\sc ssd1} (both
840	>	with an active reaction field) along with experimental results
841	>	[Refs. \citen{Gillen72} and \citen{Holz00}]. The NVE calculations were
842	>	performed at the average densities observed in the 1 atm NPT
843	>	simulations for both of the models. {\sc ssd/rf} simulates the diffusion of
844	>	water throughout this temperature range very well. The rapidly
845	>	increasing diffusion constants at high temperatures for both models
846	>	can be attributed to lower calculated densities than those observed in
847	>	experiment.}
848		\label{ssdrfdiffuse}
849	+	\end{center}
850		\end{figure}
851
852	<	In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
853	<	compared with SSD with an active reaction field. In the upper plot,
854	<	SSD/RF tracks with the experimental results incredibly well, identical
855	<	within error throughout the temperature range and only showing a
856	<	slight increasing trend at higher temperatures. SSD also tracks
857	<	experiment well, only it tends to diffuse a little more slowly at low
858	<	temperatures and deviates to diffuse too rapidly at high
859	<	temperatures. As was stated in the SSD/E comparisons, this deviation
860	<	away from the ideal trend is due to a rapid decrease in density at
861	<	higher temperatures. SSD/RF doesn't suffer from this problem as much
862	<	as SSD, because the calculated densities are more true to
863	<	experiment. This is again emphasized in the lower plot, where SSD/RF
783	<	tracks the experimental diffusion exactly while SSD's diffusion
784	<	constants are slightly too low due to its need for a lower density at
785	<	the specified temperature.
852	>	In figure \ref{ssdrfdiffuse}, the diffusion constants for {\sc ssd/rf} are
853	>	compared to {\sc ssd1} with an active reaction field. Note that {\sc ssd/rf}
854	>	tracks the experimental results quantitatively, identical within error
855	>	throughout most of the temperature range shown and exhibiting only a
856	>	slight increasing trend at higher temperatures. {\sc ssd1} tends to diffuse
857	>	more slowly at low temperatures and deviates to diffuse too rapidly at
858	>	temperatures greater than 330 K.  As stated above, this deviation away
859	>	from the ideal trend is due to a rapid decrease in density at higher
860	>	temperatures. {\sc ssd/rf} does not suffer from this problem as much as {\sc ssd1}
861	>	because the calculated densities are closer to the experimental
862	>	values. These results again emphasize the importance of careful
863	>	reparameterization when using an altered long-range correction.
864
865	<	\subsection{Additional Observations}
865	>	\begin{table}
866	>	\begin{minipage}{\linewidth}
867	>	\renewcommand{\thefootnote}{\thempfootnote}
868	>	\begin{center}
869	>	\caption{Properties of the single-point water models compared with
870	>	experimental data at ambient conditions}
871	>	\begin{tabular}{ l  c  c  c  c  c }
872	>	\hline \\[-3mm]
873	>	\ \ \ \ \ \  & \ \ \ {\sc ssd1} \ \ \ & \ {\sc ssd/e} \ \ \ & \ {\sc ssd1} (RF) \ \
874	>	\ & \ {\sc ssd/rf} \ \ \ & \ Expt. \\
875	>	\hline \\[-3mm]
876	>	\ \ \ $\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 \\
877	>	\ \ \ $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 \\
878	>	\ \ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78 $\pm$0.07 & 2.51 $\pm$0.18 &
879	>	2.00 $\pm$0.17 & 2.32 $\pm$0.06 & 2.299\cite{Mills73} \\
880	>	\ \ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 &
881	>	4.7\footnote{Calculated by integrating $g_{\text{OO}}(r)$ in
882	>	Ref. \citen{Head-Gordon00_1}} \\
883	>	\ \ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 &
884	>	3.5\footnote{Calculated by integrating $g_{\text{OH}}(r)$ in
885	>	Ref. \citen{Soper86}}  \\
886	>	\ \ \ $\tau_1$ (ps) & 10.9 $\pm$0.6 & 7.3 $\pm$0.4 & 7.5 $\pm$0.7 &
887	>	7.2 $\pm$0.4 & 5.7\footnote{Calculated for 298 K from data in Ref. \citen{Eisenberg69}} \\
888	>	\ \ \ $\tau_2$ (ps) & 4.7 $\pm$0.4 & 3.1 $\pm$0.2 & 3.5 $\pm$0.3 & 3.2
889	>	$\pm$0.2 & 2.3\footnote{Calculated for 298 K from data in
890	>	Ref. \citen{Krynicki66}}
891	>	\end{tabular}
892	>	\label{liquidproperties}
893	>	\end{center}
894	>	\end{minipage}
895	>	\end{table}
896
897	<	While performing the melting sequences of SSD/E, some interesting
898	<	observations were made. After melting at 230 K, two of the systems
899	<	underwent crystallization events near 245 K. As the heating process
900	<	continued, the two systems remained crystalline until finally melting
901	<	between 320 and 330 K. These simulations were excluded from the data
902	<	set shown in figure \ref{ssdedense} and replaced with two additional
903	<	melting sequences that did not undergo this anomalous phase
904	<	transition, while this crystallization event was investigated
905	<	separately.
897	>	Table \ref{liquidproperties} gives a synopsis of the liquid state
898	>	properties of the water models compared in this study along with the
899	>	experimental values for liquid water at ambient conditions. The
900	>	coordination number ($n_C$) and number of hydrogen bonds per particle
901	>	($n_H$) were calculated by integrating the following relations:
902	>	\begin{equation}
903	>	n_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr,
904	>	\end{equation}
905	>	\begin{equation}
906	>	n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr,
907	>	\end{equation}
908	>	where $\rho$ is the number density of the specified pair interactions,
909	>	$a$ and $b$ are the radial locations of the minima following the first
910	>	peak in g$_\text{OO}(r)$ or g$_\text{OH}(r)$ respectively. The number
911	>	of hydrogen bonds stays relatively constant across all of the models,
912	>	but the coordination numbers of {\sc ssd/e} and {\sc ssd/rf} show an improvement
913	>	over {\sc ssd1}.  This improvement is primarily due to extension of the
914	>	first solvation shell in the new parameter sets.  Because $n_H$ and
915	>	$n_C$ are nearly identical in {\sc ssd1}, it appears that all molecules in
916	>	the first solvation shell are involved in hydrogen bonds.  Since $n_H$
917	>	and $n_C$ differ in the newly parameterized models, the orientations
918	>	in the first solvation shell are a bit more ``fluid''.  Therefore {\sc ssd1}
919	>	overstructures the first solvation shell and our reparameterizations
920	>	have returned this shell to more realistic liquid-like behavior.
921
922	<	\begin{figure}
923	<	\includegraphics[width=85mm]{povIce.ps}
924	<	\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
925	<	\label{weirdice}
926	<	\end{figure}
922	>	The time constants for the orientational autocorrelation functions
923	>	are also displayed in Table \ref{liquidproperties}. The dipolar
924	>	orientational time correlation functions ($C_{l}$) are described
925	>	by:
926	>	\begin{equation}
927	>	C_{l}(t) = \langle P_l[\hat{\mathbf{u}}_j(0)\cdot\hat{\mathbf{u}}_j(t)]\rangle,
928	>	\end{equation}
929	>	where $P_l$ are Legendre polynomials of order $l$ and
930	>	$\hat{\mathbf{u}}_j$ is the unit vector pointing along the molecular
931	>	dipole.\cite{Rahman71} From these correlation functions, the
932	>	orientational relaxation time of the dipole vector can be calculated
933	>	from an exponential fit in the long-time regime ($t >
934	>	\tau_l$).\cite{Rothschild84} Calculation of these time constants were
935	>	averaged over five detailed NVE simulations performed at the ambient
936	>	conditions for each of the respective models. It should be noted that
937	>	the commonly cited value of 1.9 ps for $\tau_2$ was determined from
938	>	the NMR data in Ref. \citen{Krynicki66} at a temperature near
939	>	34$^\circ$C.\cite{Rahman71} Because of the strong temperature
940	>	dependence of $\tau_2$, it is necessary to recalculate it at 298 K to
941	>	make proper comparisons. The value shown in Table
942	>	\ref{liquidproperties} was calculated from the same NMR data in the
943	>	fashion described in Ref. \citen{Krynicki66}. Similarly, $\tau_1$ was
944	>	recomputed for 298 K from the data in Ref. \citen{Eisenberg69}.
945	>	Again, {\sc ssd/e} and {\sc ssd/rf} show improved behavior over {\sc ssd1}, both with
946	>	and without an active reaction field. Turning on the reaction field
947	>	leads to much improved time constants for {\sc ssd1}; however, these results
948	>	also include a corresponding decrease in system density.
949	>	Orientational relaxation times published in the original {\sc ssd} dynamics
950	>	papers are smaller than the values observed here, and this difference
951	>	can be attributed to the use of the Ewald sum.\cite{Ichiye99}
952
953	<	The final configurations of these two melting sequences shows an
806	<	expanded zeolite-like crystal structure that does not correspond to
807	<	any known form of ice. For convenience and to help distinguish it from
808	<	the experimentally observed forms of ice, this crystal structure will
809	<	henceforth be referred to as ice-zero (ice 0). The crystallinity was
810	<	extensive enough than a near ideal crystal structure could be
811	<	obtained. Figure \ref{weirdice} shows the repeating crystal structure
812	<	of a typical crystal at 5 K. The unit cell contains eight molecules,
813	<	and figure \ref{unitcell} shows a unit cell built from the water
814	<	particle center of masses that can be used to construct a repeating
815	<	lattice, similar to figure \ref{weirdice}. Each molecule is hydrogen
816	<	bonded to four other water molecules; however, the hydrogen bonds are
817	<	flexed rather than perfectly straight. This results in a skewed
818	<	tetrahedral geometry about the central molecule. Looking back at
819	<	figure \ref{isosurface}, it is easy to see how these flexed hydrogen
820	<	bonds are allowed in that the attractive regions are conical in shape,
821	<	with the greatest attraction in the central region. Though not ideal,
822	<	these flexed hydrogen bonds are favorable enough to stabilize an
823	<	entire crystal generated around them. In fact, the imperfect ice 0
824	<	crystals were so stable that they melted at greater than room
825	<	temperature.
953	>	\subsection{Additional Observations}
954
955		\begin{figure}
956	<	\includegraphics[width=65mm]{ice0cell.eps}
957	<	\caption{Simple unit cell for constructing ice 0. In this cell, $c$ is
958	<	equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
959	<	\label{unitcell}
956	>	\begin{center}
957	>	\epsfxsize=6in
958	>	\epsfbox{icei_bw.eps}
959	>	\caption{The most stable crystal structure assumed by the {\sc ssd} family
960	>	of water models.  We refer to this structure as Ice-{\it i} to
961	>	indicate its origins in computer simulation.  This image was taken of
962	>	the (001) face of the crystal.}
963	>	\label{weirdice}
964	>	\end{center}
965		\end{figure}
966
967	<	The initial simulations indicated that ice 0 is the preferred ice
968	<	structure for at least SSD/E. To verify this, a comparison was made
969	<	between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at
970	<	constant pressure with SSD/E, SSD/RF, and SSD. Near ideal versions of
971	<	the three types of crystals were cooled to ~1 K, and the potential
972	<	energies of each were compared using all three water models. With
973	<	every water model, ice 0 turned out to have the lowest potential
974	<	energy: 5\% lower than $I_h$ with SSD, 6.5\% lower with SSD/E, and
975	<	7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
976	<	was observed to be ~2\% less stable than ice $I_h$. In addition to
977	<	having the lowest potential energy, ice 0 was the most expanded of the
978	<	three ice crystals, ~5\% less dense than ice $I_h$ with all of the
979	<	water models. In all three water models, ice $I_c$ was observed to be
980	<	~2\% more dense than ice $I_h$.
967	>	While performing a series of melting simulations on an early iteration
968	>	of {\sc ssd/e} not discussed in this paper, we observed recrystallization
969	>	into a novel structure not previously known for water.  After melting
970	>	at 235 K, two of five systems underwent crystallization events near
971	>	245 K.  The two systems remained crystalline up to 320 and 330 K,
972	>	respectively.  The crystal exhibits an expanded zeolite-like structure
973	>	that does not correspond to any known form of ice.  This appears to be
974	>	an artifact of the point dipolar models, so to distinguish it from the
975	>	experimentally observed forms of ice, we have denoted the structure
976	>	Ice-$\sqrt{\smash[b]{-\text{I}}}$ (Ice-{\it i}).  A large enough
977	>	portion of the sample crystallized that we have been able to obtain a
978	>	near ideal crystal structure of Ice-{\it i}. Figure \ref{weirdice}
979	>	shows the repeating crystal structure of a typical crystal at 5
980	>	K. Each water molecule is hydrogen bonded to four others; however, the
981	>	hydrogen bonds are bent rather than perfectly straight. This results
982	>	in a skewed tetrahedral geometry about the central molecule.  In
983	>	figure \ref{isosurface}, it is apparent that these flexed hydrogen
984	>	bonds are allowed due to the conical shape of the attractive regions,
985	>	with the greatest attraction along the direct hydrogen bond
986	>	configuration. Though not ideal, these flexed hydrogen bonds are
987	>	favorable enough to stabilize an entire crystal generated around them.
988
989	<	In addition to the low temperature comparisons, melting sequences were
990	<	performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
991	<	and SSD both with and without a reaction field. The melting
992	<	transitions for both SSD/E and SSD without a reaction field occurred
993	<	at temperature in excess of 375 K. SSD/RF and SSD with a reaction
994	<	field had more reasonable melting transitions, down near 325 K. These
995	<	melting point observations emphasize how preferred this crystal
996	<	structure is over the most common types of ice when using these single
997	<	point water models.
989	>	Initial simulations indicated that Ice-{\it i} is the preferred ice
990	>	structure for at least the {\sc ssd/e} model. To verify this, a
991	>	comparison was made between near ideal crystals of ice~$I_h$,
992	>	ice~$I_c$, and Ice-{\it i} at constant pressure with {\sc ssd/e}, {\sc
993	>	ssd/rf}, and {\sc ssd1}. Near-ideal versions of the three types of
994	>	crystals were cooled to 1 K, and enthalpies of formation of each were
995	>	compared using all three water models.  Enthalpies were estimated from
996	>	the isobaric-isothermal simulations using $H=U+P_{\text ext}V$ where
997	>	$P_{\text ext}$ is the applied pressure.  A constant value of
998	>	-60.158 kcal / mol has been added to place our zero for the
999	>	enthalpies of formation for these systems at the traditional state
1000	>	(elemental forms at standard temperature and pressure).  With every
1001	>	model in the {\sc ssd} family, Ice-{\it i} had the lowest calculated
1002	>	enthalpy of formation.
1003
1004	<	Recognizing that the above tests show ice 0 to be both the most stable
1005	<	and lowest density crystal structure for these single point water
1006	<	models, it is interesting to speculate on the favorability of this
1007	<	crystal structure with the different charge based models. As a quick
1008	<	test, these 3 crystal types were converted from SSD type particles to
1009	<	TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
1010	<	minimizations were performed on all of these crystals to compare the
1011	<	system energies. Again, ice 0 was observed to have the lowest total
1012	<	system energy. The total energy of ice 0 was ~2\% lower than ice
1013	<	$I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial
1014	<	results, we would not be surprised if results from the other common
1015	<	water models show ice 0 to be the lowest energy crystal structure. A
1016	<	continuation on work studing ice 0 with multipoint water models will
1017	<	be published in a coming article.
1004	>	\begin{table}
1005	>	\begin{center}
1006	>	\caption{Enthalpies of Formation (in kcal / mol) of the three crystal
1007	>	structures (at 1 K) exhibited by the {\sc ssd} family of water models}
1008	>	\begin{tabular}{ l  c  c  c  }
1009	>	\hline \\[-3mm]
1010	>	\ \ \ Water Model \ \ \  & \ \ \ Ice-$I_h$ \ \ \ & \ Ice-$I_c$\ \  & \
1011	>	Ice-{\it i} \\
1012	>	\hline \\[-3mm]
1013	>	\ \ \ {\sc ssd/e} & -12.286 & -12.292 & -13.590 \\
1014	>	\ \ \ {\sc ssd/rf} & -12.935 & -12.917 & -14.022 \\
1015	>	\ \ \ {\sc ssd1} & -12.496 & -12.411 & -13.417 \\
1016	>	\ \ \ {\sc ssd1} (RF) & -12.504 & -12.411 & -13.134 \\
1017	>	\end{tabular}
1018	>	\label{iceenthalpy}
1019	>	\end{center}
1020	>	\end{table}
1021
1022	+	In addition to these energetic comparisons, melting simulations were
1023	+	performed with ice-{\it i} as the initial configuration using {\sc ssd/e},
1024	+	{\sc ssd/rf}, and {\sc ssd1} both with and without a reaction field. The melting
1025	+	transitions for both {\sc ssd/e} and {\sc ssd1} without reaction field occurred at
1026	+	temperature in excess of 375~K.  {\sc ssd/rf} and {\sc ssd1} with a reaction field
1027	+	showed more reasonable melting transitions near 325~K.  These melting
1028	+	point observations clearly show that all of the {\sc ssd}-derived models
1029	+	prefer the ice-{\it i} structure.
1030	+
1031		\section{Conclusions}
1032	<	The density maximum and temperature dependent transport for the SSD
1033	<	water model, both with and without the use of reaction field, were
1034	<	studied via a series of NPT and NVE simulations. The constant pressure
1035	<	simulations of the melting of both $I_h$ and $I_c$ ice showed a
1036	<	density maximum near 260 K. In most cases, the calculated densities
1037	<	were significantly lower than the densities calculated in simulations
1038	<	of other water models. Analysis of particle diffusion showed SSD to
1039	<	capture the transport properties of experimental very well in both the
1040	<	normal and super-cooled liquid regimes. In order to correct the
1041	<	density behavior, SSD was reparameterized for use both with and
1042	<	without a long-range interaction correction, SSD/RF and SSD/E
1043	<	respectively. In addition to correcting the abnormally low densities,
1044	<	these new versions were show to maintain or improve upon the transport
1045	<	and structural features of the original water model, all while
1046	<	maintaining the fast performance of the SSD water model. This work
1047	<	shows these simple water models, and in particular SSD/E and SSD/RF,
1048	<	to be excellent choices to represent explicit water in future
1032	>
1033	>	The density maximum and temperature dependence of the self-diffusion
1034	>	constant were studied for the {\sc ssd} water model, both with and without
1035	>	the use of reaction field, via a series of NPT and NVE
1036	>	simulations. The constant pressure simulations showed a density
1037	>	maximum near 260 K. In most cases, the calculated densities were
1038	>	significantly lower than the densities obtained from other water
1039	>	models (and experiment). Analysis of self-diffusion showed {\sc ssd} to
1040	>	capture the transport properties of water well in both the liquid and
1041	>	supercooled liquid regimes.
1042	>
1043	>	In order to correct the density behavior, the original {\sc ssd} model was
1044	>	reparameterized for use both with and without a reaction field ({\sc ssd/rf}
1045	>	and {\sc ssd/e}), and comparisons were made with {\sc ssd1}, Ichiye's density
1046	>	corrected version of {\sc ssd}. Both models improve the liquid structure,
1047	>	densities, and diffusive properties under their respective simulation
1048	>	conditions, indicating the necessity of reparameterization when
1049	>	changing the method of calculating long-range electrostatic
1050	>	interactions.  In general, however, these simple water models are
1051	>	excellent choices for representing explicit water in large scale
1052		simulations of biochemical systems.
1053
1054	+	The existence of a novel low-density ice structure that is preferred
1055	+	by the {\sc ssd} family of water models is somewhat troubling, since liquid
1056	+	simulations on this family of water models at room temperature are
1057	+	effectively simulations of supercooled or metastable liquids.  One
1058	+	way to destabilize this unphysical ice structure would be to make the
1059	+	range of angles preferred by the attractive part of the sticky
1060	+	potential much narrower.  This would require extensive
1061	+	reparameterization to maintain the same level of agreement with the
1062	+	experiments.
1063	+
1064	+	Additionally, our initial calculations show that the Ice-{\it i}
1065	+	structure may also be a preferred crystal structure for at least one
1066	+	other popular multi-point water model (TIP3P), and that much of the
1067	+	simulation work being done using this popular model could also be at
1068	+	risk for crystallization into this unphysical structure.  A future
1069	+	publication will detail the relative stability of the known ice
1070	+	structures for a wide range of popular water models.
1071	+
1072		\section{Acknowledgments}
1073	<	The authors would like to thank the National Science Foundation for
1074	<	funding under grant CHE-0134881. Computation time was provided by the
1075	<	Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR
1076	<	00 79647.
1073	>	Support for this project was provided by the National Science
1074	>	Foundation under grant CHE-0134881. Computation time was provided by
1075	>	the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
1076	>	DMR-0079647.
1077
1078	<	\bibliographystyle{jcp}
1078	>	\newpage
1079
1080	+	\bibliographystyle{jcp}
1081		\bibliography{nptSSD}
1082
1083		%\pagebreak

Diff Legend

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