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

Diff Legend

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