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 777 by chrisfen, Fri Sep 19 19:29:24 2003 UTC vs.
Revision 1036 by gezelter, Fri Feb 6 21:43:00 2004 UTC

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

Diff Legend

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