| 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 alternatives to the Ewald Sum. II. Comparison of Methods} |
| 51 |
|
|
| 52 |
|
\author{Madan Lamichhane} |
| 53 |
< |
\affiliation{Department of Physics, University |
| 55 |
< |
of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\ |
| 53 |
> |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
| 54 |
|
|
| 55 |
|
\author{Kathie E. Newman} |
| 56 |
< |
\affiliation{Department of Physics, University |
| 59 |
< |
of Notre Dame, Notre Dame, IN 46556} |
| 56 |
> |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
| 57 |
|
|
| 58 |
|
\author{J. Daniel Gezelter}% |
| 59 |
|
\email{gezelter@nd.edu.} |
| 60 |
< |
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash |
| 61 |
< |
}% |
| 60 |
> |
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556 |
| 61 |
> |
} |
| 62 |
|
|
| 63 |
< |
\date{\today}% It is always \today, today, |
| 67 |
< |
% but any date may be explicitly specified |
| 63 |
> |
\date{\today} |
| 64 |
|
|
| 65 |
|
\begin{abstract} |
| 66 |
< |
We have tested the real-space shifted potential (SP), |
| 67 |
< |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
| 68 |
< |
for multipoles that were developed in the first paper in this series |
| 69 |
< |
against the multipolar Ewald sum as a reference method. The tests |
| 70 |
< |
were carried out in a variety of condensed-phase environments which |
| 71 |
< |
were designed to test all levels of the multipole-multipole |
| 72 |
< |
interactions. Comparisons of the energy differences between |
| 73 |
< |
configurations, molecular forces, and torques were used to analyze |
| 74 |
< |
how well the real-space models perform relative to the more |
| 75 |
< |
computationally expensive Ewald treatment. We have also investigated the |
| 76 |
< |
energy conservation properties of the new methods in molecular |
| 81 |
< |
dynamics simulations using all of these methods. The SP method shows |
| 66 |
> |
We report on tests of the shifted potential (SP), gradient shifted |
| 67 |
> |
force (GSF), and Taylor shifted force (TSF) real-space methods for |
| 68 |
> |
multipole interactions developed in the first paper in this series, |
| 69 |
> |
using the multipolar Ewald sum as a reference method. The tests were |
| 70 |
> |
carried out in a variety of condensed-phase environments designed to |
| 71 |
> |
test up to quadrupole-quadrupole interactions. Comparisons of the |
| 72 |
> |
energy differences between configurations, molecular forces, and |
| 73 |
> |
torques were used to analyze how well the real-space models perform |
| 74 |
> |
relative to the more computationally expensive Ewald treatment. We |
| 75 |
> |
have also investigated the energy conservation properties of the new |
| 76 |
> |
methods in molecular dynamics simulations. The SP method shows |
| 77 |
|
excellent agreement with configurational energy differences, forces, |
| 78 |
|
and torques, and would be suitable for use in Monte Carlo |
| 79 |
|
calculations. Of the two new shifted-force methods, the GSF |
| 80 |
|
approach shows the best agreement with Ewald-derived energies, |
| 81 |
< |
forces, and torques and exhibits energy conservation properties that |
| 82 |
< |
make it an excellent choice for efficiently computing electrostatic |
| 83 |
< |
interactions in molecular dynamics simulations. |
| 81 |
> |
forces, and torques and also exhibits energy conservation properties |
| 82 |
> |
that make it an excellent choice for efficient computation of |
| 83 |
> |
electrostatic interactions in molecular dynamics simulations. |
| 84 |
|
\end{abstract} |
| 85 |
|
|
| 86 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 89 |
|
|
| 90 |
|
\maketitle |
| 91 |
|
|
| 97 |
– |
|
| 92 |
|
\section{\label{sec:intro}Introduction} |
| 93 |
|
Computing the interactions between electrostatic sites is one of the |
| 94 |
< |
most expensive aspects of molecular simulations, which is why there |
| 95 |
< |
have been significant efforts to develop practical, efficient and |
| 96 |
< |
convergent methods for handling these interactions. Ewald's method is |
| 97 |
< |
perhaps the best known and most accurate method for evaluating |
| 98 |
< |
energies, forces, and torques in explicitly-periodic simulation |
| 99 |
< |
cells. In this approach, the conditionally convergent electrostatic |
| 100 |
< |
energy is converted into two absolutely convergent contributions, one |
| 101 |
< |
which is carried out in real space with a cutoff radius, and one in |
| 102 |
< |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
| 94 |
> |
most expensive aspects of molecular simulations. There have been |
| 95 |
> |
significant efforts to develop practical, efficient and convergent |
| 96 |
> |
methods for handling these interactions. Ewald's method is perhaps the |
| 97 |
> |
best known and most accurate method for evaluating energies, forces, |
| 98 |
> |
and torques in explicitly-periodic simulation cells. In this approach, |
| 99 |
> |
the conditionally convergent electrostatic energy is converted into |
| 100 |
> |
two absolutely convergent contributions, one which is carried out in |
| 101 |
> |
real space with a cutoff radius, and one in reciprocal |
| 102 |
> |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
| 103 |
|
|
| 104 |
|
When carried out as originally formulated, the reciprocal-space |
| 105 |
|
portion of the Ewald sum exhibits relatively poor computational |
| 106 |
< |
scaling, making it prohibitive for large systems. By utilizing |
| 107 |
< |
particle meshes and three dimensional fast Fourier transforms (FFT), |
| 108 |
< |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 109 |
< |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
| 110 |
< |
the computational cost from $O(N^2)$ down to $O(N \log |
| 111 |
< |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
| 106 |
> |
scaling, making it prohibitive for large systems. By utilizing a |
| 107 |
> |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
| 108 |
> |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 109 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
| 110 |
> |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
| 111 |
> |
\log |
| 112 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
| 113 |
|
|
| 114 |
< |
Because of the artificial periodicity required for the Ewald sum, the |
| 120 |
< |
method may require modification to compute interactions for |
| 114 |
> |
Because of the artificial periodicity required for the Ewald sum, |
| 115 |
|
interfacial molecular systems such as membranes and liquid-vapor |
| 116 |
< |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 117 |
< |
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
| 118 |
< |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
| 119 |
< |
periodicity in the Ewald method can also be problematic for molecular |
| 120 |
< |
interfaces.\cite{Fennell:2006lq} Modified Ewald methods that were |
| 121 |
< |
developed to handle two-dimensional (2D) electrostatic interactions in |
| 122 |
< |
interfacial systems have not seen similar particle-mesh |
| 123 |
< |
treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 124 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
| 125 |
< |
with system size. |
| 116 |
> |
interfaces require modifications to the method. Parry's extension of |
| 117 |
> |
the three dimensional Ewald sum is appropriate for slab |
| 118 |
> |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
| 119 |
> |
developed to handle two-dimensional (2-D) electrostatic |
| 120 |
> |
interactions,\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 121 |
> |
but these methods were originally quite computationally |
| 122 |
> |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
| 123 |
> |
efforts that reduced the computational cost of 2-D lattice |
| 124 |
> |
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 125 |
> |
bringing them more in line with the scaling for the full 3-D |
| 126 |
> |
treatments. The inherent periodicity in the Ewald’s method can also |
| 127 |
> |
be problematic for interfacial molecular |
| 128 |
> |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
| 129 |
|
|
| 130 |
|
\subsection{Real-space methods} |
| 131 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
| 132 |
|
method for calculating electrostatic interactions between point |
| 133 |
< |
charges. They argued that the effective Coulomb interaction in |
| 134 |
< |
condensed phase systems is actually short ranged.\cite{Wolf92,Wolf95} |
| 135 |
< |
For an ordered lattice (e.g., when computing the Madelung constant of |
| 136 |
< |
an ionic solid), the material can be considered as a set of ions |
| 137 |
< |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
| 138 |
< |
an effective distance dependence for the electrostatic interactions of |
| 139 |
< |
$r^{-5}$ (see figure \ref{fig:schematic}). For this reason, careful |
| 140 |
< |
application of Wolf's method can obtain accurate estimates of Madelung |
| 141 |
< |
constants using relatively short cutoff radii. Recently, Fukuda used |
| 142 |
< |
neutralization of the higher order moments for calculation of the |
| 143 |
< |
electrostatic interactions in point charge |
| 144 |
< |
systems.\cite{Fukuda:2013sf} |
| 133 |
> |
charges. They argued that the effective Coulomb interaction in most |
| 134 |
> |
condensed phase systems is effectively short |
| 135 |
> |
ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when |
| 136 |
> |
computing the Madelung constant of an ionic solid), the material can |
| 137 |
> |
be considered as a set of ions interacting with neutral dipolar or |
| 138 |
> |
quadrupolar ``molecules'' giving an effective distance dependence for |
| 139 |
> |
the electrostatic interactions of $r^{-5}$ (see figure |
| 140 |
> |
\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple |
| 141 |
> |
cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the |
| 142 |
> |
electrostatic energy per ion converges more rapidly to the Madelung |
| 143 |
> |
energy than the dipolar approximation.\cite{Wolf92} To find the |
| 144 |
> |
correct Madelung constant, Lacman suggested that the NaCl structure |
| 145 |
> |
could be constructed in a way that the finite crystal terminates with |
| 146 |
> |
complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees |
| 147 |
> |
what is effectively a set of octupoles at large distances. These facts |
| 148 |
> |
suggest that the Madelung constants are relatively short ranged for |
| 149 |
> |
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
| 150 |
> |
application of Wolf's method can provide accurate estimates of |
| 151 |
> |
Madelung constants using relatively short cutoff radii. |
| 152 |
|
|
| 153 |
+ |
Direct truncation of interactions at a cutoff radius creates numerical |
| 154 |
+ |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
| 155 |
+ |
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
| 156 |
+ |
neutralize this charge they proposed placing an image charge on the |
| 157 |
+ |
surface of the cutoff sphere for every real charge inside the cutoff. |
| 158 |
+ |
These charges are present for the evaluation of both the pair |
| 159 |
+ |
interaction energy and the force, although the force expression |
| 160 |
+ |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
| 161 |
+ |
formulation, the total energy for the charge and image were not equal |
| 162 |
+ |
to the integral of the force expression, and as a result, the total |
| 163 |
+ |
energy would not be conserved in molecular dynamics (MD) |
| 164 |
+ |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
| 165 |
+ |
Gezelter later proposed shifted force variants of the Wolf method with |
| 166 |
+ |
commensurate force and energy expressions that do not exhibit this |
| 167 |
+ |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
| 168 |
+ |
were also proposed by Chen \textit{et |
| 169 |
+ |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
| 170 |
+ |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
| 171 |
+ |
used additional neutralization of higher order moments for systems of |
| 172 |
+ |
point charges.\cite{Fukuda:2013sf} |
| 173 |
+ |
|
| 174 |
|
\begin{figure} |
| 175 |
|
\centering |
| 176 |
|
\includegraphics[width=\linewidth]{schematic.pdf} |
| 188 |
|
\label{fig:schematic} |
| 189 |
|
\end{figure} |
| 190 |
|
|
| 191 |
< |
The direct truncation of interactions at a cutoff radius creates |
| 192 |
< |
truncation defects. Wolf \textit{et al.} argued that |
| 193 |
< |
truncation errors are due to net charge remaining inside the cutoff |
| 194 |
< |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
| 195 |
< |
placing an image charge on the surface of the cutoff sphere for every |
| 196 |
< |
real charge inside the cutoff. These charges are present for the |
| 197 |
< |
evaluation of both the pair interaction energy and the force, although |
| 198 |
< |
the force expression maintained a discontinuity at the cutoff sphere. |
| 199 |
< |
In the original Wolf formulation, the total energy for the charge and |
| 200 |
< |
image were not equal to the integral of their force expression, and as |
| 176 |
< |
a result, the total energy would not be conserved in molecular |
| 177 |
< |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
| 178 |
< |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
| 179 |
< |
method with commensurate force and energy expressions that do not |
| 180 |
< |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
| 181 |
< |
methods were also proposed by Chen \textit{et |
| 182 |
< |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
| 183 |
< |
and by Wu and Brooks.\cite{Wu:044107} |
| 184 |
< |
|
| 185 |
< |
Considering the interaction of one central ion in an ionic crystal |
| 186 |
< |
with a portion of the crystal at some distance, the effective Columbic |
| 187 |
< |
potential is found to decrease as $r^{-5}$. If one views the \ce{NaCl} |
| 188 |
< |
crystal as a simple cubic (SC) structure with an octupolar |
| 189 |
< |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
| 190 |
< |
rapidly to the Madelung energy than the dipolar |
| 191 |
< |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
| 192 |
< |
Lacman suggested that the NaCl structure could be constructed in a way |
| 193 |
< |
that the finite crystal terminates with complete \ce{(NaCl)4} |
| 194 |
< |
molecules.\cite{Lacman65} The central ion sees what is effectively a |
| 195 |
< |
set of octupoles at large distances. These facts suggest that the |
| 196 |
< |
Madelung constants are relatively short ranged for perfect ionic |
| 197 |
< |
crystals.\cite{Wolf:1999dn} |
| 198 |
< |
|
| 199 |
< |
One can make a similar argument for crystals of point multipoles. The |
| 200 |
< |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
| 201 |
< |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
| 202 |
< |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
| 203 |
< |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
| 204 |
< |
zero net dipole and retain contributions only from higher order |
| 205 |
< |
multipoles. The effective interaction between a dipole at the center |
| 191 |
> |
One can make a similar effective range argument for crystals of point |
| 192 |
> |
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
| 193 |
> |
constants for dipolar lattices utilizes 24 basis vectors that contain |
| 194 |
> |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
| 195 |
> |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
| 196 |
> |
moments, while the rest have zero net dipole and retain contributions |
| 197 |
> |
only from higher order multipoles. The lowest-energy crystalline |
| 198 |
> |
structures are built out of basis vectors that have only residual |
| 199 |
> |
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
| 200 |
> |
structures, the effective interaction between a dipole at the center |
| 201 |
|
of a crystal and a group of eight dipoles farther away is |
| 202 |
|
significantly shorter ranged than the $r^{-3}$ that one would expect |
| 203 |
|
for raw dipole-dipole interactions. Only in crystals which retain a |
| 217 |
|
|
| 218 |
|
The shorter effective range of electrostatic interactions is not |
| 219 |
|
limited to perfect crystals, but can also apply in disordered fluids. |
| 220 |
< |
Even at elevated temperatures, there is, on average, local charge |
| 221 |
< |
balance in an ionic liquid, where each positive ion has surroundings |
| 222 |
< |
dominated by negaitve ions and vice versa. The reversed-charge images |
| 223 |
< |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
| 224 |
< |
retain the effective multipolar interactions as the charges traverse |
| 225 |
< |
the cutoff boundary. |
| 220 |
> |
Even at elevated temperatures, there is local charge balance in an |
| 221 |
> |
ionic liquid, where each positive ion has surroundings dominated by |
| 222 |
> |
negaitve ions and vice versa. The reversed-charge images on the |
| 223 |
> |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
| 224 |
> |
the effective multipolar interactions as the charges traverse the |
| 225 |
> |
cutoff boundary. |
| 226 |
|
|
| 227 |
|
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
| 228 |
|
significant orientational averaging that additionally reduces the |
| 241 |
|
% to the non-neutralized value of the higher order moments within the |
| 242 |
|
% cutoff sphere. |
| 243 |
|
|
| 244 |
< |
The forces and torques acting on atomic sites are the fundamental |
| 245 |
< |
factors driving dynamics in molecular simulations. Fennell and |
| 246 |
< |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 247 |
< |
obtain consistent energies and forces on the atoms within the cutoff |
| 248 |
< |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 249 |
< |
aproaches the cutoff radius. The comparisons of the accuracy these |
| 250 |
< |
quantities between the DSF kernel and SPME was surprisingly |
| 251 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
| 252 |
< |
calculating electrostatic interactions in molecular systems with |
| 253 |
< |
relatively uniform charge |
| 259 |
< |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 244 |
> |
Forces and torques acting on atomic sites are fundamental in driving |
| 245 |
> |
dynamics in molecular simulations, and the damped shifted force (DSF) |
| 246 |
> |
energy kernel provides consistent energies and forces on charged atoms |
| 247 |
> |
within the cutoff sphere. Both the energy and the force go smoothly to |
| 248 |
> |
zero as an atom aproaches the cutoff radius. The comparisons of the |
| 249 |
> |
accuracy these quantities between the DSF kernel and SPME was |
| 250 |
> |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
| 251 |
> |
has seen increasing use in molecular systems with relatively uniform |
| 252 |
> |
charge |
| 253 |
> |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
| 254 |
|
|
| 255 |
|
\subsection{The damping function} |
| 256 |
< |
The damping function used in our research has been discussed in detail |
| 257 |
< |
in the first paper of this series.\cite{PaperI} The radial kernel |
| 258 |
< |
$1/r$ for the interactions between point charges can be replaced by |
| 259 |
< |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
| 260 |
< |
accelerate the rate of convergence, where $\alpha$ is a damping |
| 261 |
< |
parameter with units of inverse distance. Altering the value of |
| 262 |
< |
$\alpha$ is equivalent to changing the width of Gaussian charge |
| 263 |
< |
distributions that replace each point charge -- Gaussian overlap |
| 264 |
< |
integrals yield complementary error functions when truncated at a |
| 265 |
< |
finite distance. |
| 256 |
> |
The damping function has been discussed in detail in the first paper |
| 257 |
> |
of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the |
| 258 |
> |
interactions between point charges can be replaced by the |
| 259 |
> |
complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate |
| 260 |
> |
convergence, where $\alpha$ is a damping parameter with units of |
| 261 |
> |
inverse distance. Altering the value of $\alpha$ is equivalent to |
| 262 |
> |
changing the width of Gaussian charge distributions that replace each |
| 263 |
> |
point charge, as Coulomb integrals with Gaussian charge distributions |
| 264 |
> |
produce complementary error functions when truncated at a finite |
| 265 |
> |
distance. |
| 266 |
|
|
| 267 |
< |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
| 268 |
< |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
| 269 |
< |
agreement with SPME for the interaction energies, forces and torques |
| 270 |
< |
for charge-charge interactions.\cite{Fennell:2006lq} |
| 267 |
> |
With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method |
| 268 |
> |
produced very good agreement with SPME for interaction energies, |
| 269 |
> |
forces and torques for charge-charge |
| 270 |
> |
interactions.\cite{Fennell:2006lq} |
| 271 |
|
|
| 272 |
|
\subsection{Point multipoles in molecular modeling} |
| 273 |
|
Coarse-graining approaches which treat entire molecular subsystems as |
| 274 |
|
a single rigid body are now widely used. A common feature of many |
| 275 |
|
coarse-graining approaches is simplification of the electrostatic |
| 276 |
|
interactions between bodies so that fewer site-site interactions are |
| 277 |
< |
required to compute configurational energies. Many coarse-grained |
| 278 |
< |
molecular structures would normally consist of equal positive and |
| 285 |
< |
negative charges, and rather than use multiple site-site interactions, |
| 286 |
< |
the interaction between higher order multipoles can also be used to |
| 287 |
< |
evaluate a single molecule-molecule |
| 288 |
< |
interaction.\cite{Ren06,Essex10,Essex11} |
| 277 |
> |
required to compute configurational |
| 278 |
> |
energies.\cite{Ren06,Essex10,Essex11} |
| 279 |
|
|
| 280 |
< |
Because electrons in a molecule are not localized at specific points, |
| 281 |
< |
the assignment of partial charges to atomic centers is a relatively |
| 282 |
< |
rough approximation. Atomic sites can also be assigned point |
| 283 |
< |
multipoles and polarizabilities to increase the accuracy of the |
| 284 |
< |
molecular model. Recently, water has been modeled with point |
| 285 |
< |
multipoles up to octupolar order using the soft sticky |
| 296 |
< |
dipole-quadrupole-octupole (SSDQO) |
| 280 |
> |
Additionally, because electrons in a molecule are not localized at |
| 281 |
> |
specific points, the assignment of partial charges to atomic centers |
| 282 |
> |
is always an approximation. For increased accuracy, atomic sites can |
| 283 |
> |
also be assigned point multipoles and polarizabilities. Recently, |
| 284 |
> |
water has been modeled with point multipoles up to octupolar order |
| 285 |
> |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
| 286 |
|
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 287 |
|
multipoles up to quadrupolar order have also been coupled with point |
| 288 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
| 289 |
< |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But |
| 290 |
< |
using point multipole with the real space truncation without |
| 291 |
< |
accounting for multipolar neutrality will create energy conservation |
| 292 |
< |
issues in molecular dynamics (MD) simulations. |
| 289 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
| 290 |
> |
truncating point multipoles without smoothing the forces and torques |
| 291 |
> |
can create energy conservation issues in molecular dynamics |
| 292 |
> |
simulations. |
| 293 |
|
|
| 294 |
|
In this paper we test a set of real-space methods that were developed |
| 295 |
|
for point multipolar interactions. These methods extend the damped |
| 296 |
|
shifted force (DSF) and Wolf methods originally developed for |
| 297 |
|
charge-charge interactions and generalize them for higher order |
| 298 |
< |
multipoles. The detailed mathematical development of these methods has |
| 299 |
< |
been presented in the first paper in this series, while this work |
| 300 |
< |
covers the testing the energies, forces, torques, and energy |
| 298 |
> |
multipoles. The detailed mathematical development of these methods |
| 299 |
> |
has been presented in the first paper in this series, while this work |
| 300 |
> |
covers the testing of energies, forces, torques, and energy |
| 301 |
|
conservation properties of the methods in realistic simulation |
| 302 |
|
environments. In all cases, the methods are compared with the |
| 303 |
< |
reference method, a full multipolar Ewald treatment. |
| 303 |
> |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
| 304 |
|
|
| 305 |
|
|
| 306 |
|
%\subsection{Conservation of total energy } |
| 585 |
|
\subsection{Implementation} |
| 586 |
|
The real-space methods developed in the first paper in this series |
| 587 |
|
have been implemented in our group's open source molecular simulation |
| 588 |
< |
program, OpenMD,\cite{openmd} which was used for all calculations in |
| 588 |
> |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
| 589 |
|
this work. The complementary error function can be a relatively slow |
| 590 |
|
function on some processors, so all of the radial functions are |
| 591 |
|
precomputed on a fine grid and are spline-interpolated to provide |
