ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/ssdePaper/nptSSD.tex

Comparing trunk/ssdePaper/nptSSD.tex (file contents):
Revision 759 by chrisfen, Wed Sep 10 22:42:57 2003 UTC vs.
Revision 863 by chrisfen, Thu Nov 13 15:55:20 2003 UTC

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

Diff Legend

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