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 759 by chrisfen, Wed Sep 10 22:42:57 2003 UTC vs.
Revision 1024 by chrisfen, Wed Feb 4 22:42:52 2004 UTC

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

Diff Legend

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