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 1033 by chrisfen, Fri Feb 6 19:15:44 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[preprint,aps]{revtex4}
3	>	%\documentclass[11pt]{article}
4	>	%\usepackage{endfloat}
5		\usepackage{amsmath}
6	+	\usepackage{epsf}
7	+	\usepackage{berkeley}
8	+	\usepackage{setspace}
9	+	\usepackage{tabularx}
10		\usepackage{graphicx}
5	–
6	–	%\usepackage{endfloat}
7	–	%\usepackage{berkeley}
8	–	%\usepackage{epsf}
11		%\usepackage[ref]{overcite}
10	–	%\usepackage{setspace}
11	–	%\usepackage{tabularx}
12	–	%\usepackage{graphicx}
12		%\usepackage{curves}
14	–	%\usepackage{amsmath}
13		%\pagestyle{plain}
14		%\pagenumbering{arabic}
15		%\oddsidemargin 0.0cm \evensidemargin 0.0cm
#	Line 20 | Line 18
18		%\brokenpenalty=10000
19		%\renewcommand{\baselinestretch}{1.2}
20		%\renewcommand\citemid{\ } % no comma in optional reference note
21	+	\newcounter{captions}
22	+	\newcounter{figs}
23
24		\begin{document}
25
26	<	\title{On the temperature dependent structural and transport properties of the soft sticky dipole (SSD) and related single point water models}
26	>	\title{On the structural and transport properties of the soft sticky
27	>	dipole (SSD) and related single point water models}
28
29	<	\author{Christopher J. Fennell and J. Daniel Gezelter{\thefootnote}
29	<	\footnote[1]{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
29	>	\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu}}
30
31	<	\address{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
31	>	\affiliation{Department of Chemistry and Biochemistry\\ University of Notre Dame\\
32		Notre Dame, Indiana 46556}
33
34		\date{\today}
35
36	+
37		\begin{abstract}
38	<	NVE and NPT molecular dynamics simulations were performed in order to
39	<	investigate the density maximum and temperature dependent transport
40	<	for the SSD water model, both with and without the use of reaction
41	<	field. The constant pressure simulations of the melting of both $I_h$
42	<	and $I_c$ ice showed a density maximum near 260 K. In most cases, the
43	<	calculated densities were significantly lower than the densities
44	<	calculated in simulations of other water models. Analysis of particle
45	<	diffusion showed SSD to capture the transport properties of
46	<	experimental very well in both the normal and super-cooled liquid
47	<	regimes. In order to correct the density behavior, SSD was
48	<	reparameterized for use both with and without a long-range interaction
49	<	correction, SSD/RF and SSD/E respectively. In addition to correcting
50	<	the abnormally low densities, these new versions were show to maintain
51	<	or improve upon the transport and structural features of the original
52	<	water model.
38	>	The density maximum and temperature dependence of the self-diffusion
39	>	constant were investigated for the soft sticky dipole (SSD) water
40	>	model and two related reparameterizations of this single-point model.
41	>	A combination of microcanonical and isobaric-isothermal molecular
42	>	dynamics simulations were used to calculate these properties, both
43	>	with and without the use of reaction field to handle long-range
44	>	electrostatics.  The isobaric-isothermal (NPT) simulations of the
45	>	melting of both ice-$I_h$ and ice-$I_c$ showed a density maximum near
46	>	260~K.  In most cases, the use of the reaction field resulted in
47	>	calculated densities which were were significantly lower than
48	>	experimental densities.  Analysis of self-diffusion constants shows
49	>	that the original SSD model captures the transport properties of
50	>	experimental water very well in both the normal and super-cooled
51	>	liquid regimes.  We also present our reparameterized versions of SSD
52	>	for use both with the reaction field or without any long-range
53	>	electrostatic corrections.  These are called the SSD/RF and SSD/E
54	>	models respectively.  These modified models were shown to maintain or
55	>	improve upon the experimental agreement with the structural and
56	>	transport properties that can be obtained with either the original SSD
57	>	or the density corrected version of the original model (SSD1).
58	>	Additionally, a novel low-density ice structure is presented
59	>	which appears to be the most stable ice structure for the entire SSD
60	>	family.
61		\end{abstract}
62
63		\maketitle
64
65	+	\newpage
66	+
67		%\narrowtext
68
69
#	Line 62 | Line 73 | One of the most important tasks in simulations of bioc
73
74		\section{Introduction}
75
76	<	One of the most important tasks in simulations of biochemical systems
77	<	is the proper depiction of water and water solvation. In fact, the
78	<	bulk of the calculations performed in solvated simulations are of
79	<	interactions with or between solvent molecules. Thus, the outcomes of
80	<	these types of simulations are highly dependent on the physical
81	<	properties of water, both as individual molecules and in
82	<	groups/bulk. Due to the fact that explicit solvent accounts for a
83	<	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}
76	>	One of the most important tasks in the simulation of biochemical
77	>	systems is the proper depiction of the aqueous environment of the
78	>	molecules of interest.  In some cases (such as in the simulation of
79	>	phospholipid bilayers), the majority of the calculations that are
80	>	performed involve interactions with or between solvent molecules.
81	>	Thus, the properties one may observe in biochemical simulations are
82	>	going to be highly dependent on the physical properties of the water
83	>	model that is chosen.
84
85	<	The Soft Sticky Dipole (SSD)\ water model was developed by Ichiye
86	<	\emph{et al.} as a modified form of the hard-sphere water model
87	<	proposed by Bratko, Blum, and Luzar.\cite{Bratko85,Bratko95} SSD
88	<	consists of a single point dipole with a Lennard-Jones core and a
89	<	sticky potential that directs the particles to assume the proper
90	<	hydrogen bond orientation in the first solvation shell. Thus, the
91	<	interaction between two SSD water molecules \emph{i} and \emph{j} is
92	<	given by the potential
85	>	There is an especially delicate balance between computational
86	>	efficiency and the ability of the water model to accurately predict
87	>	the properties of bulk
88	>	water.\cite{Jorgensen83,Berendsen87,Jorgensen00} For example, the
89	>	TIP5P model improves on the structural and transport properties of
90	>	water relative to the previous TIP models, yet this comes at a greater
91	>	than 50\% increase in computational
92	>	cost.\cite{Jorgensen01,Jorgensen00}
93	>
94	>	One recently developed model that largely succeeds in retaining the
95	>	accuracy of bulk properties while greatly reducing the computational
96	>	cost is the Soft Sticky Dipole (SSD) water
97	>	model.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03} The SSD model
98	>	was developed by Ichiye \emph{et al.} as a modified form of the
99	>	hard-sphere water model proposed by Bratko, Blum, and
100	>	Luzar.\cite{Bratko85,Bratko95} SSD is a {\it single point} model
101	>	which has an interaction site that is both a point dipole and a
102	>	Lennard-Jones core.  However, since the normal aligned and
103	>	anti-aligned geometries favored by point dipoles are poor mimics of
104	>	local structure in liquid water, a short ranged ``sticky'' potential
105	>	is also added.  The sticky potential directs the molecules to assume
106	>	the proper hydrogen bond orientation in the first solvation shell.
107	>
108	>	The interaction between two SSD water molecules \emph{i} and \emph{j}
109	>	is given by the potential
110		\begin{equation}
111		u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
112	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
112	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)\ +
113		u_{ij}^{sp}
114	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
114	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j),
115		\end{equation}
116	<	where the $\mathbf{r}_{ij}$ is the position vector between molecules
117	<	\emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and
118	<	$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
119	<	orientations of the respective molecules. The Lennard-Jones, dipole,
120	<	and sticky parts of the potential are giving by the following
108	<	equations,
116	>	where the ${\bf r}_{ij}$ is the position vector between molecules
117	>	\emph{i} and \emph{j} with magnitude $r_{ij}$, and
118	>	${\bf \Omega}_i$ and ${\bf \Omega}_j$ represent the orientations of
119	>	the two molecules. The Lennard-Jones and dipole interactions are given
120	>	by the following familiar forms:
121		\begin{equation}
122	<	u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
122	>	u_{ij}^{LJ}(r_{ij}) = 4\epsilon
123	>	\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]
124	>	\ ,
125		\end{equation}
126	+	and
127		\begin{equation}
128	<	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}\ ,
128	>	u_{ij}^{dp} = \frac{|\mu_i||\mu_j|}{4 \pi \epsilon_0 r_{ij}^3} \left(
129	>	\hat{\bf u}_i \cdot \hat{\bf u}_j - 3(\hat{\bf u}_i\cdot\hat{\bf
130	>	r}_{ij})(\hat{\bf u}_j\cdot\hat{\bf r}_{ij}) \right)\ ,
131		\end{equation}
132	+	where $\hat{\bf u}_i$ and $\hat{\bf u}_j$ are the unit vectors along
133	+	the dipoles of molecules $i$ and $j$ respectively. $|\mu_i|$ and
134	+	$|\mu_j|$ are the strengths of the dipole moments, and $\hat{\bf
135	+	r}_{ij}$ is the unit vector pointing from molecule $j$ to molecule
136	+	$i$.
137	+
138	+	The sticky potential is somewhat less familiar:
139		\begin{equation}
116	–	\begin{split}
140		u_{ij}^{sp}
141	<	(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)
142	<	&=
143	<	\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
144	<	& \quad \ +
145	<	s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
123	<	\end{split}
141	>	({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j) =
142	>	\frac{\nu_0}{2}[s(r_{ij})w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)
143	>	+ s^\prime(r_{ij})w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf
144	>	\Omega}_j)]\ .
145	>	\label{stickyfunction}
146		\end{equation}
147	<	where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
148	<	unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
149	<	$\nu_0$ scales the strength of the overall sticky potential, $s$ and
150	<	$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
151	<	functions take the following forms,
147	>	Here, $\nu_0$ is a strength parameter for the sticky potential, and
148	>	$s$ and $s^\prime$ are cubic switching functions which turn off the
149	>	sticky interaction beyond the first solvation shell. The $w$ function
150	>	can be thought of as an attractive potential with tetrahedral
151	>	geometry:
152		\begin{equation}
153	<	w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
153	>	w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
154		\end{equation}
155	+	while the $w^\prime$ function counters the normal aligned and
156	+	anti-aligned structures favored by point dipoles:
157		\begin{equation}
158	<	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,
158	>	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,
159		\end{equation}
160	<	where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
161	<	term that promotes hydrogen bonding orientations within the first
162	<	solvation shell, and $w^\prime$ is a dipolar repulsion term that
163	<	repels unrealistic dipolar arrangements within the first solvation
164	<	shell. A more detailed description of the functional parts and
165	<	variables in this potential can be found in other
166	<	articles.\cite{Ichiye96,Ichiye99}
160	>	It should be noted that $w$ is proportional to the sum of the $Y_3^2$
161	>	and $Y_3^{-2}$ spherical harmonics (a linear combination which
162	>	enhances the tetrahedral geometry for hydrogen bonded structures),
163	>	while $w^\prime$ is a purely empirical function.  A more detailed
164	>	description of the functional parts and variables in this potential
165	>	can be found in the original SSD
166	>	articles.\cite{Ichiye96,Ichiye96b,Ichiye99,Ichiye03}
167
168	<	Being that this is a one-site point dipole model, the actual force
169	<	calculations are simplified significantly. In the original Monte Carlo
170	<	simulations using this model, Ichiye \emph{et al.} reported a
171	<	calculation speed up of up to an order of magnitude over other
172	<	comparable models while maintaining the structural behavior of
173	<	water.\cite{Ichiye96} In the original molecular dynamics studies of
174	<	SSD, it was shown that it actually improves upon the prediction of
175	<	water's dynamical properties 3 and 4-point models.\cite{Ichiye99} This
176	<	attractive combination of speed and accurate depiction of solvent
177	<	properties makes SSD a model of interest for the simulation of large
178	<	scale biological systems, such as membrane phase behavior, a specific
179	<	interest within our group.
168	>	Since SSD is a single-point {\it dipolar} model, the force
169	>	calculations are simplified significantly relative to the standard
170	>	{\it charged} multi-point models. In the original Monte Carlo
171	>	simulations using this model, Liu and Ichiye reported that using SSD
172	>	decreased computer time by a factor of 6-7 compared to other
173	>	models.\cite{Ichiye96} What is most impressive is that this savings
174	>	did not come at the expense of accurate depiction of the liquid state
175	>	properties.  Indeed, SSD maintains reasonable agreement with the Soper
176	>	data for the structural features of liquid
177	>	water.\cite{Soper86,Ichiye96} Additionally, the dynamical properties
178	>	exhibited by SSD agree with experiment better than those of more
179	>	computationally expensive models (like TIP3P and
180	>	SPC/E).\cite{Ichiye99} The combination of speed and accurate depiction
181	>	of solvent properties makes SSD a very attractive model for the
182	>	simulation of large scale biochemical simulations.
183
184	<	Up to this point, a detailed look at the model's structure and ion
185	<	solvation abilities has been performed.\cite{Ichiye96} In addition, a
186	<	thorough investigation of the dynamic properties of SSD was performed
187	<	by Chandra and Ichiye focusing on translational and orientational
188	<	properties at 298 K.\cite{Ichiye99} This study focuses on determining
189	<	the density maximum for SSD utilizing both microcanonical and
190	<	isobaric-isothermal ensemble molecular dynamics, while using the
191	<	reaction field method for handling long-ranged dipolar interactions. A
192	<	reaction field method has been previously implemented in Monte Carlo
193	<	simulations by Liu and Ichiye in order to study the static dielectric
194	<	constant for the model.\cite{Ichiye96b} This paper will expand the
195	<	scope of these original simulations to look on how the reaction field
196	<	affects the physical and dynamic properties of SSD systems.
184	>	One feature of the SSD model is that it was parameterized for
185	>	use with the Ewald sum to handle long-range interactions.  This would
186	>	normally be the best way of handling long-range interactions in
187	>	systems that contain other point charges.  However, our group has
188	>	recently become interested in systems with point dipoles as mimics for
189	>	neutral, but polarized regions on molecules (e.g. the zwitterionic
190	>	head group regions of phospholipids).  If the system of interest does
191	>	not contain point charges, the Ewald sum and even particle-mesh Ewald
192	>	become computational bottlenecks.  Their respective ideal
193	>	$N^\frac{3}{2}$ and $N\log N$ calculation scaling orders for $N$
194	>	particles can become prohibitive when $N$ becomes
195	>	large.\cite{Darden99} In applying this water model in these types of
196	>	systems, it would be useful to know its properties and behavior under
197	>	the more computationally efficient reaction field (RF) technique, or
198	>	even with a simple cutoff. This study addresses these issues by
199	>	looking at the structural and transport behavior of SSD over a
200	>	variety of temperatures with the purpose of utilizing the RF
201	>	correction technique.  We then suggest modifications to the parameters
202	>	that result in more realistic bulk phase behavior.  It should be noted
203	>	that in a recent publication, some of the original investigators of
204	>	the SSD water model have suggested adjustments to the SSD
205	>	water model to address abnormal density behavior (also observed here),
206	>	calling the corrected model SSD1.\cite{Ichiye03} In what
207	>	follows, we compare our reparamaterization of SSD with both the
208	>	original SSD and SSD1 models with the goal of improving
209	>	the bulk phase behavior of an SSD-derived model in simulations
210	>	utilizing the reaction field.
211
212		\section{Methods}
213
214	<	As stated previously, in this study the long-range dipole-dipole
215	<	interactions were accounted for using the reaction field method. The
216	<	magnitude of the reaction field acting on dipole \emph{i} is given by
214	>	Long-range dipole-dipole interactions were accounted for in this study
215	>	by using either the reaction field technique or by resorting to a
216	>	simple cubic switching function at a cutoff radius.  One of the early
217	>	applications of a reaction field was actually in Monte Carlo
218	>	simulations of liquid water.\cite{Barker73} Under this method, the
219	>	magnitude of the reaction field acting on dipole $i$ is
220		\begin{equation}
221		\mathcal{E}_{i} = \frac{2(\varepsilon_{s} - 1)}{2\varepsilon_{s} + 1}
222	<	\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} \boldsymbol{\mu}_{j} f(r_{ij})\  ,
222	>	\frac{1}{r_{c}^{3}} \sum_{j\in{\mathcal{R}}} {\bf \mu}_{j} s(r_{ij}),
223		\label{rfequation}
224		\end{equation}
225		where $\mathcal{R}$ is the cavity defined by the cutoff radius
226		($r_{c}$), $\varepsilon_{s}$ is the dielectric constant imposed on the
227	<	system (80 in this case), $\boldsymbol{\mu}_{j}$ is the dipole moment
228	<	vector of particle \emph{j}, and $f(r_{ij})$ is a cubic switching
227	>	system (80 in the case of liquid water), ${\bf \mu}_{j}$ is the dipole
228	>	moment vector of particle $j$, and $s(r_{ij})$ is a cubic switching
229		function.\cite{AllenTildesley} The reaction field contribution to the
230	<	total energy by particle \emph{i} is given by
231	<	$-\frac{1}{2}\boldsymbol{\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque
232	<	on dipole \emph{i} by
233	<	$\boldsymbol{\mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley} Use
234	<	of reaction field is known to alter the orientational dynamic
235	<	properties, such as the dielectric relaxation time, based on changes
236	<	in the length of the cutoff radius.\cite{Berendsen98} This variable
237	<	behavior makes reaction field a less attractive method than other
238	<	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.
230	>	total energy by particle $i$ is given by $-\frac{1}{2}{\bf
231	>	\mu}_{i}\cdot\mathcal{E}_{i}$ and the torque on dipole $i$ by ${\bf
232	>	\mu}_{i}\times\mathcal{E}_{i}$.\cite{AllenTildesley}  Use of the reaction
233	>	field is known to alter the bulk orientational properties of simulated
234	>	water, and there is particular sensitivity of these properties on
235	>	changes in the length of the cutoff radius.\cite{Berendsen98} This
236	>	variable behavior makes reaction field a less attractive method than
237	>	the Ewald sum.  However, for very large systems, the computational
238	>	benefit of reaction field is dramatic.
239
240	<	Simulations were performed in both the isobaric-isothermal and
241	<	microcanonical ensembles. The constant pressure simulations were
240	>	We have also performed a companion set of simulations {\it without} a
241	>	surrounding dielectric (i.e. using a simple cubic switching function
242	>	at the cutoff radius), and as a result we have two reparamaterizations
243	>	of SSD which could be used either with or without the reaction
244	>	field turned on.
245	>
246	>	Simulations to obtain the preferred densities of the models were
247	>	performed in the isobaric-isothermal (NPT) ensemble, while all
248	>	dynamical properties were obtained from microcanonical (NVE)
249	>	simulations done at densities matching the NPT density for a
250	>	particular target temperature.  The constant pressure simulations were
251		implemented using an integral thermostat and barostat as outlined by
252	<	Hoover.\cite{Hoover85,Hoover86} For the constant pressure
253	<	simulations, the \emph{Q} parameter for the was set to 5.0 amu
254	<	\(\cdot\)\AA\(^{2}\), and the relaxation time (\(\tau\))\ was set at
255	<	100 ps.
252	>	Hoover.\cite{Hoover85,Hoover86} All molecules were treated as
253	>	non-linear rigid bodies. Vibrational constraints are not necessary in
254	>	simulations of SSD, because there are no explicit hydrogen
255	>	atoms, and thus no molecular vibrational modes need to be considered.
256
257		Integration of the equations of motion was carried out using the
258	<	symplectic splitting method proposed by Dullweber \emph{et
259	<	al.}.\cite{Dullweber1997} The reason for this integrator selection
260	<	deals with poor energy conservation of rigid body systems using
261	<	quaternions. While quaternions work well for orientational motion in
262	<	alternate ensembles, the microcanonical ensemble has a constant energy
263	<	requirement that is actually quite sensitive to errors in the
264	<	equations of motion. The original implementation of this code utilized
265	<	quaternions for rotational motion propagation; however, a detailed
266	<	investigation showed that they resulted in a steady drift in the total
219	<	energy, something that has been observed by others.\cite{Laird97}
258	>	symplectic splitting method proposed by Dullweber, Leimkuhler, and
259	>	McLachlan ({\sc dlm}).\cite{Dullweber1997} Our reason for selecting
260	>	this integrator centers on poor energy conservation of rigid body
261	>	dynamics using traditional quaternion
262	>	integration.\cite{Evans77,Evans77b} In typical microcanonical ensemble
263	>	simulations, the energy drift when using quaternions was substantially
264	>	greater than when using the {\sc dlm} method (fig. \ref{timestep}).
265	>	This steady drift in the total energy has also been observed by Kol
266	>	{\it et al.}\cite{Laird97}
267
268		The key difference in the integration method proposed by Dullweber
269		\emph{et al.} is that the entire rotation matrix is propagated from
270	<	one time step to the next. In the past, this would not have been as
271	<	feasible a option, being that the rotation matrix for a single body is
272	<	nine elements long as opposed to 3 or 4 elements for Euler angles and
273	<	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.
270	>	one time step to the next.  The additional memory required by the
271	>	algorithm is inconsequential on modern computers, and translating the
272	>	rotation matrix into quaternions for storage purposes makes trajectory
273	>	data quite compact.
274
275	<	The symplectic splitting method allows for Verlet style integration of
276	<	both linear and angular motion of rigid bodies. In the integration
277	<	method, the orientational propagation involves a sequence of matrix
278	<	evaluations to update the rotation matrix.\cite{Dullweber1997} These
279	<	matrix rotations end up being more costly computationally than the
280	<	simpler arithmetic quaternion propagation. On average, a 1000 SSD
281	<	particle simulation shows a 7\% increase in simulation time using the
282	<	symplectic step method in place of quaternions. This cost is more than
283	<	justified when comparing the energy conservation of the two methods as
284	<	illustrated in figure \ref{timestep}.
275	>	The {\sc dlm} method allows for Verlet style integration of both
276	>	translational and orientational motion of rigid bodies. In this
277	>	integration method, the orientational propagation involves a sequence
278	>	of matrix evaluations to update the rotation
279	>	matrix.\cite{Dullweber1997} These matrix rotations are more costly
280	>	than the simpler arithmetic quaternion propagation. With the same time
281	>	step, a 1000 SSD particle simulation shows an average 7\%
282	>	increase in computation time using the {\sc dlm} method in place of
283	>	quaternions. The additional expense per step is justified when one
284	>	considers the ability to use time steps that are nearly twice as large
285	>	under {\sc dlm} than would be usable under quaternion dynamics.  The
286	>	energy conservation of the two methods using a number of different
287	>	time steps is illustrated in figure
288	>	\ref{timestep}.
289
290	<	\begin{figure}
291	<	\includegraphics[width=61mm, angle=-90]{timeStep.epsi}
292	<	\caption{Energy conservation using quaternion based integration versus
293	<	the symplectic step method proposed by Dullweber \emph{et al.} with
294	<	increasing time step. For each time step, the dotted line is total
295	<	energy using the symplectic step integrator, and the solid line comes
296	<	from the quaternion integrator. The larger time step plots are shifted
297	<	up from the true energy baseline for clarity.}
298	<	\label{timestep}
299	<	\end{figure}
290	>	%\begin{figure}
291	>	%\begin{center}
292	>	%\epsfxsize=6in
293	>	%\epsfbox{timeStep.epsi}
294	>	%\caption{Energy conservation using both quaternion-based integration and
295	>	%the {\sc dlm} method with increasing time step. The larger time step
296	>	%plots are shifted from the true energy baseline (that of $\Delta t$ =
297	>	%0.1~fs) for clarity.}
298	>	%\label{timestep}
299	>	%\end{center}
300	>	%\end{figure}
301
302		In figure \ref{timestep}, the resulting energy drift at various time
303	<	steps for both the symplectic step and quaternion integration schemes
304	<	is compared. All of the 1000 SSD particle simulations started with the
305	<	same configuration, and the only difference was the method for
306	<	handling rotational motion. At time steps of 0.1 and 0.5 fs, both
307	<	methods for propagating particle rotation conserve energy fairly well,
308	<	with the quaternion method showing a slight energy drift over time in
309	<	the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
310	<	energy conservation benefits of the symplectic step method are clearly
311	<	demonstrated.
303	>	steps for both the {\sc dlm} and quaternion integration schemes is
304	>	compared.  All of the 1000 SSD particle simulations started with
305	>	the same configuration, and the only difference was the method used to
306	>	handle orientational motion. At time steps of 0.1 and 0.5~fs, both
307	>	methods for propagating the orientational degrees of freedom conserve
308	>	energy fairly well, with the quaternion method showing a slight energy
309	>	drift over time in the 0.5~fs time step simulation. At time steps of 1
310	>	and 2~fs, the energy conservation benefits of the {\sc dlm} method are
311	>	clearly demonstrated. Thus, while maintaining the same degree of
312	>	energy conservation, one can take considerably longer time steps,
313	>	leading to an overall reduction in computation time.
314
315	<	Energy drift in these SSD particle simulations was unnoticeable for
316	<	time steps up to three femtoseconds. A slight energy drift on the
317	<	order of 0.012 kcal/mol per nanosecond was observed at a time step of
318	<	four femtoseconds, and as expected, this drift increases dramatically
319	<	with increasing time step. To insure accuracy in the constant energy
320	<	simulations, time steps were set at 2 fs and kept at this value for
315	>	Energy drift in the simulations using {\sc dlm} integration was
316	>	unnoticeable for time steps up to 3~fs. A slight energy drift on the
317	>	order of 0.012~kcal/mol per nanosecond was observed at a time step of
318	>	4~fs, and as expected, this drift increases dramatically with
319	>	increasing time step. To insure accuracy in our microcanonical
320	>	simulations, time steps were set at 2~fs and kept at this value for
321		constant pressure simulations as well.
322
323	<	Ice crystals in both the $I_h$ and $I_c$ lattices were generated as
324	<	starting points for all the simulations. The $I_h$ crystals were
325	<	formed by first arranging the center of masses of the SSD particles
326	<	into a ``hexagonal'' ice lattice of 1024 particles. Because of the
327	<	crystal structure of $I_h$ ice, the simulation box assumed a
328	<	rectangular shape with a edge length ratio of approximately
323	>	Proton-disordered ice crystals in both the $I_h$ and $I_c$ lattices
324	>	were generated as starting points for all simulations. The $I_h$
325	>	crystals were formed by first arranging the centers of mass of the SSD
326	>	particles into a ``hexagonal'' ice lattice of 1024 particles. Because
327	>	of the crystal structure of $I_h$ ice, the simulation box assumed an
328	>	orthorhombic shape with an edge length ratio of approximately
329		1.00$\times$1.06$\times$1.23. The particles were then allowed to
330		orient freely about fixed positions with angular momenta randomized at
331	<	400 K for varying times. The rotational temperature was then scaled
332	<	down in stages to slowly cool the crystals down to 25 K. The particles
333	<	were then allowed translate with fixed orientations at a constant
334	<	pressure of 1 atm for 50 ps at 25 K. Finally, all constraints were
335	<	removed and the ice crystals were allowed to equilibrate for 50 ps at
336	<	25 K and a constant pressure of 1 atm.  This procedure resulted in
331	>	400~K for varying times. The rotational temperature was then scaled
332	>	down in stages to slowly cool the crystals to 25~K. The particles were
333	>	then allowed to translate with fixed orientations at a constant
334	>	pressure of 1 atm for 50~ps at 25~K. Finally, all constraints were
335	>	removed and the ice crystals were allowed to equilibrate for 50~ps at
336	>	25~K and a constant pressure of 1~atm.  This procedure resulted in
337		structurally stable $I_h$ ice crystals that obey the Bernal-Fowler
338	<	rules\cite{Bernal33,Rahman72}.  This method was also utilized in the
338	>	rules.\cite{Bernal33,Rahman72} This method was also utilized in the
339		making of diamond lattice $I_c$ ice crystals, with each cubic
340		simulation box consisting of either 512 or 1000 particles. Only
341		isotropic volume fluctuations were performed under constant pressure,
#	Line 297 | Line 345 | constant pressure and temperature dynamics. This invol
345		\section{Results and discussion}
346
347		Melting studies were performed on the randomized ice crystals using
348	<	constant pressure and temperature dynamics. This involved an initial
349	<	randomization of velocities about the starting temperature of 25 K for
350	<	varying amounts of time. The systems were all equilibrated for 100 ps
351	<	prior to a 200 ps data collection run at each temperature setting,
352	<	ranging from 25 to 400 K, with a maximum degree increment of 25 K. For
353	<	regions of interest along this stepwise progression, the temperature
354	<	increment was decreased from 25 K to 10 and then 5 K. The above
355	<	equilibration and production times were sufficient in that the system
356	<	volume fluctuations dampened out in all but the very cold simulations
357	<	(below 225 K). In order to further improve statistics, five separate
358	<	simulation progressions were performed, and the averaged results from
359	<	the $I_h$ melting simulations are shown in figure \ref{dense1}.
348	>	isobaric-isothermal (NPT) dynamics. During melting simulations, the
349	>	melting transition and the density maximum can both be observed,
350	>	provided that the density maximum occurs in the liquid and not the
351	>	supercooled regime. An ensemble average from five separate melting
352	>	simulations was acquired, each starting from different ice crystals
353	>	generated as described previously. All simulations were equilibrated
354	>	for 100~ps prior to a 200~ps data collection run at each temperature
355	>	setting. The temperature range of study spanned from 25 to 400~K, with
356	>	a maximum degree increment of 25~K. For regions of interest along this
357	>	stepwise progression, the temperature increment was decreased from
358	>	25~K to 10 and 5~K.  The above equilibration and production times were
359	>	sufficient in that fluctuations in the volume autocorrelation function
360	>	were damped out in all simulations in under 20~ps.
361
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	–
362		\subsection{Density Behavior}
363	<	In the initial average density versus temperature plot, the density
364	<	maximum clearly appears between 255 and 265 K. The calculated
365	<	densities within this range were nearly indistinguishable, as can be
366	<	seen in the zoom of this region of interest, shown in figure
367	<	\ref{dense1}. The greater certainty of the average value at 260 K makes
368	<	a good argument for the actual density maximum residing at this
369	<	midpoint value. Figure \ref{dense1} was constructed using ice $I_h$
370	<	crystals for the initial configuration; and though not pictured, the
371	<	simulations starting from ice $I_c$ crystal configurations showed
372	<	similar results, with a liquid-phase density maximum in this same
373	<	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.
363	>
364	>	Our initial simulations focused on the original SSD water model,
365	>	and an average density versus temperature plot is shown in figure
366	>	\ref{dense1}. Note that the density maximum when using a reaction
367	>	field appears between 255 and 265~K.  There were smaller fluctuations
368	>	in the density at 260~K than at either 255 or 265~K, so we report this
369	>	value as the location of the density maximum. Figure \ref{dense1} was
370	>	constructed using ice $I_h$ crystals for the initial configuration;
371	>	though not pictured, the simulations starting from ice $I_c$ crystal
372	>	configurations showed similar results, with a liquid-phase density
373	>	maximum in this same region (between 255 and 260~K).
374
375	<	\begin{figure}
376	<	\includegraphics[width=65mm,angle=-90]{dense2.eps}
377	<	\caption{Density versus temperature for TIP4P\cite{Jorgensen98b},
378	<	TIP3P\cite{Jorgensen98b}, SPC/E\cite{Clancy94}, SSD without Reaction
379	<	Field, SSD, and Experiment\cite{CRC80}. }
380	<	\label{dense2}
381	<	\end{figure}
375	>	%\begin{figure}
376	>	%\begin{center}
377	>	%\epsfxsize=6in
378	>	%\epsfbox{denseSSDnew.eps}
379	>	%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
380	>	% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
381	>	% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
382	>	% arrows indicate the change in densities observed when turning off the
383	>	% reaction field. The the lower than expected densities for the SSD
384	>	% model were what prompted the original reparameterization of SSD1
385	>	% [Ref. \onlinecite{Ichiye03}].}
386	>	%\label{dense1}
387	>	%\end{center}
388	>	%\end{figure}
389
390	<	The density maximum for SSD actually compares quite favorably to other
391	<	simple water models. Figure \ref{dense2} shows a plot of these
392	<	findings with the density progression of several other models and
348	<	experiment obtained from other
390	>	The density maximum for SSD compares quite favorably to other
391	>	simple water models. Figure \ref{dense1} also shows calculated
392	>	densities of several other models and experiment obtained from other
393		sources.\cite{Jorgensen98b,Clancy94,CRC80} Of the listed simple water
394	<	models, SSD has results closest to the experimentally observed water
395	<	density maximum. Of the listed water models, TIP4P has a density
396	<	maximum behavior most like that seen in SSD. Though not shown, it is
397	<	useful to note that TIP5P has a water density maximum nearly identical
398	<	to experiment.
394	>	models, SSD has a temperature closest to the experimentally
395	>	observed density maximum. Of the {\it charge-based} models in
396	>	Fig. \ref{dense1}, TIP4P has a density maximum behavior most like that
397	>	seen in SSD. Though not included in this plot, it is useful to
398	>	note that TIP5P has a density maximum nearly identical to the
399	>	experimentally measured temperature.
400
401	<	Possibly of more importance is the density scaling of SSD relative to
402	<	other common models at any given temperature (Fig. \ref{dense2}). Note
403	<	that the SSD model assumes a lower density than any of the other
401	>	It has been observed that liquid state densities in water are
402	>	dependent on the cutoff radius used both with and without the use of
403	>	reaction field.\cite{Berendsen98} In order to address the possible
404	>	effect of cutoff radius, simulations were performed with a dipolar
405	>	cutoff radius of 12.0~\AA\ to complement the previous SSD
406	>	simulations, all performed with a cutoff of 9.0~\AA. All of the
407	>	resulting densities overlapped within error and showed no significant
408	>	trend toward lower or higher densities as a function of cutoff radius,
409	>	for simulations both with and without reaction field. These results
410	>	indicate that there is no major benefit in choosing a longer cutoff
411	>	radius in simulations using SSD. This is advantageous in that
412	>	the use of a longer cutoff radius results in a significant increase in
413	>	the time required to obtain a single trajectory.
414	>
415	>	The key feature to recognize in figure \ref{dense1} is the density
416	>	scaling of SSD relative to other common models at any given
417	>	temperature. SSD assumes a lower density than any of the other
418		listed models at the same pressure, behavior which is especially
419	<	apparent at temperatures greater than 300 K. Lower than expected
420	<	densities have been observed for other systems with the use of a
421	<	reaction field for long-range electrostatic interactions, so the most
422	<	likely reason for these significantly lower densities in these
423	<	simulations is the presence of the reaction field.\cite{Berendsen98}
424	<	In order to test the effect of the reaction field on the density of
425	<	the systems, the simulations were repeated for the temperature region
426	<	of interest without a reaction field present. The results of these
427	<	simulations are also displayed in figure \ref{dense2}. Without
428	<	reaction field, these densities increase considerably to more
429	<	experimentally reasonable values, especially around the freezing point
430	<	of liquid water. The shape of the curve is similar to the curve
431	<	produced from SSD simulations using reaction field, specifically the
432	<	rapidly decreasing densities at higher temperatures; however, a slight
433	<	shift in the density maximum location, down to 245 K, is
434	<	observed. This is probably a more accurate comparison to the other
435	<	listed water models in that no long range corrections were applied in
436	<	those simulations.\cite{Clancy94,Jorgensen98b}
419	>	apparent at temperatures greater than 300~K. Lower than expected
420	>	densities have been observed for other systems using a reaction field
421	>	for long-range electrostatic interactions, so the most likely reason
422	>	for the significantly lower densities seen in these simulations is the
423	>	presence of the reaction field.\cite{Berendsen98,Nezbeda02} In order
424	>	to test the effect of the reaction field on the density of the
425	>	systems, the simulations were repeated without a reaction field
426	>	present. The results of these simulations are also displayed in figure
427	>	\ref{dense1}. Without the reaction field, the densities increase
428	>	to more experimentally reasonable values, especially around the
429	>	freezing point of liquid water. The shape of the curve is similar to
430	>	the curve produced from SSD simulations using reaction field,
431	>	specifically the rapidly decreasing densities at higher temperatures;
432	>	however, a shift in the density maximum location, down to 245~K, is
433	>	observed. This is a more accurate comparison to the other listed water
434	>	models, in that no long range corrections were applied in those
435	>	simulations.\cite{Clancy94,Jorgensen98b} However, even without the
436	>	reaction field, the density around 300~K is still significantly lower
437	>	than experiment and comparable water models. This anomalous behavior
438	>	was what lead Tan {\it et al.} to recently reparameterize
439	>	SSD.\cite{Ichiye03} Throughout the remainder of the paper our
440	>	reparamaterizations of SSD will be compared with their newer SSD1
441	>	model.
442
379	–	It has been observed that densities are dependent on the cutoff radius
380	–	used for a variety of water models in simulations both with and
381	–	without the use of reaction field.\cite{Berendsen98} In order to
382	–	address the possible affect of cutoff radius, simulations were
383	–	performed with a dipolar cutoff radius of 12.0 \AA\ to compliment the
384	–	previous SSD simulations, all performed with a cutoff of 9.0 \AA. All
385	–	the resulting densities overlapped within error and showed no
386	–	significant trend in lower or higher densities as a function of cutoff
387	–	radius, both for simulations with and without reaction field. These
388	–	results indicate that there is no major benefit in choosing a longer
389	–	cutoff radius in simulations using SSD. This is comforting in that the
390	–	use of a longer cutoff radius results in a near doubling of the time
391	–	required to compute a single trajectory.
392	–
443		\subsection{Transport Behavior}
394	–	Of importance in these types of studies are the transport properties
395	–	of the particles and how they change when altering the environmental
396	–	conditions. In order to probe transport, constant energy simulations
397	–	were performed about the average density uncovered by the constant
398	–	pressure simulations. Simulations started with randomized velocities
399	–	and underwent 50 ps of temperature scaling and 50 ps of constant
400	–	energy equilibration before obtaining a 200 ps trajectory. Diffusion
401	–	constants were calculated via root-mean square deviation analysis. The
402	–	averaged results from 5 sets of these NVE simulations is displayed in
403	–	figure \ref{diffuse}, alongside experimental, SPC/E, and TIP5P
404	–	results.\cite{Gillen72,Mills73,Clancy94,Jorgensen01}
444
445	<	\begin{figure}
446	<	\includegraphics[width=65mm, angle=-90]{betterDiffuse.epsi}
447	<	\caption{Average diffusion coefficient over increasing temperature for
448	<	SSD, SPC/E\cite{Clancy94}, TIP5P\cite{Jorgensen01}, and Experimental
449	<	data from Gillen \emph{et al.}\cite{Gillen72}, and from
450	<	Mills\cite{Mills73}.}
451	<	\label{diffuse}
452	<	\end{figure}
445	>	Accurate dynamical properties of a water model are particularly
446	>	important when using the model to study permeation or transport across
447	>	biological membranes.  In order to probe transport in bulk water,
448	>	constant energy (NVE) simulations were performed at the average
449	>	density obtained by the NPT simulations at an identical target
450	>	temperature. Simulations started with randomized velocities and
451	>	underwent 50~ps of temperature scaling and 50~ps of constant energy
452	>	equilibration before a 200~ps data collection run. Diffusion constants
453	>	were calculated via linear fits to the long-time behavior of the
454	>	mean-square displacement as a function of time. The averaged results
455	>	from five sets of NVE simulations are displayed in figure
456	>	\ref{diffuse}, alongside experimental, SPC/E, and TIP5P
457	>	results.\cite{Gillen72,Holz00,Clancy94,Jorgensen01}
458
459	+	%\begin{figure}
460	+	%\begin{center}
461	+	%\epsfxsize=6in
462	+	%\epsfbox{betterDiffuse.epsi}
463	+	%\caption{Average self-diffusion constant as a function of temperature for
464	+	%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
465	+	%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
466	+	%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
467	+	%shown, SSD has the least deviation from the experimental values. The
468	+	%rapidly increasing diffusion constants for TIP5P and SSD correspond to
469	+	%significant decreases in density at the higher temperatures.}
470	+	%\label{diffuse}
471	+	%\end{center}
472	+	%\end{figure}
473	+
474		The observed values for the diffusion constant point out one of the
475	<	strengths of the SSD model. Of the three experimental models shown,
476	<	the SSD model has the most accurate depiction of the diffusion trend
477	<	seen in experiment in both the supercooled and normal regimes. SPC/E
478	<	does a respectable job by getting similar values as SSD and experiment
479	<	around 290 K; however, it deviates at both higher and lower
480	<	temperatures, failing to predict the experimental trend. TIP5P and SSD
481	<	both start off low at the colder temperatures and tend to diffuse too
482	<	rapidly at the higher temperatures. This type of trend at the higher
483	<	temperatures is not surprising in that the densities of both TIP5P and
484	<	SSD are lower than experimental water at temperatures higher than room
485	<	temperature. When calculating the diffusion coefficients for SSD at
486	<	experimental densities, the resulting values fall more in line with
487	<	experiment at these temperatures, albeit not at standard
429	<	pressure. Results under these conditions can be found later in this
430	<	paper.
475	>	strengths of the SSD model. Of the three models shown, the SSD model
476	>	has the most accurate depiction of self-diffusion in both the
477	>	supercooled and liquid regimes.  SPC/E does a respectable job by
478	>	reproducing values similar to experiment around 290~K; however, it
479	>	deviates at both higher and lower temperatures, failing to predict the
480	>	correct thermal trend. TIP5P and SSD both start off low at colder
481	>	temperatures and tend to diffuse too rapidly at higher temperatures.
482	>	This behavior at higher temperatures is not particularly surprising
483	>	since the densities of both TIP5P and SSD are lower than experimental
484	>	water densities at higher temperatures.  When calculating the
485	>	diffusion coefficients for SSD at experimental densities
486	>	(instead of the densities from the NPT simulations), the resulting
487	>	values fall more in line with experiment at these temperatures.
488
489		\subsection{Structural Changes and Characterization}
490	+
491		By starting the simulations from the crystalline state, the melting
492	<	transition and the ice structure can be studied along with the liquid
493	<	phase behavior beyond the melting point. To locate the melting
494	<	transition, the constant pressure heat capacity (C$_\text{p}$) was
495	<	monitored in each of the simulations. In the melting simulations of
496	<	the 1024 particle ice $I_h$ simulations, a large spike in C$_\text{p}$
497	<	occurs at 245 K, indicating a first order phase transition for the
498	<	melting of these ice crystals. When the reaction field is turned off,
499	<	the melting transition occurs at 235 K.  These melting transitions are
500	<	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.
492	>	transition and the ice structure can be obtained along with the liquid
493	>	phase behavior beyond the melting point. The constant pressure heat
494	>	capacity (C$_\text{p}$) was monitored to locate the melting transition
495	>	in each of the simulations. In the melting simulations of the 1024
496	>	particle ice $I_h$ simulations, a large spike in C$_\text{p}$ occurs
497	>	at 245~K, indicating a first order phase transition for the melting of
498	>	these ice crystals. When the reaction field is turned off, the melting
499	>	transition occurs at 235~K.  These melting transitions are
500	>	considerably lower than the experimental value.
501
502	<	\begin{figure}
503	<	\includegraphics[width=85mm]{fullContours.eps}
504	<	\caption{Contour plots of 2D angular g($r$)'s for 512 SSD systems at
505	<	100 K (A \& B) and 300 K (C \& D). Contour colors are inverted for
506	<	clarity: dark areas signify peaks while light areas signify
507	<	depressions. White areas have g(\emph{r}) values below 0.5 and black
508	<	areas have values above 1.5.}
509	<	\label{contour}
454	<	\end{figure}
502	>	%\begin{figure}
503	>	%\begin{center}
504	>	%\epsfxsize=6in
505	>	%\epsfbox{corrDiag.eps}
506	>	%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
507	>	%\label{corrAngle}
508	>	%\end{center}
509	>	%\end{figure}
510
511	<	Additional analyses for understanding the melting phase-transition
512	<	process were performed via two-dimensional structure and dipole angle
513	<	correlations. Expressions for the correlations are as follows:
511	>	%\begin{figure}
512	>	%\begin{center}
513	>	%\epsfxsize=6in
514	>	%\epsfbox{fullContours.eps}
515	>	%\caption{Contour plots of 2D angular pair correlation functions for
516	>	%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
517	>	%signify regions of enhanced density while light areas signify
518	>	%depletion relative to the bulk density. White areas have pair
519	>	%correlation values below 0.5 and black areas have values above 1.5.}
520	>	%\label{contour}
521	>	%\end{center}
522	>	%\end{figure}
523
524	<	\begin{figure}
525	<	\includegraphics[width=45mm]{corrDiag.eps}
526	<	\caption{Two dimensional illustration of the angles involved in the
463	<	correlations observed in figure \ref{contour}.}
464	<	\label{corrAngle}
465	<	\end{figure}
524	>	Additional analysis of the melting process was performed using
525	>	two-dimensional structure and dipole angle correlations. Expressions
526	>	for these correlations are as follows:
527
528	<	\begin{multline}
529	<	g_{\text{AB}}(r,\cos\theta) = \\
530	<	\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\ ,
531	<	\end{multline}
532	<	\begin{multline}
533	<	g_{\text{AB}}(r,\cos\omega) = \\
534	<	\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\ ,
535	<	\end{multline}
536	<	where $\theta$ and $\omega$ refer to the angles shown in the above
537	<	illustration. By binning over both distance and the cosine of the
538	<	desired angle between the two dipoles, the g(\emph{r}) can be
539	<	dissected to determine the common dipole arrangements that constitute
540	<	the peaks and troughs. Frames A and B of figure \ref{contour} show a
541	<	relatively crystalline state of an ice $I_c$ simulation. The first
542	<	peak of the g(\emph{r}) primarily consists of the preferred hydrogen
543	<	bonding arrangements as dictated by the tetrahedral sticky potential,
544	<	one peak for the donating and the other for the accepting hydrogen
545	<	bonds. Due to the high degree of crystallinity of the sample, the
485	<	second and third solvation shells show a repeated peak arrangement
528	>	\begin{equation}
529	>	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\ ,
530	>	\end{equation}
531	>	\begin{equation}
532	>	g_{\text{AB}}(r,\cos\omega) =
533	>	\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\ ,
534	>	\end{equation}
535	>	where $\theta$ and $\omega$ refer to the angles shown in figure
536	>	\ref{corrAngle}. By binning over both distance and the cosine of the
537	>	desired angle between the two dipoles, the $g(r)$ can be analyzed to
538	>	determine the common dipole arrangements that constitute the peaks and
539	>	troughs in the standard one-dimensional $g(r)$ plots. Frames A and B
540	>	of figure \ref{contour} show results from an ice $I_c$ simulation. The
541	>	first peak in the $g(r)$ consists primarily of the preferred hydrogen
542	>	bonding arrangements as dictated by the tetrahedral sticky potential -
543	>	one peak for the hydrogen bond donor and the other for the hydrogen
544	>	bond acceptor.  Due to the high degree of crystallinity of the sample,
545	>	the second and third solvation shells show a repeated peak arrangement
546		which decays at distances around the fourth solvation shell, near the
547		imposed cutoff for the Lennard-Jones and dipole-dipole interactions.
548	<	In the higher temperature simulation shown in frames C and D, the
549	<	repeated peak features are significantly blurred. The first solvation
550	<	shell still shows the strong effect of the sticky-potential, although
551	<	it covers a larger area, extending to include a fraction of aligned
548	>	In the higher temperature simulation shown in frames C and D, these
549	>	long-range features deteriorate rapidly. The first solvation shell
550	>	still shows the strong effect of the sticky-potential, although it
551	>	covers a larger area, extending to include a fraction of aligned
552		dipole peaks within the first solvation shell. The latter peaks lose
553	<	definition as thermal motion and the competing dipole force overcomes
554	<	the sticky potential's tight tetrahedral structuring of the fluid.
553	>	due to thermal motion and as the competing dipole force overcomes the
554	>	sticky potential's tight tetrahedral structuring of the crystal.
555
556		This complex interplay between dipole and sticky interactions was
557		remarked upon as a possible reason for the split second peak in the
558	<	oxygen-oxygen g(\emph{r}).\cite{Ichiye96} At low temperatures, the
559	<	second solvation shell peak appears to have two distinct parts that
558	>	oxygen-oxygen pair correlation function,
559	>	$g_\mathrm{OO}(r)$.\cite{Ichiye96} At low temperatures, the second
560	>	solvation shell peak appears to have two distinct components that
561		blend together to form one observable peak. At higher temperatures,
562	<	this split character alters to show the leading 4 \AA\ peak dominated
563	<	by equatorial anti-parallel dipole orientations, and there is tightly
564	<	bunched group of axially arranged dipoles that most likely consist of
565	<	the smaller fraction aligned dipole pairs. The trailing part of the
566	<	split peak at 5 \AA\ is dominated by aligned dipoles that range
567	<	primarily within the axial to the chief hydrogen bond arrangements
568	<	similar to those seen in the first solvation shell. This evidence
569	<	indicates that the dipole pair interaction begins to dominate outside
570	<	of the range of the dipolar repulsion term, with the primary
571	<	energetically favorable dipole arrangements populating the region
572	<	immediately outside of it's range (around 4 \AA), and arrangements
573	<	that seek to ideally satisfy both the sticky and dipole forces locate
574	<	themselves just beyond this region (around 5 \AA).
562	>	this split character alters to show the leading 4~\AA\ peak dominated
563	>	by equatorial anti-parallel dipole orientations. There is also a
564	>	tightly bunched group of axially arranged dipoles that most likely
565	>	consist of the smaller fraction of aligned dipole pairs. The trailing
566	>	component of the split peak at 5~\AA\ is dominated by aligned dipoles
567	>	that assume hydrogen bond arrangements similar to those seen in the
568	>	first solvation shell. This evidence indicates that the dipole pair
569	>	interaction begins to dominate outside of the range of the dipolar
570	>	repulsion term.  The energetically favorable dipole arrangements
571	>	populate the region immediately outside this repulsion region (around
572	>	4~\AA), while arrangements that seek to satisfy both the sticky and
573	>	dipole forces locate themselves just beyond this initial buildup
574	>	(around 5~\AA).
575
576		From these findings, the split second peak is primarily the product of
577	<	the dipolar repulsion term of the sticky potential. In fact, the
578	<	leading of the two peaks can be pushed out and merged with the outer
579	<	split peak just by extending the switching function cutoff
580	<	($s^\prime(r_{ij})$) from its normal 4.0 \AA\ to values of 4.5 or even
581	<	5 \AA. This type of correction is not recommended for improving the
582	<	liquid structure, because the second solvation shell will still be
583	<	shifted too far out. In addition, this would have an even more
584	<	detrimental effect on the system densities, leading to a liquid with a
585	<	more open structure and a density considerably lower than the normal
586	<	SSD behavior shown previously. A better correction would be to include
587	<	the quadrupole-quadrupole interactions for the water particles outside
588	<	of the first solvation shell, but this reduces the simplicity and
589	<	speed advantage of SSD, so it is not the most desirable path to take.
577	>	the dipolar repulsion term of the sticky potential. In fact, the inner
578	>	peak can be pushed out and merged with the outer split peak just by
579	>	extending the switching function ($s^\prime(r_{ij})$) from its normal
580	>	4.0~\AA\ cutoff to values of 4.5 or even 5~\AA. This type of
581	>	correction is not recommended for improving the liquid structure,
582	>	since the second solvation shell would still be shifted too far
583	>	out. In addition, this would have an even more detrimental effect on
584	>	the system densities, leading to a liquid with a more open structure
585	>	and a density considerably lower than the already low SSD
586	>	density.  A better correction would be to include the
587	>	quadrupole-quadrupole interactions for the water particles outside of
588	>	the first solvation shell, but this would remove the simplicity and
589	>	speed advantage of SSD.
590
591	<	\subsection{Adjusted Potentials: SSD/E and SSD/RF}
591	>	\subsection{Adjusted Potentials: SSD/RF and SSD/E}
592	>
593		The propensity of SSD to adopt lower than expected densities under
594		varying conditions is troubling, especially at higher temperatures. In
595	<	order to correct this behavior, it's necessary to adjust the force
596	<	field parameters for the primary intermolecular interactions. In
597	<	undergoing a reparameterization, it is important not to focus on just
598	<	one property and neglect the other important properties. In this case,
599	<	it would be ideal to correct the densities while maintaining the
600	<	accurate transport properties.
595	>	order to correct this model for use with a reaction field, it is
596	>	necessary to adjust the force field parameters for the primary
597	>	intermolecular interactions. In undergoing a reparameterization, it is
598	>	important not to focus on just one property and neglect the other
599	>	important properties. In this case, it would be ideal to correct the
600	>	densities while maintaining the accurate transport behavior.
601
602	<	The possible parameters for tuning include the $\sigma$ and $\epsilon$
603	<	Lennard-Jones parameters, the dipole strength ($\mu$), and the sticky
604	<	attractive and dipole repulsive terms with their respective
605	<	cutoffs. To alter the attractive and repulsive terms of the sticky
606	<	potential independently, it is necessary to separate the terms as
607	<	follows:
608	<	\begin{equation}
609	<	\begin{split}
610	<	u_{ij}^{sp}
611	<	(\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}
602	>	The parameters available for tuning include the $\sigma$ and
603	>	$\epsilon$ Lennard-Jones parameters, the dipole strength ($\mu$), the
604	>	strength of the sticky potential ($\nu_0$), and the cutoff distances
605	>	for the sticky attractive and dipole repulsive cubic switching
606	>	function cutoffs ($r_l$, $r_u$ and $r_l^\prime$, $r_u^\prime$
607	>	respectively). The results of the reparameterizations are shown in
608	>	table \ref{params}. We are calling these reparameterizations the Soft
609	>	Sticky Dipole / Reaction Field (SSD/RF - for use with a reaction
610	>	field) and Soft Sticky Dipole Extended (SSD/E - an attempt to improve
611	>	the liquid structure in simulations without a long-range correction).
612
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	–
613		\begin{table}
614	+	\begin{center}
615		\caption{Parameters for the original and adjusted models}
616	<	\begin{tabular}{ l  c  c  c }
616	>	\begin{tabular}{ l  c  c  c  c }
617		\hline \\[-3mm]
618	<	\ Parameters & \ \ \  SSD$^\dagger$\ \ \ \  & \ SSD/E\ \  & \ SSD/RF\ \ \\
618	>	\ \ \ Parameters\ \ \  & \ \ \ SSD [Ref. \onlinecite{Ichiye96}] \ \ \
619	>	& \ SSD1 [Ref. \onlinecite{Ichiye03}]\ \  & \ SSD/E\ \  & \ \ SSD/RF \\
620		\hline \\[-3mm]
621	<	\ \ \ $\sigma$ (\AA)  & 3.051 & 3.035 & 3.019\\
622	<	\ \ \ $\epsilon$ (kcal/mol)\ \ & 0.152 & 0.152 & 0.152\\
623	<	\ \ \ $\mu$ (D) & 2.35 & 2.418 & 2.480\\
624	<	\ \ \ $\nu_0$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
625	<	\ \ \ $r_l$ (\AA) & 2.75 & 2.40 & 2.40\\
626	<	\ \ \ $r_u$ (\AA) & 3.35 & 3.80 & 3.80\\
627	<	\ \ \ $\nu_0^\prime$ (kcal/mol)\ \ & 3.7284 & 3.90 & 3.90\\
628	<	\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75\\
629	<	\ \ \ $r_u^\prime$ (\AA) & 4.00 & 3.35 & 3.35\\
582	<	\\[-2mm]$^\dagger$ ref. \onlinecite{Ichiye96}
621	>	\ \ \ $\sigma$ (\AA)  & 3.051 & 3.016 & 3.035 & 3.019\\
622	>	\ \ \ $\epsilon$ (kcal/mol) & 0.152 & 0.152 & 0.152 & 0.152\\
623	>	\ \ \ $\mu$ (D) & 2.35 & 2.35 & 2.42 & 2.48\\
624	>	\ \ \ $\nu_0$ (kcal/mol) & 3.7284 & 3.6613 & 3.90 & 3.90\\
625	>	\ \ \ $\omega^\circ$ & 0.07715 & 0.07715 & 0.07715 & 0.07715\\
626	>	\ \ \ $r_l$ (\AA) & 2.75 & 2.75 & 2.40 & 2.40\\
627	>	\ \ \ $r_u$ (\AA) & 3.35 & 3.35 & 3.80 & 3.80\\
628	>	\ \ \ $r_l^\prime$ (\AA) & 2.75 & 2.75 & 2.75 & 2.75\\
629	>	\ \ \ $r_u^\prime$ (\AA) & 4.00 & 4.00 & 3.35 & 3.35\\
630		\end{tabular}
631		\label{params}
632	+	\end{center}
633		\end{table}
634
635	<	\begin{figure}
636	<	\includegraphics[width=85mm]{gofrCompare.epsi}
637	<	\caption{Plots comparing experiment\cite{Head-Gordon00_1} with SSD/E
638	<	and SSD without reaction field (top), as well as SSD/RF and SSD with
639	<	reaction field turned on (bottom). The insets show the respective
640	<	first peaks in detail. Solid Line - experiment, dashed line - SSD/E
641	<	and SSD/RF, and dotted line - SSD (with and without reaction field).}
642	<	\label{grcompare}
643	<	\end{figure}
644	<
645	<	\begin{figure}
646	<	\includegraphics[width=85mm]{dualsticky.ps}
647	<	\caption{Isosurfaces of the sticky potential for SSD (left) and SSD/E \&
600	<	SSD/RF (right). Light areas correspond to the tetrahedral attractive
601	<	part, and the darker areas correspond to the dipolar repulsive part.}
602	<	\label{isosurface}
603	<	\end{figure}
635	>	%\begin{figure}
636	>	%\begin{center}
637	>	%\epsfxsize=5in
638	>	%\epsfbox{GofRCompare.epsi}
639	>	%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with
640	>	%SSD/E and SSD1 without reaction field (top), as well as
641	>	%SSD/RF and SSD1 with reaction field turned on
642	>	%(bottom). The insets show the respective first peaks in detail. Note
643	>	%how the changes in parameters have lowered and broadened the first
644	>	%peak of SSD/E and SSD/RF.}
645	>	%\label{grcompare}
646	>	%\end{center}
647	>	%\end{figure}
648
649	<	In the paper detailing the development of SSD, Liu and Ichiye placed
650	<	particular emphasis on an accurate description of the first solvation
651	<	shell. This resulted in a somewhat tall and sharp first peak that
652	<	integrated to give similar coordination numbers to the experimental
653	<	data obtained by Soper and Phillips.\cite{Ichiye96,Soper86} New
654	<	experimental x-ray scattering data from the Head-Gordon lab indicates
655	<	a slightly lower and shifted first peak in the g$_\mathrm{OO}(r)$, so
656	<	adjustments to SSD were made while taking into consideration the new
657	<	experimental findings.\cite{Head-Gordon00_1} Figure \ref{grcompare}
658	<	shows the relocation of the first peak of the oxygen-oxygen
659	<	g(\emph{r}) by comparing the original SSD (with and without reaction
616	<	field), SSD-E, and SSD-RF to the new experimental results. Both the
617	<	modified water models have shorter peaks that are brought in more
618	<	closely to the experimental peak (as seen in the insets of figure
619	<	\ref{grcompare}). This structural alteration was accomplished by a
620	<	reduction in the Lennard-Jones $\sigma$ variable as well as adjustment
621	<	of the sticky potential strength and cutoffs. The cutoffs for the
622	<	tetrahedral attractive and dipolar repulsive terms were nearly swapped
623	<	with each other. Isosurfaces of the original and modified sticky
624	<	potentials are shown in figure \cite{isosurface}. In these
625	<	isosurfaces, it is easy to see how altering the cutoffs changes the
626	<	repulsive and attractive character of the particles. With a reduced
627	<	repulsive surface (the darker region), the particles can move closer
628	<	to one another, increasing the density for the overall system. This
629	<	change in interaction cutoff also results in a more gradual
630	<	orientational motion by allowing the particles to maintain preferred
631	<	dipolar arrangements before they begin to feel the pull of the
632	<	tetrahedral restructuring. Upon moving closer together, the dipolar
633	<	repulsion term becomes active and excludes the unphysical
634	<	arrangements. This compares with the original SSD's excluding dipolar
635	<	before the particles feel the pull of the ``hydrogen bonds''. Aside
636	<	from improving the shape of the first peak in the g(\emph{r}), this
637	<	improves the densities considerably by allowing the persistence of
638	<	full dipolar character below the previous 4.0 \AA\ cutoff.
649	>	%\begin{figure}
650	>	%\begin{center}
651	>	%\epsfxsize=6in
652	>	%\epsfbox{dualsticky_bw.eps}
653	>	%\caption{Positive and negative isosurfaces of the sticky potential for
654	>	%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
655	>	%correspond to the tetrahedral attractive component, and darker areas
656	>	%correspond to the dipolar repulsive component.}
657	>	%\label{isosurface}
658	>	%\end{center}
659	>	%\end{figure}
660
661	<	While adjusting the location and shape of the first peak of
662	<	g(\emph{r}) improves the densities to some degree, these changes alone
663	<	are insufficient to bring the system densities up to the values
664	<	observed experimentally. To finish bringing up the densities, the
665	<	dipole moments were increased in both the adjusted models. Being a
666	<	dipole based model, the structure and transport are very sensitive to
667	<	changes in the dipole moment. The original SSD simply used the dipole
668	<	moment calculated from the TIP3P water model, which at 2.35 D is
669	<	significantly greater than the experimental gas phase value of 1.84
670	<	D. The larger dipole moment is a more realistic value and improve the
671	<	dielectric properties of the fluid. Both theoretical and experimental
672	<	measurements indicate a liquid phase dipole moment ranging from 2.4 D
673	<	to values as high as 3.11 D, so there is quite a range available for
674	<	adjusting the dipole
675	<	moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} The increasing of
676	<	the dipole moments to 2.418 and 2.48 D for SSD/E and SSD/RF
677	<	respectively is moderate in the range of the experimental values;
678	<	however, it leads to significant changes in the density and transport
679	<	of the water models.
661	>	In the original paper detailing the development of SSD, Liu and Ichiye
662	>	placed particular emphasis on an accurate description of the first
663	>	solvation shell. This resulted in a somewhat tall and narrow first
664	>	peak in $g(r)$ that integrated to give similar coordination numbers to
665	>	the experimental data obtained by Soper and
666	>	Phillips.\cite{Ichiye96,Soper86} New experimental x-ray scattering
667	>	data from the Head-Gordon lab indicates a slightly lower and shifted
668	>	first peak in the g$_\mathrm{OO}(r)$, so our adjustments to SSD were
669	>	made after taking into consideration the new experimental
670	>	findings.\cite{Head-Gordon00_1} Figure \ref{grcompare} shows the
671	>	relocation of the first peak of the oxygen-oxygen $g(r)$ by comparing
672	>	the revised SSD model (SSD1), SSD/E, and SSD/RF to the new
673	>	experimental results. Both modified water models have shorter peaks
674	>	that match more closely to the experimental peak (as seen in the
675	>	insets of figure \ref{grcompare}).  This structural alteration was
676	>	accomplished by the combined reduction in the Lennard-Jones $\sigma$
677	>	variable and adjustment of the sticky potential strength and cutoffs.
678	>	As can be seen in table \ref{params}, the cutoffs for the tetrahedral
679	>	attractive and dipolar repulsive terms were nearly swapped with each
680	>	other.  Isosurfaces of the original and modified sticky potentials are
681	>	shown in figure \ref{isosurface}. In these isosurfaces, it is easy to
682	>	see how altering the cutoffs changes the repulsive and attractive
683	>	character of the particles. With a reduced repulsive surface (darker
684	>	region), the particles can move closer to one another, increasing the
685	>	density for the overall system.  This change in interaction cutoff
686	>	also results in a more gradual orientational motion by allowing the
687	>	particles to maintain preferred dipolar arrangements before they begin
688	>	to feel the pull of the tetrahedral restructuring. As the particles
689	>	move closer together, the dipolar repulsion term becomes active and
690	>	excludes unphysical nearest-neighbor arrangements. This compares with
691	>	how SSD and SSD1 exclude preferred dipole alignments before the
692	>	particles feel the pull of the ``hydrogen bonds''. Aside from
693	>	improving the shape of the first peak in the g(\emph{r}), this
694	>	modification improves the densities considerably by allowing the
695	>	persistence of full dipolar character below the previous 4.0~\AA\
696	>	cutoff.
697
698	<	In order to demonstrate the benefits of this reparameterization, a
698	>	While adjusting the location and shape of the first peak of $g(r)$
699	>	improves the densities, these changes alone are insufficient to bring
700	>	the system densities up to the values observed experimentally.  To
701	>	further increase the densities, the dipole moments were increased in
702	>	both of our adjusted models. Since SSD is a dipole based model, the
703	>	structure and transport are very sensitive to changes in the dipole
704	>	moment. The original SSD simply used the dipole moment calculated from
705	>	the TIP3P water model, which at 2.35~D is significantly greater than
706	>	the experimental gas phase value of 1.84~D. The larger dipole moment
707	>	is a more realistic value and improves the dielectric properties of
708	>	the fluid. Both theoretical and experimental measurements indicate a
709	>	liquid phase dipole moment ranging from 2.4~D to values as high as
710	>	3.11~D, providing a substantial range of reasonable values for a
711	>	dipole moment.\cite{Sprik91,Kusalik02,Badyal00,Barriol64} Moderately
712	>	increasing the dipole moments to 2.42 and 2.48~D for SSD/E and SSD/RF,
713	>	respectively, leads to significant changes in the density and
714	>	transport of the water models.
715	>
716	>	In order to demonstrate the benefits of these reparameterizations, a
717		series of NPT and NVE simulations were performed to probe the density
718		and transport properties of the adapted models and compare the results
719		to the original SSD model. This comparison involved full NPT melting
720		sequences for both SSD/E and SSD/RF, as well as NVE transport
721	<	calculations at both self-consistent and experimental
722	<	densities. Again, the results come from five separate simulations of
723	<	1024 particle systems, and the melting sequences were started from
724	<	different ice $I_h$ crystals constructed as stated previously. Like
725	<	before, all of the NPT simulations were equilibrated for 100 ps before
726	<	a 200 ps data collection run, and they used the previous temperature's
727	<	final configuration as a starting point. All of the NVE simulations
728	<	had the same thermalization, equilibration, and data collection times
729	<	stated earlier in this paper.
721	>	calculations at the calculated self-consistent densities. Again, the
722	>	results are obtained from five separate simulations of 1024 particle
723	>	systems, and the melting sequences were started from different ice
724	>	$I_h$ crystals constructed as described previously. Each NPT
725	>	simulation was equilibrated for 100~ps before a 200~ps data collection
726	>	run at each temperature step, and the final configuration from the
727	>	previous temperature simulation was used as a starting point. All NVE
728	>	simulations had the same thermalization, equilibration, and data
729	>	collection times as stated previously.
730
731	<	\begin{figure}
732	<	\includegraphics[width=85mm]{ssdecompare.epsi}
733	<	\caption{Comparison of densities calculated with SSD/E to SSD without a
734	<	reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
735	<	SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
736	<	includes error bars, and the calculated results from the other
737	<	references were removed for clarity.}
738	<	\label{ssdedense}
739	<	\end{figure}
731	>	%\begin{figure}
732	>	%\begin{center}
733	>	%\epsfxsize=6in
734	>	%\epsfbox{ssdeDense.epsi}
735	>	%\caption{Comparison of densities calculated with SSD/E to
736	>	%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
737	>	%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
738	>	%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
739	>	%300 K with error bars included to clarify this region of
740	>	%interest. Note that both SSD1 and SSD/E show good agreement with
741	>	%experiment when the long-range correction is neglected.}
742	>	%\label{ssdedense}
743	>	%\end{center}
744	>	%\end{figure}
745
746	<	Figure \ref{ssdedense} shows the density profile for the SSD/E water
747	<	model in comparison to the original SSD without a reaction field,
748	<	experiment, and the other common water models considered
749	<	previously. The calculated densities have increased significantly over
750	<	the original SSD model and match the experimental value just below 298
751	<	K. At 298 K, the density of SSD/E is 0.995$\pm$0.001 g/cm$^3$, which
752	<	compares well with the experimental value of 0.997 g/cm$^3$ and is
753	<	considerably better than the SSD value of 0.967$\pm$0.003
754	<	g/cm$^3$. The increased dipole moment in SSD/E has helped to flatten
755	<	out the curve at higher temperatures, only the improvement is marginal
756	<	at best. This steep drop in densities is due to the dipolar rather
757	<	than charge based interactions which decay more rapidly at longer
758	<	distances.
759	<
760	<	By monitoring C$\text{p}$ throughout these simulations, the melting
761	<	transition for SSD/E was observed at 230 K, about 5 degrees lower than
762	<	SSD. The resulting density maximum is located at 240 K, again about 5
763	<	degrees lower than the SSD value of 245 K. Though there is a decrease
764	<	in both of these values, the corrected densities near room temperature
765	<	justify the modifications taken.
746	>	Fig. \ref{ssdedense} shows the density profile for the SSD/E
747	>	model in comparison to SSD1 without a reaction field, other
748	>	common water models, and experimental results. The calculated
749	>	densities for both SSD/E and SSD1 have increased
750	>	significantly over the original SSD model (see
751	>	fig. \ref{dense1}) and are in better agreement with the experimental
752	>	values. At 298 K, the densities of SSD/E and SSD1 without
753	>	a long-range correction are 0.996$\pm$0.001 g/cm$^3$ and
754	>	0.999$\pm$0.001 g/cm$^3$ respectively.  These both compare well with
755	>	the experimental value of 0.997 g/cm$^3$, and they are considerably
756	>	better than the SSD value of 0.967$\pm$0.003 g/cm$^3$. The
757	>	changes to the dipole moment and sticky switching functions have
758	>	improved the structuring of the liquid (as seen in figure
759	>	\ref{grcompare}, but they have shifted the density maximum to much
760	>	lower temperatures. This comes about via an increase in the liquid
761	>	disorder through the weakening of the sticky potential and
762	>	strengthening of the dipolar character. However, this increasing
763	>	disorder in the SSD/E model has little effect on the melting
764	>	transition. By monitoring $C_p$ throughout these simulations, the
765	>	melting transition for SSD/E was shown to occur at 235~K.  The
766	>	same transition temperature observed with SSD and SSD1.
767
768	<	\begin{figure}
769	<	\includegraphics[width=85mm]{ssdrfcompare.epsi}
770	<	\caption{Comparison of densities calculated with SSD/RF to SSD with a
771	<	reaction field, TIP4P\cite{Jorgensen98b}, TIP3P\cite{Jorgensen98b},
772	<	SPC/E\cite{Clancy94}, and Experiment\cite{CRC80}. The upper plot
773	<	includes error bars, and the calculated results from the other
774	<	references were removed for clarity.}
775	<	\label{ssdrfdense}
776	<	\end{figure}
768	>	%\begin{figure}
769	>	%\begin{center}
770	>	%\epsfxsize=6in
771	>	%\epsfbox{ssdrfDense.epsi}
772	>	%\caption{Comparison of densities calculated with SSD/RF to
773	>	%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
774	>	%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
775	>	%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
776	>	%reparameterization when utilizing a reaction field long-ranged
777	>	%correction - SSD/RF provides significantly more accurate
778	>	%densities than SSD1 when performing room temperature
779	>	%simulations.}
780	>	%\label{ssdrfdense}
781	>	%\end{center}
782	>	%\end{figure}
783
784	<	Figure \ref{ssdrfdense} shows a density comparison between SSD/RF and
785	<	SSD with an active reaction field. Like in the simulations of SSD/E,
786	<	the densities show a dramatic increase over normal SSD. At 298 K,
787	<	SSD/RF has a density of 0.997$\pm$0.001 g/cm$^3$, right in line with
788	<	experiment and considerably better than the SSD value of
789	<	0.941$\pm$0.001 g/cm$^3$. The melting point is observed at 240 K,
790	<	which is 5 degrees lower than SSD with a reaction field, and the
791	<	density maximum at 255 K, again 5 degrees lower than SSD. The density
792	<	at higher temperature still drops off more rapidly than the charge
793	<	based models but is in better agreement than SSD/E.
784	>	Including the reaction field long-range correction in the simulations
785	>	results in a more interesting comparison.  A density profile including
786	>	SSD/RF and SSD1 with an active reaction field is shown in figure
787	>	\ref{ssdrfdense}.  As observed in the simulations without a reaction
788	>	field, the densities of SSD/RF and SSD1 show a dramatic increase over
789	>	normal SSD (see figure \ref{dense1}). At 298 K, SSD/RF has a density
790	>	of 0.997$\pm$0.001 g/cm$^3$, directly in line with experiment and
791	>	considerably better than the original SSD value of 0.941$\pm$0.001
792	>	g/cm$^3$ and the SSD1 value of 0.972$\pm$0.002 g/cm$^3$. These results
793	>	further emphasize the importance of reparameterization in order to
794	>	model the density properly under different simulation conditions.
795	>	Again, these changes have only a minor effect on the melting point,
796	>	which observed at 245~K for SSD/RF, is identical to SSD and only 5~K
797	>	lower than SSD1 with a reaction field. Additionally, the difference in
798	>	density maxima is not as extreme, with SSD/RF showing a density
799	>	maximum at 255~K, fairly close to the density maxima of 260~K and
800	>	265~K, shown by SSD and SSD1 respectively.
801
802	+	%\begin{figure}
803	+	%\begin{center}
804	+	%\epsfxsize=6in
805	+	%\epsfbox{ssdeDiffuse.epsi}
806	+	%\caption{The diffusion constants calculated from SSD/E and
807	+	%SSD1 (both without a reaction field) along with experimental results
808	+	%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
809	+	%performed at the average densities observed in the 1 atm NPT
810	+	%simulations for the respective models. SSD/E is slightly more mobile
811	+	%than experiment at all of the temperatures, but it is closer to
812	+	%experiment at biologically relevant temperatures than SSD1 without a
813	+	%long-range correction.}
814	+	%\label{ssdediffuse}
815	+	%\end{center}
816	+	%\end{figure}
817	+
818		The reparameterization of the SSD water model, both for use with and
819		without an applied long-range correction, brought the densities up to
820		what is expected for simulating liquid water. In addition to improving
821	<	the densities, it is important that particle transport be maintained
822	<	or improved. Figure \ref{ssdediffuse} compares the temperature
823	<	dependence of the diffusion constant of SSD/E to SSD without an active
824	<	reaction field, both at the densities calculated at 1 atm and at the
825	<	experimentally calculated densities for super-cooled and liquid
826	<	water. In the upper plot, the diffusion constant for SSD/E is
827	<	consistently a little faster than experiment, while SSD starts off
828	<	slower than experiment and crosses to merge with SSD/E at high
829	<	temperatures. Both models follow the experimental trend well, but
830	<	diffuse too rapidly at higher temperatures. This abnormally fast
831	<	diffusion is caused by the decreased system density. Since the
832	<	densities of SSD/E don't deviate as much from experiment as those of
833	<	SSD, it follows the experimental trend more closely. This observation
834	<	is backed up by looking at the lower plot. The diffusion constants for
835	<	SSD/E track with the experimental values while SSD deviates on the low
836	<	side of the trend with increasing temperature. This is again a product
837	<	of SSD/E having densities closer to experiment, and not deviating to
838	<	lower densities with increasing temperature as rapidly.
821	>	the densities, it is important that the diffusive behavior of SSD be
822	>	maintained or improved. Figure \ref{ssdediffuse} compares the
823	>	temperature dependence of the diffusion constant of SSD/E to SSD1
824	>	without an active reaction field at the densities calculated from
825	>	their respective NPT simulations at 1 atm. The diffusion constant for
826	>	SSD/E is consistently higher than experiment, while SSD1 remains lower
827	>	than experiment until relatively high temperatures (around 360
828	>	K). Both models follow the shape of the experimental curve well below
829	>	300~K but tend to diffuse too rapidly at higher temperatures, as seen
830	>	in SSD1's crossing above 360~K.  This increasing diffusion relative to
831	>	the experimental values is caused by the rapidly decreasing system
832	>	density with increasing temperature.  Both SSD1 and SSD/E show this
833	>	deviation in particle mobility, but this trend has different
834	>	implications on the diffusive behavior of the models.  While SSD1
835	>	shows more experimentally accurate diffusive behavior in the high
836	>	temperature regimes, SSD/E shows more accurate behavior in the
837	>	supercooled and biologically relevant temperature ranges.  Thus, the
838	>	changes made to improve the liquid structure may have had an adverse
839	>	affect on the density maximum, but they improve the transport behavior
840	>	of SSD/E relative to SSD1 under the most commonly simulated
841	>	conditions.
842
843	<	\begin{figure}
844	<	\includegraphics[width=85mm]{ssdediffuse.epsi}
845	<	\caption{Plots of the diffusion constants calculated from SSD/E and SSD,
846	<	both without a reaction field along with experimental results from
847	<	Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
848	<	upper plot is at densities calculated from the NPT simulations at a
849	<	pressure of 1 atm, while the lower plot is at the experimentally
850	<	calculated densities.}
851	<	\label{ssdediffuse}
852	<	\end{figure}
843	>	%\begin{figure}
844	>	%\begin{center}
845	>	%\epsfxsize=6in
846	>	%\epsfbox{ssdrfDiffuse.epsi}
847	>	%\caption{The diffusion constants calculated from SSD/RF and
848	>	%SSD1 (both with an active reaction field) along with
849	>	%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
850	>	%NVE calculations were performed at the average densities observed in
851	>	%the 1 atm NPT simulations for both of the models. SSD/RF
852	>	%simulates the diffusion of water throughout this temperature range
853	>	%very well. The rapidly increasing diffusion constants at high
854	>	%temperatures for both models can be attributed to lower calculated
855	>	%densities than those observed in experiment.}
856	>	%\label{ssdrfdiffuse}
857	>	%\end{center}
858	>	%\end{figure}
859
760	–	\begin{figure}
761	–	\includegraphics[width=85mm]{ssdrfdiffuse.epsi}
762	–	\caption{Plots of the diffusion constants calculated from SSD/RF and SSD,
763	–	both with an active reaction field along with experimental results
764	–	from Gillen \emph{et al.}\cite{Gillen72} and Mills\cite{Mills73}. The
765	–	upper plot is at densities calculated from the NPT simulations at a
766	–	pressure of 1 atm, while the lower plot is at the experimentally
767	–	calculated densities.}
768	–	\label{ssdrfdiffuse}
769	–	\end{figure}
770	–
860		In figure \ref{ssdrfdiffuse}, the diffusion constants for SSD/RF are
861	<	compared with SSD with an active reaction field. In the upper plot,
862	<	SSD/RF tracks with the experimental results incredibly well, identical
863	<	within error throughout the temperature range and only showing a
864	<	slight increasing trend at higher temperatures. SSD also tracks
865	<	experiment well, only it tends to diffuse a little more slowly at low
866	<	temperatures and deviates to diffuse too rapidly at high
867	<	temperatures. As was stated in the SSD/E comparisons, this deviation
868	<	away from the ideal trend is due to a rapid decrease in density at
869	<	higher temperatures. SSD/RF doesn't suffer from this problem as much
870	<	as SSD, because the calculated densities are more true to
871	<	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.
861	>	compared to SSD1 with an active reaction field. Note that SSD/RF
862	>	tracks the experimental results quantitatively, identical within error
863	>	throughout most of the temperature range shown and exhibiting only a
864	>	slight increasing trend at higher temperatures. SSD1 tends to diffuse
865	>	more slowly at low temperatures and deviates to diffuse too rapidly at
866	>	temperatures greater than 330~K.  As stated above, this deviation away
867	>	from the ideal trend is due to a rapid decrease in density at higher
868	>	temperatures. SSD/RF does not suffer from this problem as much as SSD1
869	>	because the calculated densities are closer to the experimental
870	>	values. These results again emphasize the importance of careful
871	>	reparameterization when using an altered long-range correction.
872
873	+	\begin{table}
874	+	\begin{minipage}{\linewidth}
875	+	\renewcommand{\thefootnote}{\thempfootnote}
876	+	\begin{center}
877	+	\caption{Properties of the single-point water models compared with
878	+	experimental data at ambient conditions. Deviations of the of the
879	+	averages are given in parentheses.}
880	+	\begin{tabular}{ l  c  c  c  c  c }
881	+	\hline \\[-3mm]
882	+	\ \ \ \ \ \  & \ \ \ SSD1 \ \ \ & \ \ SSD/E \ \ \ & \ \ SSD1 (RF) \ \
883	+	\ & \ \ SSD/RF \ \ \ & \ \ Expt. \\
884	+	\hline \\[-3mm]
885	+	\ \ $\rho$ (g/cm$^3$) & 0.999(0.001) & 0.996(0.001) & 0.972(0.002) & 0.997(0.001) & 0.997 \\
886	+	\ \ $C_p$ (cal/mol K) & 28.80(0.11) & 25.45(0.09) & 28.28(0.06) & 23.83(0.16) & 17.98 \\
887	+	\ \ $D$ ($10^{-5}$ cm$^2$/s) & 1.78(0.7) & 2.51(0.18) & 2.00(0.17) & 2.32(0.06) & 2.299\cite{Mills73} \\
888	+	\ \ Coordination Number ($n_C$) & 3.9 & 4.3 & 3.8 & 4.4 &
889	+	4.7\footnote{Calculated by integrating $g_{\text{OO}}(r)$ in
890	+	Ref. \onlinecite{Head-Gordon00_1}} \\
891	+	\ \ H-bonds per particle ($n_H$) & 3.7 & 3.6 & 3.7 & 3.7 &
892	+	3.5\footnote{Calculated by integrating $g_{\text{OH}}(r)$ in
893	+	Ref. \onlinecite{Soper86}}  \\
894	+	\ \ $\tau_1$ (ps) & 10.9(0.6) & 7.3(0.4) & 7.5(0.7) & 7.2(0.4) & 5.7\footnote{Calculated for 298 K from data in Ref. \onlinecite{Eisenberg69}} \\
895	+	\ \ $\tau_2$ (ps) & 4.7(0.4) & 3.1(0.2) & 3.5(0.3) & 3.2(0.2) & 2.3\footnote{Calculated for 298 K from data in
896	+	Ref. \onlinecite{Krynicki66}}
897	+	\end{tabular}
898	+	\label{liquidproperties}
899	+	\end{center}
900	+	\end{minipage}
901	+	\end{table}
902	+
903	+	Table \ref{liquidproperties} gives a synopsis of the liquid state
904	+	properties of the water models compared in this study along with the
905	+	experimental values for liquid water at ambient conditions. The
906	+	coordination number ($n_C$) and number of hydrogen bonds per particle
907	+	($n_H$) were calculated by integrating the following relations:
908	+	\begin{equation}
909	+	n_C = 4\pi\rho_{\text{OO}}\int_{0}^{a}r^2\text{g}_{\text{OO}}(r)dr,
910	+	\end{equation}
911	+	\begin{equation}
912	+	n_H = 4\pi\rho_{\text{OH}}\int_{0}^{b}r^2\text{g}_{\text{OH}}(r)dr,
913	+	\end{equation}
914	+	where $\rho$ is the number density of the specified pair interactions,
915	+	$a$ and $b$ are the radial locations of the minima following the first
916	+	peak in g$_\text{OO}(r)$ or g$_\text{OH}(r)$ respectively. The number
917	+	of hydrogen bonds stays relatively constant across all of the models,
918	+	but the coordination numbers of SSD/E and SSD/RF show an
919	+	improvement over SSD1.  This improvement is primarily due to
920	+	extension of the first solvation shell in the new parameter sets.
921	+	Because $n_H$ and $n_C$ are nearly identical in SSD1, it appears
922	+	that all molecules in the first solvation shell are involved in
923	+	hydrogen bonds.  Since $n_H$ and $n_C$ differ in the newly
924	+	parameterized models, the orientations in the first solvation shell
925	+	are a bit more ``fluid''.  Therefore SSD1 overstructures the
926	+	first solvation shell and our reparameterizations have returned this
927	+	shell to more realistic liquid-like behavior.
928	+
929	+	The time constants for the orientational autocorrelation functions
930	+	are also displayed in Table \ref{liquidproperties}. The dipolar
931	+	orientational time correlation functions ($C_{l}$) are described
932	+	by:
933	+	\begin{equation}
934	+	C_{l}(t) = \langle P_l[\hat{\mathbf{u}}_j(0)\cdot\hat{\mathbf{u}}_j(t)]\rangle,
935	+	\end{equation}
936	+	where $P_l$ are Legendre polynomials of order $l$ and
937	+	$\hat{\mathbf{u}}_j$ is the unit vector pointing along the molecular
938	+	dipole.\cite{Rahman71} From these correlation functions, the
939	+	orientational relaxation time of the dipole vector can be calculated
940	+	from an exponential fit in the long-time regime ($t >
941	+	\tau_l$).\cite{Rothschild84} Calculation of these time constants were
942	+	averaged over five detailed NVE simulations performed at the ambient
943	+	conditions for each of the respective models. It should be noted that
944	+	the commonly cited value of 1.9 ps for $\tau_2$ was determined from
945	+	the NMR data in Ref. \onlinecite{Krynicki66} at a temperature near
946	+	34$^\circ$C.\cite{Rahman71} Because of the strong temperature
947	+	dependence of $\tau_2$, it is necessary to recalculate it at 298~K to
948	+	make proper comparisons. The value shown in Table
949	+	\ref{liquidproperties} was calculated from the same NMR data in the
950	+	fashion described in Ref. \onlinecite{Krynicki66}. Similarly, $\tau_1$ was
951	+	recomputed for 298~K from the data in Ref. \onlinecite{Eisenberg69}.
952	+	Again, SSD/E and SSD/RF show improved behavior over SSD1, both with
953	+	and without an active reaction field. Turning on the reaction field
954	+	leads to much improved time constants for SSD1; however, these results
955	+	also include a corresponding decrease in system density.
956	+	Orientational relaxation times published in the original SSD dynamics
957	+	papers are smaller than the values observed here, and this difference
958	+	can be attributed to the use of the Ewald sum.\cite{Ichiye99}
959	+
960		\subsection{Additional Observations}
961
962	<	While performing the melting sequences of SSD/E, some interesting
963	<	observations were made. After melting at 230 K, two of the systems
964	<	underwent crystallization events near 245 K. As the heating process
965	<	continued, the two systems remained crystalline until finally melting
966	<	between 320 and 330 K. These simulations were excluded from the data
967	<	set shown in figure \ref{ssdedense} and replaced with two additional
968	<	melting sequences that did not undergo this anomalous phase
969	<	transition, while this crystallization event was investigated
970	<	separately.
962	>	%\begin{figure}
963	>	%\begin{center}
964	>	%\epsfxsize=6in
965	>	%\epsfbox{icei_bw.eps}
966	>	%\caption{The most stable crystal structure assumed by the SSD family
967	>	%of water models.  We refer to this structure as Ice-{\it i} to
968	>	%indicate its origins in computer simulation.  This image was taken of
969	>	%the (001) face of the crystal.}
970	>	%\label{weirdice}
971	>	%\end{center}
972	>	%\end{figure}
973
974	+	While performing a series of melting simulations on an early iteration
975	+	of SSD/E not discussed in this paper, we observed
976	+	recrystallization into a novel structure not previously known for
977	+	water.  After melting at 235~K, two of five systems underwent
978	+	crystallization events near 245~K.  The two systems remained
979	+	crystalline up to 320 and 330~K, respectively.  The crystal exhibits
980	+	an expanded zeolite-like structure that does not correspond to any
981	+	known form of ice.  This appears to be an artifact of the point
982	+	dipolar models, so to distinguish it from the experimentally observed
983	+	forms of ice, we have denoted the structure
984	+	Ice-$\sqrt{\smash[b]{-\text{I}}}$ (Ice-{\it i}).  A large enough
985	+	portion of the sample crystallized that we have been able to obtain a
986	+	near ideal crystal structure of Ice-{\it i}. Figure \ref{weirdice}
987	+	shows the repeating crystal structure of a typical crystal at 5
988	+	K. Each water molecule is hydrogen bonded to four others; however, the
989	+	hydrogen bonds are bent rather than perfectly straight. This results
990	+	in a skewed tetrahedral geometry about the central molecule.  In
991	+	figure \ref{isosurface}, it is apparent that these flexed hydrogen
992	+	bonds are allowed due to the conical shape of the attractive regions,
993	+	with the greatest attraction along the direct hydrogen bond
994	+	configuration. Though not ideal, these flexed hydrogen bonds are
995	+	favorable enough to stabilize an entire crystal generated around them.
996	+
997	+	Initial simulations indicated that Ice-{\it i} is the preferred ice
998	+	structure for at least the SSD/E model. To verify this, a comparison
999	+	was made between near ideal crystals of ice~$I_h$, ice~$I_c$, and
1000	+	Ice-{\it i} at constant pressure with SSD/E, SSD/RF, and
1001	+	SSD1. Near-ideal versions of the three types of crystals were cooled
1002	+	to 1 K, and enthalpies of formation of each were compared using all
1003	+	three water models.  Enthalpies were estimated from the
1004	+	isobaric-isothermal simulations using $H=U+P_{\text ext}V$ where
1005	+	$P_{\text ext}$ is the applied pressure.  A constant value of -60.158
1006	+	kcal / mol has been added to place our zero for the enthalpies of
1007	+	formation for these systems at the traditional state (elemental forms
1008	+	at standard temperature and pressure).  With every model in the SSD
1009	+	family, Ice-{\it i} had the lowest calculated enthalpy of formation.
1010	+
1011	+	\begin{table}
1012	+	\begin{center}
1013	+	\caption{Enthalpies of Formation (in kcal / mol) of the three crystal
1014	+	structures (at 1 K) exhibited by the SSD family of water models}
1015	+	\begin{tabular}{ l  c  c  c  }
1016	+	\hline \\[-3mm]
1017	+	\ \ \ Water Model \ \ \  & \ \ \ Ice-$I_h$ \ \ \ & \ \ \ Ice-$I_c$ \ \ \  &
1018	+	\ \ \ \ Ice-{\it i} \\
1019	+	\hline \\[-3mm]
1020	+	\ \ \ SSD/E & -72.444 & -72.450 & -73.748 \\
1021	+	\ \ \ SSD/RF & -73.093 & -73.075 & -74.180 \\
1022	+	\ \ \ SSD1 & -72.654 & -72.569 & -73.575 \\
1023	+	\ \ \ SSD1 (RF) & -72.662 & -72.569 & -73.292 \\
1024	+	\end{tabular}
1025	+	\label{iceenthalpy}
1026	+	\end{center}
1027	+	\end{table}
1028	+
1029	+	In addition to these energetic comparisons, melting simulations were
1030	+	performed with Ice-{\it i} as the initial configuration using SSD/E,
1031	+	SSD/RF, and SSD1 both with and without a reaction field. The melting
1032	+	transitions for both SSD/E and SSD1 without reaction field occurred at
1033	+	temperature in excess of 375~K.  SSD/RF and SSD1 with a reaction field
1034	+	showed more reasonable melting transitions near 325~K.  These melting
1035	+	point observations clearly show that all of the SSD-derived models
1036	+	prefer the ice-{\it i} structure.
1037	+
1038	+	\section{Conclusions}
1039	+
1040	+	The density maximum and temperature dependence of the self-diffusion
1041	+	constant were studied for the SSD water model, both with and
1042	+	without the use of reaction field, via a series of NPT and NVE
1043	+	simulations. The constant pressure simulations showed a density
1044	+	maximum near 260 K. In most cases, the calculated densities were
1045	+	significantly lower than the densities obtained from other water
1046	+	models (and experiment). Analysis of self-diffusion showed SSD
1047	+	to capture the transport properties of water well in both the liquid
1048	+	and supercooled liquid regimes.
1049	+
1050	+	In order to correct the density behavior, the original SSD model was
1051	+	reparameterized for use both with and without a reaction field (SSD/RF
1052	+	and SSD/E), and comparisons were made with SSD1, Ichiye's density
1053	+	corrected version of SSD. Both models improve the liquid structure,
1054	+	densities, and diffusive properties under their respective simulation
1055	+	conditions, indicating the necessity of reparameterization when
1056	+	changing the method of calculating long-range electrostatic
1057	+	interactions.  In general, however, these simple water models are
1058	+	excellent choices for representing explicit water in large scale
1059	+	simulations of biochemical systems.
1060	+
1061	+	The existence of a novel low-density ice structure that is preferred
1062	+	by the SSD family of water models is somewhat troubling, since
1063	+	liquid simulations on this family of water models at room temperature
1064	+	are effectively simulations of supercooled or metastable liquids.  One
1065	+	way to destabilize this unphysical ice structure would be to make the
1066	+	range of angles preferred by the attractive part of the sticky
1067	+	potential much narrower.  This would require extensive
1068	+	reparameterization to maintain the same level of agreement with the
1069	+	experiments.
1070	+
1071	+	Additionally, our initial calculations show that the Ice-{\it i}
1072	+	structure may also be a preferred crystal structure for at least one
1073	+	other popular multi-point water model (TIP3P), and that much of the
1074	+	simulation work being done using this popular model could also be at
1075	+	risk for crystallization into this unphysical structure.  A future
1076	+	publication will detail the relative stability of the known ice
1077	+	structures for a wide range of popular water models.
1078	+
1079	+	\section{Acknowledgments}
1080	+	Support for this project was provided by the National Science
1081	+	Foundation under grant CHE-0134881. Computation time was provided by
1082	+	the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
1083	+	DMR-0079647.
1084	+
1085	+	\newpage
1086	+
1087	+	\bibliographystyle{jcp}
1088	+	\bibliography{nptSSD}
1089	+
1090	+	\newpage
1091	+
1092	+	\begin{list}
1093	+	{Figure \arabic{captions}: }{\usecounter{captions}
1094	+	\setlength{\rightmargin}{\leftmargin}}
1095	+
1096	+	\item Energy conservation using both quaternion-based integration and
1097	+	the {\sc dlm} method with increasing time step. The larger time step
1098	+	plots are shifted from the true energy baseline (that of $\Delta t$ =
1099	+	0.1~fs) for clarity.
1100	+
1101	+	\item Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
1102	+	TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E
1103	+	[Ref. \onlinecite{Clancy94}], SSD without Reaction Field, SSD, and
1104	+	experiment [Ref. \onlinecite{CRC80}]. The arrows indicate the change
1105	+	in densities observed when turning off the reaction field. The the
1106	+	lower than expected densities for the SSD model were what prompted the
1107	+	original reparameterization of SSD1 [Ref. \onlinecite{Ichiye03}].
1108	+
1109	+	\item Average self-diffusion constant as a function of temperature for
1110	+	SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
1111	+	[Ref. \onlinecite{Jorgensen01}] compared with experimental data
1112	+	[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three
1113	+	water models shown, SSD has the least deviation from the experimental
1114	+	values. The rapidly increasing diffusion constants for TIP5P and SSD
1115	+	correspond to significant decreases in density at the higher
1116	+	temperatures.
1117	+
1118	+	\item An illustration of angles involved in the correlations observed in
1119	+	Fig. \ref{contour}.
1120	+
1121	+	\item Contour plots of 2D angular pair correlation functions for
1122	+	512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
1123	+	signify regions of enhanced density while light areas signify
1124	+	depletion relative to the bulk density. White areas have pair
1125	+	correlation values below 0.5 and black areas have values above 1.5.
1126	+
1127	+	\item Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with
1128	+	SSD/E and SSD1 without reaction field (top), as well as SSD/RF and
1129	+	SSD1 with reaction field turned on (bottom). The insets show the
1130	+	respective first peaks in detail. Note how the changes in parameters
1131	+	have lowered and broadened the first peak of SSD/E and SSD/RF.
1132	+
1133	+	\item Positive and negative isosurfaces of the sticky potential for
1134	+	SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
1135	+	correspond to the tetrahedral attractive component, and darker areas
1136	+	correspond to the dipolar repulsive component.
1137	+
1138	+	\item Comparison of densities calculated with SSD/E to
1139	+	SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
1140	+	TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
1141	+	experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
1142	+	300 K with error bars included to clarify this region of
1143	+	interest. Note that both SSD1 and SSD/E show good agreement with
1144	+	experiment when the long-range correction is neglected.
1145	+
1146	+	\item Comparison of densities calculated with SSD/RF to
1147	+	SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
1148	+	TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
1149	+	experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
1150	+	reparameterization when utilizing a reaction field long-ranged
1151	+	correction - SSD/RF provides significantly more accurate
1152	+	densities than SSD1 when performing room temperature
1153	+	simulations.
1154	+
1155	+	\item The diffusion constants calculated from SSD/E and
1156	+	SSD1 (both without a reaction field) along with experimental results
1157	+	[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
1158	+	performed at the average densities observed in the 1 atm NPT
1159	+	simulations for the respective models. SSD/E is slightly more mobile
1160	+	than experiment at all of the temperatures, but it is closer to
1161	+	experiment at biologically relevant temperatures than SSD1 without a
1162	+	long-range correction.
1163	+
1164	+	\item The diffusion constants calculated from SSD/RF and
1165	+	SSD1 (both with an active reaction field) along with
1166	+	experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
1167	+	NVE calculations were performed at the average densities observed in
1168	+	the 1 atm NPT simulations for both of the models. SSD/RF
1169	+	simulates the diffusion of water throughout this temperature range
1170	+	very well. The rapidly increasing diffusion constants at high
1171	+	temperatures for both models can be attributed to lower calculated
1172	+	densities than those observed in experiment.
1173	+
1174	+	\item The most stable crystal structure assumed by the SSD family
1175	+	of water models.  We refer to this structure as Ice-{\it i} to
1176	+	indicate its origins in computer simulation.  This image was taken of
1177	+	the (001) face of the crystal.
1178	+	\end{list}
1179	+
1180	+	\newpage
1181	+
1182		\begin{figure}
1183	<	\includegraphics[width=85mm]{povIce.ps}
1184	<	\caption{Crystal structure of an ice 0 lattice shown from the (001) face.}
1185	<	\label{weirdice}
1183	>	\begin{center}
1184	>	\epsfxsize=6in
1185	>	\epsfbox{timeStep.epsi}
1186	>	%\caption{Energy conservation using both quaternion-based integration and
1187	>	%the {\sc dlm} method with increasing time step. The larger time step
1188	>	%plots are shifted from the true energy baseline (that of $\Delta t$ =
1189	>	%0.1~fs) for clarity.}
1190	>	\label{timestep}
1191	>	\end{center}
1192		\end{figure}
1193
1194	<	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.
1194	>	\newpage
1195
1196		\begin{figure}
1197	<	\includegraphics[width=65mm]{ice0cell.eps}
1198	<	\caption{Simple unit cell for constructing ice 0. In this cell, $c$ is
1199	<	equal to $0.4714\times a$, and a typical value for $a$ is 8.25 \AA.}
1200	<	\label{unitcell}
1197	>	\begin{center}
1198	>	\epsfxsize=6in
1199	>	\epsfbox{denseSSDnew.eps}
1200	>	%\caption{Density versus temperature for TIP4P [Ref. \onlinecite{Jorgensen98b}],
1201	>	% TIP3P [Ref. \onlinecite{Jorgensen98b}], SPC/E [Ref. \onlinecite{Clancy94}], SSD
1202	>	% without Reaction Field, SSD, and experiment [Ref. \onlinecite{CRC80}]. The
1203	>	% arrows indicate the change in densities observed when turning off the
1204	>	% reaction field. The the lower than expected densities for the SSD
1205	>	% model were what prompted the original reparameterization of SSD1
1206	>	% [Ref. \onlinecite{Ichiye03}].}
1207	>	\label{dense1}
1208	>	\end{center}
1209		\end{figure}
1210
1211	<	The initial simulations indicated that ice 0 is the preferred ice
835	<	structure for at least SSD/E. To verify this, a comparison was made
836	<	between near ideal crystals of ice $I_h$, ice $I_c$, and ice 0 at
837	<	constant pressure with SSD/E, SSD/RF, and SSD. Near ideal versions of
838	<	the three types of crystals were cooled to ~1 K, and the potential
839	<	energies of each were compared using all three water models. With
840	<	every water model, ice 0 turned out to have the lowest potential
841	<	energy: 5\% lower than $I_h$ with SSD, 6.5\% lower with SSD/E, and
842	<	7.5\% lower with SSD/RF. In all three of these water models, ice $I_c$
843	<	was observed to be ~2\% less stable than ice $I_h$. In addition to
844	<	having the lowest potential energy, ice 0 was the most expanded of the
845	<	three ice crystals, ~5\% less dense than ice $I_h$ with all of the
846	<	water models. In all three water models, ice $I_c$ was observed to be
847	<	~2\% more dense than ice $I_h$.
1211	>	\newpage
1212
1213	<	In addition to the low temperature comparisons, melting sequences were
1214	<	performed with ice 0 as the initial configuration using SSD/E, SSD/RF,
1215	<	and SSD both with and without a reaction field. The melting
1216	<	transitions for both SSD/E and SSD without a reaction field occurred
1217	<	at temperature in excess of 375 K. SSD/RF and SSD with a reaction
1218	<	field had more reasonable melting transitions, down near 325 K. These
1219	<	melting point observations emphasize how preferred this crystal
1220	<	structure is over the most common types of ice when using these single
1221	<	point water models.
1213	>	\begin{figure}
1214	>	\begin{center}
1215	>	\epsfxsize=6in
1216	>	\epsfbox{betterDiffuse.epsi}
1217	>	%\caption{Average self-diffusion constant as a function of temperature for
1218	>	%SSD, SPC/E [Ref. \onlinecite{Clancy94}], and TIP5P
1219	>	%[Ref. \onlinecite{Jorgensen01}] compared with experimental data
1220	>	%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. Of the three water models
1221	>	%shown, SSD has the least deviation from the experimental values. The
1222	>	%rapidly increasing diffusion constants for TIP5P and SSD correspond to
1223	>	%significant decreases in density at the higher temperatures.}
1224	>	\label{diffuse}
1225	>	\end{center}
1226	>	\end{figure}
1227
1228	<	Recognizing that the above tests show ice 0 to be both the most stable
860	<	and lowest density crystal structure for these single point water
861	<	models, it is interesting to speculate on the favorability of this
862	<	crystal structure with the different charge based models. As a quick
863	<	test, these 3 crystal types were converted from SSD type particles to
864	<	TIP3P waters and read into CHARMM.\cite{Karplus83} Identical energy
865	<	minimizations were performed on all of these crystals to compare the
866	<	system energies. Again, ice 0 was observed to have the lowest total
867	<	system energy. The total energy of ice 0 was ~2\% lower than ice
868	<	$I_h$, which was in turn ~3\% lower than ice $I_c$. From these initial
869	<	results, we would not be surprised if results from the other common
870	<	water models show ice 0 to be the lowest energy crystal structure. A
871	<	continuation on work studing ice 0 with multipoint water models will
872	<	be published in a coming article.
1228	>	\newpage
1229
1230	<	\section{Conclusions}
1231	<	The density maximum and temperature dependent transport for the SSD
1232	<	water model, both with and without the use of reaction field, were
1233	<	studied via a series of NPT and NVE simulations. The constant pressure
1234	<	simulations of the melting of both $I_h$ and $I_c$ ice showed a
1235	<	density maximum near 260 K. In most cases, the calculated densities
1236	<	were significantly lower than the densities calculated in simulations
1237	<	of other water models. Analysis of particle diffusion showed SSD to
882	<	capture the transport properties of experimental very well in both the
883	<	normal and super-cooled liquid regimes. In order to correct the
884	<	density behavior, SSD was reparameterized for use both with and
885	<	without a long-range interaction correction, SSD/RF and SSD/E
886	<	respectively. In addition to correcting the abnormally low densities,
887	<	these new versions were show to maintain or improve upon the transport
888	<	and structural features of the original water model, all while
889	<	maintaining the fast performance of the SSD water model. This work
890	<	shows these simple water models, and in particular SSD/E and SSD/RF,
891	<	to be excellent choices to represent explicit water in future
892	<	simulations of biochemical systems.
1230	>	\begin{figure}
1231	>	\begin{center}
1232	>	\epsfxsize=6in
1233	>	\epsfbox{corrDiag.eps}
1234	>	%\caption{An illustration of angles involved in the correlations observed in Fig. \ref{contour}.}
1235	>	\label{corrAngle}
1236	>	\end{center}
1237	>	\end{figure}
1238
1239	<	\section{Acknowledgments}
895	<	The authors would like to thank the National Science Foundation for
896	<	funding under grant CHE-0134881. Computation time was provided by the
897	<	Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR
898	<	00 79647.
1239	>	\newpage
1240
1241	<	\bibliographystyle{jcp}
1241	>	\begin{figure}
1242	>	\begin{center}
1243	>	\epsfxsize=6in
1244	>	\epsfbox{fullContours.eps}
1245	>	%\caption{Contour plots of 2D angular pair correlation functions for
1246	>	%512 SSD molecules at 100~K (A \& B) and 300~K (C \& D). Dark areas
1247	>	%signify regions of enhanced density while light areas signify
1248	>	%depletion relative to the bulk density. White areas have pair
1249	>	%correlation values below 0.5 and black areas have values above 1.5.}
1250	>	\label{contour}
1251	>	\end{center}
1252	>	\end{figure}
1253
1254	<	\bibliography{nptSSD}
1254	>	\newpage
1255
1256	<	%\pagebreak
1256	>	\begin{figure}
1257	>	\begin{center}
1258	>	\epsfxsize=6in
1259	>	\epsfbox{GofRCompare.epsi}
1260	>	%\caption{Plots comparing experiment [Ref. \onlinecite{Head-Gordon00_1}] with
1261	>	%SSD/E and SSD1 without reaction field (top), as well as
1262	>	%SSD/RF and SSD1 with reaction field turned on
1263	>	%(bottom). The insets show the respective first peaks in detail. Note
1264	>	%how the changes in parameters have lowered and broadened the first
1265	>	%peak of SSD/E and SSD/RF.}
1266	>	\label{grcompare}
1267	>	\end{center}
1268	>	\end{figure}
1269
1270	+	\newpage
1271	+
1272	+	\begin{figure}
1273	+	\begin{center}
1274	+	\epsfxsize=7in
1275	+	\epsfbox{dualsticky_bw.eps}
1276	+	%\caption{Positive and negative isosurfaces of the sticky potential for
1277	+	%SSD1 (left) and SSD/E \& SSD/RF (right). Light areas
1278	+	%correspond to the tetrahedral attractive component, and darker areas
1279	+	%correspond to the dipolar repulsive component.}
1280	+	\label{isosurface}
1281	+	\end{center}
1282	+	\end{figure}
1283	+
1284	+	\newpage
1285	+
1286	+	\begin{figure}
1287	+	\begin{center}
1288	+	\epsfxsize=6in
1289	+	\epsfbox{ssdeDense.epsi}
1290	+	%\caption{Comparison of densities calculated with SSD/E to
1291	+	%SSD1 without a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
1292	+	%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}] and
1293	+	%experiment [Ref. \onlinecite{CRC80}]. The window shows a expansion around
1294	+	%300 K with error bars included to clarify this region of
1295	+	%interest. Note that both SSD1 and SSD/E show good agreement with
1296	+	%experiment when the long-range correction is neglected.}
1297	+	\label{ssdedense}
1298	+	\end{center}
1299	+	\end{figure}
1300	+
1301	+	\newpage
1302	+
1303	+	\begin{figure}
1304	+	\begin{center}
1305	+	\epsfxsize=6in
1306	+	\epsfbox{ssdrfDense.epsi}
1307	+	%\caption{Comparison of densities calculated with SSD/RF to
1308	+	%SSD1 with a reaction field, TIP3P [Ref. \onlinecite{Jorgensen98b}],
1309	+	%TIP5P [Ref. \onlinecite{Jorgensen00}], SPC/E [Ref. \onlinecite{Clancy94}], and
1310	+	%experiment [Ref. \onlinecite{CRC80}]. The inset shows the necessity of
1311	+	%reparameterization when utilizing a reaction field long-ranged
1312	+	%correction - SSD/RF provides significantly more accurate
1313	+	%densities than SSD1 when performing room temperature
1314	+	%simulations.}
1315	+	\label{ssdrfdense}
1316	+	\end{center}
1317	+	\end{figure}
1318	+
1319	+	\newpage
1320	+
1321	+	\begin{figure}
1322	+	\begin{center}
1323	+	\epsfxsize=6in
1324	+	\epsfbox{ssdeDiffuse.epsi}
1325	+	%\caption{The diffusion constants calculated from SSD/E and
1326	+	%SSD1 (both without a reaction field) along with experimental results
1327	+	%[Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The NVE calculations were
1328	+	%performed at the average densities observed in the 1 atm NPT
1329	+	%simulations for the respective models. SSD/E is slightly more mobile
1330	+	%than experiment at all of the temperatures, but it is closer to
1331	+	%experiment at biologically relevant temperatures than SSD1 without a
1332	+	%long-range correction.}
1333	+	\label{ssdediffuse}
1334	+	\end{center}
1335	+	\end{figure}
1336	+
1337	+	\newpage
1338	+
1339	+	\begin{figure}
1340	+	\begin{center}
1341	+	\epsfxsize=6in
1342	+	\epsfbox{ssdrfDiffuse.epsi}
1343	+	%\caption{The diffusion constants calculated from SSD/RF and
1344	+	%SSD1 (both with an active reaction field) along with
1345	+	%experimental results [Refs. \onlinecite{Gillen72} and \onlinecite{Holz00}]. The
1346	+	%NVE calculations were performed at the average densities observed in
1347	+	%the 1 atm NPT simulations for both of the models. SSD/RF
1348	+	%simulates the diffusion of water throughout this temperature range
1349	+	%very well. The rapidly increasing diffusion constants at high
1350	+	%temperatures for both models can be attributed to lower calculated
1351	+	%densities than those observed in experiment.}
1352	+	\label{ssdrfdiffuse}
1353	+	\end{center}
1354	+	\end{figure}
1355	+
1356	+	\newpage
1357	+
1358	+	\begin{figure}
1359	+	\begin{center}
1360	+	\epsfxsize=6in
1361	+	\epsfbox{icei_bw.eps}
1362	+	%\caption{The most stable crystal structure assumed by the SSD family
1363	+	%of water models.  We refer to this structure as Ice-{\it i} to
1364	+	%indicate its origins in computer simulation.  This image was taken of
1365	+	%the (001) face of the crystal.}
1366	+	\label{weirdice}
1367	+	\end{center}
1368	+	\end{figure}
1369	+
1370		\end{document}

Diff Legend

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