ViewVC Help

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

root/group/trunk/multipole/multipole_2/multipole2.tex

Comparing trunk/multipole/multipole_2/multipole2.tex (file contents):
Revision 4170 by gezelter, Wed Jun 4 18:48:27 2014 UTC vs.
Revision 4207 by gezelter, Thu Aug 7 21:13:00 2014 UTC

#	Line 35 | Line 35 | preprint,
35		%\linenumbers\relax % Commence numbering lines
36		\usepackage{amsmath}
37		\usepackage{times}
38	<	\usepackage{mathptm}
38	>	\usepackage{mathptmx}
39	>	\usepackage{tabularx}
40		\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
41		\usepackage{url}
42		\usepackage[english]{babel}
43
44	+	\newcolumntype{Y}{>{\centering\arraybackslash}X}
45
46		\begin{document}
47
48	<	\preprint{AIP/123-QED}
48	>	%\preprint{AIP/123-QED}
49
50	<	\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald
51	<	Sum. II. Comparison of Simulation Methodologies} % Force line breaks with \\
50	>	\title{Real space electrostatics for multipoles. II. Comparisons with
51	>	the Ewald Sum}
52
53		\author{Madan Lamichhane}
54	<	\affiliation{Department of Physics, University
53	<	of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\
54	>	\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556}
55
56		\author{Kathie E. Newman}
57	<	\affiliation{Department of Physics, University
57	<	of Notre Dame, Notre Dame, IN 46556}
57	>	\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556}
58
59		\author{J. Daniel Gezelter}%
60		\email{gezelter@nd.edu.}
61	<	\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash
62	<	}%
61	>	\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556
62	>	}
63
64	<	\date{\today}% It is always \today, today,
65	<	%  but any date may be explicitly specified
64	>	\date{\today}
65
66		\begin{abstract}
67	<	We have tested our recently developed shifted potential, gradient-shifted force, and Taylor-shifted force methods for the higher-order multipoles against Ewald’s method in different types of liquid and crystalline system. In this paper, we have also investigated the conservation of total energy in the molecular dynamic simulation using all of these methods. The shifted potential method shows better agreement with the Ewald in the energy differences between different configurations as compared to the direct truncation. Both the gradient shifted force and Taylor-shifted force methods reproduce very good energy conservation. But the absolute energy, force and torque evaluated from the gradient shifted force method shows better result as compared to taylor-shifted force method. Hence the gradient-shifted force method suitably mimics the electrostatic interaction in the molecular dynamic simulation.
67	>	We report on tests of the shifted potential (SP), gradient shifted
68	>	force (GSF), and Taylor shifted force (TSF) real-space methods for
69	>	multipole interactions developed in the first paper in this series,
70	>	using the multipolar Ewald sum as a reference method. The tests were
71	>	carried out in a variety of condensed-phase environments designed to
72	>	test up to quadrupole-quadrupole interactions.  Comparisons of the
73	>	energy differences between configurations, molecular forces, and
74	>	torques were used to analyze how well the real-space models perform
75	>	relative to the more computationally expensive Ewald treatment.  We
76	>	have also investigated the energy conservation properties of the new
77	>	methods in molecular dynamics simulations. The SP method shows
78	>	excellent agreement with configurational energy differences, forces,
79	>	and torques, and would be suitable for use in Monte Carlo
80	>	calculations.  Of the two new shifted-force methods, the GSF
81	>	approach shows the best agreement with Ewald-derived energies,
82	>	forces, and torques and also exhibits energy conservation properties
83	>	that make it an excellent choice for efficient computation of
84	>	electrostatic interactions in molecular dynamics simulations.
85		\end{abstract}
86
87	<	\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
87	>	%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
88		% Classification Scheme.
89	<	\keywords{Electrostatics, Multipoles, Real-space}
89	>	%\keywords{Electrostatics, Multipoles, Real-space}
90
91		\maketitle
92
77	–
93		\section{\label{sec:intro}Introduction}
94		Computing the interactions between electrostatic sites is one of the
95	<	most expensive aspects of molecular simulations, which is why there
96	<	have been significant efforts to develop practical, efficient and
97	<	convergent methods for handling these interactions. Ewald's method is
98	<	perhaps the best known and most accurate method for evaluating
99	<	energies, forces, and torques in explicitly-periodic simulation
100	<	cells. In this approach, the conditionally convergent electrostatic
101	<	energy is converted into two absolutely convergent contributions, one
102	<	which is carried out in real space with a cutoff radius, and one in
103	<	reciprocal space.\cite{Clarke:1986eu,Woodcock75}
95	>	most expensive aspects of molecular simulations. There have been
96	>	significant efforts to develop practical, efficient and convergent
97	>	methods for handling these interactions. Ewald's method is perhaps the
98	>	best known and most accurate method for evaluating energies, forces,
99	>	and torques in explicitly-periodic simulation cells. In this approach,
100	>	the conditionally convergent electrostatic energy is converted into
101	>	two absolutely convergent contributions, one which is carried out in
102	>	real space with a cutoff radius, and one in reciprocal
103	>	space.\cite{Ewald21,deLeeuw80,Smith81,Allen87}
104
105		When carried out as originally formulated, the reciprocal-space
106		portion of the Ewald sum exhibits relatively poor computational
107	<	scaling, making it prohibitive for large systems. By utilizing
108	<	particle meshes and three dimensional fast Fourier transforms (FFT),
109	<	the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald
110	<	(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease
111	<	the computational cost from $O(N^2)$ down to $O(N \log
112	<	N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}.
107	>	scaling, making it prohibitive for large systems. By utilizing a
108	>	particle mesh and three dimensional fast Fourier transforms (FFT), the
109	>	particle-mesh Ewald (PME), particle-particle particle-mesh Ewald
110	>	(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME)
111	>	methods can decrease the computational cost from $O(N^2)$ down to $O(N
112	>	\log
113	>	N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}
114
115	<	Because of the artificial periodicity required for the Ewald sum, the
100	<	method may require modification to compute interactions for
115	>	Because of the artificial periodicity required for the Ewald sum,
116		interfacial molecular systems such as membranes and liquid-vapor
117	<	interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
118	<	To simulate interfacial systems, Parry’s extension of the 3D Ewald sum
119	<	is appropriate for slab geometries.\cite{Parry:1975if} The inherent
120	<	periodicity in the Ewald’s method can also be problematic for
121	<	interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald
122	<	methods that were developed to handle two-dimensional (2D)
123	<	electrostatic interactions in interfacial systems have not had similar
124	<	particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77,
125	<	Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq}
117	>	interfaces require modifications to the method.  Parry's extension of
118	>	the three dimensional Ewald sum is appropriate for slab
119	>	geometries.\cite{Parry:1975if} Modified Ewald methods that were
120	>	developed to handle two-dimensional (2-D) electrostatic
121	>	interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
122	>	These methods were originally quite computationally
123	>	expensive.\cite{Spohr97,Yeh99} There have been several successful
124	>	efforts that reduced the computational cost of 2-D lattice summations,
125	>	bringing them more in line with the scaling for the full 3-D
126	>	treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The
127	>	inherent periodicity required by the Ewald method can also be
128	>	problematic in a number of protein/solvent and ionic solution
129	>	environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq}
130
131		\subsection{Real-space methods}
132		Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$
133		method for calculating electrostatic interactions between point
134	<	charges. They argued that the effective Coulomb interaction in
135	<	condensed systems is actually short ranged.\cite{Wolf92,Wolf95}.  For
136	<	an ordered lattice (e.g. when computing the Madelung constant of an
137	<	ionic solid), the material can be considered as a set of ions
138	<	interacting with neutral dipolar or quadrupolar ``molecules'' giving
139	<	an effective distance dependence for the electrostatic interactions of
140	<	$r^{-5}$ (see figure \ref{fig:NaCl}.  For this reason, careful
141	<	applications of Wolf's method are able to obtain accurate estimates of
142	<	Madelung constants using relatively short cutoff radii.  Recently,
143	<	Fukuda used neutralization of the higher order moments for the
144	<	calculation of the electrostatic interaction of the point charges
145	<	system.\cite{Fukuda:2013sf}
134	>	charges. They argued that the effective Coulomb interaction in most
135	>	condensed phase systems is effectively short
136	>	ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when
137	>	computing the Madelung constant of an ionic solid), the material can
138	>	be considered as a set of ions interacting with neutral dipolar or
139	>	quadrupolar ``molecules'' giving an effective distance dependence for
140	>	the electrostatic interactions of $r^{-5}$ (see figure
141	>	\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple
142	>	cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the
143	>	electrostatic energy per ion converges more rapidly to the Madelung
144	>	energy than the dipolar approximation.\cite{Wolf92} To find the
145	>	correct Madelung constant, Lacman suggested that the NaCl structure
146	>	could be constructed in a way that the finite crystal terminates with
147	>	complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees
148	>	what is effectively a set of octupoles at large distances. These facts
149	>	suggest that the Madelung constants are relatively short ranged for
150	>	perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful
151	>	application of Wolf's method can provide accurate estimates of
152	>	Madelung constants using relatively short cutoff radii.
153
154	<	\begin{figure}[h!]
154	>	Direct truncation of interactions at a cutoff radius creates numerical
155	>	errors.  Wolf \textit{et al.} suggest that truncation errors are due
156	>	to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To
157	>	neutralize this charge they proposed placing an image charge on the
158	>	surface of the cutoff sphere for every real charge inside the cutoff.
159	>	These charges are present for the evaluation of both the pair
160	>	interaction energy and the force, although the force expression
161	>	maintains a discontinuity at the cutoff sphere.  In the original Wolf
162	>	formulation, the total energy for the charge and image were not equal
163	>	to the integral of the force expression, and as a result, the total
164	>	energy would not be conserved in molecular dynamics (MD)
165	>	simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and
166	>	Gezelter later proposed shifted force variants of the Wolf method with
167	>	commensurate force and energy expressions that do not exhibit this
168	>	problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods
169	>	were also proposed by Chen \textit{et
170	>	al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw}
171	>	and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly
172	>	used additional neutralization of higher order moments for systems of
173	>	point charges.\cite{Fukuda:2013sf}
174	>
175	>	\begin{figure}
176		\centering
177	<	\includegraphics[width=0.50 \textwidth]{chargesystem.pdf}
178	<	\caption{Top: NaCl crystal showing how spherical truncation can
179	<	breaking effective charge ordering, and how complete \ce{(NaCl)4}
180	<	molecules interact with the central ion.  Bottom: A dipolar
181	<	crystal exhibiting similar behavior and illustrating how the
182	<	effective dipole-octupole interactions can be disrupted by
183	<	spherical truncation.}
184	<	\label{fig:NaCl}
177	>	\includegraphics[width=\linewidth]{schematic.eps}
178	>	\caption{Top: Ionic systems exhibit local clustering of dissimilar
179	>	charges (in the smaller grey circle), so interactions are
180	>	effectively charge-multipole at longer distances.  With hard
181	>	cutoffs, motion of individual charges in and out of the cutoff
182	>	sphere can break the effective multipolar ordering.  Bottom:
183	>	dipolar crystals and fluids have a similar effective
184	>	\textit{quadrupolar} ordering (in the smaller grey circles), and
185	>	orientational averaging helps to reduce the effective range of the
186	>	interactions in the fluid.  Placement of reversed image multipoles
187	>	on the surface of the cutoff sphere recovers the effective
188	>	higher-order multipole behavior.}
189	>	\label{fig:schematic}
190		\end{figure}
191
192	<	The direct truncation of interactions at a cutoff radius creates
193	<	truncation defects. Wolf \textit{et al.} further argued that
194	<	truncation errors are due to net charge remaining inside the cutoff
195	<	sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed
196	<	placing an image charge on the surface of the cutoff sphere for every
197	<	real charge inside the cutoff.  These charges are present for the
198	<	evaluation of both the pair interaction energy and the force, although
199	<	the force expression maintained a discontinuity at the cutoff sphere.
200	<	In the original Wolf formulation, the total energy for the charge and
201	<	image were not equal to the integral of their force expression, and as
150	<	a result, the total energy would not be conserved in molecular
151	<	dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and
152	<	Fennel and Gezelter later proposed shifted force variants of the Wolf
153	<	method with commensurate force and energy expressions that do not
154	<	exhibit this problem.\cite{Fennell:2006lq}   Related real-space
155	<	methods were also proposed by Chen \textit{et
156	<	al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw}
157	<	and by Wu and Brooks.\cite{Wu:044107}
158	<
159	<	Considering the interaction of one central ion in an ionic crystal
160	<	with a portion of the crystal at some distance, the effective Columbic
161	<	potential is found to be decreasing as $r^{-5}$. If one views the
162	<	\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar
163	<	\ce{(NaCl)4} basis, the electrostatic energy per ion converges more
164	<	rapidly to the Madelung energy than the dipolar
165	<	approximation.\cite{Wolf92} To find the correct Madelung constant,
166	<	Lacman suggested that the NaCl structure could be constructed in a way
167	<	that the finite crystal terminates with complete \ce{(NaCl)4}
168	<	molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded
169	<	by opposite charges. Similarly for each pair of charges, there is an
170	<	opposite pair of charge adjacent to it.  The central ion sees what is
171	<	effectively a set of octupoles at large distances. These facts suggest
172	<	that the Madelung constants are relatively short ranged for perfect
173	<	ionic crystals.\cite{Wolf:1999dn}
174	<
175	<	One can make a similar argument for crystals of point multipoles. The
176	<	Luttinger and Tisza treatment of energy constants for dipolar lattices
177	<	utilizes 24 basis vectors that contain dipoles at the eight corners of
178	<	a unit cube.  Only three of these basis vectors, $X_1, Y_1,
179	<	\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have
180	<	zero net dipole and retain contributions only from higher order
181	<	multipoles.  The effective interaction between a dipole at the center
192	>	One can make a similar effective range argument for crystals of point
193	>	\textit{multipoles}. The Luttinger and Tisza treatment of energy
194	>	constants for dipolar lattices utilizes 24 basis vectors that contain
195	>	dipoles at the eight corners of a unit cube.\cite{LT} Only three of
196	>	these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole
197	>	moments, while the rest have zero net dipole and retain contributions
198	>	only from higher order multipoles.  The lowest-energy crystalline
199	>	structures are built out of basis vectors that have only residual
200	>	quadrupolar moments (e.g. the $Z_5$ array). In these low energy
201	>	structures, the effective interaction between a dipole at the center
202		of a crystal and a group of eight dipoles farther away is
203		significantly shorter ranged than the $r^{-3}$ that one would expect
204		for raw dipole-dipole interactions.  Only in crystals which retain a
#	Line 188 | Line 208 | multipolar arrangements (see Fig. \ref{fig:NaCl}), cau
208		unstable.
209
210		In ionic crystals, real-space truncation can break the effective
211	<	multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant
212	<	swings in the electrostatic energy as the cutoff radius is increased
213	<	(or as individual ions move back and forth across the boundary).  This
214	<	is why the image charges were necessary for the Wolf sum to exhibit
215	<	rapid convergence.  Similarly, the real-space truncation of point
216	<	multipole interactions breaks higher order multipole arrangements, and
217	<	image multipoles are required for real-space treatments of
198	<	electrostatic energies.
211	>	multipolar arrangements (see Fig. \ref{fig:schematic}), causing
212	>	significant swings in the electrostatic energy as individual ions move
213	>	back and forth across the boundary.  This is why the image charges are
214	>	necessary for the Wolf sum to exhibit rapid convergence.  Similarly,
215	>	the real-space truncation of point multipole interactions breaks
216	>	higher order multipole arrangements, and image multipoles are required
217	>	for real-space treatments of electrostatic energies.
218
219	+	The shorter effective range of electrostatic interactions is not
220	+	limited to perfect crystals, but can also apply in disordered fluids.
221	+	Even at elevated temperatures, there is local charge balance in an
222	+	ionic liquid, where each positive ion has surroundings dominated by
223	+	negative ions and vice versa.  The reversed-charge images on the
224	+	cutoff sphere that are integral to the Wolf and DSF approaches retain
225	+	the effective multipolar interactions as the charges traverse the
226	+	cutoff boundary.
227	+
228	+	In multipolar fluids (see Fig. \ref{fig:schematic}) there is
229	+	significant orientational averaging that additionally reduces the
230	+	effect of long-range multipolar interactions.  The image multipoles
231	+	that are introduced in the TSF, GSF, and SP methods mimic this effect
232	+	and reduce the effective range of the multipolar interactions as
233	+	interacting molecules traverse each other's cutoff boundaries.
234	+
235		% Because of this reason, although the nature of electrostatic
236		% interaction short ranged, the hard cutoff sphere creates very large
237		% fluctuation in the electrostatic energy for the perfect crystal. In
#	Line 207 | Line 242 | The forces and torques acting on atomic sites are the
242		% to the non-neutralized value of the higher order moments within the
243		% cutoff sphere.
244
245	<	The forces and torques acting on atomic sites are the fundamental
246	<	factors driving dynamics in molecular simulations. Fennell and
247	<	Gezelter proposed the damped shifted force (DSF) energy kernel to
248	<	obtain consistent energies and forces on the atoms within the cutoff
249	<	sphere. Both the energy and the force go smoothly to zero as an atom
250	<	aproaches the cutoff radius. The comparisons of the accuracy these
251	<	quantities between the DSF kernel and SPME was surprisingly
252	<	good.\cite{Fennell:2006lq} The DSF method has seen increasing use for
253	<	calculating electrostatic interactions in molecular systems with
254	<	relatively uniform charge
220	<	densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13}
245	>	Forces and torques acting on atomic sites are fundamental in driving
246	>	dynamics in molecular simulations, and the damped shifted force (DSF)
247	>	energy kernel provides consistent energies and forces on charged atoms
248	>	within the cutoff sphere. Both the energy and the force go smoothly to
249	>	zero as an atom aproaches the cutoff radius. The comparisons of the
250	>	accuracy these quantities between the DSF kernel and SPME was
251	>	surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method
252	>	has seen increasing use in molecular systems with relatively uniform
253	>	charge
254	>	densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13}
255
256		\subsection{The damping function}
257	<	The damping function used in our research has been discussed in detail
258	<	in the first paper of this series.\cite{PaperI} The radial kernel
259	<	$1/r$ for the interactions between point charges can be replaced by
260	<	the complementary error function $\mathrm{erfc}(\alpha r)/r$ to
261	<	accelerate the rate of convergence, where $\alpha$ is a damping
262	<	parameter with units of inverse distance.  Altering the value of
263	<	$\alpha$ is equivalent to changing the width of Gaussian charge
264	<	distributions that replace each point charge -- Gaussian overlap
265	<	integrals yield complementary error functions when truncated at a
266	<	finite distance.
257	>	The damping function has been discussed in detail in the first paper
258	>	of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the
259	>	interactions between point charges can be replaced by the
260	>	complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate
261	>	convergence, where $\alpha$ is a damping parameter with units of
262	>	inverse distance.  Altering the value of $\alpha$ is equivalent to
263	>	changing the width of Gaussian charge distributions that replace each
264	>	point charge, as Coulomb integrals with Gaussian charge distributions
265	>	produce complementary error functions when truncated at a finite
266	>	distance.
267
268	<	By using suitable value of damping alpha ($\alpha \sim 0.2$) for a
269	<	cutoff radius ($r_{­c}=9 A$), Fennel and Gezelter produced very good
270	<	agreement with SPME for the interaction energies, forces and torques
271	<	for charge-charge interactions.\cite{Fennell:2006lq}
268	>	With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method
269	>	produced very good agreement with SPME for interaction energies,
270	>	forces and torques for charge-charge
271	>	interactions.\cite{Fennell:2006lq}
272
273		\subsection{Point multipoles in molecular modeling}
274		Coarse-graining approaches which treat entire molecular subsystems as
275		a single rigid body are now widely used. A common feature of many
276		coarse-graining approaches is simplification of the electrostatic
277		interactions between bodies so that fewer site-site interactions are
278	<	required to compute configurational energies.  Many coarse-grained
279	<	molecular structures would normally consist of equal positive and
246	<	negative charges, and rather than use multiple site-site interactions,
247	<	the interaction between higher order multipoles can also be used to
248	<	evaluate a single molecule-molecule
249	<	interaction.\cite{Ren06,Essex10,Essex11}
278	>	required to compute configurational
279	>	energies.\cite{Ren06,Essex10,Essex11}
280
281	<	Because electrons in a molecule are not localized at specific points,
282	<	the assignment of partial charges to atomic centers is a relatively
283	<	rough approximation.  Atomic sites can also be assigned point
284	<	multipoles and polarizabilities to increase the accuracy of the
285	<	molecular model.  Recently, water has been modeled with point
286	<	multipoles up to octupolar
287	<	order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point
281	>	Additionally, because electrons in a molecule are not localized at
282	>	specific points, the assignment of partial charges to atomic centers
283	>	is always an approximation.  For increased accuracy, atomic sites can
284	>	also be assigned point multipoles and polarizabilities.  Recently,
285	>	water has been modeled with point multipoles up to octupolar order
286	>	using the soft sticky dipole-quadrupole-octupole (SSDQO)
287	>	model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point
288		multipoles up to quadrupolar order have also been coupled with point
289		polarizabilities in the high-quality AMOEBA and iAMOEBA water
290	<	models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}.  But
291	<	using point multipole with the real space truncation without
292	<	accounting for multipolar neutrality will create energy conservation
293	<	issues in molecular dynamics (MD) simulations.
290	>	models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However,
291	>	truncating point multipoles without smoothing the forces and torques
292	>	can create energy conservation issues in molecular dynamics
293	>	simulations.
294
295		In this paper we test a set of real-space methods that were developed
296		for point multipolar interactions.  These methods extend the damped
297		shifted force (DSF) and Wolf methods originally developed for
298		charge-charge interactions and generalize them for higher order
299	<	multipoles. The detailed mathematical development of these methods has
300	<	been presented in the first paper in this series, while this work
301	<	covers the testing the energies, forces, torques, and energy
299	>	multipoles.  The detailed mathematical development of these methods
300	>	has been presented in the first paper in this series, while this work
301	>	covers the testing of energies, forces, torques, and energy
302		conservation properties of the methods in realistic simulation
303		environments.  In all cases, the methods are compared with the
304	<	reference method, a full multipolar Ewald treatment.
304	>	reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98}
305
306
307		%\subsection{Conservation of total energy }
#	Line 295 | Line 325 | U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1
325		expressed as the product of two multipole operators and a Coulombic
326		kernel,
327		\begin{equation}
328	<	U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
328	>	U_{ab}(r)= M_{a} M_{b} \frac{1}{r}  \label{kernel}.
329		\end{equation}
330	<	where the multipole operator for site $\bf a$,
331	<	\begin{equation}
332	<	\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}}
303	<	+  Q_{{\bf a}\alpha\beta}
304	<	\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
305	<	\end{equation}
306	<	is expressed in terms of the point charge, $C_{\bf a}$, dipole,
307	<	$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for
308	<	object $\bf a$.  Note that in this work, we use the primitive
309	<	quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k
310	<	r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site.
330	>	where the multipole operator for site $a$, $M_{a}$, is
331	>	expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object
332	>	$a$, etc.
333
334	<	Interactions between multipoles can be expressed as higher derivatives
335	<	of the bare Coulomb potential, so one way of ensuring that the forces
336	<	and torques vanish at the cutoff distance is to include a larger
337	<	number of terms in the truncated Taylor expansion, e.g.,
338	<	%
339	<	\begin{equation}
340	<	f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-R_c)^m}{m!} f^{(m)} \Big \lvert  _{R_c}  .
341	<	\end{equation}
342	<	%
343	<	The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
344	<	Thus, for $f(r)=1/r$, we find
345	<	%
346	<	\begin{equation}
347	<	f_1(r)=\frac{1}{r}- \frac{1}{R_c} + (r - R_c) \frac{1}{R_c^2} - \frac{(r-R_c)^2}{R_c^3} .
348	<	\end{equation}
349	<	This function is an approximate electrostatic potential that has
350	<	vanishing second derivatives at the cutoff radius, making it suitable
351	<	for shifting the forces and torques of charge-dipole interactions.
334	>	% Interactions between multipoles can be expressed as higher derivatives
335	>	% of the bare Coulomb potential, so one way of ensuring that the forces
336	>	% and torques vanish at the cutoff distance is to include a larger
337	>	% number of terms in the truncated Taylor expansion, e.g.,
338	>	% %
339	>	% \begin{equation}
340	>	% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert  _{r_c}  .
341	>	% \end{equation}
342	>	% %
343	>	% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
344	>	% Thus, for $f(r)=1/r$, we find
345	>	% %
346	>	% \begin{equation}
347	>	% f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} .
348	>	% \end{equation}
349	>	% This function is an approximate electrostatic potential that has
350	>	% vanishing second derivatives at the cutoff radius, making it suitable
351	>	% for shifting the forces and torques of charge-dipole interactions.
352
353	<	In general, the TSF potential for any multipole-multipole interaction
354	<	can be written
353	>	The TSF potential for any multipole-multipole interaction can be
354	>	written
355		\begin{equation}
356		U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r)
357		\label{generic}
358		\end{equation}
359	<	with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for
360	<	charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and
361	<	$n=4$ for quadrupole-quadrupole.  To ensure smooth convergence of the
362	<	energy, force, and torques, the required number of terms from Taylor
363	<	series expansion in $f_n(r)$ must be performed for different
364	<	multipole-multipole interactions.
359	>	where $f_n(r)$ is a shifted kernel that is appropriate for the order
360	>	of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for
361	>	charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
362	>	and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for
363	>	quadrupole-quadrupole.  To ensure smooth convergence of the energy,
364	>	force, and torques, a Taylor expansion with $n$ terms must be
365	>	performed at cutoff radius ($r_c$) to obtain $f_n(r)$.
366
367	<	To carry out the same procedure for a damped electrostatic kernel, we
368	<	replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$.
369	<	Many of the derivatives of the damped kernel are well known from
370	<	Smith's early work on multipoles for the Ewald
371	<	summation.\cite{Smith82,Smith98}
367	>	% To carry out the same procedure for a damped electrostatic kernel, we
368	>	% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$.
369	>	% Many of the derivatives of the damped kernel are well known from
370	>	% Smith's early work on multipoles for the Ewald
371	>	% summation.\cite{Smith82,Smith98}
372
373	<	Note that increasing the value of $n$ will add additional terms to the
374	<	electrostatic potential, e.g., $f_2(r)$ includes orders up to
375	<	$(r-R_c)^3/R_c^4$, and so on.  Successive derivatives of the $f_n(r)$
376	<	functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) =
377	<	f^{\prime\prime}_2(r)$, etc.  These higher derivatives are required
378	<	for computing multipole energies, forces, and torques, and smooth
379	<	cutoffs of these quantities can be guaranteed as long as the number of
380	<	terms in the Taylor series exceeds the derivative order required.
373	>	% Note that increasing the value of $n$ will add additional terms to the
374	>	% electrostatic potential, e.g., $f_2(r)$ includes orders up to
375	>	% $(r-r_c)^3/r_c^4$, and so on.  Successive derivatives of the $f_n(r)$
376	>	% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) =
377	>	% f^{\prime\prime}_2(r)$, etc.  These higher derivatives are required
378	>	% for computing multipole energies, forces, and torques, and smooth
379	>	% cutoffs of these quantities can be guaranteed as long as the number of
380	>	% terms in the Taylor series exceeds the derivative order required.
381
382		For multipole-multipole interactions, following this procedure results
383	<	in separate radial functions for each distinct orientational
384	<	contribution to the potential, and ensures that the forces and torques
385	<	from {\it each} of these contributions will vanish at the cutoff
386	<	radius.  For example, the direct dipole dot product ($\mathbf{D}_{i}
387	<	\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance
383	>	in separate radial functions for each of the distinct orientational
384	>	contributions to the potential, and ensures that the forces and
385	>	torques from each of these contributions will vanish at the cutoff
386	>	radius.  For example, the direct dipole dot product
387	>	($\mathbf{D}_{a}
388	>	\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance
389		dot products:
390		\begin{equation}
391	<	U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot
392	<	\mathbf{D}_{j} \right) \frac{g_2(r)}{r}
393	<	-\frac{1}{4\pi \epsilon_0}
394	<	\left( \mathbf{D}_{i} \cdot \hat{r} \right)
395	<	\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) -
372	<	\frac{g_2(r)}{r} \right)
391	>	U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
392	>	\mathbf{D}_{a} \cdot
393	>	\mathbf{D}_{b} \right) v_{21}(r) +
394	>	\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right)
395	>	\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right]
396		\end{equation}
397
398	<	The electrostatic forces and torques acting on the central multipole
399	<	site due to another site within cutoff sphere are derived from
398	>	For the Taylor shifted (TSF) method with the undamped kernel,
399	>	$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} +
400	>	\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4}
401	>	- \frac{6}{r r_c^2}$.  In these functions, one can easily see the
402	>	connection to unmodified electrostatics as well as the smooth
403	>	transition to zero in both these functions as $r\rightarrow r_c$.  The
404	>	electrostatic forces and torques acting on the central multipole due
405	>	to another site within the cutoff sphere are derived from
406		Eq.~\ref{generic}, accounting for the appropriate number of
407		derivatives. Complete energy, force, and torque expressions are
408		presented in the first paper in this series (Reference
409	<	\citep{PaperI}).
409	>	\onlinecite{PaperI}).
410
411		\subsection{Gradient-shifted force (GSF)}
412
413	<	A second (and significantly simpler) method involves shifting the
414	<	gradient of the raw coulomb potential for each particular multipole
413	>	A second (and conceptually simpler) method involves shifting the
414	>	gradient of the raw Coulomb potential for each particular multipole
415		order.  For example, the raw dipole-dipole potential energy may be
416		shifted smoothly by finding the gradient for two interacting dipoles
417		which have been projected onto the surface of the cutoff sphere
418		without changing their relative orientation,
419	<	\begin{displaymath}
420	<	U_{D_{i}D_{j}}(r_{ij})  = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(R_c)
421	<	- (r_{ij}-R_c) \hat{r}_{ij} \cdot
422	<	\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{R_c}
423	<	\end{displaymath}
424	<	Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$
425	<	and $\mathbf{D}_{j}$, are retained at the cutoff distance (although
426	<	the signs are reversed for the dipole that has been projected onto the
427	<	cutoff sphere).  In many ways, this simpler approach is closer in
428	<	spirit to the original shifted force method, in that it projects a
429	<	neutralizing multipole (and the resulting forces from this multipole)
430	<	onto a cutoff sphere. The resulting functional forms for the
431	<	potentials, forces, and torques turn out to be quite similar in form
432	<	to the Taylor-shifted approach, although the radial contributions are
433	<	significantly less perturbed by the Gradient-shifted approach than
434	<	they are in the Taylor-shifted method.
419	>	\begin{equation}
420	>	U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
421	>	U_{D_{a}D_{b}}(r_c)
422	>	- (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot
423	>	\nabla U_{D_{a}D_{b}}(r_c).
424	>	\end{equation}
425	>	Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance
426	>	(although the signs are reversed for the dipole that has been
427	>	projected onto the cutoff sphere).  In many ways, this simpler
428	>	approach is closer in spirit to the original shifted force method, in
429	>	that it projects a neutralizing multipole (and the resulting forces
430	>	from this multipole) onto a cutoff sphere. The resulting functional
431	>	forms for the potentials, forces, and torques turn out to be quite
432	>	similar in form to the Taylor-shifted approach, although the radial
433	>	contributions are significantly less perturbed by the gradient-shifted
434	>	approach than they are in the Taylor-shifted method.
435
436	+	For the gradient shifted (GSF) method with the undamped kernel,
437	+	$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and
438	+	$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$.
439	+	Again, these functions go smoothly to zero as $r\rightarrow r_c$, and
440	+	because the Taylor expansion retains only one term, they are
441	+	significantly less perturbed than the TSF functions.
442	+
443		In general, the gradient shifted potential between a central multipole
444		and any multipolar site inside the cutoff radius is given by,
445		\begin{equation}
446	<	U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
447	<	U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r}
448	<	\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert  _{r_c} \right]
446	>	U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
447	>	U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c)
448	>	\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
449		\label{generic2}
450		\end{equation}
451		where the sum describes a separate force-shifting that is applied to
452	<	each orientational contribution to the energy.
452	>	each orientational contribution to the energy.  In this expression,
453	>	$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles
454	>	($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$
455	>	represent the orientations the multipoles.
456
457		The third term converges more rapidly than the first two terms as a
458		function of radius, hence the contribution of the third term is very
459		small for large cutoff radii.  The force and torque derived from
460	<	equation \ref{generic2} are consistent with the energy expression and
461	<	approach zero as $r \rightarrow R_c$.  Both the GSF and TSF methods
460	>	Eq. \ref{generic2} are consistent with the energy expression and
461	>	approach zero as $r \rightarrow r_c$.  Both the GSF and TSF methods
462		can be considered generalizations of the original DSF method for
463		higher order multipole interactions. GSF and TSF are also identical up
464		to the charge-dipole interaction but generate different expressions in
465		the energy, force and torque for higher order multipole-multipole
466		interactions. Complete energy, force, and torque expressions for the
467		GSF potential are presented in the first paper in this series
468	<	(Reference \citep{PaperI})
468	>	(Reference~\onlinecite{PaperI}).
469
470
471		\subsection{Shifted potential (SP) }
#	Line 439 | Line 478 | U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c)
478		interactions with the central multipole and the image. This
479		effectively shifts the total potential to zero at the cutoff radius,
480		\begin{equation}
481	<	U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c)
481	>	U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
482	>	U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
483		\label{eq:SP}
484		\end{equation}
485		where the sum describes separate potential shifting that is done for
486		each orientational contribution to the energy (e.g. the direct dipole
487		product contribution is shifted {\it separately} from the
488		dipole-distance terms in dipole-dipole interactions).  Note that this
489	<	is not a simple shifting of the total potential at $R_c$. Each radial
489	>	is not a simple shifting of the total potential at $r_c$. Each radial
490		contribution is shifted separately.  One consequence of this is that
491		multipoles that reorient after leaving the cutoff sphere can re-enter
492		the cutoff sphere without perturbing the total energy.
493
494	<	The potential energy between a central multipole and other multipolar
495	<	sites then goes smoothly to zero as $r \rightarrow R_c$. However, the
496	<	force and torque obtained from the shifted potential (SP) are
497	<	discontinuous at $R_c$. Therefore, MD simulations will still
498	<	experience energy drift while operating under the SP potential, but it
499	<	may be suitable for Monte Carlo approaches where the configurational
500	<	energy differences are the primary quantity of interest.
494	>	For the shifted potential (SP) method with the undamped kernel,
495	>	$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) =
496	>	\frac{3}{r^3} - \frac{3}{r_c^3}$.  The potential energy between a
497	>	central multipole and other multipolar sites goes smoothly to zero as
498	>	$r \rightarrow r_c$.  However, the force and torque obtained from the
499	>	shifted potential (SP) are discontinuous at $r_c$.  MD simulations
500	>	will still experience energy drift while operating under the SP
501	>	potential, but it may be suitable for Monte Carlo approaches where the
502	>	configurational energy differences are the primary quantity of
503	>	interest.
504
505	<	\subsection{The Self term}
505	>	\subsection{The Self Term}
506		In the TSF, GSF, and SP methods, a self-interaction is retained for
507		the central multipole interacting with its own image on the surface of
508		the cutoff sphere.  This self interaction is nearly identical with the
509		self-terms that arise in the Ewald sum for multipoles.  Complete
510		expressions for the self terms are presented in the first paper in
511	<	this series (Reference \citep{PaperI})
511	>	this series (Reference \onlinecite{PaperI}).
512
513
514		\section{\label{sec:methodology}Methodology}
#	Line 477 | Line 520 | arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} Thi
520		real-space cutoffs.  In the first paper of this series, we compared
521		the dipolar and quadrupolar energy expressions against analytic
522		expressions for ordered dipolar and quadrupolar
523	<	arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} This work uses the
524	<	multipolar Ewald sum as a reference method for comparing energies,
525	<	forces, and torques for molecular models that mimic disordered and
526	<	ordered condensed-phase systems.  These test-cases include:
523	>	arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} In this work, we
524	>	used the multipolar Ewald sum as a reference method for comparing
525	>	energies, forces, and torques for molecular models that mimic
526	>	disordered and ordered condensed-phase systems.  The parameters used
527	>	in the test cases are given in table~\ref{tab:pars}.
528
529	<	\begin{itemize}
530	<	\item Soft Dipolar fluids ($\sigma = , \epsilon = , |D| = $)
531	<	\item Soft Dipolar solids ($\sigma = , \epsilon = , |D| = $)
532	<	\item Soft Quadrupolar fluids ($\sigma = , \epsilon = , Q_{xx} = ...$)
533	<	\item Soft Quadrupolar solids  ($\sigma = , \epsilon = , Q_{xx} = ...$)
534	<	\item A mixed multipole model for water
535	<	\item A mixed multipole models for water with dissolved ions
536	<	\end{itemize}
537	<	This last test case exercises all levels of the multipole-multipole
538	<	interactions we have derived so far and represents the most complete
539	<	test of the new methods.
529	>	\begin{table}
530	>	\label{tab:pars}
531	>	\caption{The parameters used in the systems used to evaluate the new
532	>	real-space methods.  The most comprehensive test was a liquid
533	>	composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-}
534	>	ions).  This test excercises all orders of the multipolar
535	>	interactions developed in the first paper.}
536	>	\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline
537	>	& \multicolumn{2}{c|}{LJ parameters} &
538	>	\multicolumn{5}{c|}{Electrostatic moments} & & & & \\
539	>	Test system & $\sigma$& $\epsilon$ & $C$ & $D$  &
540	>	$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass  & $I_{xx}$ & $I_{yy}$ &
541	>	$I_{zz}$ \\ \cline{6-8}\cline{10-12}
542	>	& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu
543	>	\AA\textsuperscript{2})} \\ \hline
544	>	Soft Dipolar fluid & 3.051 & 0.152 &  & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\
545	>	Soft Dipolar solid & 2.837 & 1.0   &  & 2.35 & & & & $10^4$  & 17.6 &17.6 & 0 \\
546	>	Soft Quadrupolar fluid & 3.051 & 0.152 &  &  & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155  \\
547	>	Soft Quadrupolar solid & 2.837 & 1.0   &  &  & -1&-1&-2.5 & $10^4$  & 17.6&17.6&0 \\
548	>	SSDQ water  & 3.051 & 0.152 &  & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\
549	>	\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\
550	>	\ce{Cl-} & 4.445 & 0.1   & -1& & & & & 35.4527& & & \\ \hline
551	>	\end{tabularx}
552	>	\end{table}
553	>	The systems consist of pure multipolar solids (both dipole and
554	>	quadrupole), pure multipolar liquids (both dipole and quadrupole), a
555	>	fluid composed of sites containing both dipoles and quadrupoles
556	>	simultaneously, and a final test case that includes ions with point
557	>	charges in addition to the multipolar fluid.  The solid-phase
558	>	parameters were chosen so that the systems can explore some
559	>	orientational freedom for the multipolar sites, while maintaining
560	>	relatively strict translational order.  The SSDQ model used here is
561	>	not a particularly accurate water model, but it does test
562	>	dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole
563	>	interactions at roughly the same magnitudes. The last test case, SSDQ
564	>	water with dissolved ions, exercises \textit{all} levels of the
565	>	multipole-multipole interactions we have derived so far and represents
566	>	the most complete test of the new methods.
567
568		In the following section, we present results for the total
569		electrostatic energy, as well as the electrostatic contributions to
570		the force and torque on each molecule.  These quantities have been
571		computed using the SP, TSF, and GSF methods, as well as a hard cutoff,
572	<	and have been compared with the values obtaine from the multipolar
573	<	Ewald sum.  In Mote Carlo (MC) simulations, the energy differences
572	>	and have been compared with the values obtained from the multipolar
573	>	Ewald sum.  In Monte Carlo (MC) simulations, the energy differences
574		between two configurations is the primary quantity that governs how
575	<	the simulation proceeds. These differences are the most imporant
575	>	the simulation proceeds. These differences are the most important
576		indicators of the reliability of a method even if the absolute
577		energies are not exact.  For each of the multipolar systems listed
578		above, we have compared the change in electrostatic potential energy
#	Line 509 | Line 580 | contributions to the forces and torques.
580		molecular dynamics (MD) simulations, the forces and torques govern the
581		behavior of the simulation, so we also compute the electrostatic
582		contributions to the forces and torques.
583	+
584	+	\subsection{Implementation}
585	+	The real-space methods developed in the first paper in this series
586	+	have been implemented in our group's open source molecular simulation
587	+	program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in
588	+	this work.  The complementary error function can be a relatively slow
589	+	function on some processors, so all of the radial functions are
590	+	precomputed on a fine grid and are spline-interpolated to provide
591	+	values when required.
592
593	<	\subsection{Model systems}
594	<	To sample independent configurations of multipolar crystals, a body
595	<	centered cubic (BCC) crystal which is a minimum energy structure for
596	<	point dipoles was generated using 3,456 molecules.  The multipoles
597	<	were translationally locked in their respective crystal sites for
598	<	equilibration at a relatively low temperature (50K), so that dipoles
599	<	or quadrupoles could freely explore all accessible orientations.  The
600	<	translational constraints were removed, and the crystals were
601	<	simulated for 10 ps in the microcanonical (NVE) ensemble with an
602	<	average temperature of 50 K.  Configurations were sampled at equal
603	<	time intervals for the comparison of the configurational energy
524	<	differences.  The crystals were not simulated close to the melting
525	<	points in order to avoid translational deformation away of the ideal
526	<	lattice geometry.
593	>	Using the same simulation code, we compare to a multipolar Ewald sum
594	>	with a reciprocal space cutoff, $k_\mathrm{max} = 7$.  Our version of
595	>	the Ewald sum is a re-implementation of the algorithm originally
596	>	proposed by Smith that does not use the particle mesh or smoothing
597	>	approximations.\cite{Smith82,Smith98} This implementation was tested
598	>	extensively against the analytic energy constants for the multipolar
599	>	lattices that are discussed in reference \onlinecite{PaperI}.  In all
600	>	cases discussed below, the quantities being compared are the
601	>	electrostatic contributions to energies, force, and torques.  All
602	>	other contributions to these quantities (i.e. from Lennard-Jones
603	>	interactions) are removed prior to the comparisons.
604
605	<	For dipolar, quadrupolar, and mixed-multipole liquid simulations, each
606	<	system was created with 2048 molecules oriented randomly.  These were
605	>	The convergence parameter ($\alpha$) also plays a role in the balance
606	>	of the real-space and reciprocal-space portions of the Ewald
607	>	calculation.  Typical molecular mechanics packages set this to a value
608	>	that depends on the cutoff radius and a tolerance (typically less than
609	>	$1 \times 10^{-4}$ kcal/mol).  Smaller tolerances are typically
610	>	associated with increasing accuracy at the expense of computational
611	>	time spent on the reciprocal-space portion of the
612	>	summation.\cite{Perram88,Essmann:1995pb} A default tolerance of $1 \times
613	>	10^{-8}$ kcal/mol was used in all Ewald calculations, resulting in
614	>	Ewald coefficient 0.3119 \AA$^{-1}$ for a cutoff radius of 12 \AA.
615
616	<	system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K
617	<	temperature after equilibration.  We collected 250 different
618	<	configurations in equal interval of time. For the ions mixed liquid
619	<	system, we converted 48 different molecules into 24 \ce{Na+} and 24
620	<	\ce{Cl-} ions and equilibrated. After equilibration, the system was run
536	<	at the same environment for 1ns and 250 configurations were
537	<	collected. While comparing energies, forces, and torques with Ewald
538	<	method, Lennard-Jones potentials were turned off and purely
539	<	electrostatic interaction had been compared.
616	>	The real-space models have self-interactions that provide
617	>	contributions to the energies only.  Although the self interaction is
618	>	a rapid calculation, we note that in systems with fluctuating charges
619	>	or point polarizabilities, the self-term is not static and must be
620	>	recomputed at each time step.
621
622	+	\subsection{Model systems}
623	+	To sample independent configurations of the multipolar crystals, body
624	+	centered cubic (bcc) crystals, which exhibit the minimum energy
625	+	structures for point dipoles, were generated using 3,456 molecules.
626	+	The multipoles were translationally locked in their respective crystal
627	+	sites for equilibration at a relatively low temperature (50K) so that
628	+	dipoles or quadrupoles could freely explore all accessible
629	+	orientations.  The translational constraints were then removed, the
630	+	systems were re-equilibrated, and the crystals were simulated for an
631	+	additional 10 ps in the microcanonical (NVE) ensemble with an average
632	+	temperature of 50 K.  The balance between moments of inertia and
633	+	particle mass were chosen to allow orientational sampling without
634	+	significant translational motion.  Configurations were sampled at
635	+	equal time intervals in order to compare configurational energy
636	+	differences.  The crystals were simulated far from the melting point
637	+	in order to avoid translational deformation away of the ideal lattice
638	+	geometry.
639	+
640	+	For dipolar, quadrupolar, and mixed-multipole \textit{liquid}
641	+	simulations, each system was created with 2,048 randomly-oriented
642	+	molecules.  These were equilibrated at a temperature of 300K for 1 ns.
643	+	Each system was then simulated for 1 ns in the microcanonical (NVE)
644	+	ensemble.  We collected 250 different configurations at equal time
645	+	intervals. For the liquid system that included ionic species, we
646	+	converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24
647	+	\ce{Cl-} ions and re-equilibrated. After equilibration, the system was
648	+	run under the same conditions for 1 ns. A total of 250 configurations
649	+	were collected. In the following comparisons of energies, forces, and
650	+	torques, the Lennard-Jones potentials were turned off and only the
651	+	purely electrostatic quantities were compared with the same values
652	+	obtained via the Ewald sum.
653	+
654		\subsection{Accuracy of Energy Differences, Forces and Torques}
655		The pairwise summation techniques (outlined above) were evaluated for
656		use in MC simulations by studying the energy differences between
#	Line 550 | Line 663 | we used least square regressions analysiss for the six
663		should be identical for all methods.
664
665		Since none of the real-space methods provide exact energy differences,
666	<	we used least square regressions analysiss for the six different
666	>	we used least square regressions analysis for the six different
667		molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF
668		with the multipolar Ewald reference method.  Unitary results for both
669		the correlation (slope) and correlation coefficient for these
#	Line 561 | Line 674 | also been compared by using least squares regression a
674		configurations and 250 configurations were recorded for comparison.
675		Each system provided 31,125 energy differences for a total of 186,750
676		data points.  Similarly, the magnitudes of the forces and torques have
677	<	also been compared by using least squares regression analyses. In the
677	>	also been compared using least squares regression analysis. In the
678		forces and torques comparison, the magnitudes of the forces acting in
679		each molecule for each configuration were evaluated. For example, our
680		dipolar liquid simulation contains 2048 molecules and there are 250
#	Line 647 | Line 760 | model must allow for long simulation times with minima
760
761		%        \label{fig:barGraph2}
762		%      \end{figure}
763	<	%The correlation coefficient ($R^2$) and slope of the linear regression plots for the energy differences for all six different molecular systems is shown in figure 4a and 4b.The plot shows that the correlation coefficient improves for the SP cutoff method as compared to the undamped hard cutoff method in the case of SSDQC, SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar crystal and liquid, the correlation coefficient is almost unchanged and close to 1.  The correlation coefficient is smallest (0.696276 for $r_c$ = 9 $A^o$) for the SSDQC liquid because of the presence of charge-charge and charge-multipole interactions. Since the charge-charge and charge-multipole interaction is long ranged, there is huge deviation of correlation coefficient from 1. Similarly, the quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with compared to interactions in the other multipolar systems, thus the correlation coefficient very close to 1 even for hard cutoff method. The idea of placing image multipole on the surface of the cutoff sphere improves the correlation coefficient and makes it close to 1 for all types of multipolar systems. Similarly the slope is hugely deviated from the correct value for the lower order multipole-multipole interaction and slightly deviated for higher order multipole – multipole interaction. The SP method improves both correlation coefficient ($R^2$) and slope significantly in SSDQC and dipolar systems.  The Slope is found to be deviated more in dipolar crystal as compared to liquid which is associated with the large fluctuation in the electrostatic energy in crystal. The GSF also produced better values of correlation coefficient and slope with the proper selection of the damping alpha (Interested reader can consult accompanying supporting material). The TSF method gives good value of correlation coefficient for the dipolar crystal, dipolar liquid, SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the regression slopes are significantly deviated.
763	>	%The correlation coefficient ($R^2$) and slope of the linear
764	>	%regression plots for the energy differences for all six different
765	>	%molecular systems is shown in figure 4a and 4b.The plot shows that
766	>	%the correlation coefficient improves for the SP cutoff method as
767	>	%compared to the undamped hard cutoff method in the case of SSDQC,
768	>	%SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar
769	>	%crystal and liquid, the correlation coefficient is almost unchanged
770	>	%and close to 1.  The correlation coefficient is smallest (0.696276
771	>	%for $r_c$ = 9 $A^\circ$) for the SSDQC liquid because of the presence of
772	>	%charge-charge and charge-multipole interactions. Since the
773	>	%charge-charge and charge-multipole interaction is long ranged, there
774	>	%is huge deviation of correlation coefficient from 1. Similarly, the
775	>	%quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with
776	>	%compared to interactions in the other multipolar systems, thus the
777	>	%correlation coefficient very close to 1 even for hard cutoff
778	>	%method. The idea of placing image multipole on the surface of the
779	>	%cutoff sphere improves the correlation coefficient and makes it close
780	>	%to 1 for all types of multipolar systems. Similarly the slope is
781	>	%hugely deviated from the correct value for the lower order
782	>	%multipole-multipole interaction and slightly deviated for higher
783	>	%order multipole – multipole interaction. The SP method improves both
784	>	%correlation coefficient ($R^2$) and slope significantly in SSDQC and
785	>	%dipolar systems.  The Slope is found to be deviated more in dipolar
786	>	%crystal as compared to liquid which is associated with the large
787	>	%fluctuation in the electrostatic energy in crystal. The GSF also
788	>	%produced better values of correlation coefficient and slope with the
789	>	%proper selection of the damping alpha (Interested reader can consult
790	>	%accompanying supporting material). The TSF method gives good value of
791	>	%correlation coefficient for the dipolar crystal, dipolar liquid,
792	>	%SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the
793	>	%regression slopes are significantly deviated.
794	>
795		\begin{figure}
796	<	\centering
797	<	\includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf}
798	<	\caption{The correlation coefficient and regression slope of configurational energy differences for a given method with compared with the reference Ewald method. The value of result equal to 1(dashed line) indicates energy difference is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\  = circle, 12 \AA\  = square 15 \AA\  = inverted triangle)}
799	<	\label{fig:slopeCorr_energy}
800	<	\end{figure}
801	<	The combined correlation coefficient and slope for all six systems is shown in Figure ~\ref{fig:slopeCorr_energy}. The correlation coefficient for the undamped hard cutoff method is does not have good agreement with the Ewald because of the fluctuation of the electrostatic energy in the direct truncation method. This deviation in correlation coefficient is improved by using SP, GSF, and TSF method. But the TSF method worsens the regression slope stating that this method produces statistically more biased result as compared to Ewald. Also the GSF method slightly deviate slope but it can be alleviated by using proper value of damping alpha and cutoff radius. The SP method shows good agreement with Ewald method for all values of damping alpha and radii.
802	<	\subsection{Magnitude of the force and torque vectors}
803	<	The comparison of the magnitude of the combined forces and torques for the data accumulated from all system types are shown in Figure ~\ref{fig:slopeCorr_force}. The correlation and slope for the forces agree with the Ewald even for the hard cutoff method. For the system of molecules with higher order multipoles, the interaction is short ranged. Moreover, the force decays more rapidly than the electrostatic energy hence the hard cutoff method also produces good results. Although the pure cutoff gives the good match of the electrostatic force, the discontinuity in the force at the cutoff radius causes problem in the total energy conservation in MD simulations, which will be discussed in detail in subsection D. The correlation coefficient for GSF method also perfectly matches with Ewald but the slope is slightly deviated (due to extra term obtained from the angular differentiation). This deviation in the slope can be alleviated with proper selection of the damping alpha and radii ($\alpha = 0.2$ and $r_c = 12 A^o$ are good choice). The TSF method shows good agreement in the correlation coefficient but the slope is not good as compared to the Ewald.
804	<	\begin{figure}
805	<	\centering
806	<	\includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf}
807	<	\caption{The correlation coefficient and regression slope of the magnitude of the force for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\  = circle, 12 \AA\  = square 15 \AA\  = inverted triangle). }
808	<	\label{fig:slopeCorr_force}
809	<	\end{figure}
810	<	The torques appears to be very influenced because of extra term generated when the potential energy is modified to get consistent force and torque.  The result shows that the torque from the hard cutoff method has good agreement with Ewald. As the potential is modified to make it consistent with the force and torque, the correlation and slope is deviated as shown in Figure~\ref{fig:slopeCorr_torque} for SP, GSF and TSF cutoff methods.  But the proper value of the damping alpha and radius can improve the agreement of the GSF with the Ewald method. The TSF method shows worst agreement in the slope as compared to Ewald even for larger cutoff radii.
811	<	\begin{figure}
812	<	\centering
813	<	\includegraphics[width=0.5 \textwidth]{torquePlot_slopeCorrelation_combined-crop.pdf}
814	<	\caption{The correlation coefficient and regression slope of the magnitude of the torque for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle).}
815	<	\label{fig:slopeCorr_torque}
816	<	\end{figure}
817	<	\subsection{Directionality of the force and torque vectors}
818	<	The accurate evaluation of the direction of the force and torques are also important for the dynamic simulation.In our research, the direction data sets were computed from the purposed method and compared with Ewald using Fisher statistics and results are expressed in terms of circular variance ($Var(\theta$).The force and torque vectors from the purposed method followed Fisher probability distribution function expressed in equation~\ref{eq:pdf}. The circular variance for the force and torque vectors of each molecule in the 250 configurations for all system types is shown in Figure~\ref{fig:slopeCorr_circularVariance}. The direction of the force and torque vectors from hard and SP cutoff methods showed best directional agreement with the Ewald. The force and torque vectors from GSF method also showed good agreement with the Ewald method, which can also be improved by varying damping alpha and cutoff radius.For $\alpha = 0.2$ and $r_c = 12 A^o$, $ Var(\theta) $ for direction of the force was found to be 0.002061 and corresponding value of $\kappa $ was 485.20. Integration of equation ~\ref{eq:pdf} for that corresponding value of $\kappa$ showed that 95\% of force vectors are with in $6.37^o$. The TSF method is the poorest in evaluating accurate direction with compared to Hard, SP, and GSF methods. The circular variance for the direction of the torques is larger as compared to force. For same $\alpha = 0.2, r_c = 12 A^o$ and GSF method, the circular variance was 0.01415, which showed 95\% of torque vectors are within $16.75^o$.The direction of the force and torque vectors can be improved by varying $\alpha$ and $r_c$.
796	>	\centering
797	>	\includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps}
798	>	\caption{Statistical analysis of the quality of configurational
799	>	energy differences for the real-space electrostatic methods
800	>	compared with the reference Ewald sum.  Results with a value equal
801	>	to 1 (dashed line) indicate $\Delta E$ values indistinguishable
802	>	from those obtained using the multipolar Ewald sum.  Different
803	>	values of the cutoff radius are indicated with different symbols
804	>	(9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted
805	>	triangles).}
806	>	\label{fig:slopeCorr_energy}
807	>	\end{figure}
808	>
809	>	The combined correlation coefficient and slope for all six systems is
810	>	shown in Figure ~\ref{fig:slopeCorr_energy}.  Most of the methods
811	>	reproduce the Ewald configurational energy differences with remarkable
812	>	fidelity.  Undamped hard cutoffs introduce a significant amount of
813	>	random scatter in the energy differences which is apparent in the
814	>	reduced value of the correlation coefficient for this method.  This
815	>	can be easily understood as configurations which exhibit small
816	>	traversals of a few dipoles or quadrupoles out of the cutoff sphere
817	>	will see large energy jumps when hard cutoffs are used.  The
818	>	orientations of the multipoles (particularly in the ordered crystals)
819	>	mean that these energy jumps can go in either direction, producing a
820	>	significant amount of random scatter, but no systematic error.
821
822	+	The TSF method produces energy differences that are highly correlated
823	+	with the Ewald results, but it also introduces a significant
824	+	systematic bias in the values of the energies, particularly for
825	+	smaller cutoff values. The TSF method alters the distance dependence
826	+	of different orientational contributions to the energy in a
827	+	non-uniform way, so the size of the cutoff sphere can have a large
828	+	effect, particularly for the crystalline systems.
829	+
830	+	Both the SP and GSF methods appear to reproduce the Ewald results with
831	+	excellent fidelity, particularly for moderate damping ($\alpha =
832	+	0.1-0.2$\AA$^{-1}$) and with a commonly-used cutoff value ($r_c =
833	+	12$\AA).  With the exception of the undamped hard cutoff, and the TSF
834	+	method with short cutoffs, all of the methods would be appropriate for
835	+	use in Monte Carlo simulations.
836	+
837	+	\subsection{Magnitude of the force and torque vectors}
838	+
839	+	The comparisons of the magnitudes of the forces and torques for the
840	+	data accumulated from all six systems are shown in Figures
841	+	~\ref{fig:slopeCorr_force} and \ref{fig:slopeCorr_torque}. The
842	+	correlation and slope for the forces agree well with the Ewald sum
843	+	even for the hard cutoffs.
844	+
845	+	For systems of molecules with only multipolar interactions, the pair
846	+	energy contributions are quite short ranged.  Moreover, the force
847	+	decays more rapidly than the electrostatic energy, hence the hard
848	+	cutoff method can also produce reasonable agreement for this quantity.
849	+	Although the pure cutoff gives reasonably good electrostatic forces
850	+	for pairs of molecules included within each other's cutoff spheres,
851	+	the discontinuity in the force at the cutoff radius can potentially
852	+	cause energy conservation problems as molecules enter and leave the
853	+	cutoff spheres.  This is discussed in detail in section
854	+	\ref{sec:conservation}.
855	+
856	+	The two shifted-force methods (GSF and TSF) exhibit a small amount of
857	+	systematic variation and scatter compared with the Ewald forces.  The
858	+	shifted-force models intentionally perturb the forces between pairs of
859	+	molecules inside each other's cutoff spheres in order to correct the
860	+	energy conservation issues, and this perturbation is evident in the
861	+	statistics accumulated for the molecular forces.  The GSF
862	+	perturbations are minimal, particularly for moderate damping and
863	+	commonly-used cutoff values ($r_c = 12$\AA).  The TSF method shows
864	+	reasonable agreement in the correlation coefficient but again the
865	+	systematic error in the forces is concerning if replication of Ewald
866	+	forces is desired.
867	+
868		\begin{figure}
869	<	\centering
870	<	\includegraphics[width=0.5 \textwidth]{Variance_forceNtorque_modified-crop.pdf}
871	<	\caption{The circular variance of the data sets of the direction of the  force and torque vectors obtained from a given method about reference Ewald method. The result equal to 0 (dashed line) indicates direction of the vectors are indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle)}
872	<	\label{fig:slopeCorr_circularVariance}
873	<	\end{figure}
874	<	\subsection{Total energy conservation}
875	<	We have tested the conservation of energy in the SSDQC liquid system by running system for 1ns in the Hard, SP, GSF and TSF method. The Hard cutoff method shows very high energy drifts 433.53 KCal/Mol/ns/particle and energy fluctuation 32574.53 Kcal/Mol (measured by the SD from the slope) for the undamped case, which makes it completely unusable in MD simulations. The SP method also shows large value of energy drift 1.289 Kcal/Mol/ns/particle and energy fluctuation 0.01953 Kcal/Mol. The energy fluctuation in the SP method is due to the non-vanishing nature of the torque and force at the cutoff radius. We can improve the energy conservation in some extent by the proper selection of the damping alpha but the improvement is not good enough, which can be observed in Figure 9a and 9b .The GSF and TSF shows very low value of energy drift 0.09016, 0.07371 KCal/Mol/ns/particle and fluctuation 3.016E-6, 1.217E-6 Kcal/Mol respectively for the undamped case. Since the absolute value of the evaluated electrostatic energy, force and torque from TSF method are deviated from the Ewald, it does not mimic MD simulations appropriately. The electrostatic energy, force and torque from the GSF method have very good agreement with the Ewald. In addition, the energy drift and energy fluctuation from the GSF method is much better than Ewald’s method for reciprocal space vector value ($k_f$) equal to 7 as shown in Figure~\ref{fig:energyDrift} and ~\ref{fig:fluctuation}. We can improve the total energy fluctuation and drift for the Ewald’s method by increasing size of the reciprocal space, which extremely increseses the simulation time. In our current simulation, the simulation time for the Hard, SP, and GSF methods are about 5.5 times faster than the Ewald method.
869	>	\centering
870	>	\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps}
871	>	\caption{Statistical analysis of the quality of the force vector
872	>	magnitudes for the real-space electrostatic methods compared with
873	>	the reference Ewald sum. Results with a value equal to 1 (dashed
874	>	line) indicate force magnitude values indistinguishable from those
875	>	obtained using the multipolar Ewald sum.  Different values of the
876	>	cutoff radius are indicated with different symbols (9\AA\ =
877	>	circles, 12\AA\ = squares, and 15\AA\ = inverted triangles). }
878	>	\label{fig:slopeCorr_force}
879	>	\end{figure}
880	>
881	>
882		\begin{figure}
883	<	\centering
884	<	\includegraphics[width=0.5 \textwidth]{log(energyDrift)-crop.pdf}
885	<	\label{fig:energyDrift}
886	<	\end{figure}
883	>	\centering
884	>	\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps}
885	>	\caption{Statistical analysis of the quality of the torque vector
886	>	magnitudes for the real-space electrostatic methods compared with
887	>	the reference Ewald sum. Results with a value equal to 1 (dashed
888	>	line) indicate force magnitude values indistinguishable from those
889	>	obtained using the multipolar Ewald sum.  Different values of the
890	>	cutoff radius are indicated with different symbols (9\AA\ =
891	>	circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
892	>	\label{fig:slopeCorr_torque}
893	>	\end{figure}
894	>
895	>	The torques (Fig. \ref{fig:slopeCorr_torque}) appear to be
896	>	significantly influenced by the choice of real-space method.  The
897	>	torque expressions have the same distance dependence as the energies,
898	>	which are naturally longer-ranged expressions than the inter-site
899	>	forces.  Torques are also quite sensitive to orientations of
900	>	neighboring molecules, even those that are near the cutoff distance.
901	>
902	>	The results shows that the torque from the hard cutoff method
903	>	reproduces the torques in quite good agreement with the Ewald sum.
904	>	The other real-space methods can cause some deviations, but excellent
905	>	agreement with the Ewald sum torques is recovered at moderate values
906	>	of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1}$) and cutoff
907	>	radius ($r_c \ge 12$\AA).  The TSF method exhibits only fair agreement
908	>	in the slope when compared with the Ewald torques even for larger
909	>	cutoff radii.  It appears that the severity of the perturbations in
910	>	the TSF method are most in evidence for the torques.
911	>
912	>	\subsection{Directionality of the force and torque vectors}
913	>
914	>	The accurate evaluation of force and torque directions is just as
915	>	important for molecular dynamics simulations as the magnitudes of
916	>	these quantities. Force and torque vectors for all six systems were
917	>	analyzed using Fisher statistics, and the quality of the vector
918	>	directionality is shown in terms of circular variance
919	>	($\mathrm{Var}(\theta)$) in figure
920	>	\ref{fig:slopeCorr_circularVariance}. The force and torque vectors
921	>	from the new real-space methods exhibit nearly-ideal Fisher probability
922	>	distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods
923	>	exhibit the best vectorial agreement with the Ewald sum. The force and
924	>	torque vectors from GSF method also show good agreement with the Ewald
925	>	method, which can also be systematically improved by using moderate
926	>	damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c =
927	>	12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds
928	>	to a distribution with 95\% of force vectors within $6.37^\circ$ of
929	>	the corresponding Ewald forces. The TSF method produces the poorest
930	>	agreement with the Ewald force directions.
931	>
932	>	Torques are again more perturbed than the forces by the new real-space
933	>	methods, but even here the variance is reasonably small.  For the same
934	>	method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA),
935	>	the circular variance was 0.01415, corresponds to a distribution which
936	>	has 95\% of torque vectors are within $16.75^\circ$ of the Ewald
937	>	results. Again, the direction of the force and torque vectors can be
938	>	systematically improved by varying $\alpha$ and $r_c$.
939	>
940		\begin{figure}
941	<	\centering
942	<	\includegraphics[width=0.5 \textwidth]{logSD-crop.pdf}
943	<	\caption{The plot showing (a) standard deviation, and (b) total energy drift in the total energy conservation plot for different values of the damping alpha for different cut off methods. }
944	<	\label{fig:fluctuation}
945	<	\end{figure}
941	>	\centering
942	>	\includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps}
943	>	\caption{The circular variance of the direction of the force and
944	>	torque vectors obtained from the real-space methods around the
945	>	reference Ewald vectors. A variance equal to 0 (dashed line)
946	>	indicates direction of the force or torque vectors are
947	>	indistinguishable from those obtained from the Ewald sum. Here
948	>	different symbols represent different values of the cutoff radius
949	>	(9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)}
950	>	\label{fig:slopeCorr_circularVariance}
951	>	\end{figure}
952	>
953	>	\subsection{Energy conservation\label{sec:conservation}}
954	>
955	>	We have tested the conservation of energy one can expect to see with
956	>	the new real-space methods using the SSDQ water model with a small
957	>	fraction of solvated ions. This is a test system which exercises all
958	>	orders of multipole-multipole interactions derived in the first paper
959	>	in this series and provides the most comprehensive test of the new
960	>	methods.  A liquid-phase system was created with 2000 water molecules
961	>	and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a
962	>	temperature of 300K.  After equilibration in the canonical (NVT)
963	>	ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system
964	>	was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald,
965	>	Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA.  The
966	>	value of the damping coefficient was also varied from the undamped
967	>	case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$
968	>	\AA$^{-1}$) for all of the real space methods.  A sample was also run
969	>	using the multipolar Ewald sum with the same real-space cutoff.
970	>
971	>	In figure~\ref{fig:energyDrift} we show the both the linear drift in
972	>	energy over time, $\delta E_1$, and the standard deviation of energy
973	>	fluctuations around this drift $\delta E_0$.  Both of the
974	>	shifted-force methods (GSF and TSF) provide excellent energy
975	>	conservation (drift less than $10^{-5}$ kcal / mol / ns / particle),
976	>	while the hard cutoff is essentially unusable for molecular dynamics.
977	>	SP provides some benefit over the hard cutoff because the energetic
978	>	jumps that happen as particles leave and enter the cutoff sphere are
979	>	somewhat reduced, but like the Wolf method for charges, the SP method
980	>	would not be as useful for molecular dynamics as either of the
981	>	shifted-force methods.
982	>
983	>	We note that for all tested values of the cutoff radius, the new
984	>	real-space methods can provide better energy conservation behavior
985	>	than the multipolar Ewald sum, even when utilizing a relatively large
986	>	$k$-space cutoff values.
987	>
988	>	\begin{figure}
989	>	\centering
990	>	\includegraphics[width=\textwidth]{newDrift_12.eps}
991	>	\label{fig:energyDrift}
992	>	\caption{Analysis of the energy conservation of the real-space
993	>	methods. $\delta \mathrm{E}_1$ is the linear drift in energy over
994	>	time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is
995	>	the standard deviation of energy fluctuations around this drift (in
996	>	kcal / mol / particle).  Points that appear below the dashed grey
997	>	(Ewald) lines exhibit better energy conservation than commonly-used
998	>	parameters for Ewald-based electrostatics.  All simulations were of
999	>	a 2000-molecule simulation of SSDQ water with 48 ionic charges at
1000	>	300 K starting from the same initial configuration. All runs
1001	>	utilized the same real-space cutoff, $r_c = 12$\AA.}
1002	>	\end{figure}
1003	>
1004	>	\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}}
1005	>	The most important test of the modified interaction potentials is the
1006	>	fidelity with which they can reproduce structural features and
1007	>	dynamical properties in a liquid.  One commonly-utilized measure of
1008	>	structural ordering is the pair distribution function, $g(r)$, which
1009	>	measures local density deviations in relation to the bulk density.  In
1010	>	the electrostatic approaches studied here, the short-range repulsion
1011	>	from the Lennard-Jones potential is identical for the various
1012	>	electrostatic methods, and since short range repulsion determines much
1013	>	of the local liquid ordering, one would not expect to see many
1014	>	differences in $g(r)$.  Indeed, the pair distributions are essentially
1015	>	identical for all of the electrostatic methods studied (for each of
1016	>	the different systems under investigation).  An example of this
1017	>	agreement for the SSDQ water/ion system is shown in
1018	>	Fig. \ref{fig:gofr}.
1019	>
1020	>	\begin{figure}
1021	>	\centering
1022	>	\includegraphics[width=\textwidth]{gofr_ssdqc.eps}
1023	>	\label{fig:gofr}
1024	>	\caption{The pair distribution functions, $g(r)$, for the SSDQ
1025	>	water/ion system obtained using the different real-space methods are
1026	>	essentially identical with the result from the Ewald
1027	>	treatment.}
1028	>	\end{figure}
1029	>
1030	>	There is a very slight overstructuring of the first solvation shell
1031	>	when using when using TSF at lower values of the damping coefficient
1032	>	($\alpha \le 0.1$) or when using undamped GSF.  With moderate damping,
1033	>	GSF and SP produce pair distributions that are identical (within
1034	>	numerical noise) to their Ewald counterparts.
1035	>
1036	>	A structural property that is a more demanding test of modified
1037	>	electrostatics is the mean value of the electrostatic energy $\langle
1038	>	U_\mathrm{elect} \rangle / N$ which is obtained by sampling the
1039	>	liquid-state configurations experienced by a liquid evolving entirely
1040	>	under the influence of each of the methods.  In table \ref{tab:Props}
1041	>	we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with
1042	>	the damping parameter, $\alpha$, for each of the methods.
1043	>
1044	>	As in the crystals studied in the first paper, damping is important
1045	>	for converging the mean electrostatic energy values, particularly for
1046	>	the two shifted force methods (GSF and TSF).  A value of $\alpha
1047	>	\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF
1048	>	energies with a cutoff of 12 \AA, while shorter cutoffs require more
1049	>	dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA).
1050	>	Overdamping the real-space electrostatic methods occurs with $\alpha >
1051	>	0.4$, causing the estimate of the energy to drop below the Ewald
1052	>	results.
1053	>
1054	>	These ``optimal'' values of the damping coefficient are slightly
1055	>	larger than what were observed for DSF electrostatics for purely
1056	>	point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for
1057	>	$r_c = 12$\AA appears to be an excellent compromise for mixed charge
1058	>	multipole systems.
1059	>
1060	>	To test the fidelity of the electrostatic methods at reproducing
1061	>	dynamics in a multipolar liquid, it is also useful to look at
1062	>	transport properties, particularly the diffusion constant,
1063	>	\begin{equation}
1064	>	D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left|
1065	>	\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle
1066	>	\label{eq:diff}
1067	>	\end{equation}
1068	>	which measures long-time behavior and is sensitive to the forces on
1069	>	the multipoles.  For the soft dipolar fluid and the SSDQ liquid
1070	>	systems, the self-diffusion constants (D) were calculated from linear
1071	>	fits to the long-time portion of the mean square displacement,
1072	>	$\langle r^{2}(t) \rangle$.\cite{Allen87}
1073	>
1074	>	In addition to translational diffusion, orientational relaxation times
1075	>	were calculated for comparisons with the Ewald simulations and with
1076	>	experiments. These values were determined from the same 1~ns
1077	>	microcanonical trajectories used for translational diffusion by
1078	>	calculating the orientational time correlation function,
1079	>	\begin{equation}
1080	>	C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t)
1081	>	\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle,
1082	>	\label{eq:OrientCorr}
1083	>	\end{equation}
1084	>	where $P_l$ is the Legendre polynomial of order $l$ and
1085	>	$\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis
1086	>	$\gamma$ on a molecule..  Th body-fixed reference frame used for our
1087	>	models has the $z$-axis running along the dipoles, and for the SSDQ
1088	>	water model, the $y$-axis connects the two implied hydrogen atom
1089	>	positions.  From the orientation autocorrelation functions, we can
1090	>	obtain time constants for rotational relaxation either by fitting an
1091	>	exponential function or by integrating the entire correlation
1092	>	function.  In a good water model, these decay times would be
1093	>	comparable to water orientational relaxation times from nuclear
1094	>	magnetic resonance (NMR). The relaxation constant obtained from
1095	>	$C_2^y(t)$ is normally of experimental interest because it describes
1096	>	the relaxation of the principle axis connecting the hydrogen
1097	>	atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
1098	>	of the dipole-dipole relaxation from a proton NMR signal and should
1099	>	provide an estimate of the NMR relaxation time constant.\cite{Impey82}
1100	>
1101	>	Results for the diffusion constants and orientational relaxation times
1102	>	are shown in figure \ref{tab:Props}. From this data, it is apparent
1103	>	that the values for both $D$ and $\tau_2$ using the Ewald sum are
1104	>	reproduced with reasonable fidelity by the GSF method.
1105	>
1106	>	The $\tau_2$ results in \ref{tab:Props} show a much greater difference
1107	>	between the real-space and the Ewald results.
1108	>
1109	>	\begin{table}
1110	>	\label{tab:Props}
1111	>	\caption{Comparison of the structural and dynamic properties for the
1112	>	soft dipolar liquid test for all of the real-space methods.}
1113	>	\begin{tabular}{l|c|cccc|cccc|cccc}
1114	>	& Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\
1115	>	$\alpha$ (\AA$^{-1}$) & &
1116	>	0.0 & 0.1 & 0.2 & 0.3 &
1117	>	0.0 & 0.1 & 0.2 & 0.3 &
1118	>	0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14}
1119	>
1120	>	$\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\
1121	>	D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)&
1122	>	470.2(6) &
1123	>	416.6(5) &
1124	>	379.6(5) &
1125	>	438.6(5) &
1126	>	476.0(6) &
1127	>	412.8(5) &
1128	>	421.1(5) &
1129	>	400.5(5) &
1130	>	437.5(6) &
1131	>	434.6(5) &
1132	>	411.4(5) &
1133	>	545.3(7) &
1134	>	459.6(6) \\
1135	>	$\tau_2$ (fs) &
1136	>	1.136 &
1137	>	1.041 &
1138	>	1.064 &
1139	>	1.109 &
1140	>	1.211 &
1141	>	1.119 &
1142	>	1.039 &
1143	>	1.058 &
1144	>	1.21  &
1145	>	1.15  &
1146	>	1.172 &
1147	>	1.153 &
1148	>	1.125 \\
1149	>	\end{tabular}
1150	>	\end{table}
1151	>
1152	>
1153		\section{CONCLUSION}
1154	<	We have generalized the charged neutralized potential energy originally developed by the Wolf et al.\cite{Wolf:1999dn} for the charge-charge interaction to the charge-multipole and multipole-multipole interaction in the SP method for higher order multipoles. Also, we have developed GSF and TSF methods by implementing the modification purposed by Fennel and Gezelter\cite{Fennell:2006lq} for the charge-charge interaction to the higher order multipoles to ensure consistency and smooth truncation of the electrostatic energy, force, and torque for the spherical truncation. The SP methods for multipoles proved its suitability in MC simulations. On the other hand, the results from the GSF method produced good agreement with the Ewald's energy, force, and torque. Also, it shows very good energy conservation in MD simulations.
1155	<	The direct truncation of any molecular system without multipole neutralization creates the fluctuation in the electrostatic energy. This fluctuation in the energy is very large for the case of crystal because of long range of multipole ordering (Refer paper I).\cite{PaperI} This is also significant in the case of the liquid because of the local multipole ordering in the molecules. If the net multipole within cutoff radius neutralized within cutoff sphere by placing image multiples on the surface of the sphere, this fluctuation in the energy reduced significantly. Also, the multipole neutralization in the generalized SP method showed very good agreement with the Ewald as compared to direct truncation for the evaluation of the $\triangle E$ between the configurations.
1156	<	In MD simulations, the energy conservation is very important. The
1157	<	conservation of the total energy can be ensured by  i) enforcing the
1158	<	smooth truncation of the energy, force and torque in the cutoff radius
1159	<	and ii) making the energy, force and torque consistent with each
1160	<	other. The GSF and TSF methods ensure the consistency and smooth
1161	<	truncation of the energy, force and torque at the cutoff radius, as a
1162	<	result show very good total energy conservation. But the TSF method
705	<	does not show good agreement in the absolute value of the
706	<	electrostatic energy, force and torque with the Ewald.  The GSF method
707	<	has mimicked Ewald’s force, energy and torque accurately and also
708	<	conserved energy. Therefore, the GSF method is the suitable method for
709	<	evaluating required force field in MD simulations. In addition, the
710	<	energy drift and fluctuation from the GSF method is much better than
711	<	Ewald’s method for finite-sized reciprocal space.
1154	>	In the first paper in this series, we generalized the
1155	>	charge-neutralized electrostatic energy originally developed by Wolf
1156	>	\textit{et al.}\cite{Wolf:1999dn} to multipole-multipole interactions
1157	>	up to quadrupolar order.  The SP method is essentially a
1158	>	multipole-capable version of the Wolf model.  The SP method for
1159	>	multipoles provides excellent agreement with Ewald-derived energies,
1160	>	forces and torques, and is suitable for Monte Carlo simulations,
1161	>	although the forces and torques retain discontinuities at the cutoff
1162	>	distance that prevents its use in molecular dynamics.
1163
1164	<	Note that the TSF, GSF, and SP models are $\mathcal{O}(N)$ methods
1165	<	that can be made extremely efficient using spline interpolations of
1166	<	the radial functions.  They require no Fourier transforms or $k$-space
1167	<	sums, and guarantee the smooth handling of energies, forces, and
1168	<	torques as multipoles cross the real-space cutoff boundary.
1164	>	We also developed two natural extensions of the damped shifted-force
1165	>	(DSF) model originally proposed by Fennel and
1166	>	Gezelter.\cite{Fennell:2006lq} The GSF and TSF approaches provide
1167	>	smooth truncation of energies, forces, and torques at the real-space
1168	>	cutoff, and both converge to DSF electrostatics for point-charge
1169	>	interactions.  The TSF model is based on a high-order truncated Taylor
1170	>	expansion which can be relatively perturbative inside the cutoff
1171	>	sphere.  The GSF model takes the gradient from an images of the
1172	>	interacting multipole that has been projected onto the cutoff sphere
1173	>	to derive shifted force and torque expressions, and is a significantly
1174	>	more gentle approach.
1175
1176	+	Of the two newly-developed shifted force models, the GSF method
1177	+	produced quantitative agreement with Ewald energy, force, and torques.
1178	+	It also performs well in conserving energy in MD simulations.  The
1179	+	Taylor-shifted (TSF) model provides smooth dynamics, but these take
1180	+	place on a potential energy surface that is significantly perturbed
1181	+	from Ewald-based electrostatics.
1182	+
1183	+	% The direct truncation of any electrostatic potential energy without
1184	+	% multipole neutralization creates large fluctuations in molecular
1185	+	% simulations.  This fluctuation in the energy is very large for the case
1186	+	% of crystal because of long range of multipole ordering (Refer paper
1187	+	% I).\cite{PaperI} This is also significant in the case of the liquid
1188	+	% because of the local multipole ordering in the molecules. If the net
1189	+	% multipole within cutoff radius neutralized within cutoff sphere by
1190	+	% placing image multiples on the surface of the sphere, this fluctuation
1191	+	% in the energy reduced significantly. Also, the multipole
1192	+	% neutralization in the generalized SP method showed very good agreement
1193	+	% with the Ewald as compared to direct truncation for the evaluation of
1194	+	% the $\triangle E$ between the configurations.  In MD simulations, the
1195	+	% energy conservation is very important. The conservation of the total
1196	+	% energy can be ensured by i) enforcing the smooth truncation of the
1197	+	% energy, force and torque in the cutoff radius and ii) making the
1198	+	% energy, force and torque consistent with each other. The GSF and TSF
1199	+	% methods ensure the consistency and smooth truncation of the energy,
1200	+	% force and torque at the cutoff radius, as a result show very good
1201	+	% total energy conservation. But the TSF method does not show good
1202	+	% agreement in the absolute value of the electrostatic energy, force and
1203	+	% torque with the Ewald.  The GSF method has mimicked Ewald’s force,
1204	+	% energy and torque accurately and also conserved energy.
1205	+
1206	+	The only cases we have found where the new GSF and SP real-space
1207	+	methods can be problematic are those which retain a bulk dipole moment
1208	+	at large distances (e.g. the $Z_1$ dipolar lattice).  In ferroelectric
1209	+	materials, uniform weighting of the orientational contributions can be
1210	+	important for converging the total energy.  In these cases, the
1211	+	damping function which causes the non-uniform weighting can be
1212	+	replaced by the bare electrostatic kernel, and the energies return to
1213	+	the expected converged values.
1214	+
1215	+	Based on the results of this work, the GSF method is a suitable and
1216	+	efficient replacement for the Ewald sum for evaluating electrostatic
1217	+	interactions in MD simulations.  Both methods retain excellent
1218	+	fidelity to the Ewald energies, forces and torques.  Additionally, the
1219	+	energy drift and fluctuations from the GSF electrostatics are better
1220	+	than a multipolar Ewald sum for finite-sized reciprocal spaces.
1221	+	Because they use real-space cutoffs with moderate cutoff radii, the
1222	+	GSF and SP models approach $\mathcal{O}(N)$ scaling as the system size
1223	+	increases.  Additionally, they can be made extremely efficient using
1224	+	spline interpolations of the radial functions.  They require no
1225	+	Fourier transforms or $k$-space sums, and guarantee the smooth
1226	+	handling of energies, forces, and torques as multipoles cross the
1227	+	real-space cutoff boundary.
1228	+
1229	+	\begin{acknowledgments}
1230	+	JDG acknowledges helpful discussions with Christopher
1231	+	Fennell. Support for this project was provided by the National
1232	+	Science Foundation under grant CHE-1362211. Computational time was
1233	+	provided by the Center for Research Computing (CRC) at the
1234	+	University of Notre Dame.
1235	+	\end{acknowledgments}
1236	+
1237		%\bibliographystyle{aip}
1238		\newpage
1239		\bibliography{references}

Diff Legend

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