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} |
47 |
|
|
48 |
|
%\preprint{AIP/123-QED} |
49 |
|
|
50 |
< |
\title{Real space alternatives to the Ewald |
51 |
< |
Sum. II. Comparison of Methods} % 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 |
55 |
< |
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 |
59 |
< |
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, |
67 |
< |
% but any date may be explicitly specified |
64 |
> |
\date{\today} |
65 |
|
|
66 |
|
\begin{abstract} |
67 |
< |
We have tested the real-space shifted potential (SP), |
68 |
< |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
69 |
< |
for multipoles that were developed in the first paper in this series |
70 |
< |
against a reference method. The tests were carried out in a variety |
71 |
< |
of condensed-phase environments which were designed to test all |
72 |
< |
levels of the multipole-multipole interactions. Comparisons of the |
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 sum. We have |
76 |
< |
also investigated the energy conservation properties of the new |
77 |
< |
methods in molecular dynamics simulations using all of these |
78 |
< |
methods. 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 exhibits energy |
83 |
< |
conservation properties that make it an excellent choice for |
84 |
< |
efficiently computing electrostatic interactions in molecular |
88 |
< |
dynamics simulations. |
75 |
> |
relative to the more computationally expensive Ewald treatment. We |
76 |
> |
have also investigated the energy conservation properties of the new |
77 |
> |
methods in molecular dynamics simulations. The SP method shows |
78 |
> |
excellent agreement with configurational energy differences, forces, |
79 |
> |
and torques, and would be suitable for use in Monte Carlo |
80 |
> |
calculations. Of the two new shifted-force methods, the GSF |
81 |
> |
approach shows the best agreement with Ewald-derived energies, |
82 |
> |
forces, and torques and also exhibits energy conservation properties |
83 |
> |
that make it an excellent choice for efficient computation of |
84 |
> |
electrostatic interactions in molecular dynamics simulations. |
85 |
|
\end{abstract} |
86 |
|
|
87 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
88 |
|
% Classification Scheme. |
89 |
< |
\keywords{Electrostatics, Multipoles, Real-space} |
89 |
> |
%\keywords{Electrostatics, Multipoles, Real-space} |
90 |
|
|
91 |
|
\maketitle |
92 |
|
|
97 |
– |
|
93 |
|
\section{\label{sec:intro}Introduction} |
94 |
|
Computing the interactions between electrostatic sites is one of the |
95 |
< |
most expensive aspects of molecular simulations, which is why there |
96 |
< |
have been significant efforts to develop practical, efficient and |
97 |
< |
convergent methods for handling these interactions. Ewald's method is |
98 |
< |
perhaps the best known and most accurate method for evaluating |
99 |
< |
energies, forces, and torques in explicitly-periodic simulation |
100 |
< |
cells. In this approach, the conditionally convergent electrostatic |
101 |
< |
energy is converted into two absolutely convergent contributions, one |
102 |
< |
which is carried out in real space with a cutoff radius, and one in |
103 |
< |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
95 |
> |
most expensive aspects of molecular simulations. There have been |
96 |
> |
significant efforts to develop practical, efficient and convergent |
97 |
> |
methods for handling these interactions. Ewald's method is perhaps the |
98 |
> |
best known and most accurate method for evaluating energies, forces, |
99 |
> |
and torques in explicitly-periodic simulation cells. In this approach, |
100 |
> |
the conditionally convergent electrostatic energy is converted into |
101 |
> |
two absolutely convergent contributions, one which is carried out in |
102 |
> |
real space with a cutoff radius, and one in reciprocal |
103 |
> |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
104 |
|
|
105 |
|
When carried out as originally formulated, the reciprocal-space |
106 |
|
portion of the Ewald sum exhibits relatively poor computational |
107 |
< |
scaling, making it prohibitive for large systems. By utilizing |
108 |
< |
particle meshes and three dimensional fast Fourier transforms (FFT), |
109 |
< |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
110 |
< |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
111 |
< |
the computational cost from $O(N^2)$ down to $O(N \log |
112 |
< |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
107 |
> |
scaling, making it prohibitive for large systems. By utilizing a |
108 |
> |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
109 |
> |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
110 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
111 |
> |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
112 |
> |
\log |
113 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
114 |
|
|
115 |
< |
Because of the artificial periodicity required for the Ewald sum, the |
120 |
< |
method may require modification to compute interactions for |
115 |
> |
Because of the artificial periodicity required for the Ewald sum, |
116 |
|
interfacial molecular systems such as membranes and liquid-vapor |
117 |
< |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
118 |
< |
To simulate interfacial systems, Parry’s extension of the 3D Ewald sum |
119 |
< |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
120 |
< |
periodicity in the Ewald’s method can also be problematic for |
121 |
< |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
122 |
< |
methods that were developed to handle two-dimensional (2D) |
123 |
< |
electrostatic interactions in interfacial systems have not had similar |
124 |
< |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
125 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
117 |
> |
interfaces require modifications to the method. Parry's extension of |
118 |
> |
the three dimensional Ewald sum is appropriate for slab |
119 |
> |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
120 |
> |
developed to handle two-dimensional (2-D) electrostatic |
121 |
> |
interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
122 |
> |
These methods were originally quite computationally |
123 |
> |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
124 |
> |
efforts that reduced the computational cost of 2-D lattice summations, |
125 |
> |
bringing them more in line with the scaling for the full 3-D |
126 |
> |
treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The |
127 |
> |
inherent periodicity required by the Ewald method can also be |
128 |
> |
problematic in a number of protein/solvent and ionic solution |
129 |
> |
environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
130 |
|
|
131 |
|
\subsection{Real-space methods} |
132 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
133 |
|
method for calculating electrostatic interactions between point |
134 |
< |
charges. They argued that the effective Coulomb interaction in |
135 |
< |
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
136 |
< |
an ordered lattice (e.g. when computing the Madelung constant of an |
137 |
< |
ionic solid), the material can be considered as a set of ions |
138 |
< |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
139 |
< |
an effective distance dependence for the electrostatic interactions of |
140 |
< |
$r^{-5}$ (see figure \ref{fig:NaCl}. For this reason, careful |
141 |
< |
applications of Wolf's method are able to obtain accurate estimates of |
142 |
< |
Madelung constants using relatively short cutoff radii. Recently, |
143 |
< |
Fukuda used neutralization of the higher order moments for the |
144 |
< |
calculation of the electrostatic interaction of the point charges |
145 |
< |
system.\cite{Fukuda:2013sf} |
134 |
> |
charges. They argued that the effective Coulomb interaction in most |
135 |
> |
condensed phase systems is effectively short |
136 |
> |
ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when |
137 |
> |
computing the Madelung constant of an ionic solid), the material can |
138 |
> |
be considered as a set of ions interacting with neutral dipolar or |
139 |
> |
quadrupolar ``molecules'' giving an effective distance dependence for |
140 |
> |
the electrostatic interactions of $r^{-5}$ (see figure |
141 |
> |
\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple |
142 |
> |
cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the |
143 |
> |
electrostatic energy per ion converges more rapidly to the Madelung |
144 |
> |
energy than the dipolar approximation.\cite{Wolf92} To find the |
145 |
> |
correct Madelung constant, Lacman suggested that the NaCl structure |
146 |
> |
could be constructed in a way that the finite crystal terminates with |
147 |
> |
complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees |
148 |
> |
what is effectively a set of octupoles at large distances. These facts |
149 |
> |
suggest that the Madelung constants are relatively short ranged for |
150 |
> |
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
151 |
> |
application of Wolf's method can provide accurate estimates of |
152 |
> |
Madelung constants using relatively short cutoff radii. |
153 |
|
|
154 |
< |
\begin{figure}[h!] |
155 |
< |
\centering |
156 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
157 |
< |
\caption{Top: NaCl crystal showing how spherical truncation can |
158 |
< |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
159 |
< |
molecules interact with the central ion. Bottom: A dipolar |
160 |
< |
crystal exhibiting similar behavior and illustrating how the |
161 |
< |
effective dipole-octupole interactions can be disrupted by |
162 |
< |
spherical truncation.} |
163 |
< |
\label{fig:NaCl} |
164 |
< |
\end{figure} |
165 |
< |
|
166 |
< |
The direct truncation of interactions at a cutoff radius creates |
167 |
< |
truncation defects. Wolf \textit{et al.} further argued that |
168 |
< |
truncation errors are due to net charge remaining inside the cutoff |
169 |
< |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
164 |
< |
placing an image charge on the surface of the cutoff sphere for every |
165 |
< |
real charge inside the cutoff. These charges are present for the |
166 |
< |
evaluation of both the pair interaction energy and the force, although |
167 |
< |
the force expression maintained a discontinuity at the cutoff sphere. |
168 |
< |
In the original Wolf formulation, the total energy for the charge and |
169 |
< |
image were not equal to the integral of their force expression, and as |
170 |
< |
a result, the total energy would not be conserved in molecular |
171 |
< |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
172 |
< |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
173 |
< |
method with commensurate force and energy expressions that do not |
174 |
< |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
175 |
< |
methods were also proposed by Chen \textit{et |
154 |
> |
Direct truncation of interactions at a cutoff radius creates numerical |
155 |
> |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
156 |
> |
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
157 |
> |
neutralize this charge they proposed placing an image charge on the |
158 |
> |
surface of the cutoff sphere for every real charge inside the cutoff. |
159 |
> |
These charges are present for the evaluation of both the pair |
160 |
> |
interaction energy and the force, although the force expression |
161 |
> |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
162 |
> |
formulation, the total energy for the charge and image were not equal |
163 |
> |
to the integral of the force expression, and as a result, the total |
164 |
> |
energy would not be conserved in molecular dynamics (MD) |
165 |
> |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
166 |
> |
Gezelter later proposed shifted force variants of the Wolf method with |
167 |
> |
commensurate force and energy expressions that do not exhibit this |
168 |
> |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
169 |
> |
were also proposed by Chen \textit{et |
170 |
|
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
171 |
< |
and by Wu and Brooks.\cite{Wu:044107} |
171 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
172 |
> |
used additional neutralization of higher order moments for systems of |
173 |
> |
point charges.\cite{Fukuda:2013sf} |
174 |
|
|
175 |
< |
Considering the interaction of one central ion in an ionic crystal |
176 |
< |
with a portion of the crystal at some distance, the effective Columbic |
177 |
< |
potential is found to be decreasing as $r^{-5}$. If one views the |
178 |
< |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
179 |
< |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
180 |
< |
rapidly to the Madelung energy than the dipolar |
181 |
< |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
182 |
< |
Lacman suggested that the NaCl structure could be constructed in a way |
183 |
< |
that the finite crystal terminates with complete \ce{(NaCl)4} |
184 |
< |
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded |
185 |
< |
by opposite charges. Similarly for each pair of charges, there is an |
186 |
< |
opposite pair of charge adjacent to it. The central ion sees what is |
187 |
< |
effectively a set of octupoles at large distances. These facts suggest |
188 |
< |
that the Madelung constants are relatively short ranged for perfect |
189 |
< |
ionic crystals.\cite{Wolf:1999dn} |
175 |
> |
\begin{figure} |
176 |
> |
\centering |
177 |
> |
\includegraphics[width=\linewidth]{schematic.eps} |
178 |
> |
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
179 |
> |
charges (in the smaller grey circle), so interactions are |
180 |
> |
effectively charge-multipole at longer distances. With hard |
181 |
> |
cutoffs, motion of individual charges in and out of the cutoff |
182 |
> |
sphere can break the effective multipolar ordering. Bottom: |
183 |
> |
dipolar crystals and fluids have a similar effective |
184 |
> |
\textit{quadrupolar} ordering (in the smaller grey circles), and |
185 |
> |
orientational averaging helps to reduce the effective range of the |
186 |
> |
interactions in the fluid. Placement of reversed image multipoles |
187 |
> |
on the surface of the cutoff sphere recovers the effective |
188 |
> |
higher-order multipole behavior.} |
189 |
> |
\label{fig:schematic} |
190 |
> |
\end{figure} |
191 |
|
|
192 |
< |
One can make a similar argument for crystals of point multipoles. The |
193 |
< |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
194 |
< |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
195 |
< |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
196 |
< |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
197 |
< |
zero net dipole and retain contributions only from higher order |
198 |
< |
multipoles. The effective interaction between a dipole at the center |
192 |
> |
One can make a similar effective range argument for crystals of point |
193 |
> |
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
194 |
> |
constants for dipolar lattices utilizes 24 basis vectors that contain |
195 |
> |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
196 |
> |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
197 |
> |
moments, while the rest have zero net dipole and retain contributions |
198 |
> |
only from higher order multipoles. The lowest-energy crystalline |
199 |
> |
structures are built out of basis vectors that have only residual |
200 |
> |
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
201 |
> |
structures, the effective interaction between a dipole at the center |
202 |
|
of a crystal and a group of eight dipoles farther away is |
203 |
|
significantly shorter ranged than the $r^{-3}$ that one would expect |
204 |
|
for raw dipole-dipole interactions. Only in crystals which retain a |
208 |
|
unstable. |
209 |
|
|
210 |
|
In ionic crystals, real-space truncation can break the effective |
211 |
< |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
212 |
< |
swings in the electrostatic energy as the cutoff radius is increased |
213 |
< |
(or as individual ions move back and forth across the boundary). This |
214 |
< |
is why the image charges were necessary for the Wolf sum to exhibit |
215 |
< |
rapid convergence. Similarly, the real-space truncation of point |
216 |
< |
multipole interactions breaks higher order multipole arrangements, and |
217 |
< |
image multipoles are required for real-space treatments of |
218 |
< |
electrostatic energies. |
211 |
> |
multipolar arrangements (see Fig. \ref{fig:schematic}), causing |
212 |
> |
significant swings in the electrostatic energy as individual ions move |
213 |
> |
back and forth across the boundary. This is why the image charges are |
214 |
> |
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
215 |
> |
the real-space truncation of point multipole interactions breaks |
216 |
> |
higher order multipole arrangements, and image multipoles are required |
217 |
> |
for real-space treatments of electrostatic energies. |
218 |
|
|
219 |
+ |
The shorter effective range of electrostatic interactions is not |
220 |
+ |
limited to perfect crystals, but can also apply in disordered fluids. |
221 |
+ |
Even at elevated temperatures, there is local charge balance in an |
222 |
+ |
ionic liquid, where each positive ion has surroundings dominated by |
223 |
+ |
negaitve ions and vice versa. The reversed-charge images on the |
224 |
+ |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
225 |
+ |
the effective multipolar interactions as the charges traverse the |
226 |
+ |
cutoff boundary. |
227 |
+ |
|
228 |
+ |
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
229 |
+ |
significant orientational averaging that additionally reduces the |
230 |
+ |
effect of long-range multipolar interactions. The image multipoles |
231 |
+ |
that are introduced in the TSF, GSF, and SP methods mimic this effect |
232 |
+ |
and reduce the effective range of the multipolar interactions as |
233 |
+ |
interacting molecules traverse each other's cutoff boundaries. |
234 |
+ |
|
235 |
|
% Because of this reason, although the nature of electrostatic |
236 |
|
% interaction short ranged, the hard cutoff sphere creates very large |
237 |
|
% fluctuation in the electrostatic energy for the perfect crystal. In |
242 |
|
% to the non-neutralized value of the higher order moments within the |
243 |
|
% cutoff sphere. |
244 |
|
|
245 |
< |
The forces and torques acting on atomic sites are the fundamental |
246 |
< |
factors driving dynamics in molecular simulations. Fennell and |
247 |
< |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
248 |
< |
obtain consistent energies and forces on the atoms within the cutoff |
249 |
< |
sphere. Both the energy and the force go smoothly to zero as an atom |
250 |
< |
aproaches the cutoff radius. The comparisons of the accuracy these |
251 |
< |
quantities between the DSF kernel and SPME was surprisingly |
252 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
253 |
< |
calculating electrostatic interactions in molecular systems with |
254 |
< |
relatively uniform charge |
240 |
< |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
245 |
> |
Forces and torques acting on atomic sites are fundamental in driving |
246 |
> |
dynamics in molecular simulations, and the damped shifted force (DSF) |
247 |
> |
energy kernel provides consistent energies and forces on charged atoms |
248 |
> |
within the cutoff sphere. Both the energy and the force go smoothly to |
249 |
> |
zero as an atom aproaches the cutoff radius. The comparisons of the |
250 |
> |
accuracy these quantities between the DSF kernel and SPME was |
251 |
> |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
252 |
> |
has seen increasing use in molecular systems with relatively uniform |
253 |
> |
charge |
254 |
> |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
255 |
|
|
256 |
|
\subsection{The damping function} |
257 |
< |
The damping function used in our research has been discussed in detail |
258 |
< |
in the first paper of this series.\cite{PaperI} The radial kernel |
259 |
< |
$1/r$ for the interactions between point charges can be replaced by |
260 |
< |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
261 |
< |
accelerate the rate of convergence, where $\alpha$ is a damping |
262 |
< |
parameter with units of inverse distance. Altering the value of |
263 |
< |
$\alpha$ is equivalent to changing the width of Gaussian charge |
264 |
< |
distributions that replace each point charge -- Gaussian overlap |
265 |
< |
integrals yield complementary error functions when truncated at a |
266 |
< |
finite distance. |
257 |
> |
The damping function has been discussed in detail in the first paper |
258 |
> |
of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the |
259 |
> |
interactions between point charges can be replaced by the |
260 |
> |
complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate |
261 |
> |
convergence, where $\alpha$ is a damping parameter with units of |
262 |
> |
inverse distance. Altering the value of $\alpha$ is equivalent to |
263 |
> |
changing the width of Gaussian charge distributions that replace each |
264 |
> |
point charge, as Coulomb integrals with Gaussian charge distributions |
265 |
> |
produce complementary error functions when truncated at a finite |
266 |
> |
distance. |
267 |
|
|
268 |
< |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
269 |
< |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
270 |
< |
agreement with SPME for the interaction energies, forces and torques |
271 |
< |
for charge-charge interactions.\cite{Fennell:2006lq} |
268 |
> |
With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method |
269 |
> |
produced very good agreement with SPME for interaction energies, |
270 |
> |
forces and torques for charge-charge |
271 |
> |
interactions.\cite{Fennell:2006lq} |
272 |
|
|
273 |
|
\subsection{Point multipoles in molecular modeling} |
274 |
|
Coarse-graining approaches which treat entire molecular subsystems as |
275 |
|
a single rigid body are now widely used. A common feature of many |
276 |
|
coarse-graining approaches is simplification of the electrostatic |
277 |
|
interactions between bodies so that fewer site-site interactions are |
278 |
< |
required to compute configurational energies. Many coarse-grained |
279 |
< |
molecular structures would normally consist of equal positive and |
266 |
< |
negative charges, and rather than use multiple site-site interactions, |
267 |
< |
the interaction between higher order multipoles can also be used to |
268 |
< |
evaluate a single molecule-molecule |
269 |
< |
interaction.\cite{Ren06,Essex10,Essex11} |
278 |
> |
required to compute configurational |
279 |
> |
energies.\cite{Ren06,Essex10,Essex11} |
280 |
|
|
281 |
< |
Because electrons in a molecule are not localized at specific points, |
282 |
< |
the assignment of partial charges to atomic centers is a relatively |
283 |
< |
rough approximation. Atomic sites can also be assigned point |
284 |
< |
multipoles and polarizabilities to increase the accuracy of the |
285 |
< |
molecular model. Recently, water has been modeled with point |
286 |
< |
multipoles up to octupolar |
287 |
< |
order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
281 |
> |
Additionally, because electrons in a molecule are not localized at |
282 |
> |
specific points, the assignment of partial charges to atomic centers |
283 |
> |
is always an approximation. For increased accuracy, atomic sites can |
284 |
> |
also be assigned point multipoles and polarizabilities. Recently, |
285 |
> |
water has been modeled with point multipoles up to octupolar order |
286 |
> |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
287 |
> |
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
288 |
|
multipoles up to quadrupolar order have also been coupled with point |
289 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
290 |
< |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}. But |
291 |
< |
using point multipole with the real space truncation without |
292 |
< |
accounting for multipolar neutrality will create energy conservation |
293 |
< |
issues in molecular dynamics (MD) simulations. |
290 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
291 |
> |
truncating point multipoles without smoothing the forces and torques |
292 |
> |
can create energy conservation issues in molecular dynamics |
293 |
> |
simulations. |
294 |
|
|
295 |
|
In this paper we test a set of real-space methods that were developed |
296 |
|
for point multipolar interactions. These methods extend the damped |
297 |
|
shifted force (DSF) and Wolf methods originally developed for |
298 |
|
charge-charge interactions and generalize them for higher order |
299 |
< |
multipoles. The detailed mathematical development of these methods has |
300 |
< |
been presented in the first paper in this series, while this work |
301 |
< |
covers the testing the energies, forces, torques, and energy |
299 |
> |
multipoles. The detailed mathematical development of these methods |
300 |
> |
has been presented in the first paper in this series, while this work |
301 |
> |
covers the testing of energies, forces, torques, and energy |
302 |
|
conservation properties of the methods in realistic simulation |
303 |
|
environments. In all cases, the methods are compared with the |
304 |
< |
reference method, a full multipolar Ewald treatment. |
304 |
> |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
305 |
|
|
306 |
|
|
307 |
|
%\subsection{Conservation of total energy } |
325 |
|
expressed as the product of two multipole operators and a Coulombic |
326 |
|
kernel, |
327 |
|
\begin{equation} |
328 |
< |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
328 |
> |
U_{ab}(r)= M_{a} M_{b} \frac{1}{r} \label{kernel}. |
329 |
|
\end{equation} |
330 |
< |
where the multipole operator for site $\bf a$, |
331 |
< |
\begin{equation} |
332 |
< |
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} |
323 |
< |
+ Q_{{\bf a}\alpha\beta} |
324 |
< |
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
325 |
< |
\end{equation} |
326 |
< |
is expressed in terms of the point charge, $C_{\bf a}$, dipole, |
327 |
< |
$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for |
328 |
< |
object $\bf a$. Note that in this work, we use the primitive |
329 |
< |
quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k |
330 |
< |
r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site. |
330 |
> |
where the multipole operator for site $a$, $M_{a}$, is |
331 |
> |
expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object |
332 |
> |
$a$, etc. |
333 |
|
|
334 |
< |
Interactions between multipoles can be expressed as higher derivatives |
335 |
< |
of the bare Coulomb potential, so one way of ensuring that the forces |
336 |
< |
and torques vanish at the cutoff distance is to include a larger |
337 |
< |
number of terms in the truncated Taylor expansion, e.g., |
338 |
< |
% |
339 |
< |
\begin{equation} |
340 |
< |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
341 |
< |
\end{equation} |
342 |
< |
% |
343 |
< |
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
344 |
< |
Thus, for $f(r)=1/r$, we find |
345 |
< |
% |
346 |
< |
\begin{equation} |
347 |
< |
f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . |
348 |
< |
\end{equation} |
349 |
< |
This function is an approximate electrostatic potential that has |
350 |
< |
vanishing second derivatives at the cutoff radius, making it suitable |
351 |
< |
for shifting the forces and torques of charge-dipole interactions. |
334 |
> |
% Interactions between multipoles can be expressed as higher derivatives |
335 |
> |
% of the bare Coulomb potential, so one way of ensuring that the forces |
336 |
> |
% and torques vanish at the cutoff distance is to include a larger |
337 |
> |
% number of terms in the truncated Taylor expansion, e.g., |
338 |
> |
% % |
339 |
> |
% \begin{equation} |
340 |
> |
% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
341 |
> |
% \end{equation} |
342 |
> |
% % |
343 |
> |
% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
344 |
> |
% Thus, for $f(r)=1/r$, we find |
345 |
> |
% % |
346 |
> |
% \begin{equation} |
347 |
> |
% f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . |
348 |
> |
% \end{equation} |
349 |
> |
% This function is an approximate electrostatic potential that has |
350 |
> |
% vanishing second derivatives at the cutoff radius, making it suitable |
351 |
> |
% for shifting the forces and torques of charge-dipole interactions. |
352 |
|
|
353 |
< |
In general, the TSF potential for any multipole-multipole interaction |
354 |
< |
can be written |
353 |
> |
The TSF potential for any multipole-multipole interaction can be |
354 |
> |
written |
355 |
|
\begin{equation} |
356 |
|
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
357 |
|
\label{generic} |
358 |
|
\end{equation} |
359 |
< |
with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for |
360 |
< |
charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and |
361 |
< |
$n=4$ for quadrupole-quadrupole. To ensure smooth convergence of the |
362 |
< |
energy, force, and torques, the required number of terms from Taylor |
363 |
< |
series expansion in $f_n(r)$ must be performed for different |
364 |
< |
multipole-multipole interactions. |
365 |
< |
|
364 |
< |
To carry out the same procedure for a damped electrostatic kernel, we |
365 |
< |
replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
366 |
< |
Many of the derivatives of the damped kernel are well known from |
367 |
< |
Smith's early work on multipoles for the Ewald |
368 |
< |
summation.\cite{Smith82,Smith98} |
359 |
> |
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
360 |
> |
of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for |
361 |
> |
charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole |
362 |
> |
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for |
363 |
> |
quadrupole-quadrupole. To ensure smooth convergence of the energy, |
364 |
> |
force, and torques, a Taylor expansion with $n$ terms must be |
365 |
> |
performed at cutoff radius ($r_c$) to obtain $f_n(r)$. |
366 |
|
|
367 |
< |
Note that increasing the value of $n$ will add additional terms to the |
368 |
< |
electrostatic potential, e.g., $f_2(r)$ includes orders up to |
369 |
< |
$(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
370 |
< |
functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
371 |
< |
f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
375 |
< |
for computing multipole energies, forces, and torques, and smooth |
376 |
< |
cutoffs of these quantities can be guaranteed as long as the number of |
377 |
< |
terms in the Taylor series exceeds the derivative order required. |
367 |
> |
% To carry out the same procedure for a damped electrostatic kernel, we |
368 |
> |
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
369 |
> |
% Many of the derivatives of the damped kernel are well known from |
370 |
> |
% Smith's early work on multipoles for the Ewald |
371 |
> |
% summation.\cite{Smith82,Smith98} |
372 |
|
|
373 |
+ |
% Note that increasing the value of $n$ will add additional terms to the |
374 |
+ |
% electrostatic potential, e.g., $f_2(r)$ includes orders up to |
375 |
+ |
% $(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
376 |
+ |
% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
377 |
+ |
% f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
378 |
+ |
% for computing multipole energies, forces, and torques, and smooth |
379 |
+ |
% cutoffs of these quantities can be guaranteed as long as the number of |
380 |
+ |
% terms in the Taylor series exceeds the derivative order required. |
381 |
+ |
|
382 |
|
For multipole-multipole interactions, following this procedure results |
383 |
< |
in separate radial functions for each distinct orientational |
384 |
< |
contribution to the potential, and ensures that the forces and torques |
385 |
< |
from {\it each} of these contributions will vanish at the cutoff |
386 |
< |
radius. For example, the direct dipole dot product ($\mathbf{D}_{i} |
387 |
< |
\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance |
383 |
> |
in separate radial functions for each of the distinct orientational |
384 |
> |
contributions to the potential, and ensures that the forces and |
385 |
> |
torques from each of these contributions will vanish at the cutoff |
386 |
> |
radius. For example, the direct dipole dot product |
387 |
> |
($\mathbf{D}_{a} |
388 |
> |
\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance |
389 |
|
dot products: |
390 |
|
\begin{equation} |
391 |
< |
U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot |
392 |
< |
\mathbf{D}_{j} \right) \frac{g_2(r)}{r} |
393 |
< |
-\frac{1}{4\pi \epsilon_0} |
394 |
< |
\left( \mathbf{D}_{i} \cdot \hat{r} \right) |
395 |
< |
\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) - |
392 |
< |
\frac{g_2(r)}{r} \right) |
391 |
> |
U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
392 |
> |
\mathbf{D}_{a} \cdot |
393 |
> |
\mathbf{D}_{b} \right) v_{21}(r) + |
394 |
> |
\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right) |
395 |
> |
\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right] |
396 |
|
\end{equation} |
397 |
|
|
398 |
< |
The electrostatic forces and torques acting on the central multipole |
399 |
< |
site due to another site within cutoff sphere are derived from |
398 |
> |
For the Taylor shifted (TSF) method with the undamped kernel, |
399 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} + |
400 |
> |
\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4} |
401 |
> |
- \frac{6}{r r_c^2}$. In these functions, one can easily see the |
402 |
> |
connection to unmodified electrostatics as well as the smooth |
403 |
> |
transition to zero in both these functions as $r\rightarrow r_c$. The |
404 |
> |
electrostatic forces and torques acting on the central multipole due |
405 |
> |
to another site within the cutoff sphere are derived from |
406 |
|
Eq.~\ref{generic}, accounting for the appropriate number of |
407 |
|
derivatives. Complete energy, force, and torque expressions are |
408 |
|
presented in the first paper in this series (Reference |
410 |
|
|
411 |
|
\subsection{Gradient-shifted force (GSF)} |
412 |
|
|
413 |
< |
A second (and significantly simpler) method involves shifting the |
414 |
< |
gradient of the raw coulomb potential for each particular multipole |
413 |
> |
A second (and conceptually simpler) method involves shifting the |
414 |
> |
gradient of the raw Coulomb potential for each particular multipole |
415 |
|
order. For example, the raw dipole-dipole potential energy may be |
416 |
|
shifted smoothly by finding the gradient for two interacting dipoles |
417 |
|
which have been projected onto the surface of the cutoff sphere |
418 |
|
without changing their relative orientation, |
419 |
< |
\begin{displaymath} |
420 |
< |
U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(r_c) |
421 |
< |
- (r_{ij}-r_c) \hat{r}_{ij} \cdot |
422 |
< |
\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{r_c} |
423 |
< |
\end{displaymath} |
424 |
< |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$ |
425 |
< |
and $\mathbf{D}_{j}$, are retained at the cutoff distance (although |
426 |
< |
the signs are reversed for the dipole that has been projected onto the |
427 |
< |
cutoff sphere). In many ways, this simpler approach is closer in |
428 |
< |
spirit to the original shifted force method, in that it projects a |
429 |
< |
neutralizing multipole (and the resulting forces from this multipole) |
430 |
< |
onto a cutoff sphere. The resulting functional forms for the |
431 |
< |
potentials, forces, and torques turn out to be quite similar in form |
432 |
< |
to the Taylor-shifted approach, although the radial contributions are |
433 |
< |
significantly less perturbed by the Gradient-shifted approach than |
434 |
< |
they are in the Taylor-shifted method. |
419 |
> |
\begin{equation} |
420 |
> |
U_{D_{a}D_{b}}(r) = U_{D_{a}D_{b}}(r) - |
421 |
> |
U_{D_{a}D_{b}}(r_c) |
422 |
> |
- (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot |
423 |
> |
\nabla U_{D_{a}D_{b}}(r_c). |
424 |
> |
\end{equation} |
425 |
> |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance |
426 |
> |
(although the signs are reversed for the dipole that has been |
427 |
> |
projected onto the cutoff sphere). In many ways, this simpler |
428 |
> |
approach is closer in spirit to the original shifted force method, in |
429 |
> |
that it projects a neutralizing multipole (and the resulting forces |
430 |
> |
from this multipole) onto a cutoff sphere. The resulting functional |
431 |
> |
forms for the potentials, forces, and torques turn out to be quite |
432 |
> |
similar in form to the Taylor-shifted approach, although the radial |
433 |
> |
contributions are significantly less perturbed by the gradient-shifted |
434 |
> |
approach than they are in the Taylor-shifted method. |
435 |
|
|
436 |
+ |
For the gradient shifted (GSF) method with the undamped kernel, |
437 |
+ |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
438 |
+ |
$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. |
439 |
+ |
Again, these functions go smoothly to zero as $r\rightarrow r_c$, and |
440 |
+ |
because the Taylor expansion retains only one term, they are |
441 |
+ |
significantly less perturbed than the TSF functions. |
442 |
+ |
|
443 |
|
In general, the gradient shifted potential between a central multipole |
444 |
|
and any multipolar site inside the cutoff radius is given by, |
445 |
|
\begin{equation} |
446 |
< |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
447 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
448 |
< |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
446 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
447 |
> |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c) |
448 |
> |
\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
449 |
|
\label{generic2} |
450 |
|
\end{equation} |
451 |
|
where the sum describes a separate force-shifting that is applied to |
452 |
< |
each orientational contribution to the energy. |
452 |
> |
each orientational contribution to the energy. In this expression, |
453 |
> |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
454 |
> |
($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$ |
455 |
> |
represent the orientations the multipoles. |
456 |
|
|
457 |
|
The third term converges more rapidly than the first two terms as a |
458 |
|
function of radius, hence the contribution of the third term is very |
459 |
|
small for large cutoff radii. The force and torque derived from |
460 |
< |
equation \ref{generic2} are consistent with the energy expression and |
460 |
> |
Eq. \ref{generic2} are consistent with the energy expression and |
461 |
|
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
462 |
|
can be considered generalizations of the original DSF method for |
463 |
|
higher order multipole interactions. GSF and TSF are also identical up |
465 |
|
the energy, force and torque for higher order multipole-multipole |
466 |
|
interactions. Complete energy, force, and torque expressions for the |
467 |
|
GSF potential are presented in the first paper in this series |
468 |
< |
(Reference~\onlinecite{PaperI}) |
468 |
> |
(Reference~\onlinecite{PaperI}). |
469 |
|
|
470 |
|
|
471 |
|
\subsection{Shifted potential (SP) } |
478 |
|
interactions with the central multipole and the image. This |
479 |
|
effectively shifts the total potential to zero at the cutoff radius, |
480 |
|
\begin{equation} |
481 |
< |
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) |
481 |
> |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
482 |
> |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
483 |
|
\label{eq:SP} |
484 |
|
\end{equation} |
485 |
|
where the sum describes separate potential shifting that is done for |
491 |
|
multipoles that reorient after leaving the cutoff sphere can re-enter |
492 |
|
the cutoff sphere without perturbing the total energy. |
493 |
|
|
494 |
< |
The potential energy between a central multipole and other multipolar |
495 |
< |
sites then goes smoothly to zero as $r \rightarrow r_c$. However, the |
496 |
< |
force and torque obtained from the shifted potential (SP) are |
497 |
< |
discontinuous at $r_c$. Therefore, MD simulations will still |
498 |
< |
experience energy drift while operating under the SP potential, but it |
499 |
< |
may be suitable for Monte Carlo approaches where the configurational |
500 |
< |
energy differences are the primary quantity of interest. |
494 |
> |
For the shifted potential (SP) method with the undamped kernel, |
495 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) = |
496 |
> |
\frac{3}{r^3} - \frac{3}{r_c^3}$. The potential energy between a |
497 |
> |
central multipole and other multipolar sites goes smoothly to zero as |
498 |
> |
$r \rightarrow r_c$. However, the force and torque obtained from the |
499 |
> |
shifted potential (SP) are discontinuous at $r_c$. MD simulations |
500 |
> |
will still experience energy drift while operating under the SP |
501 |
> |
potential, but it may be suitable for Monte Carlo approaches where the |
502 |
> |
configurational energy differences are the primary quantity of |
503 |
> |
interest. |
504 |
|
|
505 |
< |
\subsection{The Self term} |
505 |
> |
\subsection{The Self Term} |
506 |
|
In the TSF, GSF, and SP methods, a self-interaction is retained for |
507 |
|
the central multipole interacting with its own image on the surface of |
508 |
|
the cutoff sphere. This self interaction is nearly identical with the |
524 |
|
used the multipolar Ewald sum as a reference method for comparing |
525 |
|
energies, forces, and torques for molecular models that mimic |
526 |
|
disordered and ordered condensed-phase systems. The parameters used |
527 |
< |
in the test-cases are given in table~\ref{tab:pars}. |
527 |
> |
in the test cases are given in table~\ref{tab:pars}. |
528 |
|
|
529 |
|
\begin{table} |
530 |
|
\label{tab:pars} |
542 |
|
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
543 |
|
\AA\textsuperscript{2})} \\ \hline |
544 |
|
Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ |
545 |
< |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & 10,000 & 17.6 &17.6 & 0 \\ |
545 |
> |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ |
546 |
|
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
547 |
< |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & 10,000 & 17.6&17.6&0 \\ |
547 |
> |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ |
548 |
|
SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
549 |
|
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
550 |
|
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
569 |
|
electrostatic energy, as well as the electrostatic contributions to |
570 |
|
the force and torque on each molecule. These quantities have been |
571 |
|
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
572 |
< |
and have been compared with the values obtaine from the multipolar |
573 |
< |
Ewald sum. In Mote Carlo (MC) simulations, the energy differences |
572 |
> |
and have been compared with the values obtained from the multipolar |
573 |
> |
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
574 |
|
between two configurations is the primary quantity that governs how |
575 |
< |
the simulation proceeds. These differences are the most imporant |
575 |
> |
the simulation proceeds. These differences are the most important |
576 |
|
indicators of the reliability of a method even if the absolute |
577 |
|
energies are not exact. For each of the multipolar systems listed |
578 |
|
above, we have compared the change in electrostatic potential energy |
584 |
|
\subsection{Implementation} |
585 |
|
The real-space methods developed in the first paper in this series |
586 |
|
have been implemented in our group's open source molecular simulation |
587 |
< |
program, OpenMD,\cite{openmd} which was used for all calculations in |
587 |
> |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
588 |
|
this work. The complementary error function can be a relatively slow |
589 |
|
function on some processors, so all of the radial functions are |
590 |
|
precomputed on a fine grid and are spline-interpolated to provide |
617 |
|
recomputed at each time step. |
618 |
|
|
619 |
|
\subsection{Model systems} |
620 |
< |
To sample independent configurations of multipolar crystals, a body |
621 |
< |
centered cubic (bcc) crystal which is a minimum energy structure for |
622 |
< |
point dipoles was generated using 3,456 molecules. The multipoles |
623 |
< |
were translationally locked in their respective crystal sites for |
624 |
< |
equilibration at a relatively low temperature (50K), so that dipoles |
625 |
< |
or quadrupoles could freely explore all accessible orientations. The |
626 |
< |
translational constraints were removed, and the crystals were |
627 |
< |
simulated for 10 ps in the microcanonical (NVE) ensemble with an |
628 |
< |
average temperature of 50 K. Configurations were sampled at equal |
629 |
< |
time intervals for the comparison of the configurational energy |
630 |
< |
differences. The crystals were not simulated close to the melting |
631 |
< |
points in order to avoid translational deformation away of the ideal |
632 |
< |
lattice geometry. |
620 |
> |
To sample independent configurations of the multipolar crystals, body |
621 |
> |
centered cubic (bcc) crystals, which exhibit the minimum energy |
622 |
> |
structures for point dipoles, were generated using 3,456 molecules. |
623 |
> |
The multipoles were translationally locked in their respective crystal |
624 |
> |
sites for equilibration at a relatively low temperature (50K) so that |
625 |
> |
dipoles or quadrupoles could freely explore all accessible |
626 |
> |
orientations. The translational constraints were then removed, the |
627 |
> |
systems were re-equilibrated, and the crystals were simulated for an |
628 |
> |
additional 10 ps in the microcanonical (NVE) ensemble with an average |
629 |
> |
temperature of 50 K. The balance between moments of inertia and |
630 |
> |
particle mass were chosen to allow orientational sampling without |
631 |
> |
significant translational motion. Configurations were sampled at |
632 |
> |
equal time intervals in order to compare configurational energy |
633 |
> |
differences. The crystals were simulated far from the melting point |
634 |
> |
in order to avoid translational deformation away of the ideal lattice |
635 |
> |
geometry. |
636 |
|
|
637 |
< |
For dipolar, quadrupolar, and mixed-multipole liquid simulations, each |
638 |
< |
system was created with 2048 molecules oriented randomly. These were |
639 |
< |
|
640 |
< |
system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K |
641 |
< |
temperature after equilibration. We collected 250 different |
642 |
< |
configurations in equal interval of time. For the ions mixed liquid |
643 |
< |
system, we converted 48 different molecules into 24 \ce{Na+} and 24 |
644 |
< |
\ce{Cl-} ions and equilibrated. After equilibration, the system was run |
645 |
< |
at the same environment for 1ns and 250 configurations were |
646 |
< |
collected. While comparing energies, forces, and torques with Ewald |
647 |
< |
method, Lennard-Jones potentials were turned off and purely |
648 |
< |
electrostatic interaction had been compared. |
637 |
> |
For dipolar, quadrupolar, and mixed-multipole \textit{liquid} |
638 |
> |
simulations, each system was created with 2,048 randomly-oriented |
639 |
> |
molecules. These were equilibrated at a temperature of 300K for 1 ns. |
640 |
> |
Each system was then simulated for 1 ns in the microcanonical (NVE) |
641 |
> |
ensemble. We collected 250 different configurations at equal time |
642 |
> |
intervals. For the liquid system that included ionic species, we |
643 |
> |
converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24 |
644 |
> |
\ce{Cl-} ions and re-equilibrated. After equilibration, the system was |
645 |
> |
run under the same conditions for 1 ns. A total of 250 configurations |
646 |
> |
were collected. In the following comparisons of energies, forces, and |
647 |
> |
torques, the Lennard-Jones potentials were turned off and only the |
648 |
> |
purely electrostatic quantities were compared with the same values |
649 |
> |
obtained via the Ewald sum. |
650 |
|
|
651 |
|
\subsection{Accuracy of Energy Differences, Forces and Torques} |
652 |
|
The pairwise summation techniques (outlined above) were evaluated for |
660 |
|
should be identical for all methods. |
661 |
|
|
662 |
|
Since none of the real-space methods provide exact energy differences, |
663 |
< |
we used least square regressions analysiss for the six different |
663 |
> |
we used least square regressions analysis for the six different |
664 |
|
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
665 |
|
with the multipolar Ewald reference method. Unitary results for both |
666 |
|
the correlation (slope) and correlation coefficient for these |
671 |
|
configurations and 250 configurations were recorded for comparison. |
672 |
|
Each system provided 31,125 energy differences for a total of 186,750 |
673 |
|
data points. Similarly, the magnitudes of the forces and torques have |
674 |
< |
also been compared by using least squares regression analyses. In the |
674 |
> |
also been compared using least squares regression analysis. In the |
675 |
|
forces and torques comparison, the magnitudes of the forces acting in |
676 |
|
each molecule for each configuration were evaluated. For example, our |
677 |
|
dipolar liquid simulation contains 2048 molecules and there are 250 |
791 |
|
|
792 |
|
\begin{figure} |
793 |
|
\centering |
794 |
< |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
794 |
> |
\includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
795 |
|
\caption{Statistical analysis of the quality of configurational |
796 |
|
energy differences for the real-space electrostatic methods |
797 |
|
compared with the reference Ewald sum. Results with a value equal |
856 |
|
molecules inside each other's cutoff spheres in order to correct the |
857 |
|
energy conservation issues, and this perturbation is evident in the |
858 |
|
statistics accumulated for the molecular forces. The GSF |
859 |
< |
perturbations are minimal, particularly for moderate damping and and |
859 |
> |
perturbations are minimal, particularly for moderate damping and |
860 |
|
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
861 |
|
reasonable agreement in the correlation coefficient but again the |
862 |
|
systematic error in the forces is concerning if replication of Ewald |
864 |
|
|
865 |
|
\begin{figure} |
866 |
|
\centering |
867 |
< |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
867 |
> |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} |
868 |
|
\caption{Statistical analysis of the quality of the force vector |
869 |
|
magnitudes for the real-space electrostatic methods compared with |
870 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
878 |
|
|
879 |
|
\begin{figure} |
880 |
|
\centering |
881 |
< |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
881 |
> |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps} |
882 |
|
\caption{Statistical analysis of the quality of the torque vector |
883 |
|
magnitudes for the real-space electrostatic methods compared with |
884 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
913 |
|
these quantities. Force and torque vectors for all six systems were |
914 |
|
analyzed using Fisher statistics, and the quality of the vector |
915 |
|
directionality is shown in terms of circular variance |
916 |
< |
($\mathrm{Var}(\theta$) in figure |
916 |
> |
($\mathrm{Var}(\theta)$) in figure |
917 |
|
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
918 |
|
from the new real-space methods exhibit nearly-ideal Fisher probability |
919 |
|
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
936 |
|
|
937 |
|
\begin{figure} |
938 |
|
\centering |
939 |
< |
\includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} |
939 |
> |
\includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps} |
940 |
|
\caption{The circular variance of the direction of the force and |
941 |
|
torque vectors obtained from the real-space methods around the |
942 |
|
reference Ewald vectors. A variance equal to 0 (dashed line) |
956 |
|
in this series and provides the most comprehensive test of the new |
957 |
|
methods. A liquid-phase system was created with 2000 water molecules |
958 |
|
and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
959 |
< |
temperature of 300K. After equilibration, this liquid-phase system |
960 |
< |
was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with |
961 |
< |
a cutoff radius of 12\AA. The value of the damping coefficient was |
962 |
< |
also varied from the undamped case ($\alpha = 0$) to a heavily damped |
963 |
< |
case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods. A |
964 |
< |
sample was also run using the multipolar Ewald sum with the same |
965 |
< |
real-space cutoff. |
959 |
> |
temperature of 300K. After equilibration in the canonical (NVT) |
960 |
> |
ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system |
961 |
> |
was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald, |
962 |
> |
Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA. The |
963 |
> |
value of the damping coefficient was also varied from the undamped |
964 |
> |
case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$ |
965 |
> |
\AA$^{-1}$) for all of the real space methods. A sample was also run |
966 |
> |
using the multipolar Ewald sum with the same real-space cutoff. |
967 |
|
|
968 |
|
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
969 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
970 |
|
fluctuations around this drift $\delta E_0$. Both of the |
971 |
|
shifted-force methods (GSF and TSF) provide excellent energy |
972 |
< |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
972 |
> |
conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), |
973 |
|
while the hard cutoff is essentially unusable for molecular dynamics. |
974 |
|
SP provides some benefit over the hard cutoff because the energetic |
975 |
|
jumps that happen as particles leave and enter the cutoff sphere are |
984 |
|
|
985 |
|
\begin{figure} |
986 |
|
\centering |
987 |
< |
\includegraphics[width=\textwidth]{newDrift.pdf} |
987 |
> |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
988 |
|
\label{fig:energyDrift} |
989 |
|
\caption{Analysis of the energy conservation of the real-space |
990 |
< |
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
991 |
< |
energy over time and $\delta \mathrm{E}_0$ is the standard deviation |
992 |
< |
of energy fluctuations around this drift. All simulations were of a |
993 |
< |
2000-molecule simulation of SSDQ water with 48 ionic charges at 300 |
994 |
< |
K starting from the same initial configuration. All runs utilized |
995 |
< |
the same real-space cutoff, $r_c = 12$\AA.} |
990 |
> |
methods. $\delta \mathrm{E}_1$ is the linear drift in energy over |
991 |
> |
time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is |
992 |
> |
the standard deviation of energy fluctuations around this drift (in |
993 |
> |
kcal / mol / particle). Points that appear below the dashed grey |
994 |
> |
(Ewald) lines exhibit better energy conservation than commonly-used |
995 |
> |
parameters for Ewald-based electrostatics. All simulations were of |
996 |
> |
a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
997 |
> |
300 K starting from the same initial configuration. All runs |
998 |
> |
utilized the same real-space cutoff, $r_c = 12$\AA.} |
999 |
|
\end{figure} |
1000 |
+ |
|
1001 |
+ |
\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}} |
1002 |
+ |
The most important test of the modified interaction potentials is the |
1003 |
+ |
fidelity with which they can reproduce structural features and |
1004 |
+ |
dynamical properties in a liquid. One commonly-utilized measure of |
1005 |
+ |
structural ordering is the pair distribution function, $g(r)$, which |
1006 |
+ |
measures local density deviations in relation to the bulk density. In |
1007 |
+ |
the electrostatic approaches studied here, the short-range repulsion |
1008 |
+ |
from the Lennard-Jones potential is identical for the various |
1009 |
+ |
electrostatic methods, and since short range repulsion determines much |
1010 |
+ |
of the local liquid ordering, one would not expect to see many |
1011 |
+ |
differences in $g(r)$. Indeed, the pair distributions are essentially |
1012 |
+ |
identical for all of the electrostatic methods studied (for each of |
1013 |
+ |
the different systems under investigation). An example of this |
1014 |
+ |
agreement for the SSDQ water/ion system is shown in |
1015 |
+ |
Fig. \ref{fig:gofr}. |
1016 |
+ |
|
1017 |
+ |
\begin{figure} |
1018 |
+ |
\centering |
1019 |
+ |
\includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
1020 |
+ |
\label{fig:gofr} |
1021 |
+ |
\caption{The pair distribution functions, $g(r)$, for the SSDQ |
1022 |
+ |
water/ion system obtained using the different real-space methods are |
1023 |
+ |
essentially identical with the result from the Ewald |
1024 |
+ |
treatment.} |
1025 |
+ |
\end{figure} |
1026 |
+ |
|
1027 |
+ |
There is a very slight overstructuring of the first solvation shell |
1028 |
+ |
when using when using TSF at lower values of the damping coefficient |
1029 |
+ |
($\alpha \le 0.1$) or when using undamped GSF. With moderate damping, |
1030 |
+ |
GSF and SP produce pair distributions that are identical (within |
1031 |
+ |
numerical noise) to their Ewald counterparts. |
1032 |
+ |
|
1033 |
+ |
A structural property that is a more demanding test of modified |
1034 |
+ |
electrostatics is the mean value of the electrostatic energy $\langle |
1035 |
+ |
U_\mathrm{elect} \rangle / N$ which is obtained by sampling the |
1036 |
+ |
liquid-state configurations experienced by a liquid evolving entirely |
1037 |
+ |
under the influence of each of the methods. In table \ref{tab:Props} |
1038 |
+ |
we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with |
1039 |
+ |
the damping parameter, $\alpha$, for each of the methods. |
1040 |
+ |
|
1041 |
+ |
As in the crystals studied in the first paper, damping is important |
1042 |
+ |
for converging the mean electrostatic energy values, particularly for |
1043 |
+ |
the two shifted force methods (GSF and TSF). A value of $\alpha |
1044 |
+ |
\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF |
1045 |
+ |
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
1046 |
+ |
dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA). |
1047 |
+ |
Overdamping the real-space electrostatic methods occurs with $\alpha > |
1048 |
+ |
0.4$, causing the estimate of the energy to drop below the Ewald |
1049 |
+ |
results. |
1050 |
+ |
|
1051 |
+ |
These ``optimal'' values of the damping coefficient are slightly |
1052 |
+ |
larger than what were observed for DSF electrostatics for purely |
1053 |
+ |
point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for |
1054 |
+ |
$r_c = 12$\AA appears to be an excellent compromise for mixed charge |
1055 |
+ |
multipole systems. |
1056 |
|
|
1057 |
+ |
To test the fidelity of the electrostatic methods at reproducing |
1058 |
+ |
dynamics in a multipolar liquid, it is also useful to look at |
1059 |
+ |
transport properties, particularly the diffusion constant, |
1060 |
+ |
\begin{equation} |
1061 |
+ |
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
1062 |
+ |
\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle |
1063 |
+ |
\label{eq:diff} |
1064 |
+ |
\end{equation} |
1065 |
+ |
which measures long-time behavior and is sensitive to the forces on |
1066 |
+ |
the multipoles. For the soft dipolar fluid and the SSDQ liquid |
1067 |
+ |
systems, the self-diffusion constants (D) were calculated from linear |
1068 |
+ |
fits to the long-time portion of the mean square displacement, |
1069 |
+ |
$\langle r^{2}(t) \rangle$.\cite{Allen87} |
1070 |
|
|
1071 |
+ |
In addition to translational diffusion, orientational relaxation times |
1072 |
+ |
were calculated for comparisons with the Ewald simulations and with |
1073 |
+ |
experiments. These values were determined from the same 1~ns |
1074 |
+ |
microcanonical trajectories used for translational diffusion by |
1075 |
+ |
calculating the orientational time correlation function, |
1076 |
+ |
\begin{equation} |
1077 |
+ |
C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t) |
1078 |
+ |
\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle, |
1079 |
+ |
\label{eq:OrientCorr} |
1080 |
+ |
\end{equation} |
1081 |
+ |
where $P_l$ is the Legendre polynomial of order $l$ and |
1082 |
+ |
$\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis |
1083 |
+ |
$\gamma$ on a molecule.. Th body-fixed reference frame used for our |
1084 |
+ |
models has the $z$-axis running along the dipoles, and for the SSDQ |
1085 |
+ |
water model, the $y$-axis connects the two implied hydrogen atom |
1086 |
+ |
positions. From the orientation autocorrelation functions, we can |
1087 |
+ |
obtain time constants for rotational relaxation either by fitting an |
1088 |
+ |
exponential function or by integrating the entire correlation |
1089 |
+ |
function. In a good water model, these decay times would be |
1090 |
+ |
comparable to water orientational relaxation times from nuclear |
1091 |
+ |
magnetic resonance (NMR). The relaxation constant obtained from |
1092 |
+ |
$C_2^y(t)$ is normally of experimental interest because it describes |
1093 |
+ |
the relaxation of the principle axis connecting the hydrogen |
1094 |
+ |
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
1095 |
+ |
of the dipole-dipole relaxation from a proton NMR signal and should |
1096 |
+ |
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
1097 |
+ |
|
1098 |
+ |
Results for the diffusion constants and orientational relaxation times |
1099 |
+ |
are shown in figure \ref{tab:Props}. From this data, it is apparent |
1100 |
+ |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
1101 |
+ |
reproduced with reasonable fidelity by the GSF method. |
1102 |
+ |
|
1103 |
+ |
The $\tau_2$ results in \ref{tab:Props} show a much greater difference |
1104 |
+ |
between the real-space and the Ewald results. |
1105 |
+ |
|
1106 |
+ |
\begin{table} |
1107 |
+ |
\label{tab:Props} |
1108 |
+ |
\caption{Comparison of the structural and dynamic properties for the |
1109 |
+ |
soft dipolar liquid test for all of the real-space methods.} |
1110 |
+ |
\begin{tabular}{l|c|cccc|cccc|cccc} |
1111 |
+ |
& Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\ |
1112 |
+ |
$\alpha$ (\AA$^{-1}$) & & |
1113 |
+ |
0.0 & 0.1 & 0.2 & 0.3 & |
1114 |
+ |
0.0 & 0.1 & 0.2 & 0.3 & |
1115 |
+ |
0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14} |
1116 |
+ |
|
1117 |
+ |
$\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\ |
1118 |
+ |
D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)& |
1119 |
+ |
470.2(6) & |
1120 |
+ |
416.6(5) & |
1121 |
+ |
379.6(5) & |
1122 |
+ |
438.6(5) & |
1123 |
+ |
476.0(6) & |
1124 |
+ |
412.8(5) & |
1125 |
+ |
421.1(5) & |
1126 |
+ |
400.5(5) & |
1127 |
+ |
437.5(6) & |
1128 |
+ |
434.6(5) & |
1129 |
+ |
411.4(5) & |
1130 |
+ |
545.3(7) & |
1131 |
+ |
459.6(6) \\ |
1132 |
+ |
$\tau_2$ (fs) & |
1133 |
+ |
1.136 & |
1134 |
+ |
1.041 & |
1135 |
+ |
1.064 & |
1136 |
+ |
1.109 & |
1137 |
+ |
1.211 & |
1138 |
+ |
1.119 & |
1139 |
+ |
1.039 & |
1140 |
+ |
1.058 & |
1141 |
+ |
1.21 & |
1142 |
+ |
1.15 & |
1143 |
+ |
1.172 & |
1144 |
+ |
1.153 & |
1145 |
+ |
1.125 \\ |
1146 |
+ |
\end{tabular} |
1147 |
+ |
\end{table} |
1148 |
+ |
|
1149 |
+ |
|
1150 |
|
\section{CONCLUSION} |
1151 |
|
In the first paper in this series, we generalized the |
1152 |
|
charge-neutralized electrostatic energy originally developed by Wolf |
1223 |
|
handling of energies, forces, and torques as multipoles cross the |
1224 |
|
real-space cutoff boundary. |
1225 |
|
|
1226 |
+ |
\begin{acknowledgments} |
1227 |
+ |
JDG acknowledges helpful discussions with Christopher |
1228 |
+ |
Fennell. Support for this project was provided by the National |
1229 |
+ |
Science Foundation under grant CHE-1362211. Computational time was |
1230 |
+ |
provided by the Center for Research Computing (CRC) at the |
1231 |
+ |
University of Notre Dame. |
1232 |
+ |
\end{acknowledgments} |
1233 |
+ |
|
1234 |
|
%\bibliographystyle{aip} |
1235 |
|
\newpage |
1236 |
|
\bibliography{references} |