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 861 by chrisfen, Wed Nov 12 13:37:15 2003 UTC vs.
Revision 1017 by chrisfen, Wed Feb 4 18:51:43 2004 UTC

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

Diff Legend

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