ViewVC Help

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

root/group/trunk/fennellDissertation/dissertation.tex

Comparing trunk/fennellDissertation/dissertation.tex (file contents):
Revision 2927 by chrisfen, Sat Jul 8 13:13:27 2006 UTC vs.
Revision 2973 by chrisfen, Fri Aug 18 15:49:29 2006 UTC

#	Line 1 | Line 1
1	<	\documentclass[12pt]{ndthesis}
1	>	\documentclass[11pt]{ndthesis}
2
3		% some packages for things like equations and graphics
4	+	\usepackage[tbtags]{amsmath}
5		\usepackage{amsmath,bm}
6		\usepackage{amssymb}
7		\usepackage{mathrsfs}
8		\usepackage{tabularx}
9		\usepackage{graphicx}
10		\usepackage{booktabs}
11	+	\usepackage{cite}
12	+	\usepackage{enumitem}
13
14		\begin{document}
15
#	Line 39 | Line 42 | STUDY OF WATER}
42		\mainmatter
43
44		\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
42	–
43	–	\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION
44	–	TECHNIQUES}
45	–
46	–	In molecular simulations, proper accumulation of the electrostatic
47	–	interactions is essential and is one of the most
48	–	computationally-demanding tasks.  The common molecular mechanics force
49	–	fields represent atomic sites with full or partial charges protected
50	–	by Lennard-Jones (short range) interactions.  This means that nearly
51	–	every pair interaction involves a calculation of charge-charge forces.
52	–	Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
53	–	interactions quickly become the most expensive part of molecular
54	–	simulations.  Historically, the electrostatic pair interaction would
55	–	not have decayed appreciably within the typical box lengths that could
56	–	be feasibly simulated.  In the larger systems that are more typical of
57	–	modern simulations, large cutoffs should be used to incorporate
58	–	electrostatics correctly.
59	–
60	–	There have been many efforts to address the proper and practical
61	–	handling of electrostatic interactions, and these have resulted in a
62	–	variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are
63	–	typically classified as implicit methods (i.e., continuum dielectrics,
64	–	static dipolar fields),\cite{Born20,Grossfield00} explicit methods
65	–	(i.e., Ewald summations, interaction shifting or
66	–	truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
67	–	reaction field type methods, fast multipole
68	–	methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
69	–	often preferred because they physically incorporate solvent molecules
70	–	in the system of interest, but these methods are sometimes difficult
71	–	to utilize because of their high computational cost.\cite{Roux99} In
72	–	addition to the computational cost, there have been some questions
73	–	regarding possible artifacts caused by the inherent periodicity of the
74	–	explicit Ewald summation.\cite{Tobias01}
75	–
76	–	In this chapter, we focus on a new set of pairwise methods devised by
77	–	Wolf {\it et al.},\cite{Wolf99} which we further extend.  These
78	–	methods along with a few other mixed methods (i.e. reaction field) are
79	–	compared with the smooth particle mesh Ewald
80	–	sum,\cite{Onsager36,Essmann99} which is our reference method for
81	–	handling long-range electrostatic interactions. The new methods for
82	–	handling electrostatics have the potential to scale linearly with
83	–	increasing system size since they involve only a simple modification
84	–	to the direct pairwise sum.  They also lack the added periodicity of
85	–	the Ewald sum, so they can be used for systems which are non-periodic
86	–	or which have one- or two-dimensional periodicity.  Below, these
87	–	methods are evaluated using a variety of model systems to
88	–	establish their usability in molecular simulations.
89	–
90	–	\section{The Ewald Sum}
91	–
92	–	The complete accumulation of the electrostatic interactions in a system with
93	–	periodic boundary conditions (PBC) requires the consideration of the
94	–	effect of all charges within a (cubic) simulation box as well as those
95	–	in the periodic replicas,
96	–	\begin{equation}
97	–	V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime
98	–	\left[ \sum_{i=1}^N\sum_{j=1}^N \phi
99	–	\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right)
100	–	\right],
101	–	\label{eq:PBCSum}
102	–	\end{equation}
103	–	where the sum over $\mathbf{n}$ is a sum over all periodic box
104	–	replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the
105	–	prime indicates $i = j$ are neglected for $\mathbf{n} =
106	–	0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
107	–	particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
108	–	the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
109	–	$j$, and $\phi$ is the solution to Poisson's equation
110	–	($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for
111	–	charge-charge interactions). In the case of monopole electrostatics,
112	–	eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for
113	–	non-neutral systems.
114	–
115	–	The electrostatic summation problem was originally studied by Ewald
116	–	for the case of an infinite crystal.\cite{Ewald21}. The approach he
117	–	took was to convert this conditionally convergent sum into two
118	–	absolutely convergent summations: a short-ranged real-space summation
119	–	and a long-ranged reciprocal-space summation,
120	–	\begin{equation}
121	–	\begin{split}
122	–	V_\textrm{elec} = \frac{1}{2}&
123	–	\sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime
124	–	\frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}
125	–	{|\mathbf{r}_{ij}+\mathbf{n}|} \\
126	–	&+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2}
127	–	\exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right)
128	–	\cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\
129	–	&- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2
130	–	+ \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}
131	–	\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2,
132	–	\end{split}
133	–	\label{eq:EwaldSum}
134	–	\end{equation}
135	–	where $\alpha$ is the damping or convergence parameter with units of
136	–	\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
137	–	$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
138	–	constant of the surrounding medium. The final two terms of
139	–	eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
140	–	for interacting with a surrounding dielectric.\cite{Allen87} This
141	–	dipolar term was neglected in early applications in molecular
142	–	simulations,\cite{Brush66,Woodcock71} until it was introduced by de
143	–	Leeuw {\it et al.} to address situations where the unit cell has a
144	–	dipole moment which is magnified through replication of the periodic
145	–	images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
146	–	system is said to be using conducting (or ``tin-foil'') boundary
147	–	conditions, $\epsilon_{\rm S} = \infty$. Figure
148	–	\ref{fig:ewaldTime} shows how the Ewald sum has been applied over
149	–	time.  Initially, due to the small system sizes that could be
150	–	simulated feasibly, the entire simulation box was replicated to
151	–	convergence.  In more modern simulations, the systems have grown large
152	–	enough that a real-space cutoff could potentially give convergent
153	–	behavior.  Indeed, it has been observed that with the choice of a
154	–	small $\alpha$, the reciprocal-space portion of the Ewald sum can be
155	–	rapidly convergent and small relative to the real-space
156	–	portion.\cite{Karasawa89,Kolafa92}
157	–
158	–	\begin{figure}
159	–	\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
160	–	\caption{The change in the need for the Ewald sum with
161	–	increasing computational power.  A:~Initially, only small systems
162	–	could be studied, and the Ewald sum replicated the simulation box to
163	–	convergence.  B:~Now, radial cutoff methods should be able to reach
164	–	convergence for the larger systems of charges that are common today.}
165	–	\label{fig:ewaldTime}
166	–	\end{figure}
167	–
168	–	The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
169	–	convergence parameter $(\alpha)$ plays an important role in balancing
170	–	the computational cost between the direct and reciprocal-space
171	–	portions of the summation.  The choice of this value allows one to
172	–	select whether the real-space or reciprocal space portion of the
173	–	summation is an $\mathscr{O}(N^2)$ calculation (with the other being
174	–	$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of
175	–	$\alpha$ and thoughtful algorithm development, this cost can be
176	–	reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route
177	–	taken to reduce the cost of the Ewald summation even further is to set
178	–	$\alpha$ such that the real-space interactions decay rapidly, allowing
179	–	for a short spherical cutoff. Then the reciprocal space summation is
180	–	optimized.  These optimizations usually involve utilization of the
181	–	fast Fourier transform (FFT),\cite{Hockney81} leading to the
182	–	particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
183	–	methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
184	–	methods, the cost of the reciprocal-space portion of the Ewald
185	–	summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
186	–	\log N)$.
187	–
188	–	These developments and optimizations have made the use of the Ewald
189	–	summation routine in simulations with periodic boundary
190	–	conditions. However, in certain systems, such as vapor-liquid
191	–	interfaces and membranes, the intrinsic three-dimensional periodicity
192	–	can prove problematic.  The Ewald sum has been reformulated to handle
193	–	2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these
194	–	methods are computationally expensive.\cite{Spohr97,Yeh99} More
195	–	recently, there have been several successful efforts toward reducing
196	–	the computational cost of 2-D lattice
197	–	summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
198	–	bringing them more in line with the cost of the full 3-D summation.
199	–
200	–	Several studies have recognized that the inherent periodicity in the
201	–	Ewald sum can also have an effect on three-dimensional
202	–	systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
203	–	Solvated proteins are essentially kept at high concentration due to
204	–	the periodicity of the electrostatic summation method.  In these
205	–	systems, the more compact folded states of a protein can be
206	–	artificially stabilized by the periodic replicas introduced by the
207	–	Ewald summation.\cite{Weber00} Thus, care must be taken when
208	–	considering the use of the Ewald summation where the assumed
209	–	periodicity would introduce spurious effects in the system dynamics.
210	–
211	–
212	–	\section{The Wolf and Zahn Methods}
213	–
214	–	In a recent paper by Wolf \textit{et al.}, a procedure was outlined
215	–	for the accurate accumulation of electrostatic interactions in an
216	–	efficient pairwise fashion.  This procedure lacks the inherent
217	–	periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.}
218	–	observed that the electrostatic interaction is effectively
219	–	short-ranged in condensed phase systems and that neutralization of the
220	–	charge contained within the cutoff radius is crucial for potential
221	–	stability. They devised a pairwise summation method that ensures
222	–	charge neutrality and gives results similar to those obtained with the
223	–	Ewald summation.  The resulting shifted Coulomb potential includes
224	–	image-charges subtracted out through placement on the cutoff sphere
225	–	and a distance-dependent damping function (identical to that seen in
226	–	the real-space portion of the Ewald sum) to aid convergence
227	–	\begin{equation}
228	–	V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}
229	–	- \lim_{r_{ij}\rightarrow R_\textrm{c}}
230	–	\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}.
231	–	\label{eq:WolfPot}
232	–	\end{equation}
233	–	Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
234	–	potential.  However, neutralizing the charge contained within each
235	–	cutoff sphere requires the placement of a self-image charge on the
236	–	surface of the cutoff sphere.  This additional self-term in the total
237	–	potential enabled Wolf {\it et al.}  to obtain excellent estimates of
238	–	Madelung energies for many crystals.
239	–
240	–	In order to use their charge-neutralized potential in molecular
241	–	dynamics simulations, Wolf \textit{et al.} suggested taking the
242	–	derivative of this potential prior to evaluation of the limit.  This
243	–	procedure gives an expression for the forces,
244	–	\begin{equation}
245	–	\begin{split}
246	–	F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{&
247	–	\Biggr[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}}
248	–	+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}
249	–	\Biggr]\\
250	–	&-\Biggr[
251	–	\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}
252	–	+ \frac{2\alpha}{\pi^{1/2}}
253	–	\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}
254	–	\Biggr]\Biggr\},
255	–	\end{split}
256	–	\label{eq:WolfForces}
257	–	\end{equation}
258	–	that incorporates both image charges and damping of the electrostatic
259	–	interaction.
260	–
261	–	More recently, Zahn \textit{et al.} investigated these potential and
262	–	force expressions for use in simulations involving water.\cite{Zahn02}
263	–	In their work, they pointed out that the forces and derivative of
264	–	the potential are not commensurate.  Attempts to use both
265	–	eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead
266	–	to poor energy conservation.  They correctly observed that taking the
267	–	limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the
268	–	derivatives gives forces for a different potential energy function
269	–	than the one shown in eq. (\ref{eq:WolfPot}).
270	–
271	–	Zahn \textit{et al.} introduced a modified form of this summation
272	–	method as a way to use the technique in Molecular Dynamics
273	–	simulations.  They proposed a new damped Coulomb potential,
274	–	\begin{equation}
275	–	\begin{split}
276	–	V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{&
277	–	\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} \\
278	–	&- \left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
279	–	+ \frac{2\alpha}{\pi^{1/2}}
280	–	\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
281	–	\right]\left(r_{ij}-R_\mathrm{c}\right)\Biggr\},
282	–	\end{split}
283	–	\label{eq:ZahnPot}
284	–	\end{equation}
285	–	and showed that this potential does fairly well at capturing the
286	–	structural and dynamic properties of water compared the same
287	–	properties obtained using the Ewald sum.
288	–
289	–	\section{Simple Forms for Pairwise Electrostatics}
290	–
291	–	The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
292	–	al.} are constructed using two different (and separable) computational
293	–	tricks:
294	–
295	–	\begin{enumerate}
296	–	\item shifting through the use of image charges, and
297	–	\item damping the electrostatic interaction.
298	–	\end{enumerate}
299	–	Wolf \textit{et al.} treated the
300	–	development of their summation method as a progressive application of
301	–	these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
302	–	their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
303	–	post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
304	–	both techniques.  It is possible, however, to separate these
305	–	tricks and study their effects independently.
306	–
307	–	Starting with the original observation that the effective range of the
308	–	electrostatic interaction in condensed phases is considerably less
309	–	than $r^{-1}$, either the cutoff sphere neutralization or the
310	–	distance-dependent damping technique could be used as a foundation for
311	–	a new pairwise summation method.  Wolf \textit{et al.} made the
312	–	observation that charge neutralization within the cutoff sphere plays
313	–	a significant role in energy convergence; therefore we will begin our
314	–	analysis with the various shifted forms that maintain this charge
315	–	neutralization.  We can evaluate the methods of Wolf
316	–	\textit{et al.}  and Zahn \textit{et al.} by considering the standard
317	–	shifted potential,
318	–	\begin{equation}
319	–	V_\textrm{SP}(r) =      \begin{cases}
320	–	v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
321	–	R_\textrm{c}
322	–	\end{cases},
323	–	\label{eq:shiftingPotForm}
324	–	\end{equation}
325	–	and shifted force,
326	–	\begin{equation}
327	–	V_\textrm{SF}(r) = \begin{cases}
328	–	v(r) - v_\textrm{c}
329	–	- \left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
330	–	&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c}
331	–	\end{cases},
332	–	\label{eq:shiftingForm}
333	–	\end{equation}
334	–	functions where $v(r)$ is the unshifted form of the potential, and
335	–	$v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
336	–	that both the potential and the forces goes to zero at the cutoff
337	–	radius, while the Shifted Potential ({\sc sp}) form only ensures the
338	–	potential is smooth at the cutoff radius
339	–	($R_\textrm{c}$).\cite{Allen87}
340	–
341	–	The forces associated with the shifted potential are simply the forces
342	–	of the unshifted potential itself (when inside the cutoff sphere),
343	–	\begin{equation}
344	–	F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
345	–	\end{equation}
346	–	and are zero outside.  Inside the cutoff sphere, the forces associated
347	–	with the shifted force form can be written,
348	–	\begin{equation}
349	–	F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
350	–	v(r)}{dr} \right)_{r=R_\textrm{c}}.
351	–	\end{equation}
352	–
353	–	If the potential, $v(r)$, is taken to be the normal Coulomb potential,
354	–	\begin{equation}
355	–	v(r) = \frac{q_i q_j}{r},
356	–	\label{eq:Coulomb}
357	–	\end{equation}
358	–	then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
359	–	al.}'s undamped prescription:
360	–	\begin{equation}
361	–	V_\textrm{SP}(r) =
362	–	q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
363	–	r\leqslant R_\textrm{c},
364	–	\label{eq:SPPot}
365	–	\end{equation}
366	–	with associated forces,
367	–	\begin{equation}
368	–	F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right)
369	–	\quad r\leqslant R_\textrm{c}.
370	–	\label{eq:SPForces}
371	–	\end{equation}
372	–	These forces are identical to the forces of the standard Coulomb
373	–	interaction, and cutting these off at $R_c$ was addressed by Wolf
374	–	\textit{et al.} as undesirable.  They pointed out that the effect of
375	–	the image charges is neglected in the forces when this form is
376	–	used,\cite{Wolf99} thereby eliminating any benefit from the method in
377	–	molecular dynamics.  Additionally, there is a discontinuity in the
378	–	forces at the cutoff radius which results in energy drift during MD
379	–	simulations.
380	–
381	–	The shifted force ({\sc sf}) form using the normal Coulomb potential
382	–	will give,
383	–	\begin{equation}
384	–	V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}
385	–	+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right]
386	–	\quad r\leqslant R_\textrm{c}.
387	–	\label{eq:SFPot}
388	–	\end{equation}
389	–	with associated forces,
390	–	\begin{equation}
391	–	F_\textrm{SF}(r) =  q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right)
392	–	\quad r\leqslant R_\textrm{c}.
393	–	\label{eq:SFForces}
394	–	\end{equation}
395	–	This formulation has the benefits that there are no discontinuities at
396	–	the cutoff radius, while the neutralizing image charges are present in
397	–	both the energy and force expressions.  It would be simple to add the
398	–	self-neutralizing term back when computing the total energy of the
399	–	system, thereby maintaining the agreement with the Madelung energies.
400	–	A side effect of this treatment is the alteration in the shape of the
401	–	potential that comes from the derivative term.  Thus, a degree of
402	–	clarity about agreement with the empirical potential is lost in order
403	–	to gain functionality in dynamics simulations.
404	–
405	–	Wolf \textit{et al.} originally discussed the energetics of the
406	–	shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was
407	–	insufficient for accurate determination of the energy with reasonable
408	–	cutoff distances.  The calculated Madelung energies fluctuated around
409	–	the expected value as the cutoff radius was increased, but the
410	–	oscillations converged toward the correct value.\cite{Wolf99} A
411	–	damping function was incorporated to accelerate the convergence; and
412	–	though alternative forms for the damping function could be
413	–	used,\cite{Jones56,Heyes81} the complimentary error function was
414	–	chosen to mirror the effective screening used in the Ewald summation.
415	–	Incorporating this error function damping into the simple Coulomb
416	–	potential,
417	–	\begin{equation}
418	–	v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
419	–	\label{eq:dampCoulomb}
420	–	\end{equation}
421	–	the shifted potential (eq. (\ref{eq:SPPot})) becomes
422	–	\begin{equation}
423	–	V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}
424	–	- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right)
425	–	\quad r\leqslant R_\textrm{c},
426	–	\label{eq:DSPPot}
427	–	\end{equation}
428	–	with associated forces,
429	–	\begin{equation}
430	–	F_{\textrm{DSP}}(r) = q_iq_j
431	–	\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}
432	–	+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right)
433	–	\quad r\leqslant R_\textrm{c}.
434	–	\label{eq:DSPForces}
435	–	\end{equation}
436	–	Again, this damped shifted potential suffers from a
437	–	force-discontinuity at the cutoff radius, and the image charges play
438	–	no role in the forces.  To remedy these concerns, one may derive a
439	–	{\sc sf} variant by including the derivative term in
440	–	eq. (\ref{eq:shiftingForm}),
441	–	\begin{equation}
442	–	\begin{split}
443	–	V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
444	–	\frac{\mathrm{erfc}\left(\alpha r\right)}{r}
445	–	- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\
446	–	&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
447	–	+ \frac{2\alpha}{\pi^{1/2}}
448	–	\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
449	–	\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]}
450	–	\quad r\leqslant R_\textrm{c}.
451	–	\label{eq:DSFPot}
452	–	\end{split}
453	–	\end{equation}
454	–	The derivative of the above potential will lead to the following forces,
455	–	\begin{equation}
456	–	\begin{split}
457	–	F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
458	–	\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}
459	–	+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\
460	–	&- \left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}
461	–	{R_{\textrm{c}}^2}
462	–	+ \frac{2\alpha}{\pi^{1/2}}
463	–	\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}
464	–	\right)\Biggr{]}
465	–	\quad r\leqslant R_\textrm{c}.
466	–	\label{eq:DSFForces}
467	–	\end{split}
468	–	\end{equation}
469	–	If the damping parameter $(\alpha)$ is set to zero, the undamped case,
470	–	eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
471	–	recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}).
472	–
473	–	This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
474	–	derived by Zahn \textit{et al.}; however, there are two important
475	–	differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from
476	–	eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
477	–	with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
478	–	in the Zahn potential, resulting in a potential discontinuity as
479	–	particles cross $R_\textrm{c}$.  Second, the sign of the derivative
480	–	portion is different.  The missing $v_\textrm{c}$ term would not
481	–	affect molecular dynamics simulations (although the computed energy
482	–	would be expected to have sudden jumps as particle distances crossed
483	–	$R_c$).  The sign problem is a potential source of errors, however.
484	–	In fact, it introduces a discontinuity in the forces at the cutoff,
485	–	because the force function is shifted in the wrong direction and
486	–	doesn't cross zero at $R_\textrm{c}$.
487	–
488	–	Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
489	–	electrostatic summation method in which the potential and forces are
490	–	continuous at the cutoff radius and which incorporates the damping
491	–	function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
492	–	this paper, we will evaluate exactly how good these methods ({\sc sp},
493	–	{\sc sf}, damping) are at reproducing the correct electrostatic
494	–	summation performed by the Ewald sum.
495	–
496	–
497	–	\section{Evaluating Pairwise Summation Techniques}
498	–
499	–	In classical molecular mechanics simulations, there are two primary
500	–	techniques utilized to obtain information about the system of
501	–	interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
502	–	techniques utilize pairwise summations of interactions between
503	–	particle sites, but they use these summations in different ways.
504	–
505	–	In MC, the potential energy difference between configurations dictates
506	–	the progression of MC sampling.  Going back to the origins of this
507	–	method, the acceptance criterion for the canonical ensemble laid out
508	–	by Metropolis \textit{et al.} states that a subsequent configuration
509	–	is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where
510	–	$\xi$ is a random number between 0 and 1.\cite{Metropolis53}
511	–	Maintaining the correct $\Delta E$ when using an alternate method for
512	–	handling the long-range electrostatics will ensure proper sampling
513	–	from the ensemble.
514	–
515	–	In MD, the derivative of the potential governs how the system will
516	–	progress in time.  Consequently, the force and torque vectors on each
517	–	body in the system dictate how the system evolves.  If the magnitude
518	–	and direction of these vectors are similar when using alternate
519	–	electrostatic summation techniques, the dynamics in the short term
520	–	will be indistinguishable.  Because error in MD calculations is
521	–	cumulative, one should expect greater deviation at longer times,
522	–	although methods which have large differences in the force and torque
523	–	vectors will diverge from each other more rapidly.
524	–
525	–	\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
526	–
527	–	The pairwise summation techniques (outlined in section
528	–	\ref{sec:ESMethods}) were evaluated for use in MC simulations by
529	–	studying the energy differences between conformations.  We took the
530	–	{\sc spme}-computed energy difference between two conformations to be the
531	–	correct behavior. An ideal performance by an alternative method would
532	–	reproduce these energy differences exactly (even if the absolute
533	–	energies calculated by the methods are different).  Since none of the
534	–	methods provide exact energy differences, we used linear least squares
535	–	regressions of energy gap data to evaluate how closely the methods
536	–	mimicked the Ewald energy gaps.  Unitary results for both the
537	–	correlation (slope) and correlation coefficient for these regressions
538	–	indicate perfect agreement between the alternative method and {\sc spme}.
539	–	Sample correlation plots for two alternate methods are shown in
540	–	Fig. \ref{fig:linearFit}.
45
46	<	\begin{figure}
543	<	\centering
544	<	\includegraphics[width = \linewidth]{./figures/dualLinear.pdf}
545	<	\caption{Example least squares regressions of the configuration energy
546	<	differences for SPC/E water systems. The upper plot shows a data set
547	<	with a poor correlation coefficient ($R^2$), while the lower plot
548	<	shows a data set with a good correlation coefficient.}
549	<	\label{fig:linearFit}
550	<	\end{figure}
551	<
552	<	Each of the seven system types (detailed in section \ref{sec:RepSims})
553	<	were represented using 500 independent configurations.  Thus, each of
554	<	the alternative (non-Ewald) electrostatic summation methods was
555	<	evaluated using an accumulated 873,250 configurational energy
556	<	differences.
557	<
558	<	Results and discussion for the individual analysis of each of the
559	<	system types appear in sections \ref{sec:IndividualResults}, while the
560	<	cumulative results over all the investigated systems appear below in
561	<	sections \ref{sec:EnergyResults}.
562	<
563	<	\subsection{Molecular Dynamics and the Force and Torque
564	<	Vectors}\label{sec:MDMethods} We evaluated the pairwise methods
565	<	(outlined in section \ref{sec:ESMethods}) for use in MD simulations by
566	<	comparing the force and torque vectors with those obtained using the
567	<	reference Ewald summation ({\sc spme}).  Both the magnitude and the
568	<	direction of these vectors on each of the bodies in the system were
569	<	analyzed.  For the magnitude of these vectors, linear least squares
570	<	regression analyses were performed as described previously for
571	<	comparing $\Delta E$ values.  Instead of a single energy difference
572	<	between two system configurations, we compared the magnitudes of the
573	<	forces (and torques) on each molecule in each configuration.  For a
574	<	system of 1000 water molecules and 40 ions, there are 1040 force
575	<	vectors and 1000 torque vectors.  With 500 configurations, this
576	<	results in 520,000 force and 500,000 torque vector comparisons.
577	<	Additionally, data from seven different system types was aggregated
578	<	before the comparison was made.
579	<
580	<	The {\it directionality} of the force and torque vectors was
581	<	investigated through measurement of the angle ($\theta$) formed
582	<	between those computed from the particular method and those from {\sc spme},
583	<	\begin{equation}
584	<	\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME}
585	<	\cdot \hat{F}_\textrm{M}\right),
586	<	\end{equation}
587	<	where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force
588	<	vector computed using method M.  Each of these $\theta$ values was
589	<	accumulated in a distribution function and weighted by the area on the
590	<	unit sphere.  Since this distribution is a measure of angular error
591	<	between two different electrostatic summation methods, there is no
592	<	{\it a priori} reason for the profile to adhere to any specific
593	<	shape. Thus, gaussian fits were used to measure the width of the
594	<	resulting distributions. The variance ($\sigma^2$) was extracted from
595	<	each of these fits and was used to compare distribution widths.
596	<	Values of $\sigma^2$ near zero indicate vector directions
597	<	indistinguishable from those calculated when using the reference
598	<	method ({\sc spme}).
599	<
600	<	\subsection{Short-time Dynamics}
601	<
602	<	The effects of the alternative electrostatic summation methods on the
603	<	short-time dynamics of charged systems were evaluated by considering a
604	<	NaCl crystal at a temperature of 1000 K.  A subset of the best
605	<	performing pairwise methods was used in this comparison.  The NaCl
606	<	crystal was chosen to avoid possible complications from the treatment
607	<	of orientational motion in molecular systems.  All systems were
608	<	started with the same initial positions and velocities.  Simulations
609	<	were performed under the microcanonical ensemble, and velocity
610	<	autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each
611	<	of the trajectories,
612	<	\begin{equation}
613	<	C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}.
614	<	\label{eq:vCorr}
615	<	\end{equation}
616	<	Velocity autocorrelation functions require detailed short time data,
617	<	thus velocity information was saved every 2 fs over 10 ps
618	<	trajectories. Because the NaCl crystal is composed of two different
619	<	atom types, the average of the two resulting velocity autocorrelation
620	<	functions was used for comparisons.
621	<
622	<	\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
623	<
624	<	The effects of the same subset of alternative electrostatic methods on
625	<	the {\it long-time} dynamics of charged systems were evaluated using
626	<	the same model system (NaCl crystals at 1000K).  The power spectrum
627	<	($I(\omega)$) was obtained via Fourier transform of the velocity
628	<	autocorrelation function, \begin{equation} I(\omega) =
629	<	\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
630	<	\label{eq:powerSpec}
631	<	\end{equation}
632	<	where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
633	<	NaCl crystal is composed of two different atom types, the average of
634	<	the two resulting power spectra was used for comparisons. Simulations
635	<	were performed under the microcanonical ensemble, and velocity
636	<	information was saved every 5~fs over 100~ps trajectories.
637	<
638	<	\subsection{Representative Simulations}\label{sec:RepSims}
639	<	A variety of representative molecular simulations were analyzed to
640	<	determine the relative effectiveness of the pairwise summation
641	<	techniques in reproducing the energetics and dynamics exhibited by
642	<	{\sc spme}.  We wanted to span the space of typical molecular
643	<	simulations (i.e. from liquids of neutral molecules to ionic
644	<	crystals), so the systems studied were:
645	<
646	<	\begin{enumerate}
647	<	\item liquid water (SPC/E),\cite{Berendsen87}
648	<	\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E),
649	<	\item NaCl crystals,
650	<	\item NaCl melts,
651	<	\item a low ionic strength solution of NaCl in water (0.11 M),
652	<	\item a high ionic strength solution of NaCl in water (1.1 M), and
653	<	\item a 6\AA\  radius sphere of Argon in water.
654	<	\end{enumerate}
655	<
656	<	By utilizing the pairwise techniques (outlined in section
657	<	\ref{sec:ESMethods}) in systems composed entirely of neutral groups,
658	<	charged particles, and mixtures of the two, we hope to discern under
659	<	which conditions it will be possible to use one of the alternative
660	<	summation methodologies instead of the Ewald sum.
661	<
662	<	For the solid and liquid water configurations, configurations were
663	<	taken at regular intervals from high temperature trajectories of 1000
664	<	SPC/E water molecules.  Each configuration was equilibrated
665	<	independently at a lower temperature (300K for the liquid, 200K for
666	<	the crystal).  The solid and liquid NaCl systems consisted of 500
667	<	$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions.  Configurations for
668	<	these systems were selected and equilibrated in the same manner as the
669	<	water systems. In order to introduce measurable fluctuations in the
670	<	configuration energy differences, the crystalline simulations were
671	<	equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid
672	<	NaCl configurations needed to represent a fully disordered array of
673	<	point charges, so the high temperature of 7000K was selected for
674	<	equilibration. The ionic solutions were made by solvating 4 (or 40)
675	<	ions in a periodic box containing 1000 SPC/E water molecules.  Ion and
676	<	water positions were then randomly swapped, and the resulting
677	<	configurations were again equilibrated individually.  Finally, for the
678	<	Argon / Water ``charge void'' systems, the identities of all the SPC/E
679	<	waters within 6\AA\ of the center of the equilibrated water
680	<	configurations were converted to argon.
681	<
682	<	These procedures guaranteed us a set of representative configurations
683	<	from chemically-relevant systems sampled from appropriate
684	<	ensembles. Force field parameters for the ions and Argon were taken
685	<	from the force field utilized by {\sc oopse}.\cite{Meineke05}
686	<
687	<	\subsection{Comparison of Summation Methods}\label{sec:ESMethods}
688	<	We compared the following alternative summation methods with results
689	<	from the reference method ({\sc spme}):
690	<
691	<	\begin{enumerate}
692	<	\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
693	<	and 0.3\AA$^{-1}$,
694	<	\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
695	<	and 0.3\AA$^{-1}$,
696	<	\item reaction field with an infinite dielectric constant, and
697	<	\item an unmodified cutoff.
698	<	\end{enumerate}
699	<
700	<	Group-based cutoffs with a fifth-order polynomial switching function
701	<	were utilized for the reaction field simulations.  Additionally, we
702	<	investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure
703	<	cutoff.  The {\sc spme} electrostatics were performed using the {\sc tinker}
704	<	implementation of {\sc spme},\cite{Ponder87} while all other calculations
705	<	were performed using the {\sc oopse} molecular mechanics
706	<	package.\cite{Meineke05} All other portions of the energy calculation
707	<	(i.e. Lennard-Jones interactions) were handled in exactly the same
708	<	manner across all systems and configurations.
709	<
710	<	The alternative methods were also evaluated with three different
711	<	cutoff radii (9, 12, and 15\AA).  As noted previously, the
712	<	convergence parameter ($\alpha$) plays a role in the balance of the
713	<	real-space and reciprocal-space portions of the Ewald calculation.
714	<	Typical molecular mechanics packages set this to a value dependent on
715	<	the cutoff radius and a tolerance (typically less than $1 \times
716	<	10^{-4}$ kcal/mol).  Smaller tolerances are typically associated with
717	<	increasing accuracy at the expense of computational time spent on the
718	<	reciprocal-space portion of the summation.\cite{Perram88,Essmann95}
719	<	The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used
720	<	in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200,
721	<	0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\
722	<	respectively.
723	<
724	<	\section{Configuration Energy Difference Results}\label{sec:EnergyResults}
725	<	In order to evaluate the performance of the pairwise electrostatic
726	<	summation methods for Monte Carlo (MC) simulations, the energy
727	<	differences between configurations were compared to the values
728	<	obtained when using {\sc spme}.  The results for the combined
729	<	regression analysis of all of the systems are shown in figure
730	<	\ref{fig:delE}.
731	<
732	<	\begin{figure}
733	<	\centering
734	<	\includegraphics[width=4.75in]{./figures/delEplot.pdf}
735	<	\caption{Statistical analysis of the quality of configurational energy
736	<	differences for a given electrostatic method compared with the
737	<	reference Ewald sum.  Results with a value equal to 1 (dashed line)
738	<	indicate $\Delta E$ values indistinguishable from those obtained using
739	<	{\sc spme}.  Different values of the cutoff radius are indicated with
740	<	different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
741	<	inverted triangles).}
742	<	\label{fig:delE}
743	<	\end{figure}
744	<
745	<	The most striking feature of this plot is how well the Shifted Force
746	<	({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
747	<	differences.  For the undamped {\sc sf} method, and the
748	<	moderately-damped {\sc sp} methods, the results are nearly
749	<	indistinguishable from the Ewald results.  The other common methods do
750	<	significantly less well.
751	<
752	<	The unmodified cutoff method is essentially unusable.  This is not
753	<	surprising since hard cutoffs give large energy fluctuations as atoms
754	<	or molecules move in and out of the cutoff
755	<	radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to
756	<	some degree by using group based cutoffs with a switching
757	<	function.\cite{Adams79,Steinbach94,Leach01} However, we do not see
758	<	significant improvement using the group-switched cutoff because the
759	<	salt and salt solution systems contain non-neutral groups.  Section
760	<	\ref{sec:IndividualResults} includes results for systems comprised entirely
761	<	of neutral groups.
762	<
763	<	For the {\sc sp} method, inclusion of electrostatic damping improves
764	<	the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$
765	<	shows an excellent correlation and quality of fit with the {\sc spme}
766	<	results, particularly with a cutoff radius greater than 12
767	<	\AA .  Use of a larger damping parameter is more helpful for the
768	<	shortest cutoff shown, but it has a detrimental effect on simulations
769	<	with larger cutoffs.
770	<
771	<	In the {\sc sf} sets, increasing damping results in progressively {\it
772	<	worse} correlation with Ewald.  Overall, the undamped case is the best
773	<	performing set, as the correlation and quality of fits are
774	<	consistently superior regardless of the cutoff distance.  The undamped
775	<	case is also less computationally demanding (because no evaluation of
776	<	the complementary error function is required).
777	<
778	<	The reaction field results illustrates some of that method's
779	<	limitations, primarily that it was developed for use in homogenous
780	<	systems; although it does provide results that are an improvement over
781	<	those from an unmodified cutoff.
782	<
783	<	\section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults}
784	<
785	<	Evaluation of pairwise methods for use in Molecular Dynamics
786	<	simulations requires consideration of effects on the forces and
787	<	torques.  Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the
788	<	regression results for the force and torque vector magnitudes,
789	<	respectively.  The data in these figures was generated from an
790	<	accumulation of the statistics from all of the system types.
791	<
792	<	\begin{figure}
793	<	\centering
794	<	\includegraphics[width=4.75in]{./figures/frcMagplot.pdf}
795	<	\caption{Statistical analysis of the quality of the force vector
796	<	magnitudes for a given electrostatic method compared with the
797	<	reference Ewald sum.  Results with a value equal to 1 (dashed line)
798	<	indicate force magnitude values indistinguishable from those obtained
799	<	using {\sc spme}.  Different values of the cutoff radius are indicated with
800	<	different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
801	<	inverted triangles).}
802	<	\label{fig:frcMag}
803	<	\end{figure}
804	<
805	<	Again, it is striking how well the Shifted Potential and Shifted Force
806	<	methods are doing at reproducing the {\sc spme} forces.  The undamped and
807	<	weakly-damped {\sc sf} method gives the best agreement with Ewald.
808	<	This is perhaps expected because this method explicitly incorporates a
809	<	smooth transition in the forces at the cutoff radius as well as the
810	<	neutralizing image charges.
811	<
812	<	Figure \ref{fig:frcMag}, for the most part, parallels the results seen
813	<	in the previous $\Delta E$ section.  The unmodified cutoff results are
814	<	poor, but using group based cutoffs and a switching function provides
815	<	an improvement much more significant than what was seen with $\Delta
816	<	E$.
817	<
818	<	With moderate damping and a large enough cutoff radius, the {\sc sp}
819	<	method is generating usable forces.  Further increases in damping,
820	<	while beneficial for simulations with a cutoff radius of 9\AA\ , is
821	<	detrimental to simulations with larger cutoff radii.
822	<
823	<	The reaction field results are surprisingly good, considering the poor
824	<	quality of the fits for the $\Delta E$ results.  There is still a
825	<	considerable degree of scatter in the data, but the forces correlate
826	<	well with the Ewald forces in general.  We note that the reaction
827	<	field calculations do not include the pure NaCl systems, so these
828	<	results are partly biased towards conditions in which the method
829	<	performs more favorably.
830	<
831	<	\begin{figure}
832	<	\centering
833	<	\includegraphics[width=4.75in]{./figures/trqMagplot.pdf}
834	<	\caption{Statistical analysis of the quality of the torque vector
835	<	magnitudes for a given electrostatic method compared with the
836	<	reference Ewald sum.  Results with a value equal to 1 (dashed line)
837	<	indicate torque magnitude values indistinguishable from those obtained
838	<	using {\sc spme}.  Different values of the cutoff radius are indicated with
839	<	different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
840	<	inverted triangles).}
841	<	\label{fig:trqMag}
842	<	\end{figure}
843	<
844	<	Molecular torques were only available from the systems which contained
845	<	rigid molecules (i.e. the systems containing water).  The data in
846	<	fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
847	<
848	<	Torques appear to be much more sensitive to charges at a longer
849	<	distance.   The striking feature in comparing the new electrostatic
850	<	methods with {\sc spme} is how much the agreement improves with increasing
851	<	cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
852	<	appears to be reproducing the {\sc spme} torques most accurately.
853	<
854	<	Water molecules are dipolar, and the reaction field method reproduces
855	<	the effect of the surrounding polarized medium on each of the
856	<	molecular bodies. Therefore it is not surprising that reaction field
857	<	performs best of all of the methods on molecular torques.
858	<
859	<	\section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}
860	<
861	<	It is clearly important that a new electrostatic method can reproduce
862	<	the magnitudes of the force and torque vectors obtained via the Ewald
863	<	sum. However, the {\it directionality} of these vectors will also be
864	<	vital in calculating dynamical quantities accurately.  Force and
865	<	torque directionalities were investigated by measuring the angles
866	<	formed between these vectors and the same vectors calculated using
867	<	{\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
868	<	variance ($\sigma^2$) of the Gaussian fits of the angle error
869	<	distributions of the combined set over all system types.
870	<
871	<	\begin{figure}
872	<	\centering
873	<	\includegraphics[width=4.75in]{./figures/frcTrqAngplot.pdf}
874	<	\caption{Statistical analysis of the width of the angular distribution
875	<	that the force and torque vectors from a given electrostatic method
876	<	make with their counterparts obtained using the reference Ewald sum.
877	<	Results with a variance ($\sigma^2$) equal to zero (dashed line)
878	<	indicate force and torque directions indistinguishable from those
879	<	obtained using {\sc spme}.  Different values of the cutoff radius are
880	<	indicated with different symbols (9\AA\ = circles, 12\AA\ = squares,
881	<	and 15\AA\ = inverted triangles).}
882	<	\label{fig:frcTrqAng}
883	<	\end{figure}
884	<
885	<	Both the force and torque $\sigma^2$ results from the analysis of the
886	<	total accumulated system data are tabulated in figure
887	<	\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
888	<	sp}) method would be essentially unusable for molecular dynamics
889	<	unless the damping function is added.  The Shifted Force ({\sc sf})
890	<	method, however, is generating force and torque vectors which are
891	<	within a few degrees of the Ewald results even with weak (or no)
892	<	damping.
893	<
894	<	All of the sets (aside from the over-damped case) show the improvement
895	<	afforded by choosing a larger cutoff radius.  Increasing the cutoff
896	<	from 9 to 12\AA\ typically results in a halving of the width of the
897	<	distribution, with a similar improvement when going from 12 to 15
898	<	\AA .
899	<
900	<	The undamped {\sc sf}, group-based cutoff, and reaction field methods
901	<	all do equivalently well at capturing the direction of both the force
902	<	and torque vectors.  Using the electrostatic damping improves the
903	<	angular behavior significantly for the {\sc sp} and moderately for the
904	<	{\sc sf} methods.  Overdamping is detrimental to both methods.  Again
905	<	it is important to recognize that the force vectors cover all
906	<	particles in all seven systems, while torque vectors are only
907	<	available for neutral molecular groups.  Damping is more beneficial to
908	<	charged bodies, and this observation is investigated further in
909	<	section \ref{IndividualResults}.
910	<
911	<	Although not discussed previously, group based cutoffs can be applied
912	<	to both the {\sc sp} and {\sc sf} methods.  The group-based cutoffs
913	<	will reintroduce small discontinuities at the cutoff radius, but the
914	<	effects of these can be minimized by utilizing a switching function.
915	<	Though there are no significant benefits or drawbacks observed in
916	<	$\Delta E$ and the force and torque magnitudes when doing this, there
917	<	is a measurable improvement in the directionality of the forces and
918	<	torques. Table \ref{tab:groupAngle} shows the angular variances
919	<	obtained both without (N) and with (Y) group based cutoffs and a
920	<	switching function.  Note that the $\alpha$ values have units of
921	<	\AA$^{-1}$ and the variance values have units of degrees$^2$.  The
922	<	{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much
923	<	narrower angular distributions when using group-based cutoffs. The
924	<	{\sc sf} method likewise shows improvement in the undamped and lightly
925	<	damped cases.
926	<
927	<	\begin{table}
928	<	\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$)
929	<	THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A
930	<	GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING
931	<	THE REFERENCE EWALD SUMMATION}
932	<
933	<	\footnotesize
934	<	\begin{center}
935	<	\begin{tabular}{@{} ccrrrrrrrr @{}} \\
936	<	\toprule
937	<	\toprule
938	<
939	<	& &\multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted
940	<	Force} \\
941	<	\cmidrule(lr){3-6} \cmidrule(l){7-10} $R_\textrm{c}$ & Groups &
942	<	$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ &
943	<	$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\
944	<
945	<	\midrule
946	<
947	<	9\AA   & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\
948	<	& \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\
949	<	12\AA  & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\
950	<	& \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\
951	<	15\AA  & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\
952	<	& \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\
953	<
954	<	\midrule
955	<
956	<	9\AA   & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\
957	<	& \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\
958	<	12\AA  & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\
959	<	& \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\
960	<	15\AA  & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\
961	<	& \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\
962	<
963	<	\bottomrule
964	<	\end{tabular}
965	<	\end{center}
966	<	\label{tab:groupAngle}
967	<	\end{table}
968	<
969	<	One additional trend in table \ref{tab:groupAngle} is that the
970	<	$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$
971	<	increases, something that is more obvious with group-based cutoffs.
972	<	The complimentary error function inserted into the potential weakens
973	<	the electrostatic interaction as the value of $\alpha$ is increased.
974	<	However, at larger values of $\alpha$, it is possible to overdamp the
975	<	electrostatic interaction and to remove it completely.  Kast
976	<	\textit{et al.}  developed a method for choosing appropriate $\alpha$
977	<	values for these types of electrostatic summation methods by fitting
978	<	to $g(r)$ data, and their methods indicate optimal values of 0.34,
979	<	0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\
980	<	respectively.\cite{Kast03} These appear to be reasonable choices to
981	<	obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on
982	<	these findings, choices this high would introduce error in the
983	<	molecular torques, particularly for the shorter cutoffs.  Based on our
984	<	observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial,
985	<	but damping may be unnecessary when using the {\sc sf} method.
986	<
987	<	\section{Individual System Analysis Results}\label{sec:IndividualResults}
988	<
989	<	The combined results of the previous sections show how the pairwise
990	<	methods compare to the Ewald summation in the general sense over all
991	<	of the system types.  It is also useful to consider each of the
992	<	studied systems in an individual fashion, so that we can identify
993	<	conditions that are particularly difficult for a selected pairwise
994	<	method to address. This allows us to further establish the limitations
995	<	of these pairwise techniques.  Below, the energy difference, force
996	<	vector, and torque vector analyses are presented on an individual
997	<	system basis.
998	<
999	<	\subsection{SPC/E Water Results}\label{sec:WaterResults}
1000	<
1001	<	The first system considered was liquid water at 300K using the SPC/E
1002	<	model of water.\cite{Berendsen87} The results for the energy gap
1003	<	comparisons and the force and torque vector magnitude comparisons are
1004	<	shown in table \ref{tab:spce}.  The force and torque vector
1005	<	directionality results are displayed separately in table
1006	<	\ref{tab:spceAng}, where the effect of group-based cutoffs and
1007	<	switching functions on the {\sc sp} and {\sc sf} potentials are also
1008	<	investigated.  In all of the individual results table, the method
1009	<	abbreviations are as follows:
1010	<
1011	<	\begin{itemize}
1012	<	\item PC = Pure Cutoff,
1013	<	\item SP = Shifted Potential,
1014	<	\item SF = Shifted Force,
1015	<	\item GSC = Group Switched Cutoff,
1016	<	\item RF = Reaction Field (where $\varepsilon \approx\infty$),
1017	<	\item GSSP = Group Switched Shifted Potential, and
1018	<	\item GSSF = Group Switched Shifted Force.
1019	<	\end{itemize}
1020	<
1021	<	\begin{table}[htbp]
1022	<	\centering
1023	<	\caption{REGRESSION RESULTS OF THE LIQUID WATER SYSTEM FOR THE
1024	<	$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it middle})
1025	<	AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1026	<
1027	<	\footnotesize
1028	<	\begin{tabular}{@{} ccrrrrrr @{}}
1029	<	\\
1030	<	\toprule
1031	<	\toprule
1032	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1033	<	\cmidrule(lr){3-4}
1034	<	\cmidrule(lr){5-6}
1035	<	\cmidrule(l){7-8}
1036	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1037	<	\midrule
1038	<	PC  &     & 3.046 & 0.002 & -3.018 & 0.002 & 4.719 & 0.005 \\
1039	<	SP  & 0.0 & 1.035 & 0.218 & 0.908 & 0.313 & 1.037 & 0.470 \\
1040	<	& 0.1 & 1.021 & 0.387 & 0.965 & 0.752 & 1.006 & 0.947 \\
1041	<	& 0.2 & 0.997 & 0.962 & 1.001 & 0.994 & 0.994 & 0.996 \\
1042	<	& 0.3 & 0.984 & 0.980 & 0.997 & 0.985 & 0.982 & 0.987 \\
1043	<	SF  & 0.0 & 0.977 & 0.974 & 0.996 & 0.992 & 0.991 & 0.997 \\
1044	<	& 0.1 & 0.983 & 0.974 & 1.001 & 0.994 & 0.996 & 0.998 \\
1045	<	& 0.2 & 0.992 & 0.989 & 1.001 & 0.995 & 0.994 & 0.996 \\
1046	<	& 0.3 & 0.984 & 0.980 & 0.996 & 0.985 & 0.982 & 0.987 \\
1047	<	GSC &     & 0.918 & 0.862 & 0.852 & 0.756 & 0.801 & 0.700 \\
1048	<	RF  &     & 0.971 & 0.958 & 0.975 & 0.987 & 0.959 & 0.983 \\
1049	<	\midrule
1050	<	PC  &     & -1.647 & 0.000 & -0.127 & 0.000 & -0.979 & 0.000 \\
1051	<	SP  & 0.0 & 0.735 & 0.368 & 0.813 & 0.537 & 0.865 & 0.659 \\
1052	<	& 0.1 & 0.850 & 0.612 & 0.956 & 0.887 & 0.992 & 0.979 \\
1053	<	& 0.2 & 0.996 & 0.989 & 1.000 & 1.000 & 1.000 & 1.000 \\
1054	<	& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\
1055	<	SF  & 0.0 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 0.999 \\
1056	<	& 0.1 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\
1057	<	& 0.2 & 0.999 & 0.998 & 1.000 & 1.000 & 1.000 & 1.000 \\
1058	<	& 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\
1059	<	GSC &     & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\
1060	<	RF  &     & 0.999 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\
1061	<	\midrule
1062	<	PC  &     & 2.387 & 0.000 & 0.183 & 0.000 & 1.282 & 0.000 \\
1063	<	SP  & 0.0 & 0.847 & 0.543 & 0.904 & 0.694 & 0.935 & 0.786 \\
1064	<	& 0.1 & 0.922 & 0.749 & 0.980 & 0.934 & 0.996 & 0.988 \\
1065	<	& 0.2 & 0.987 & 0.985 & 0.989 & 0.992 & 0.990 & 0.993 \\
1066	<	& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\
1067	<	SF  & 0.0 & 0.978 & 0.990 & 0.988 & 0.997 & 0.993 & 0.999 \\
1068	<	& 0.1 & 0.983 & 0.991 & 0.993 & 0.997 & 0.997 & 0.999 \\
1069	<	& 0.2 & 0.986 & 0.989 & 0.989 & 0.992 & 0.990 & 0.993 \\
1070	<	& 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\
1071	<	GSC &     & 0.995 & 0.981 & 0.999 & 0.991 & 1.001 & 0.994 \\
1072	<	RF  &     & 0.993 & 0.989 & 0.998 & 0.996 & 1.000 & 0.999 \\
1073	<	\bottomrule
1074	<	\end{tabular}
1075	<	\label{tab:spce}
1076	<	\end{table}
1077	<
1078	<	\begin{table}[htbp]
1079	<	\centering
1080	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1081	<	DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE LIQUID WATER
1082	<	SYSTEM}
1083	<
1084	<	\footnotesize
1085	<	\begin{tabular}{@{} ccrrrrrr @{}}
1086	<	\\
1087	<	\toprule
1088	<	\toprule
1089	<	& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\
1090	<	\cmidrule(lr){3-5}
1091	<	\cmidrule(l){6-8}
1092	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\
1093	<	\midrule
1094	<	PC  &     & 783.759 & 481.353 & 332.677 & 248.674 & 144.382 & 98.535 \\
1095	<	SP  & 0.0 & 659.440 & 380.699 & 250.002 & 235.151 & 134.661 & 88.135 \\
1096	<	& 0.1 & 293.849 & 67.772 & 11.609 & 105.090 & 23.813 & 4.369 \\
1097	<	& 0.2 & 5.975 & 0.136 & 0.094 & 5.553 & 1.784 & 1.536 \\
1098	<	& 0.3 & 0.725 & 0.707 & 0.693 & 7.293 & 6.933 & 6.748 \\
1099	<	SF  & 0.0 & 2.238 & 0.713 & 0.292 & 3.290 & 1.090 & 0.416 \\
1100	<	& 0.1 & 2.238 & 0.524 & 0.115 & 3.184 & 0.945 & 0.326 \\
1101	<	& 0.2 & 0.374 & 0.102 & 0.094 & 2.598 & 1.755 & 1.537 \\
1102	<	& 0.3 & 0.721 & 0.707 & 0.693 & 7.322 & 6.933 & 6.748 \\
1103	<	GSC &     & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\
1104	<	RF  &     & 2.091 & 0.403 & 0.113 & 3.583 & 1.071 & 0.399 \\
1105	<	\midrule
1106	<	GSSP  & 0.0 & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\
1107	<	& 0.1 & 1.879 & 0.291 & 0.057 & 3.983 & 1.117 & 0.370 \\
1108	<	& 0.2 & 0.443 & 0.103 & 0.093 & 2.821 & 1.794 & 1.532 \\
1109	<	& 0.3 & 0.728 & 0.694 & 0.692 & 7.387 & 6.942 & 6.748 \\
1110	<	GSSF  & 0.0 & 1.298 & 0.270 & 0.083 & 3.098 & 0.992 & 0.375 \\
1111	<	& 0.1 & 1.296 & 0.210 & 0.044 & 3.055 & 0.922 & 0.330 \\
1112	<	& 0.2 & 0.433 & 0.104 & 0.093 & 2.895 & 1.797 & 1.532 \\
1113	<	& 0.3 & 0.728 & 0.694 & 0.692 & 7.410 & 6.942 & 6.748 \\
1114	<	\bottomrule
1115	<	\end{tabular}
1116	<	\label{tab:spceAng}
1117	<	\end{table}
1118	<
1119	<	The water results parallel the combined results seen in sections
1120	<	\ref{sec:EnergyResults} through \ref{sec:FTDirResults}.  There is good
1121	<	agreement with {\sc spme} in both energetic and dynamic behavior when
1122	<	using the {\sc sf} method with and without damping. The {\sc sp}
1123	<	method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly
1124	<	with cutoff radii greater than 12\AA. Overdamping the electrostatics
1125	<	reduces the agreement between both these methods and {\sc spme}.
1126	<
1127	<	The pure cutoff ({\sc pc}) method performs poorly, again mirroring the
1128	<	observations from the combined results.  In contrast to these results, however, the use of a switching function and group
1129	<	based cutoffs greatly improves the results for these neutral water
1130	<	molecules.  The group switched cutoff ({\sc gsc}) does not mimic the
1131	<	energetics of {\sc spme} as well as the {\sc sp} (with moderate
1132	<	damping) and {\sc sf} methods, but the dynamics are quite good.  The
1133	<	switching functions correct discontinuities in the potential and
1134	<	forces, leading to these improved results.  Such improvements with the
1135	<	use of a switching function have been recognized in previous
1136	<	studies,\cite{Andrea83,Steinbach94} and this proves to be a useful
1137	<	tactic for stably incorporating local area electrostatic effects.
1138	<
1139	<	The reaction field ({\sc rf}) method simply extends upon the results
1140	<	observed in the {\sc gsc} case.  Both methods are similar in form
1141	<	(i.e. neutral groups, switching function), but {\sc rf} incorporates
1142	<	an added effect from the external dielectric. This similarity
1143	<	translates into the same good dynamic results and improved energetic
1144	<	agreement with {\sc spme}.  Though this agreement is not to the level
1145	<	of the moderately damped {\sc sp} and {\sc sf} methods, these results
1146	<	show how incorporating some implicit properties of the surroundings
1147	<	(i.e. $\epsilon_\textrm{S}$) can improve the solvent depiction.
1148	<
1149	<	As a final note for the liquid water system, use of group cutoffs and a
1150	<	switching function leads to noticeable improvements in the {\sc sp}
1151	<	and {\sc sf} methods, primarily in directionality of the force and
1152	<	torque vectors (table \ref{tab:spceAng}). The {\sc sp} method shows
1153	<	significant narrowing of the angle distribution when using little to
1154	<	no damping and only modest improvement for the recommended conditions
1155	<	($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA).  The
1156	<	{\sc sf} method shows modest narrowing across all damping and cutoff
1157	<	ranges of interest.  When overdamping these methods, group cutoffs and
1158	<	the switching function do not improve the force and torque
1159	<	directionalities.
1160	<
1161	<	\subsection{SPC/E Ice I$_\textrm{c}$ Results}\label{sec:IceResults}
1162	<
1163	<	In addition to the disordered molecular system above, the ordered
1164	<	molecular system of ice I$_\textrm{c}$ was also considered.  Ice
1165	<	polymorph could have been used to fit this role; however, ice
1166	<	I$_\textrm{c}$ was chosen because it can form an ideal periodic
1167	<	lattice with the same number of water molecules used in the disordered
1168	<	liquid state case.  The results for the energy gap comparisons and the
1169	<	force and torque vector magnitude comparisons are shown in table
1170	<	\ref{tab:ice}.  The force and torque vector directionality results are
1171	<	displayed separately in table \ref{tab:iceAng}, where the effect of
1172	<	group-based cutoffs and switching functions on the {\sc sp} and {\sc
1173	<	sf} potentials are also displayed.
1174	<
1175	<	\begin{table}[htbp]
1176	<	\centering
1177	<	\caption{REGRESSION RESULTS OF THE ICE I$_\textrm{c}$ SYSTEM FOR
1178	<	$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it
1179	<	middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1180	<
1181	<	\footnotesize
1182	<	\begin{tabular}{@{} ccrrrrrr @{}}
1183	<	\\
1184	<	\toprule
1185	<	\toprule
1186	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1187	<	\cmidrule(lr){3-4}
1188	<	\cmidrule(lr){5-6}
1189	<	\cmidrule(l){7-8}
1190	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1191	<	\midrule
1192	<	PC  &     & 19.897 & 0.047 & -29.214 & 0.048 & -3.771 & 0.001 \\
1193	<	SP  & 0.0 & -0.014 & 0.000 & 2.135 & 0.347 & 0.457 & 0.045 \\
1194	<	& 0.1 & 0.321 & 0.017 & 1.490 & 0.584 & 0.886 & 0.796 \\
1195	<	& 0.2 & 0.896 & 0.872 & 1.011 & 0.998 & 0.997 & 0.999 \\
1196	<	& 0.3 & 0.983 & 0.997 & 0.992 & 0.997 & 0.991 & 0.997 \\
1197	<	SF  & 0.0 & 0.943 & 0.979 & 1.048 & 0.978 & 0.995 & 0.999 \\
1198	<	& 0.1 & 0.948 & 0.979 & 1.044 & 0.983 & 1.000 & 0.999 \\
1199	<	& 0.2 & 0.982 & 0.997 & 0.969 & 0.960 & 0.997 & 0.999 \\
1200	<	& 0.3 & 0.985 & 0.997 & 0.961 & 0.961 & 0.991 & 0.997 \\
1201	<	GSC &     & 0.983 & 0.985 & 0.966 & 0.994 & 1.003 & 0.999 \\
1202	<	RF  &     & 0.924 & 0.944 & 0.990 & 0.996 & 0.991 & 0.998 \\
1203	<	\midrule
1204	<	PC  &     & -4.375 & 0.000 & 6.781 & 0.000 & -3.369 & 0.000 \\
1205	<	SP  & 0.0 & 0.515 & 0.164 & 0.856 & 0.426 & 0.743 & 0.478 \\
1206	<	& 0.1 & 0.696 & 0.405 & 0.977 & 0.817 & 0.974 & 0.964 \\
1207	<	& 0.2 & 0.981 & 0.980 & 1.001 & 1.000 & 1.000 & 1.000 \\
1208	<	& 0.3 & 0.996 & 0.998 & 0.997 & 0.999 & 0.997 & 0.999 \\
1209	<	SF  & 0.0 & 0.991 & 0.995 & 1.003 & 0.998 & 0.999 & 1.000 \\
1210	<	& 0.1 & 0.992 & 0.995 & 1.003 & 0.998 & 1.000 & 1.000 \\
1211	<	& 0.2 & 0.998 & 0.998 & 0.981 & 0.962 & 1.000 & 1.000 \\
1212	<	& 0.3 & 0.996 & 0.998 & 0.976 & 0.957 & 0.997 & 0.999 \\
1213	<	GSC &     & 0.997 & 0.996 & 0.998 & 0.999 & 1.000 & 1.000 \\
1214	<	RF  &     & 0.988 & 0.989 & 1.000 & 0.999 & 1.000 & 1.000 \\
1215	<	\midrule
1216	<	PC  &     & -6.367 & 0.000 & -3.552 & 0.000 & -3.447 & 0.000 \\
1217	<	SP  & 0.0 & 0.643 & 0.409 & 0.833 & 0.607 & 0.961 & 0.805 \\
1218	<	& 0.1 & 0.791 & 0.683 & 0.957 & 0.914 & 1.000 & 0.989 \\
1219	<	& 0.2 & 0.974 & 0.991 & 0.993 & 0.998 & 0.993 & 0.998 \\
1220	<	& 0.3 & 0.976 & 0.992 & 0.977 & 0.992 & 0.977 & 0.992 \\
1221	<	SF  & 0.0 & 0.979 & 0.997 & 0.992 & 0.999 & 0.994 & 1.000 \\
1222	<	& 0.1 & 0.984 & 0.997 & 0.996 & 0.999 & 0.998 & 1.000 \\
1223	<	& 0.2 & 0.991 & 0.997 & 0.974 & 0.958 & 0.993 & 0.998 \\
1224	<	& 0.3 & 0.977 & 0.992 & 0.956 & 0.948 & 0.977 & 0.992 \\
1225	<	GSC &     & 0.999 & 0.997 & 0.996 & 0.999 & 1.002 & 1.000 \\
1226	<	RF  &     & 0.994 & 0.997 & 0.997 & 0.999 & 1.000 & 1.000 \\
1227	<	\bottomrule
1228	<	\end{tabular}
1229	<	\label{tab:ice}
1230	<	\end{table}
1231	<
1232	<	\begin{table}[htbp]
1233	<	\centering
1234	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1235	<	OF THE FORCE AND TORQUE VECTORS IN THE ICE I$_\textrm{c}$ SYSTEM}
1236	<
1237	<	\footnotesize
1238	<	\begin{tabular}{@{} ccrrrrrr @{}}
1239	<	\\
1240	<	\toprule
1241	<	\toprule
1242	<	& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque
1243	<	$\sigma^2$} \\
1244	<	\cmidrule(lr){3-5}
1245	<	\cmidrule(l){6-8}
1246	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\
1247	<	\midrule
1248	<	PC  &     & 2128.921 & 603.197 & 715.579 & 329.056 & 221.397 & 81.042 \\
1249	<	SP  & 0.0 & 1429.341 & 470.320 & 447.557 & 301.678 & 197.437 & 73.840 \\
1250	<	& 0.1 & 590.008 & 107.510 & 18.883 & 118.201 & 32.472 & 3.599 \\
1251	<	& 0.2 & 10.057 & 0.105 & 0.038 & 2.875 & 0.572 & 0.518 \\
1252	<	& 0.3 & 0.245 & 0.260 & 0.262 & 2.365 & 2.396 & 2.327 \\
1253	<	SF  & 0.0 & 1.745 & 1.161 & 0.212 & 1.135 & 0.426 & 0.155 \\
1254	<	& 0.1 & 1.721 & 0.868 & 0.082 & 1.118 & 0.358 & 0.118 \\
1255	<	& 0.2 & 0.201 & 0.040 & 0.038 & 0.786 & 0.555 & 0.518 \\
1256	<	& 0.3 & 0.241 & 0.260 & 0.262 & 2.368 & 2.400 & 2.327 \\
1257	<	GSC &     & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\
1258	<	RF  &     & 2.887 & 0.217 & 0.107 & 1.006 & 0.281 & 0.085 \\
1259	<	\midrule
1260	<	GSSP  & 0.0 & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\
1261	<	& 0.1 & 1.341 & 0.123 & 0.037 & 0.835 & 0.234 & 0.085 \\
1262	<	& 0.2 & 0.558 & 0.040 & 0.037 & 0.823 & 0.557 & 0.519 \\
1263	<	& 0.3 & 0.250 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\
1264	<	GSSF  & 0.0 & 2.124 & 0.132 & 0.069 & 0.919 & 0.263 & 0.099 \\
1265	<	& 0.1 & 2.165 & 0.101 & 0.035 & 0.895 & 0.244 & 0.096 \\
1266	<	& 0.2 & 0.706 & 0.040 & 0.037 & 0.870 & 0.559 & 0.519 \\
1267	<	& 0.3 & 0.251 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\
1268	<	\bottomrule
1269	<	\end{tabular}
1270	<	\label{tab:iceAng}
1271	<	\end{table}
1272	<
1273	<	Highly ordered systems are a difficult test for the pairwise methods
1274	<	in that they lack the implicit periodicity of the Ewald summation.  As
1275	<	expected, the energy gap agreement with {\sc spme} is reduced for the
1276	<	{\sc sp} and {\sc sf} methods with parameters that were ideal for the
1277	<	disordered liquid system.  Moving to higher $R_\textrm{c}$ helps
1278	<	improve the agreement, though at an increase in computational cost.
1279	<	The dynamics of this crystalline system (both in magnitude and
1280	<	direction) are little affected. Both methods still reproduce the Ewald
1281	<	behavior with the same parameter recommendations from the previous
1282	<	section.
1283	<
1284	<	It is also worth noting that {\sc rf} exhibits improved energy gap
1285	<	results over the liquid water system.  One possible explanation is
1286	<	that the ice I$_\textrm{c}$ crystal is ordered such that the net
1287	<	dipole moment of the crystal is zero.  With $\epsilon_\textrm{S} =
1288	<	\infty$, the reaction field incorporates this structural organization
1289	<	by actively enforcing a zeroed dipole moment within each cutoff
1290	<	sphere.
1291	<
1292	<	\subsection{NaCl Melt Results}\label{sec:SaltMeltResults}
1293	<
1294	<	A high temperature NaCl melt was tested to gauge the accuracy of the
1295	<	pairwise summation methods in a disordered system of charges. The
1296	<	results for the energy gap comparisons and the force vector magnitude
1297	<	comparisons are shown in table \ref{tab:melt}.  The force vector
1298	<	directionality results are displayed separately in table
1299	<	\ref{tab:meltAng}.
1300	<
1301	<	\begin{table}[htbp]
1302	<	\centering
1303	<	\caption{REGRESSION RESULTS OF THE MOLTEN SODIUM CHLORIDE SYSTEM FOR
1304	<	$\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES ({\it
1305	<	lower})}
1306	<
1307	<	\footnotesize
1308	<	\begin{tabular}{@{} ccrrrrrr @{}}
1309	<	\\
1310	<	\toprule
1311	<	\toprule
1312	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1313	<	\cmidrule(lr){3-4}
1314	<	\cmidrule(lr){5-6}
1315	<	\cmidrule(l){7-8}
1316	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1317	<	\midrule
1318	<	PC  &     & -0.008 & 0.000 & -0.049 & 0.005 & -0.136 & 0.020 \\
1319	<	SP  & 0.0 & 0.928 & 0.996 & 0.931 & 0.998 & 0.950 & 0.999 \\
1320	<	& 0.1 & 0.977 & 0.998 & 0.998 & 1.000 & 0.997 & 1.000 \\
1321	<	& 0.2 & 0.960 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\
1322	<	& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\
1323	<	SF  & 0.0 & 0.996 & 1.000 & 0.995 & 1.000 & 0.997 & 1.000 \\
1324	<	& 0.1 & 1.021 & 1.000 & 1.024 & 1.000 & 1.007 & 1.000 \\
1325	<	& 0.2 & 0.966 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\
1326	<	& 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\
1327	<	\midrule
1328	<	PC  &     & 1.103 & 0.000 & 0.989 & 0.000 & 0.802 & 0.000 \\
1329	<	SP  & 0.0 & 0.973 & 0.981 & 0.975 & 0.988 & 0.979 & 0.992 \\
1330	<	& 0.1 & 0.987 & 0.992 & 0.993 & 0.998 & 0.997 & 0.999 \\
1331	<	& 0.2 & 0.993 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\
1332	<	& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\
1333	<	SF  & 0.0 & 0.996 & 0.997 & 0.997 & 0.999 & 0.998 & 1.000 \\
1334	<	& 0.1 & 1.000 & 0.997 & 1.001 & 0.999 & 1.000 & 1.000 \\
1335	<	& 0.2 & 0.994 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\
1336	<	& 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\
1337	<	\bottomrule
1338	<	\end{tabular}
1339	<	\label{tab:melt}
1340	<	\end{table}
1341	<
1342	<	\begin{table}[htbp]
1343	<	\centering
1344	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1345	<	OF THE FORCE VECTORS IN THE MOLTEN SODIUM CHLORIDE SYSTEM}
1346	<
1347	<	\footnotesize
1348	<	\begin{tabular}{@{} ccrrrrrr @{}}
1349	<	\\
1350	<	\toprule
1351	<	\toprule
1352	<	& & \multicolumn{3}{c}{Force $\sigma^2$} \\
1353	<	\cmidrule(lr){3-5}
1354	<	\cmidrule(l){6-8}
1355	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\
1356	<	\midrule
1357	<	PC  &     & 13.294 & 8.035 & 5.366 \\
1358	<	SP  & 0.0 & 13.316 & 8.037 & 5.385 \\
1359	<	& 0.1 & 5.705 & 1.391 & 0.360 \\
1360	<	& 0.2 & 2.415 & 7.534 & 13.927 \\
1361	<	& 0.3 & 23.769 & 67.306 & 57.252 \\
1362	<	SF  & 0.0 & 1.693 & 0.603 & 0.256 \\
1363	<	& 0.1 & 1.687 & 0.653 & 0.272 \\
1364	<	& 0.2 & 2.598 & 7.523 & 13.930 \\
1365	<	& 0.3 & 23.734 & 67.305 & 57.252 \\
1366	<	\bottomrule
1367	<	\end{tabular}
1368	<	\label{tab:meltAng}
1369	<	\end{table}
1370	<
1371	<	The molten NaCl system shows more sensitivity to the electrostatic
1372	<	damping than the water systems. The most noticeable point is that the
1373	<	undamped {\sc sf} method does very well at replicating the {\sc spme}
1374	<	configurational energy differences and forces. Light damping appears
1375	<	to minimally improve the dynamics, but this comes with a deterioration
1376	<	of the energy gap results. In contrast, this light damping improves
1377	<	the {\sc sp} energy gaps and forces. Moderate and heavy electrostatic
1378	<	damping reduce the agreement with {\sc spme} for both methods. From
1379	<	these observations, the undamped {\sc sf} method is the best choice
1380	<	for disordered systems of charges.
1381	<
1382	<	\subsection{NaCl Crystal Results}\label{sec:SaltCrystalResults}
1383	<
1384	<	Similar to the use of ice I$_\textrm{c}$ to investigate the role of
1385	<	order in molecular systems on the effectiveness of the pairwise
1386	<	methods, the 1000K NaCl crystal system was used to investigate the
1387	<	accuracy of the pairwise summation methods in an ordered system of
1388	<	charged particles. The results for the energy gap comparisons and the
1389	<	force vector magnitude comparisons are shown in table \ref{tab:salt}.
1390	<	The force vector directionality results are displayed separately in
1391	<	table \ref{tab:saltAng}.
1392	<
1393	<	\begin{table}[htbp]
1394	<	\centering
1395	<	\caption{REGRESSION RESULTS OF THE CRYSTALLINE SODIUM CHLORIDE
1396	<	SYSTEM FOR $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES
1397	<	({\it lower})}
1398	<
1399	<	\footnotesize
1400	<	\begin{tabular}{@{} ccrrrrrr @{}}
1401	<	\\
1402	<	\toprule
1403	<	\toprule
1404	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1405	<	\cmidrule(lr){3-4}
1406	<	\cmidrule(lr){5-6}
1407	<	\cmidrule(l){7-8}
1408	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1409	<	\midrule
1410	<	PC  &     & -20.241 & 0.228 & -20.248 & 0.229 & -20.239 & 0.228 \\
1411	<	SP  & 0.0 & 1.039 & 0.733 & 2.037 & 0.565 & 1.225 & 0.743 \\
1412	<	& 0.1 & 1.049 & 0.865 & 1.424 & 0.784 & 1.029 & 0.980 \\
1413	<	& 0.2 & 0.982 & 0.976 & 0.969 & 0.980 & 0.960 & 0.980 \\
1414	<	& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.945 \\
1415	<	SF  & 0.0 & 1.041 & 0.967 & 0.994 & 0.989 & 0.957 & 0.993 \\
1416	<	& 0.1 & 1.050 & 0.968 & 0.996 & 0.991 & 0.972 & 0.995 \\
1417	<	& 0.2 & 0.982 & 0.975 & 0.959 & 0.980 & 0.960 & 0.980 \\
1418	<	& 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.944 \\
1419	<	\midrule
1420	<	PC  &     & 0.795 & 0.000 & 0.792 & 0.000 & 0.793 & 0.000 \\
1421	<	SP  & 0.0 & 0.916 & 0.829 & 1.086 & 0.791 & 1.010 & 0.936 \\
1422	<	& 0.1 & 0.958 & 0.917 & 1.049 & 0.943 & 1.001 & 0.995 \\
1423	<	& 0.2 & 0.981 & 0.981 & 0.982 & 0.984 & 0.981 & 0.984 \\
1424	<	& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\
1425	<	SF  & 0.0 & 1.002 & 0.983 & 0.997 & 0.994 & 0.991 & 0.997 \\
1426	<	& 0.1 & 1.003 & 0.984 & 0.996 & 0.995 & 0.993 & 0.997 \\
1427	<	& 0.2 & 0.983 & 0.980 & 0.981 & 0.984 & 0.981 & 0.984 \\
1428	<	& 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\
1429	<	\bottomrule
1430	<	\end{tabular}
1431	<	\label{tab:salt}
1432	<	\end{table}
1433	<
1434	<	\begin{table}[htbp]
1435	<	\centering
1436	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1437	<	DISTRIBUTIONS OF THE FORCE VECTORS IN THE CRYSTALLINE SODIUM CHLORIDE
1438	<	SYSTEM}
1439	<
1440	<	\footnotesize
1441	<	\begin{tabular}{@{} ccrrrrrr @{}}
1442	<	\\
1443	<	\toprule
1444	<	\toprule
1445	<	& & \multicolumn{3}{c}{Force $\sigma^2$} \\
1446	<	\cmidrule(lr){3-5}
1447	<	\cmidrule(l){6-8}
1448	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\
1449	<	\midrule
1450	<	PC  &     & 111.945 & 111.824 & 111.866 \\
1451	<	SP  & 0.0 & 112.414 & 152.215 & 38.087 \\
1452	<	& 0.1 & 52.361 & 42.574 & 2.819 \\
1453	<	& 0.2 & 10.847 & 9.709 & 9.686 \\
1454	<	& 0.3 & 31.128 & 31.104 & 31.029 \\
1455	<	SF  & 0.0 & 10.025 & 3.555 & 1.648 \\
1456	<	& 0.1 & 9.462 & 3.303 & 1.721 \\
1457	<	& 0.2 & 11.454 & 9.813 & 9.701 \\
1458	<	& 0.3 & 31.120 & 31.105 & 31.029 \\
1459	<	\bottomrule
1460	<	\end{tabular}
1461	<	\label{tab:saltAng}
1462	<	\end{table}
1463	<
1464	<	The crystalline NaCl system is the most challenging test case for the
1465	<	pairwise summation methods, as evidenced by the results in tables
1466	<	\ref{tab:salt} and \ref{tab:saltAng}. The undamped and weakly damped
1467	<	{\sc sf} methods seem to be the best choices. These methods match well
1468	<	with {\sc spme} across the energy gap, force magnitude, and force
1469	<	directionality tests.  The {\sc sp} method struggles in all cases,
1470	<	with the exception of good dynamics reproduction when using weak
1471	<	electrostatic damping with a large cutoff radius.
1472	<
1473	<	The moderate electrostatic damping case is not as good as we would
1474	<	expect given the long-time dynamics results observed for this system
1475	<	(see section \ref{sec:LongTimeDynamics}). Since the data tabulated in
1476	<	tables \ref{tab:salt} and \ref{tab:saltAng} are a test of
1477	<	instantaneous dynamics, this indicates that good long-time dynamics
1478	<	comes in part at the expense of short-time dynamics.
1479	<
1480	<	\subsection{0.11M NaCl Solution Results}
1481	<
1482	<	In an effort to bridge the charged atomic and neutral molecular
1483	<	systems, Na$^+$ and Cl$^-$ ion charge defects were incorporated into
1484	<	the liquid water system. This low ionic strength system consists of 4
1485	<	ions in the 1000 SPC/E water solvent ($\approx$0.11 M). The results
1486	<	for the energy gap comparisons and the force and torque vector
1487	<	magnitude comparisons are shown in table \ref{tab:solnWeak}.  The
1488	<	force and torque vector directionality results are displayed
1489	<	separately in table \ref{tab:solnWeakAng}, where the effect of
1490	<	group-based cutoffs and switching functions on the {\sc sp} and {\sc
1491	<	sf} potentials are investigated.
1492	<
1493	<	\begin{table}[htbp]
1494	<	\centering
1495	<	\caption{REGRESSION RESULTS OF THE WEAK SODIUM CHLORIDE SOLUTION
1496	<	SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES
1497	<	({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1498	<
1499	<	\footnotesize
1500	<	\begin{tabular}{@{} ccrrrrrr @{}}
1501	<	\\
1502	<	\toprule
1503	<	\toprule
1504	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1505	<	\cmidrule(lr){3-4}
1506	<	\cmidrule(lr){5-6}
1507	<	\cmidrule(l){7-8}
1508	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1509	<	\midrule
1510	<	PC  &     & 0.247 & 0.000 & -1.103 & 0.001 & 5.480 & 0.015 \\
1511	<	SP  & 0.0 & 0.935 & 0.388 & 0.984 & 0.541 & 1.010 & 0.685 \\
1512	<	& 0.1 & 0.951 & 0.603 & 0.993 & 0.875 & 1.001 & 0.979 \\
1513	<	& 0.2 & 0.969 & 0.968 & 0.996 & 0.997 & 0.994 & 0.997 \\
1514	<	& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\
1515	<	SF  & 0.0 & 0.963 & 0.971 & 0.989 & 0.996 & 0.991 & 0.998 \\
1516	<	& 0.1 & 0.970 & 0.971 & 0.995 & 0.997 & 0.997 & 0.999 \\
1517	<	& 0.2 & 0.972 & 0.975 & 0.996 & 0.997 & 0.994 & 0.997 \\
1518	<	& 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\
1519	<	GSC &     & 0.964 & 0.731 & 0.984 & 0.704 & 1.005 & 0.770 \\
1520	<	RF  &     & 0.968 & 0.605 & 0.974 & 0.541 & 1.014 & 0.614 \\
1521	<	\midrule
1522	<	PC  &     & 1.354 & 0.000 & -1.190 & 0.000 & -0.314 & 0.000 \\
1523	<	SP  & 0.0 & 0.720 & 0.338 & 0.808 & 0.523 & 0.860 & 0.643 \\
1524	<	& 0.1 & 0.839 & 0.583 & 0.955 & 0.882 & 0.992 & 0.978 \\
1525	<	& 0.2 & 0.995 & 0.987 & 0.999 & 1.000 & 0.999 & 1.000 \\
1526	<	& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\
1527	<	SF  & 0.0 & 0.998 & 0.994 & 1.000 & 0.998 & 1.000 & 0.999 \\
1528	<	& 0.1 & 0.997 & 0.994 & 1.000 & 0.999 & 1.000 & 1.000 \\
1529	<	& 0.2 & 0.999 & 0.998 & 0.999 & 1.000 & 0.999 & 1.000 \\
1530	<	& 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\
1531	<	GSC &     & 0.995 & 0.990 & 0.998 & 0.997 & 0.998 & 0.996 \\
1532	<	RF  &     & 0.998 & 0.993 & 0.999 & 0.998 & 0.999 & 0.996 \\
1533	<	\midrule
1534	<	PC  &     & 2.437 & 0.000 & -1.872 & 0.000 & 2.138 & 0.000 \\
1535	<	SP  & 0.0 & 0.838 & 0.525 & 0.901 & 0.686 & 0.932 & 0.779 \\
1536	<	& 0.1 & 0.914 & 0.733 & 0.979 & 0.932 & 0.995 & 0.987 \\
1537	<	& 0.2 & 0.977 & 0.969 & 0.988 & 0.990 & 0.989 & 0.990 \\
1538	<	& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\
1539	<	SF  & 0.0 & 0.969 & 0.977 & 0.987 & 0.996 & 0.993 & 0.998 \\
1540	<	& 0.1 & 0.975 & 0.978 & 0.993 & 0.996 & 0.997 & 0.998 \\
1541	<	& 0.2 & 0.976 & 0.973 & 0.988 & 0.990 & 0.989 & 0.990 \\
1542	<	& 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\
1543	<	GSC &     & 0.980 & 0.959 & 0.990 & 0.983 & 0.992 & 0.989 \\
1544	<	RF  &     & 0.984 & 0.975 & 0.996 & 0.995 & 0.998 & 0.998 \\
1545	<	\bottomrule
1546	<	\end{tabular}
1547	<	\label{tab:solnWeak}
1548	<	\end{table}
1549	<
1550	<	\begin{table}[htbp]
1551	<	\centering
1552	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1553	<	DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE WEAK SODIUM
1554	<	CHLORIDE SOLUTION SYSTEM}
1555	<
1556	<	\footnotesize
1557	<	\begin{tabular}{@{} ccrrrrrr @{}}
1558	<	\\
1559	<	\toprule
1560	<	\toprule
1561	<	& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\
1562	<	\cmidrule(lr){3-5}
1563	<	\cmidrule(l){6-8}
1564	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\
1565	<	\midrule
1566	<	PC  &     & 882.863 & 510.435 & 344.201 & 277.691 & 154.231 & 100.131 \\
1567	<	SP  & 0.0 & 732.569 & 405.704 & 257.756 & 261.445 & 142.245 & 91.497 \\
1568	<	& 0.1 & 329.031 & 70.746 & 12.014 & 118.496 & 25.218 & 4.711 \\
1569	<	& 0.2 & 6.772 & 0.153 & 0.118 & 9.780 & 2.101 & 2.102 \\
1570	<	& 0.3 & 0.951 & 0.774 & 0.784 & 12.108 & 7.673 & 7.851 \\
1571	<	SF  & 0.0 & 2.555 & 0.762 & 0.313 & 6.590 & 1.328 & 0.558 \\
1572	<	& 0.1 & 2.561 & 0.560 & 0.123 & 6.464 & 1.162 & 0.457 \\
1573	<	& 0.2 & 0.501 & 0.118 & 0.118 & 5.698 & 2.074 & 2.099 \\
1574	<	& 0.3 & 0.943 & 0.774 & 0.784 & 12.118 & 7.674 & 7.851 \\
1575	<	GSC &     & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\
1576	<	RF  &     & 2.415 & 0.452 & 0.130 & 6.915 & 1.423 & 0.507 \\
1577	<	\midrule
1578	<	GSSP  & 0.0 & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\
1579	<	& 0.1 & 2.251 & 0.324 & 0.064 & 7.628 & 1.639 & 0.497 \\
1580	<	& 0.2 & 0.590 & 0.118 & 0.116 & 6.080 & 2.096 & 2.103 \\
1581	<	& 0.3 & 0.953 & 0.759 & 0.780 & 12.347 & 7.683 & 7.849 \\
1582	<	GSSF  & 0.0 & 1.541 & 0.301 & 0.096 & 6.407 & 1.316 & 0.496 \\
1583	<	& 0.1 & 1.541 & 0.237 & 0.050 & 6.356 & 1.202 & 0.457 \\
1584	<	& 0.2 & 0.568 & 0.118 & 0.116 & 6.166 & 2.105 & 2.105 \\
1585	<	& 0.3 & 0.954 & 0.759 & 0.780 & 12.337 & 7.684 & 7.849 \\
1586	<	\bottomrule
1587	<	\end{tabular}
1588	<	\label{tab:solnWeakAng}
1589	<	\end{table}
1590	<
1591	<	Because this system is a perturbation of the pure liquid water system,
1592	<	comparisons are best drawn between these two sets. The {\sc sp} and
1593	<	{\sc sf} methods are not significantly affected by the inclusion of a
1594	<	few ions. The aspect of cutoff sphere neutralization aids in the
1595	<	smooth incorporation of these ions; thus, all of the observations
1596	<	regarding these methods carry over from section
1597	<	\ref{sec:WaterResults}. The differences between these systems are more
1598	<	visible for the {\sc rf} method. Though good force agreement is still
1599	<	maintained, the energy gaps show a significant increase in the scatter
1600	<	of the data.
1601	<
1602	<	\subsection{1.1M NaCl Solution Results}
1603	<
1604	<	The bridging of the charged atomic and neutral molecular systems was
1605	<	further developed by considering a high ionic strength system
1606	<	consisting of 40 ions in the 1000 SPC/E water solvent ($\approx$1.1
1607	<	M). The results for the energy gap comparisons and the force and
1608	<	torque vector magnitude comparisons are shown in table
1609	<	\ref{tab:solnStr}.  The force and torque vector directionality
1610	<	results are displayed separately in table \ref{tab:solnStrAng}, where
1611	<	the effect of group-based cutoffs and switching functions on the {\sc
1612	<	sp} and {\sc sf} potentials are investigated.
1613	<
1614	<	\begin{table}[htbp]
1615	<	\centering
1616	<	\caption{REGRESSION RESULTS OF THE STRONG SODIUM CHLORIDE SOLUTION
1617	<	SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES
1618	<	({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1619	<
1620	<	\footnotesize
1621	<	\begin{tabular}{@{} ccrrrrrr @{}}
1622	<	\\
1623	<	\toprule
1624	<	\toprule
1625	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1626	<	\cmidrule(lr){3-4}
1627	<	\cmidrule(lr){5-6}
1628	<	\cmidrule(l){7-8}
1629	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1630	<	\midrule
1631	<	PC  &     & -0.081 & 0.000 & 0.945 & 0.001 & 0.073 & 0.000 \\
1632	<	SP  & 0.0 & 0.978 & 0.469 & 0.996 & 0.672 & 0.975 & 0.668 \\
1633	<	& 0.1 & 0.944 & 0.645 & 0.997 & 0.886 & 0.991 & 0.978 \\
1634	<	& 0.2 & 0.873 & 0.896 & 0.985 & 0.993 & 0.980 & 0.993 \\
1635	<	& 0.3 & 0.831 & 0.860 & 0.960 & 0.979 & 0.955 & 0.977 \\
1636	<	SF  & 0.0 & 0.858 & 0.905 & 0.985 & 0.970 & 0.990 & 0.998 \\
1637	<	& 0.1 & 0.865 & 0.907 & 0.992 & 0.974 & 0.994 & 0.999 \\
1638	<	& 0.2 & 0.862 & 0.894 & 0.985 & 0.993 & 0.980 & 0.993 \\
1639	<	& 0.3 & 0.831 & 0.859 & 0.960 & 0.979 & 0.955 & 0.977 \\
1640	<	GSC &     & 1.985 & 0.152 & 0.760 & 0.031 & 1.106 & 0.062 \\
1641	<	RF  &     & 2.414 & 0.116 & 0.813 & 0.017 & 1.434 & 0.047 \\
1642	<	\midrule
1643	<	PC  &     & -7.028 & 0.000 & -9.364 & 0.000 & 0.925 & 0.865 \\
1644	<	SP  & 0.0 & 0.701 & 0.319 & 0.909 & 0.773 & 0.861 & 0.665 \\
1645	<	& 0.1 & 0.824 & 0.565 & 0.970 & 0.930 & 0.990 & 0.979 \\
1646	<	& 0.2 & 0.988 & 0.981 & 0.995 & 0.998 & 0.991 & 0.998 \\
1647	<	& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\
1648	<	SF  & 0.0 & 0.993 & 0.988 & 0.992 & 0.984 & 0.998 & 0.999 \\
1649	<	& 0.1 & 0.993 & 0.989 & 0.993 & 0.986 & 0.998 & 1.000 \\
1650	<	& 0.2 & 0.993 & 0.992 & 0.995 & 0.998 & 0.991 & 0.998 \\
1651	<	& 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\
1652	<	GSC &     & 0.964 & 0.897 & 0.970 & 0.917 & 0.925 & 0.865 \\
1653	<	RF  &     & 0.994 & 0.864 & 0.988 & 0.865 & 0.980 & 0.784 \\
1654	<	\midrule
1655	<	PC  &     & -2.212 & 0.000 & -0.588 & 0.000 & 0.953 & 0.925 \\
1656	<	SP  & 0.0 & 0.800 & 0.479 & 0.930 & 0.804 & 0.924 & 0.759 \\
1657	<	& 0.1 & 0.883 & 0.694 & 0.976 & 0.942 & 0.993 & 0.986 \\
1658	<	& 0.2 & 0.952 & 0.943 & 0.980 & 0.984 & 0.980 & 0.983 \\
1659	<	& 0.3 & 0.914 & 0.909 & 0.943 & 0.948 & 0.944 & 0.946 \\
1660	<	SF  & 0.0 & 0.945 & 0.953 & 0.980 & 0.984 & 0.991 & 0.998 \\
1661	<	& 0.1 & 0.951 & 0.954 & 0.987 & 0.986 & 0.995 & 0.998 \\
1662	<	& 0.2 & 0.951 & 0.946 & 0.980 & 0.984 & 0.980 & 0.983 \\
1663	<	& 0.3 & 0.914 & 0.908 & 0.943 & 0.948 & 0.944 & 0.946 \\
1664	<	GSC &     & 0.882 & 0.818 & 0.939 & 0.902 & 0.953 & 0.925 \\
1665	<	RF  &     & 0.949 & 0.939 & 0.988 & 0.988 & 0.992 & 0.993 \\
1666	<	\bottomrule
1667	<	\end{tabular}
1668	<	\label{tab:solnStr}
1669	<	\end{table}
1670	<
1671	<	\begin{table}[htbp]
1672	<	\centering
1673	<	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1674	<	OF THE FORCE AND TORQUE VECTORS IN THE STRONG SODIUM CHLORIDE SOLUTION
1675	<	SYSTEM}
1676	<
1677	<	\footnotesize
1678	<	\begin{tabular}{@{} ccrrrrrr @{}}
1679	<	\\
1680	<	\toprule
1681	<	\toprule
1682	<	& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\
1683	<	\cmidrule(lr){3-5}
1684	<	\cmidrule(l){6-8}
1685	<	Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\
1686	<	\midrule
1687	<	PC  &     & 957.784 & 513.373 & 2.260 & 340.043 & 179.443 & 13.079 \\
1688	<	SP  & 0.0 & 786.244 & 139.985 & 259.289 & 311.519 & 90.280 & 105.187 \\
1689	<	& 0.1 & 354.697 & 38.614 & 12.274 & 144.531 & 23.787 & 5.401 \\
1690	<	& 0.2 & 7.674 & 0.363 & 0.215 & 16.655 & 3.601 & 3.634 \\
1691	<	& 0.3 & 1.745 & 1.456 & 1.449 & 23.669 & 14.376 & 14.240 \\
1692	<	SF  & 0.0 & 3.282 & 8.567 & 0.369 & 11.904 & 6.589 & 0.717 \\
1693	<	& 0.1 & 3.263 & 7.479 & 0.142 & 11.634 & 5.750 & 0.591 \\
1694	<	& 0.2 & 0.686 & 0.324 & 0.215 & 10.809 & 3.580 & 3.635 \\
1695	<	& 0.3 & 1.749 & 1.456 & 1.449 & 23.635 & 14.375 & 14.240 \\
1696	<	GSC &     & 6.181 & 2.904 & 2.263 & 44.349 & 19.442 & 12.873 \\
1697	<	RF  &     & 3.891 & 0.847 & 0.323 & 18.628 & 3.995 & 2.072 \\
1698	<	\midrule
1699	<	GSSP  & 0.0 & 6.197 & 2.929 & 2.290 & 44.441 & 19.442 & 12.873 \\
1700	<	& 0.1 & 4.688 & 1.064 & 0.260 & 31.208 & 6.967 & 2.303 \\
1701	<	& 0.2 & 1.021 & 0.218 & 0.213 & 14.425 & 3.629 & 3.649 \\
1702	<	& 0.3 & 1.752 & 1.454 & 1.451 & 23.540 & 14.390 & 14.245 \\
1703	<	GSSF  & 0.0 & 2.494 & 0.546 & 0.217 & 16.391 & 3.230 & 1.613 \\
1704	<	& 0.1 & 2.448 & 0.429 & 0.106 & 16.390 & 2.827 & 1.159 \\
1705	<	& 0.2 & 0.899 & 0.214 & 0.213 & 13.542 & 3.583 & 3.645 \\
1706	<	& 0.3 & 1.752 & 1.454 & 1.451 & 23.587 & 14.390 & 14.245 \\
1707	<	\bottomrule
1708	<	\end{tabular}
1709	<	\label{tab:solnStrAng}
1710	<	\end{table}
1711	<
1712	<	The {\sc rf} method struggles with the jump in ionic strength. The
1713	<	configuration energy differences degrade to unusable levels while the
1714	<	forces and torques show a more modest reduction in the agreement with
1715	<	{\sc spme}. The {\sc rf} method was designed for homogeneous systems,
1716	<	and this attribute is apparent in these results.
1717	<
1718	<	The {\sc sp} and {\sc sf} methods require larger cutoffs to maintain
1719	<	their agreement with {\sc spme}. With these results, we still
1720	<	recommend undamped to moderate damping for the {\sc sf} method and
1721	<	moderate damping for the {\sc sp} method, both with cutoffs greater
1722	<	than 12\AA.
1723	<
1724	<	\subsection{6\AA\ Argon Sphere in SPC/E Water Results}
1725	<
1726	<	The final model system studied was a 6\AA\ sphere of Argon solvated
1727	<	by SPC/E water. This serves as a test case of a specifically sized
1728	<	electrostatic defect in a disordered molecular system. The results for
1729	<	the energy gap comparisons and the force and torque vector magnitude
1730	<	comparisons are shown in table \ref{tab:argon}.  The force and torque
1731	<	vector directionality results are displayed separately in table
1732	<	\ref{tab:argonAng}, where the effect of group-based cutoffs and
1733	<	switching functions on the {\sc sp} and {\sc sf} potentials are
1734	<	investigated.
1735	<
1736	<	\begin{table}[htbp]
1737	<	\centering
1738	<	\caption{REGRESSION RESULTS OF THE 6\AA\ ARGON SPHERE IN LIQUID
1739	<	WATER SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR
1740	<	MAGNITUDES ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1741	<
1742	<	\footnotesize
1743	<	\begin{tabular}{@{} ccrrrrrr @{}}
1744	<	\\
1745	<	\toprule
1746	<	\toprule
1747	<	& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\
1748	<	\cmidrule(lr){3-4}
1749	<	\cmidrule(lr){5-6}
1750	<	\cmidrule(l){7-8}
1751	<	Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1752	<	\midrule
1753	<	PC  &     & 2.320 & 0.008 & -0.650 & 0.001 & 3.848 & 0.029 \\
1754	<	SP  & 0.0 & 1.053 & 0.711 & 0.977 & 0.820 & 0.974 & 0.882 \\
1755	<	& 0.1 & 1.032 & 0.846 & 0.989 & 0.965 & 0.992 & 0.994 \\
1756	<	& 0.2 & 0.993 & 0.995 & 0.982 & 0.998 & 0.986 & 0.998 \\
1757	<	& 0.3 & 0.968 & 0.995 & 0.954 & 0.992 & 0.961 & 0.994 \\
1758	<	SF  & 0.0 & 0.982 & 0.996 & 0.992 & 0.999 & 0.993 & 1.000 \\
1759	<	& 0.1 & 0.987 & 0.996 & 0.996 & 0.999 & 0.997 & 1.000 \\
1760	<	& 0.2 & 0.989 & 0.998 & 0.984 & 0.998 & 0.989 & 0.998 \\
1761	<	& 0.3 & 0.971 & 0.995 & 0.957 & 0.992 & 0.965 & 0.994 \\
1762	<	GSC &     & 1.002 & 0.983 & 0.992 & 0.973 & 0.996 & 0.971 \\
1763	<	RF  &     & 0.998 & 0.995 & 0.999 & 0.998 & 0.998 & 0.998 \\
1764	<	\midrule
1765	<	PC  &     & -36.559 & 0.002 & -44.917 & 0.004 & -52.945 & 0.006 \\
1766	<	SP  & 0.0 & 0.890 & 0.786 & 0.927 & 0.867 & 0.949 & 0.909 \\
1767	<	& 0.1 & 0.942 & 0.895 & 0.984 & 0.974 & 0.997 & 0.995 \\
1768	<	& 0.2 & 0.999 & 0.997 & 1.000 & 1.000 & 1.000 & 1.000 \\
1769	<	& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\
1770	<	SF  & 0.0 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1771	<	& 0.1 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1772	<	& 0.2 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
1773	<	& 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\
1774	<	GSC &     & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1775	<	RF  &     & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1776	<	\midrule
1777	<	PC  &     & 1.984 & 0.000 & 0.012 & 0.000 & 1.357 & 0.000 \\
1778	<	SP  & 0.0 & 0.850 & 0.552 & 0.907 & 0.703 & 0.938 & 0.793 \\
1779	<	& 0.1 & 0.924 & 0.755 & 0.980 & 0.936 & 0.995 & 0.988 \\
1780	<	& 0.2 & 0.985 & 0.983 & 0.986 & 0.988 & 0.987 & 0.988 \\
1781	<	& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\
1782	<	SF  & 0.0 & 0.977 & 0.989 & 0.987 & 0.995 & 0.992 & 0.998 \\
1783	<	& 0.1 & 0.982 & 0.989 & 0.992 & 0.996 & 0.997 & 0.998 \\
1784	<	& 0.2 & 0.984 & 0.987 & 0.986 & 0.987 & 0.987 & 0.988 \\
1785	<	& 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\
1786	<	GSC &     & 0.995 & 0.981 & 0.999 & 0.990 & 1.000 & 0.993 \\
1787	<	RF  &     & 0.993 & 0.988 & 0.997 & 0.995 & 0.999 & 0.998 \\
1788	<	\bottomrule
1789	<	\end{tabular}
1790	<	\label{tab:argon}
1791	<	\end{table}
46	>	\input{electrostaticsChapter}
47
1793	–	\begin{table}[htbp]
1794	–	\centering
1795	–	\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1796	–	DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6\AA\ SPHERE OF
1797	–	ARGON IN LIQUID WATER SYSTEM}
1798	–
1799	–	\footnotesize
1800	–	\begin{tabular}{@{} ccrrrrrr @{}}
1801	–	\\
1802	–	\toprule
1803	–	\toprule
1804	–	& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\
1805	–	\cmidrule(lr){3-5}
1806	–	\cmidrule(l){6-8}
1807	–	Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\
1808	–	\midrule
1809	–	PC  &     & 568.025 & 265.993 & 195.099 & 246.626 & 138.600 & 91.654 \\
1810	–	SP  & 0.0 & 504.578 & 251.694 & 179.932 & 231.568 & 131.444 & 85.119 \\
1811	–	& 0.1 & 224.886 & 49.746 & 9.346 & 104.482 & 23.683 & 4.480 \\
1812	–	& 0.2 & 4.889 & 0.197 & 0.155 & 6.029 & 2.507 & 2.269 \\
1813	–	& 0.3 & 0.817 & 0.833 & 0.812 & 8.286 & 8.436 & 8.135 \\
1814	–	SF  & 0.0 & 1.924 & 0.675 & 0.304 & 3.658 & 1.448 & 0.600 \\
1815	–	& 0.1 & 1.937 & 0.515 & 0.143 & 3.565 & 1.308 & 0.546 \\
1816	–	& 0.2 & 0.407 & 0.166 & 0.156 & 3.086 & 2.501 & 2.274 \\
1817	–	& 0.3 & 0.815 & 0.833 & 0.812 & 8.330 & 8.437 & 8.135 \\
1818	–	GSC &     & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\
1819	–	RF  &     & 1.822 & 0.408 & 0.142 & 3.799 & 1.362 & 0.550 \\
1820	–	\midrule
1821	–	GSSP  & 0.0 & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\
1822	–	& 0.1 & 1.652 & 0.309 & 0.087 & 4.197 & 1.401 & 0.590 \\
1823	–	& 0.2 & 0.465 & 0.165 & 0.153 & 3.323 & 2.529 & 2.273 \\
1824	–	& 0.3 & 0.813 & 0.825 & 0.816 & 8.316 & 8.447 & 8.132 \\
1825	–	GSSF  & 0.0 & 1.173 & 0.292 & 0.113 & 3.452 & 1.347 & 0.583 \\
1826	–	& 0.1 & 1.166 & 0.240 & 0.076 & 3.381 & 1.281 & 0.575 \\
1827	–	& 0.2 & 0.459 & 0.165 & 0.153 & 3.430 & 2.542 & 2.273 \\
1828	–	& 0.3 & 0.814 & 0.825 & 0.816 & 8.325 & 8.447 & 8.132 \\
1829	–	\bottomrule
1830	–	\end{tabular}
1831	–	\label{tab:argonAng}
1832	–	\end{table}
1833	–
1834	–	This system does not appear to show any significant deviations from
1835	–	the previously observed results. The {\sc sp} and {\sc sf} methods
1836	–	have aggrements similar to those observed in section
1837	–	\ref{sec:WaterResults}. The only significant difference is the
1838	–	improvement in the configuration energy differences for the {\sc rf}
1839	–	method. This is surprising in that we are introducing an inhomogeneity
1840	–	to the system; however, this inhomogeneity is charge-neutral and does
1841	–	not result in charged cutoff spheres. The charge-neutrality of the
1842	–	cutoff spheres, which the {\sc sp} and {\sc sf} methods explicitly
1843	–	enforce, seems to play a greater role in the stability of the {\sc rf}
1844	–	method than the required homogeneity of the environment.
1845	–
1846	–
1847	–	\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics}
1848	–
1849	–	Zahn {\it et al.} investigated the structure and dynamics of water
1850	–	using eqs. (\ref{eq:ZahnPot}) and
1851	–	(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated
1852	–	that a method similar (but not identical with) the damped {\sc sf}
1853	–	method resulted in properties very similar to those obtained when
1854	–	using the Ewald summation.  The properties they studied (pair
1855	–	distribution functions, diffusion constants, and velocity and
1856	–	orientational correlation functions) may not be particularly sensitive
1857	–	to the long-range and collective behavior that governs the
1858	–	low-frequency behavior in crystalline systems.  Additionally, the
1859	–	ionic crystals are the worst case scenario for the pairwise methods
1860	–	because they lack the reciprocal space contribution contained in the
1861	–	Ewald summation.
1862	–
1863	–	We are using two separate measures to probe the effects of these
1864	–	alternative electrostatic methods on the dynamics in crystalline
1865	–	materials.  For short- and intermediate-time dynamics, we are
1866	–	computing the velocity autocorrelation function, and for long-time
1867	–	and large length-scale collective motions, we are looking at the
1868	–	low-frequency portion of the power spectrum.
1869	–
1870	–	\begin{figure}
1871	–	\centering
1872	–	\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf}
1873	–	\caption{Velocity autocorrelation functions of NaCl crystals at
1874	–	1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \&
1875	–	0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2\AA$^{-1}$). The inset is
1876	–	a magnification of the area around the first minimum.  The times to
1877	–	first collision are nearly identical, but differences can be seen in
1878	–	the peaks and troughs, where the undamped and weakly damped methods
1879	–	are stiffer than the moderately damped and {\sc spme} methods.}
1880	–	\label{fig:vCorrPlot}
1881	–	\end{figure}
1882	–
1883	–	The short-time decay of the velocity autocorrelation function through
1884	–	the first collision are nearly identical in figure
1885	–	\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
1886	–	how the methods differ.  The undamped {\sc sf} method has deeper
1887	–	troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than
1888	–	any of the other methods.  As the damping parameter ($\alpha$) is
1889	–	increased, these peaks are smoothed out, and the {\sc sf} method
1890	–	approaches the {\sc spme} results.  With $\alpha$ values of 0.2\AA$^{-1}$,
1891	–	the {\sc sf} and {\sc sp} functions are nearly identical and track the
1892	–	{\sc spme} features quite well.  This is not surprising because the {\sc sf}
1893	–	and {\sc sp} potentials become nearly identical with increased
1894	–	damping.  However, this appears to indicate that once damping is
1895	–	utilized, the details of the form of the potential (and forces)
1896	–	constructed out of the damped electrostatic interaction are less
1897	–	important.
1898	–
1899	–	\section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}
1900	–
1901	–	To evaluate how the differences between the methods affect the
1902	–	collective long-time motion, we computed power spectra from long-time
1903	–	traces of the velocity autocorrelation function. The power spectra for
1904	–	the best-performing alternative methods are shown in
1905	–	fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1906	–	a cubic switching function between 40 and 50 ps was used to reduce the
1907	–	ringing resulting from data truncation.  This procedure had no
1908	–	noticeable effect on peak location or magnitude.
1909	–
1910	–	\begin{figure}
1911	–	\centering
1912	–	\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
1913	–	\caption{Power spectra obtained from the velocity auto-correlation
1914	–	functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf}
1915	–	($\alpha$ = 0, 0.1, \& 0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ =
1916	–	0.2\AA$^{-1}$).  The inset shows the frequency region below 100
1917	–	cm$^{-1}$ to highlight where the spectra differ.}
1918	–	\label{fig:methodPS}
1919	–	\end{figure}
1920	–
1921	–	While the high frequency regions of the power spectra for the
1922	–	alternative methods are quantitatively identical with Ewald spectrum,
1923	–	the low frequency region shows how the summation methods differ.
1924	–	Considering the low-frequency inset (expanded in the upper frame of
1925	–	figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the
1926	–	correlated motions are blue-shifted when using undamped or weakly
1927	–	damped {\sc sf}.  When using moderate damping ($\alpha = 0.2$
1928	–	\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical
1929	–	correlated motion to the Ewald method (which has a convergence
1930	–	parameter of 0.3119\AA$^{-1}$).  This weakening of the electrostatic
1931	–	interaction with increased damping explains why the long-ranged
1932	–	correlated motions are at lower frequencies for the moderately damped
1933	–	methods than for undamped or weakly damped methods.
1934	–
1935	–	To isolate the role of the damping constant, we have computed the
1936	–	spectra for a single method ({\sc sf}) with a range of damping
1937	–	constants and compared this with the {\sc spme} spectrum.
1938	–	Fig. \ref{fig:dampInc} shows more clearly that increasing the
1939	–	electrostatic damping red-shifts the lowest frequency phonon modes.
1940	–	However, even without any electrostatic damping, the {\sc sf} method
1941	–	has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
1942	–	Without the {\sc sf} modifications, an undamped (pure cutoff) method
1943	–	would predict the lowest frequency peak near 325 cm$^{-1}$.  {\it
1944	–	Most} of the collective behavior in the crystal is accurately captured
1945	–	using the {\sc sf} method.  Quantitative agreement with Ewald can be
1946	–	obtained using moderate damping in addition to the shifting at the
1947	–	cutoff distance.
1948	–
1949	–	\begin{figure}
1950	–	\centering
1951	–	\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
1952	–	\caption{Effect of damping on the two lowest-frequency phonon modes in
1953	–	the NaCl crystal at 1000K.  The undamped shifted force ({\sc sf})
1954	–	method is off by less than 10 cm$^{-1}$, and increasing the
1955	–	electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement
1956	–	with the power spectrum obtained using the Ewald sum.  Overdamping can
1957	–	result in underestimates of frequencies of the long-wavelength
1958	–	motions.}
1959	–	\label{fig:dampInc}
1960	–	\end{figure}
1961	–
1962	–	\section{Synopsis of the Pairwise Method Evaluation}\label{sec:PairwiseSynopsis}
1963	–
1964	–	The above investigation of pairwise electrostatic summation techniques
1965	–	shows that there are viable and computationally efficient alternatives
1966	–	to the Ewald summation.  These methods are derived from the damped and
1967	–	cutoff-neutralized Coulombic sum originally proposed by Wolf
1968	–	\textit{et al.}\cite{Wolf99} In particular, the {\sc sf}
1969	–	method, reformulated above as eqs. (\ref{eq:DSFPot}) and
1970	–	(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the
1971	–	energetic and dynamic characteristics exhibited by simulations
1972	–	employing lattice summation techniques.  The cumulative energy
1973	–	difference results showed the undamped {\sc sf} and moderately damped
1974	–	{\sc sp} methods produced results nearly identical to {\sc spme}.
1975	–	Similarly for the dynamic features, the undamped or moderately damped
1976	–	{\sc sf} and moderately damped {\sc sp} methods produce force and
1977	–	torque vector magnitude and directions very similar to the expected
1978	–	values.  These results translate into long-time dynamic behavior
1979	–	equivalent to that produced in simulations using {\sc spme}.
1980	–
1981	–	As in all purely-pairwise cutoff methods, these methods are expected
1982	–	to scale approximately {\it linearly} with system size, and they are
1983	–	easily parallelizable.  This should result in substantial reductions
1984	–	in the computational cost of performing large simulations.
1985	–
1986	–	Aside from the computational cost benefit, these techniques have
1987	–	applicability in situations where the use of the Ewald sum can prove
1988	–	problematic.  Of greatest interest is their potential use in
1989	–	interfacial systems, where the unmodified lattice sum techniques
1990	–	artificially accentuate the periodicity of the system in an
1991	–	undesirable manner.  There have been alterations to the standard Ewald
1992	–	techniques, via corrections and reformulations, to compensate for
1993	–	these systems; but the pairwise techniques discussed here require no
1994	–	modifications, making them natural tools to tackle these problems.
1995	–	Additionally, this transferability gives them benefits over other
1996	–	pairwise methods, like reaction field, because estimations of physical
1997	–	properties (e.g. the dielectric constant) are unnecessary.
1998	–
1999	–	If a researcher is using Monte Carlo simulations of large chemical
2000	–	systems containing point charges, most structural features will be
2001	–	accurately captured using the undamped {\sc sf} method or the {\sc sp}
2002	–	method with an electrostatic damping of 0.2\AA$^{-1}$.  These methods
2003	–	would also be appropriate for molecular dynamics simulations where the
2004	–	data of interest is either structural or short-time dynamical
2005	–	quantities.  For long-time dynamics and collective motions, the safest
2006	–	pairwise method we have evaluated is the {\sc sf} method with an
2007	–	electrostatic damping between 0.2 and 0.25\AA$^{-1}$.
2008	–
2009	–	We are not suggesting that there is any flaw with the Ewald sum; in
2010	–	fact, it is the standard by which these simple pairwise sums have been
2011	–	judged.  However, these results do suggest that in the typical
2012	–	simulations performed today, the Ewald summation may no longer be
2013	–	required to obtain the level of accuracy most researchers have come to
2014	–	expect.
2015	–
2016	–	\section{An Application: TIP5P-E Water}
2017	–
2018	–
48		\chapter{\label{chap:water}SIMPLE MODELS FOR WATER}
49
50		\chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}

Diff Legend

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