| 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 multipole interactions that were developed in the first paper in |
| 70 |
< |
this series, using the multipolar Ewald sum as a reference |
| 71 |
< |
method. The tests were carried out in a variety of condensed-phase |
| 72 |
< |
environments which were designed to test all levels of the |
| 73 |
< |
multipole-multipole interactions. Comparisons of the energy |
| 74 |
< |
differences between configurations, molecular forces, and torques |
| 75 |
< |
were used to analyze how well the real-space models perform relative |
| 76 |
< |
to the more computationally expensive Ewald treatment. We have also |
| 77 |
< |
investigated the energy conservation properties of the new methods |
| 78 |
< |
in molecular dynamics simulations. The SP method shows excellent |
| 79 |
< |
agreement with configurational energy differences, forces, and |
| 80 |
< |
torques, and would be suitable for use in Monte Carlo calculations. |
| 81 |
< |
Of the two new shifted-force methods, the GSF approach shows the |
| 82 |
< |
best agreement with Ewald-derived energies, forces, and torques and |
| 83 |
< |
exhibits energy conservation properties that make it an excellent |
| 84 |
< |
choice for efficient computation of electrostatic interactions in |
| 85 |
< |
molecular dynamics simulations. |
| 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 |
| 93 |
|
|
| 94 |
|
\maketitle |
| 95 |
|
|
| 97 |
– |
|
| 96 |
|
\section{\label{sec:intro}Introduction} |
| 97 |
|
Computing the interactions between electrostatic sites is one of the |
| 98 |
|
most expensive aspects of molecular simulations. There have been |
| 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. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
| 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, |
| 119 |
|
interfacial molecular systems such as membranes and liquid-vapor |
| 120 |
< |
interfaces require modifications to the |
| 121 |
< |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 122 |
< |
Parry's extension of the three dimensional Ewald sum is appropriate |
| 123 |
< |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
| 124 |
< |
were developed to handle two-dimensional (2D) electrostatic |
| 125 |
< |
interactions in interfacial systems have not seen similar |
| 126 |
< |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 127 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
| 128 |
< |
with system size. The inherent periodicity in the Ewald’s method can |
| 129 |
< |
also be problematic for interfacial molecular |
| 130 |
< |
systems.\cite{Fennell:2006lq} |
| 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)$ |
| 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 are able to obtain accurate estimates of |
| 154 |
> |
application of Wolf's method can provide accurate estimates of |
| 155 |
|
Madelung constants using relatively short cutoff radii. |
| 156 |
|
|
| 157 |
|
Direct truncation of interactions at a cutoff radius creates numerical |
| 158 |
< |
errors. Wolf \textit{et al.} argued that truncation errors are due |
| 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 |
< |
maintained a discontinuity at the cutoff sphere. In the original Wolf |
| 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 their force expression, and as a result, the total |
| 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{Fennell:2006lq} Related real-space methods were also |
| 172 |
< |
proposed by Chen \textit{et |
| 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 used |
| 175 |
< |
neutralization of the higher order moments for the calculation of the |
| 176 |
< |
electrostatic interaction of the point charge |
| 176 |
< |
systems.\cite{Fukuda:2013sf} |
| 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=\linewidth]{schematic.pdf} |
| 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 |
| 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.} |
| 192 |
< |
\label{fig:schematic} |
| 191 |
> |
higher-order multipole behavior. \label{fig:schematic}} |
| 192 |
|
\end{figure} |
| 193 |
|
|
| 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. Only three of these |
| 198 |
< |
basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
| 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 |
| 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 |
| 220 |
|
|
| 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, on average, local charge |
| 224 |
< |
balance in an ionic liquid, where each positive ion has surroundings |
| 225 |
< |
dominated by negaitve ions and vice versa. The reversed-charge images |
| 226 |
< |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
| 227 |
< |
retain the effective multipolar interactions as the charges traverse |
| 228 |
< |
the cutoff boundary. |
| 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 |
|
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
| 231 |
|
significant orientational averaging that additionally reduces the |
| 234 |
|
and reduce the effective range of the multipolar interactions as |
| 235 |
|
interacting molecules traverse each other's cutoff boundaries. |
| 236 |
|
|
| 237 |
< |
% Because of this reason, although the nature of electrostatic |
| 238 |
< |
% interaction short ranged, the hard cutoff sphere creates very large |
| 239 |
< |
% fluctuation in the electrostatic energy for the perfect crystal. In |
| 240 |
< |
% addition, the charge neutralized potential proposed by Wolf et |
| 241 |
< |
% al. converged to correct Madelung constant but still holds oscillation |
| 242 |
< |
% in the energy about correct Madelung energy.\cite{Wolf:1999dn}. This |
| 243 |
< |
% oscillation in the energy around its fully converged value can be due |
| 244 |
< |
% to the non-neutralized value of the higher order moments within the |
| 245 |
< |
% cutoff sphere. |
| 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 |
– |
The forces and torques acting on atomic sites are the fundamental |
| 249 |
– |
factors driving dynamics in molecular simulations. Fennell and |
| 250 |
– |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 251 |
– |
obtain consistent energies and forces on the atoms within the cutoff |
| 252 |
– |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 253 |
– |
aproaches the cutoff radius. The comparisons of the accuracy these |
| 254 |
– |
quantities between the DSF kernel and SPME was surprisingly |
| 255 |
– |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
| 256 |
– |
calculating electrostatic interactions in molecular systems with |
| 257 |
– |
relatively uniform charge |
| 258 |
– |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 259 |
– |
|
| 248 |
|
\subsection{The damping function} |
| 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 |
| 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 always an |
| 275 |
< |
approximation. Atomic sites can also be assigned point multipoles and |
| 276 |
< |
polarizabilities to increase the accuracy of the molecular model. |
| 277 |
< |
Recently, water has been modeled with point multipoles up to octupolar |
| 278 |
< |
order using the soft sticky dipole-quadrupole-octupole (SSDQO) |
| 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} However, |
| 283 |
|
truncating point multipoles without smoothing the forces and torques |
| 284 |
< |
will create energy conservation issues in molecular dynamics simulations. |
| 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 |
|
|
| 310 |
– |
%\subsection{Conservation of total energy } |
| 311 |
– |
%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. |
| 312 |
– |
|
| 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) |
| 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$, $\hat{M}_{\bf a}$, is |
| 320 |
< |
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
| 321 |
< |
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
| 336 |
< |
$\bf a$, etc. |
| 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 |
|
|
| 338 |
– |
% Interactions between multipoles can be expressed as higher derivatives |
| 339 |
– |
% of the bare Coulomb potential, so one way of ensuring that the forces |
| 340 |
– |
% and torques vanish at the cutoff distance is to include a larger |
| 341 |
– |
% number of terms in the truncated Taylor expansion, e.g., |
| 342 |
– |
% % |
| 343 |
– |
% \begin{equation} |
| 344 |
– |
% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
| 345 |
– |
% \end{equation} |
| 346 |
– |
% % |
| 347 |
– |
% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 348 |
– |
% Thus, for $f(r)=1/r$, we find |
| 349 |
– |
% % |
| 350 |
– |
% \begin{equation} |
| 351 |
– |
% 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} . |
| 352 |
– |
% \end{equation} |
| 353 |
– |
% This function is an approximate electrostatic potential that has |
| 354 |
– |
% vanishing second derivatives at the cutoff radius, making it suitable |
| 355 |
– |
% for shifting the forces and torques of charge-dipole interactions. |
| 356 |
– |
|
| 323 |
|
The TSF potential for any multipole-multipole interaction can be |
| 324 |
|
written |
| 325 |
|
\begin{equation} |
| 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 |
|
|
| 371 |
– |
% To carry out the same procedure for a damped electrostatic kernel, we |
| 372 |
– |
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
| 373 |
– |
% Many of the derivatives of the damped kernel are well known from |
| 374 |
– |
% Smith's early work on multipoles for the Ewald |
| 375 |
– |
% summation.\cite{Smith82,Smith98} |
| 376 |
– |
|
| 377 |
– |
% Note that increasing the value of $n$ will add additional terms to the |
| 378 |
– |
% electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 379 |
– |
% $(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 380 |
– |
% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 381 |
– |
% f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 382 |
– |
% for computing multipole energies, forces, and torques, and smooth |
| 383 |
– |
% cutoffs of these quantities can be guaranteed as long as the number of |
| 384 |
– |
% terms in the Taylor series exceeds the derivative order required. |
| 385 |
– |
|
| 337 |
|
For multipole-multipole interactions, following this procedure results |
| 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}_{\bf a} |
| 343 |
< |
\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance |
| 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_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 347 |
< |
\mathbf{D}_{\bf a} \cdot |
| 348 |
< |
\mathbf{D}_{\bf b} \right) v_{21}(r) + |
| 349 |
< |
\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right) |
| 350 |
< |
\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(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 |
|
For the Taylor shifted (TSF) method with the undamped kernel, |
| 372 |
|
which have been projected onto the surface of the cutoff sphere |
| 373 |
|
without changing their relative orientation, |
| 374 |
|
\begin{equation} |
| 375 |
< |
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
| 376 |
< |
U_{D_{\bf a} D_{\bf b}}(r_c) |
| 377 |
< |
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
| 378 |
< |
\nabla U_{D_{\bf a}D_{\bf b}}(r_c). |
| 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}_{\bf |
| 430 |
< |
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
| 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 |
| 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(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} |
| 403 |
< |
\cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \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. In this expression, |
| 408 |
|
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
| 409 |
< |
($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ |
| 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 |
| 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^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 437 |
< |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 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 |
| 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. The parameters used |
| 482 |
< |
in the test cases are given in table~\ref{tab:pars}. |
| 482 |
> |
in the test cases are given in table~\ref{tab:pars}. |
| 483 |
|
|
| 484 |
|
\begin{table} |
| 535 |
– |
\label{tab:pars} |
| 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.} |
| 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} & & & & \\ |
| 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 imporant |
| 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 |
| 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{openmd} which was used for all calculations in |
| 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 |
| 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} In all cases, the quantities |
| 552 |
< |
being compared are the electrostatic contributions to energies, force, |
| 553 |
< |
and torques. All other contributions to these quantities (i.e. from |
| 554 |
< |
Lennard-Jones interactions) are removed prior to the comparisons. |
| 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 |
|
The convergence parameter ($\alpha$) also plays a role in the balance |
| 560 |
|
of the real-space and reciprocal-space portions of the Ewald |
| 693 |
|
|
| 694 |
|
\section{\label{sec:result}RESULTS} |
| 695 |
|
\subsection{Configurational energy differences} |
| 744 |
– |
%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. |
| 745 |
– |
%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. |
| 746 |
– |
%\begin{figure}[h!] |
| 747 |
– |
% \centering |
| 748 |
– |
% \includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf} |
| 749 |
– |
% \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. } |
| 750 |
– |
% \label{fig:rcutConvergence_dipolarLiquid} |
| 751 |
– |
% \end{figure} |
| 752 |
– |
%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. |
| 753 |
– |
%\begin{figure} |
| 754 |
– |
% \centering |
| 755 |
– |
% \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 756 |
– |
% \label{fig:barGraph1} |
| 757 |
– |
% \end{figure} |
| 758 |
– |
% \begin{figure} |
| 759 |
– |
% \centering |
| 760 |
– |
% \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 761 |
– |
% \caption{} |
| 762 |
– |
|
| 763 |
– |
% \label{fig:barGraph2} |
| 764 |
– |
% \end{figure} |
| 765 |
– |
%The correlation coefficient ($R^2$) and slope of the linear |
| 766 |
– |
%regression plots for the energy differences for all six different |
| 767 |
– |
%molecular systems is shown in figure 4a and 4b.The plot shows that |
| 768 |
– |
%the correlation coefficient improves for the SP cutoff method as |
| 769 |
– |
%compared to the undamped hard cutoff method in the case of SSDQC, |
| 770 |
– |
%SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar |
| 771 |
– |
%crystal and liquid, the correlation coefficient is almost unchanged |
| 772 |
– |
%and close to 1. The correlation coefficient is smallest (0.696276 |
| 773 |
– |
%for $r_c$ = 9 $A^\circ$) for the SSDQC liquid because of the presence of |
| 774 |
– |
%charge-charge and charge-multipole interactions. Since the |
| 775 |
– |
%charge-charge and charge-multipole interaction is long ranged, there |
| 776 |
– |
%is huge deviation of correlation coefficient from 1. Similarly, the |
| 777 |
– |
%quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with |
| 778 |
– |
%compared to interactions in the other multipolar systems, thus the |
| 779 |
– |
%correlation coefficient very close to 1 even for hard cutoff |
| 780 |
– |
%method. The idea of placing image multipole on the surface of the |
| 781 |
– |
%cutoff sphere improves the correlation coefficient and makes it close |
| 782 |
– |
%to 1 for all types of multipolar systems. Similarly the slope is |
| 783 |
– |
%hugely deviated from the correct value for the lower order |
| 784 |
– |
%multipole-multipole interaction and slightly deviated for higher |
| 785 |
– |
%order multipole – multipole interaction. The SP method improves both |
| 786 |
– |
%correlation coefficient ($R^2$) and slope significantly in SSDQC and |
| 787 |
– |
%dipolar systems. The Slope is found to be deviated more in dipolar |
| 788 |
– |
%crystal as compared to liquid which is associated with the large |
| 789 |
– |
%fluctuation in the electrostatic energy in crystal. The GSF also |
| 790 |
– |
%produced better values of correlation coefficient and slope with the |
| 791 |
– |
%proper selection of the damping alpha (Interested reader can consult |
| 792 |
– |
%accompanying supporting material). The TSF method gives good value of |
| 793 |
– |
%correlation coefficient for the dipolar crystal, dipolar liquid, |
| 794 |
– |
%SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the |
| 795 |
– |
%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 |
| 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).} |
| 808 |
< |
\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 |
| 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 |
| 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)} |
| 952 |
< |
\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\label{sec:conservation}} |
| 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 12\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 all of the real space methods. A |
| 875 |
< |
sample was also run using the multipolar Ewald sum with the same |
| 876 |
< |
real-space cutoff. |
| 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 |
| 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_12.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 |
| 995 |
< |
energy over time (in kcal / mol / particle / ns) and $\delta |
| 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). All simulations were |
| 904 |
< |
of a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 905 |
< |
300 K starting from the same initial configuration. All runs |
| 906 |
< |
utilized the same real-space cutoff, $r_c = 12$\AA.} |
| 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 |
|
In the first paper in this series, we generalized the |
| 1058 |
|
charge-neutralized electrostatic energy originally developed by Wolf |
| 1065 |
|
distance that prevents its use in molecular dynamics. |
| 1066 |
|
|
| 1067 |
|
We also developed two natural extensions of the damped shifted-force |
| 1068 |
< |
(DSF) model originally proposed by Fennel and |
| 1069 |
< |
Gezelter.\cite{Fennell:2006lq} The GSF and TSF approaches provide |
| 1070 |
< |
smooth truncation of energies, forces, and torques at the real-space |
| 1071 |
< |
cutoff, and both converge to DSF electrostatics for point-charge |
| 1072 |
< |
interactions. The TSF model is based on a high-order truncated Taylor |
| 1073 |
< |
expansion which can be relatively perturbative inside the cutoff |
| 1074 |
< |
sphere. The GSF model takes the gradient from an images of the |
| 1075 |
< |
interacting multipole that has been projected onto the cutoff sphere |
| 1076 |
< |
to derive shifted force and torque expressions, and is a significantly |
| 1077 |
< |
more gentle approach. |
| 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 |
< |
Of the two newly-developed shifted force models, the GSF method |
| 1080 |
< |
produced quantitative agreement with Ewald energy, force, and torques. |
| 1081 |
< |
It also performs well in conserving energy in MD simulations. The |
| 1082 |
< |
Taylor-shifted (TSF) model provides smooth dynamics, but these take |
| 1083 |
< |
place on a potential energy surface that is significantly perturbed |
| 1084 |
< |
from Ewald-based electrostatics. |
| 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 |
< |
% The direct truncation of any electrostatic potential energy without |
| 1095 |
< |
% multipole neutralization creates large fluctuations in molecular |
| 1096 |
< |
% simulations. This fluctuation in the energy is very large for the case |
| 1097 |
< |
% of crystal because of long range of multipole ordering (Refer paper |
| 1098 |
< |
% I).\cite{PaperI} This is also significant in the case of the liquid |
| 1099 |
< |
% because of the local multipole ordering in the molecules. If the net |
| 1100 |
< |
% multipole within cutoff radius neutralized within cutoff sphere by |
| 1101 |
< |
% placing image multiples on the surface of the sphere, this fluctuation |
| 1102 |
< |
% in the energy reduced significantly. Also, the multipole |
| 1103 |
< |
% neutralization in the generalized SP method showed very good agreement |
| 1104 |
< |
% with the Ewald as compared to direct truncation for the evaluation of |
| 1045 |
< |
% the $\triangle E$ between the configurations. In MD simulations, the |
| 1046 |
< |
% energy conservation is very important. The conservation of the total |
| 1047 |
< |
% energy can be ensured by i) enforcing the smooth truncation of the |
| 1048 |
< |
% energy, force and torque in the cutoff radius and ii) making the |
| 1049 |
< |
% energy, force and torque consistent with each other. The GSF and TSF |
| 1050 |
< |
% methods ensure the consistency and smooth truncation of the energy, |
| 1051 |
< |
% force and torque at the cutoff radius, as a result show very good |
| 1052 |
< |
% total energy conservation. But the TSF method does not show good |
| 1053 |
< |
% agreement in the absolute value of the electrostatic energy, force and |
| 1054 |
< |
% torque with the Ewald. The GSF method has mimicked Ewald’s force, |
| 1055 |
< |
% energy and torque accurately and also conserved energy. |
| 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 |
| 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, the GSF method is a suitable and |
| 1116 |
< |
efficient replacement for the Ewald sum for evaluating electrostatic |
| 1117 |
< |
interactions in MD simulations. Both methods retain excellent |
| 1118 |
< |
fidelity to the Ewald energies, forces and torques. Additionally, the |
| 1119 |
< |
energy drift and fluctuations from the GSF electrostatics are better |
| 1120 |
< |
than a multipolar Ewald sum for finite-sized reciprocal spaces. |
| 1121 |
< |
Because they use real-space cutoffs with moderate cutoff radii, the |
| 1122 |
< |
GSF and SP models approach $\mathcal{O}(N)$ scaling as the system size |
| 1123 |
< |
increases. Additionally, they can be made extremely efficient using |
| 1075 |
< |
spline interpolations of the radial functions. They require no |
| 1076 |
< |
Fourier transforms or $k$-space sums, and guarantee the smooth |
| 1077 |
< |
handling of energies, forces, and torques as multipoles cross the |
| 1078 |
< |
real-space cutoff boundary. |
| 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 |
