| 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 |
| 88 |
< |
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 properties of the new |
| 77 |
> |
methods in molecular dynamics simulations. The SP method shows |
| 78 |
> |
excellent agreement with configurational energy differences, forces, |
| 79 |
> |
and torques, and would be suitable for use in Monte Carlo |
| 80 |
> |
calculations. Of the two new shifted-force methods, the GSF |
| 81 |
> |
approach shows the best agreement with Ewald-derived energies, |
| 82 |
> |
forces, and torques and also exhibits energy conservation properties |
| 83 |
> |
that make it an excellent choice for efficient computation of |
| 84 |
> |
electrostatic interactions in molecular dynamics simulations. |
| 85 |
|
\end{abstract} |
| 86 |
|
|
| 87 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 90 |
|
|
| 91 |
|
\maketitle |
| 92 |
|
|
| 97 |
– |
|
| 93 |
|
\section{\label{sec:intro}Introduction} |
| 94 |
|
Computing the interactions between electrostatic sites is one of the |
| 95 |
|
most expensive aspects of molecular simulations. There have been |
| 100 |
|
the conditionally convergent electrostatic energy is converted into |
| 101 |
|
two absolutely convergent contributions, one which is carried out in |
| 102 |
|
real space with a cutoff radius, and one in reciprocal |
| 103 |
< |
space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
| 103 |
> |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
| 104 |
|
|
| 105 |
|
When carried out as originally formulated, the reciprocal-space |
| 106 |
|
portion of the Ewald sum exhibits relatively poor computational |
| 107 |
< |
scaling, making it prohibitive for large systems. By utilizing |
| 108 |
< |
particle meshes and three dimensional fast Fourier transforms (FFT), |
| 109 |
< |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 110 |
< |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
| 111 |
< |
the computational cost from $O(N^2)$ down to $O(N \log |
| 112 |
< |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
| 107 |
> |
scaling, making it prohibitive for large systems. By utilizing a |
| 108 |
> |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
| 109 |
> |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 110 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
| 111 |
> |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
| 112 |
> |
\log |
| 113 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
| 114 |
|
|
| 115 |
|
Because of the artificial periodicity required for the Ewald sum, |
| 116 |
|
interfacial molecular systems such as membranes and liquid-vapor |
| 117 |
< |
interfaces require modifications to the |
| 118 |
< |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 119 |
< |
Parry's extension of the three dimensional Ewald sum is appropriate |
| 120 |
< |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
| 121 |
< |
were developed to handle two-dimensional (2D) electrostatic |
| 122 |
< |
interactions in interfacial systems have not seen similar |
| 123 |
< |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 124 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
| 125 |
< |
with system size. The inherent periodicity in the Ewald’s method can |
| 126 |
< |
also be problematic for interfacial molecular |
| 127 |
< |
systems.\cite{Fennell:2006lq} |
| 117 |
> |
interfaces require modifications to the method. Parry's extension of |
| 118 |
> |
the three dimensional Ewald sum is appropriate for slab |
| 119 |
> |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
| 120 |
> |
developed to handle two-dimensional (2-D) electrostatic |
| 121 |
> |
interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 122 |
> |
These methods were originally quite computationally |
| 123 |
> |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
| 124 |
> |
efforts that reduced the computational cost of 2-D lattice summations, |
| 125 |
> |
bringing them more in line with the scaling for the full 3-D |
| 126 |
> |
treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The |
| 127 |
> |
inherent periodicity required by the Ewald method can also be |
| 128 |
> |
problematic in a number of protein/solvent and ionic solution |
| 129 |
> |
environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
| 130 |
|
|
| 131 |
|
\subsection{Real-space methods} |
| 132 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
| 148 |
|
what is effectively a set of octupoles at large distances. These facts |
| 149 |
|
suggest that the Madelung constants are relatively short ranged for |
| 150 |
|
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
| 151 |
< |
application of Wolf's method are able to obtain accurate estimates of |
| 151 |
> |
application of Wolf's method can provide accurate estimates of |
| 152 |
|
Madelung constants using relatively short cutoff radii. |
| 153 |
|
|
| 154 |
|
Direct truncation of interactions at a cutoff radius creates numerical |
| 155 |
< |
errors. Wolf \textit{et al.} argued that truncation errors are due |
| 155 |
> |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
| 156 |
|
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
| 157 |
|
neutralize this charge they proposed placing an image charge on the |
| 158 |
|
surface of the cutoff sphere for every real charge inside the cutoff. |
| 159 |
|
These charges are present for the evaluation of both the pair |
| 160 |
|
interaction energy and the force, although the force expression |
| 161 |
< |
maintained a discontinuity at the cutoff sphere. In the original Wolf |
| 161 |
> |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
| 162 |
|
formulation, the total energy for the charge and image were not equal |
| 163 |
< |
to the integral of their force expression, and as a result, the total |
| 163 |
> |
to the integral of the force expression, and as a result, the total |
| 164 |
|
energy would not be conserved in molecular dynamics (MD) |
| 165 |
|
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
| 166 |
|
Gezelter later proposed shifted force variants of the Wolf method with |
| 167 |
|
commensurate force and energy expressions that do not exhibit this |
| 168 |
< |
problem.\cite{Fennell:2006lq} Related real-space methods were also |
| 169 |
< |
proposed by Chen \textit{et |
| 168 |
> |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
| 169 |
> |
were also proposed by Chen \textit{et |
| 170 |
|
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
| 171 |
< |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has used |
| 172 |
< |
neutralization of the higher order moments for the calculation of the |
| 173 |
< |
electrostatic interaction of the point charge |
| 176 |
< |
systems.\cite{Fukuda:2013sf} |
| 171 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
| 172 |
> |
used additional neutralization of higher order moments for systems of |
| 173 |
> |
point charges.\cite{Fukuda:2013sf} |
| 174 |
|
|
| 175 |
|
\begin{figure} |
| 176 |
|
\centering |
| 177 |
< |
\includegraphics[width=\linewidth]{schematic.pdf} |
| 177 |
> |
\includegraphics[width=\linewidth]{schematic.eps} |
| 178 |
|
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
| 179 |
|
charges (in the smaller grey circle), so interactions are |
| 180 |
|
effectively charge-multipole at longer distances. With hard |
| 192 |
|
One can make a similar effective range argument for crystals of point |
| 193 |
|
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
| 194 |
|
constants for dipolar lattices utilizes 24 basis vectors that contain |
| 195 |
< |
dipoles at the eight corners of a unit cube. Only three of these |
| 196 |
< |
basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
| 195 |
> |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
| 196 |
> |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
| 197 |
|
moments, while the rest have zero net dipole and retain contributions |
| 198 |
< |
only from higher order multipoles. The lowest energy crystalline |
| 198 |
> |
only from higher order multipoles. The lowest-energy crystalline |
| 199 |
|
structures are built out of basis vectors that have only residual |
| 200 |
|
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
| 201 |
|
structures, the effective interaction between a dipole at the center |
| 218 |
|
|
| 219 |
|
The shorter effective range of electrostatic interactions is not |
| 220 |
|
limited to perfect crystals, but can also apply in disordered fluids. |
| 221 |
< |
Even at elevated temperatures, there is, on average, local charge |
| 222 |
< |
balance in an ionic liquid, where each positive ion has surroundings |
| 223 |
< |
dominated by negaitve ions and vice versa. The reversed-charge images |
| 224 |
< |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
| 225 |
< |
retain the effective multipolar interactions as the charges traverse |
| 226 |
< |
the cutoff boundary. |
| 221 |
> |
Even at elevated temperatures, there is local charge balance in an |
| 222 |
> |
ionic liquid, where each positive ion has surroundings dominated by |
| 223 |
> |
negaitve ions and vice versa. The reversed-charge images on the |
| 224 |
> |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
| 225 |
> |
the effective multipolar interactions as the charges traverse the |
| 226 |
> |
cutoff boundary. |
| 227 |
|
|
| 228 |
|
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
| 229 |
|
significant orientational averaging that additionally reduces the |
| 242 |
|
% to the non-neutralized value of the higher order moments within the |
| 243 |
|
% cutoff sphere. |
| 244 |
|
|
| 245 |
< |
The forces and torques acting on atomic sites are the fundamental |
| 246 |
< |
factors driving dynamics in molecular simulations. Fennell and |
| 247 |
< |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 248 |
< |
obtain consistent energies and forces on the atoms within the cutoff |
| 249 |
< |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 250 |
< |
aproaches the cutoff radius. The comparisons of the accuracy these |
| 251 |
< |
quantities between the DSF kernel and SPME was surprisingly |
| 252 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
| 253 |
< |
calculating electrostatic interactions in molecular systems with |
| 254 |
< |
relatively uniform charge |
| 258 |
< |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 245 |
> |
Forces and torques acting on atomic sites are fundamental in driving |
| 246 |
> |
dynamics in molecular simulations, and the damped shifted force (DSF) |
| 247 |
> |
energy kernel provides consistent energies and forces on charged atoms |
| 248 |
> |
within the cutoff sphere. Both the energy and the force go smoothly to |
| 249 |
> |
zero as an atom aproaches the cutoff radius. The comparisons of the |
| 250 |
> |
accuracy these quantities between the DSF kernel and SPME was |
| 251 |
> |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
| 252 |
> |
has seen increasing use in molecular systems with relatively uniform |
| 253 |
> |
charge |
| 254 |
> |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
| 255 |
|
|
| 256 |
|
\subsection{The damping function} |
| 257 |
|
The damping function has been discussed in detail in the first paper |
| 278 |
|
required to compute configurational |
| 279 |
|
energies.\cite{Ren06,Essex10,Essex11} |
| 280 |
|
|
| 281 |
< |
Because electrons in a molecule are not localized at specific points, |
| 282 |
< |
the assignment of partial charges to atomic centers is always an |
| 283 |
< |
approximation. Atomic sites can also be assigned point multipoles and |
| 284 |
< |
polarizabilities to increase the accuracy of the molecular model. |
| 285 |
< |
Recently, water has been modeled with point multipoles up to octupolar |
| 286 |
< |
order using the soft sticky dipole-quadrupole-octupole (SSDQO) |
| 281 |
> |
Additionally, because electrons in a molecule are not localized at |
| 282 |
> |
specific points, the assignment of partial charges to atomic centers |
| 283 |
> |
is always an approximation. For increased accuracy, atomic sites can |
| 284 |
> |
also be assigned point multipoles and polarizabilities. Recently, |
| 285 |
> |
water has been modeled with point multipoles up to octupolar order |
| 286 |
> |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
| 287 |
|
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 288 |
|
multipoles up to quadrupolar order have also been coupled with point |
| 289 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
| 290 |
|
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
| 291 |
|
truncating point multipoles without smoothing the forces and torques |
| 292 |
< |
will create energy conservation issues in molecular dynamics simulations. |
| 292 |
> |
can create energy conservation issues in molecular dynamics |
| 293 |
> |
simulations. |
| 294 |
|
|
| 295 |
|
In this paper we test a set of real-space methods that were developed |
| 296 |
|
for point multipolar interactions. These methods extend the damped |
| 297 |
|
shifted force (DSF) and Wolf methods originally developed for |
| 298 |
|
charge-charge interactions and generalize them for higher order |
| 299 |
< |
multipoles. The detailed mathematical development of these methods has |
| 300 |
< |
been presented in the first paper in this series, while this work |
| 301 |
< |
covers the testing the energies, forces, torques, and energy |
| 299 |
> |
multipoles. The detailed mathematical development of these methods |
| 300 |
> |
has been presented in the first paper in this series, while this work |
| 301 |
> |
covers the testing of energies, forces, torques, and energy |
| 302 |
|
conservation properties of the methods in realistic simulation |
| 303 |
|
environments. In all cases, the methods are compared with the |
| 304 |
< |
reference method, a full multipolar Ewald treatment. |
| 304 |
> |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
| 305 |
|
|
| 306 |
|
|
| 307 |
|
%\subsection{Conservation of total energy } |
| 325 |
|
expressed as the product of two multipole operators and a Coulombic |
| 326 |
|
kernel, |
| 327 |
|
\begin{equation} |
| 328 |
< |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
| 328 |
> |
U_{ab}(r)= M_{a} M_{b} \frac{1}{r} \label{kernel}. |
| 329 |
|
\end{equation} |
| 330 |
< |
where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is |
| 331 |
< |
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
| 332 |
< |
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
| 336 |
< |
$\bf a$, etc. |
| 330 |
> |
where the multipole operator for site $a$, $M_{a}$, is |
| 331 |
> |
expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object |
| 332 |
> |
$a$, etc. |
| 333 |
|
|
| 334 |
|
% Interactions between multipoles can be expressed as higher derivatives |
| 335 |
|
% of the bare Coulomb potential, so one way of ensuring that the forces |
| 384 |
|
contributions to the potential, and ensures that the forces and |
| 385 |
|
torques from each of these contributions will vanish at the cutoff |
| 386 |
|
radius. For example, the direct dipole dot product |
| 387 |
< |
($\mathbf{D}_{\bf a} |
| 388 |
< |
\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance |
| 387 |
> |
($\mathbf{D}_{a} |
| 388 |
> |
\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance |
| 389 |
|
dot products: |
| 390 |
|
\begin{equation} |
| 391 |
< |
U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 392 |
< |
\mathbf{D}_{\bf a} \cdot |
| 393 |
< |
\mathbf{D}_{\bf b} \right) v_{21}(r) + |
| 394 |
< |
\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right) |
| 395 |
< |
\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right] |
| 391 |
> |
U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 392 |
> |
\mathbf{D}_{a} \cdot |
| 393 |
> |
\mathbf{D}_{b} \right) v_{21}(r) + |
| 394 |
> |
\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right) |
| 395 |
> |
\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right] |
| 396 |
|
\end{equation} |
| 397 |
|
|
| 398 |
|
For the Taylor shifted (TSF) method with the undamped kernel, |
| 417 |
|
which have been projected onto the surface of the cutoff sphere |
| 418 |
|
without changing their relative orientation, |
| 419 |
|
\begin{equation} |
| 420 |
< |
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
| 421 |
< |
U_{D_{\bf a} D_{\bf b}}(r_c) |
| 422 |
< |
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
| 423 |
< |
\nabla U_{D_{\bf a}D_{\bf b}}(r_c). |
| 420 |
> |
U_{D_{a}D_{b}}(r) = U_{D_{a}D_{b}}(r) - |
| 421 |
> |
U_{D_{a}D_{b}}(r_c) |
| 422 |
> |
- (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot |
| 423 |
> |
\nabla U_{D_{a}D_{b}}(r_c). |
| 424 |
|
\end{equation} |
| 425 |
< |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf |
| 430 |
< |
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
| 425 |
> |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance |
| 426 |
|
(although the signs are reversed for the dipole that has been |
| 427 |
|
projected onto the cutoff sphere). In many ways, this simpler |
| 428 |
|
approach is closer in spirit to the original shifted force method, in |
| 443 |
|
In general, the gradient shifted potential between a central multipole |
| 444 |
|
and any multipolar site inside the cutoff radius is given by, |
| 445 |
|
\begin{equation} |
| 446 |
< |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 447 |
< |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} |
| 448 |
< |
\cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 446 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
| 447 |
> |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c) |
| 448 |
> |
\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
| 449 |
|
\label{generic2} |
| 450 |
|
\end{equation} |
| 451 |
|
where the sum describes a separate force-shifting that is applied to |
| 452 |
|
each orientational contribution to the energy. In this expression, |
| 453 |
|
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
| 454 |
< |
($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ |
| 454 |
> |
($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$ |
| 455 |
|
represent the orientations the multipoles. |
| 456 |
|
|
| 457 |
|
The third term converges more rapidly than the first two terms as a |
| 478 |
|
interactions with the central multipole and the image. This |
| 479 |
|
effectively shifts the total potential to zero at the cutoff radius, |
| 480 |
|
\begin{equation} |
| 481 |
< |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 482 |
< |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 481 |
> |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
| 482 |
> |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
| 483 |
|
\label{eq:SP} |
| 484 |
|
\end{equation} |
| 485 |
|
where the sum describes separate potential shifting that is done for |
| 572 |
|
and have been compared with the values obtained from the multipolar |
| 573 |
|
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
| 574 |
|
between two configurations is the primary quantity that governs how |
| 575 |
< |
the simulation proceeds. These differences are the most imporant |
| 575 |
> |
the simulation proceeds. These differences are the most important |
| 576 |
|
indicators of the reliability of a method even if the absolute |
| 577 |
|
energies are not exact. For each of the multipolar systems listed |
| 578 |
|
above, we have compared the change in electrostatic potential energy |
| 584 |
|
\subsection{Implementation} |
| 585 |
|
The real-space methods developed in the first paper in this series |
| 586 |
|
have been implemented in our group's open source molecular simulation |
| 587 |
< |
program, OpenMD,\cite{openmd} which was used for all calculations in |
| 587 |
> |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
| 588 |
|
this work. The complementary error function can be a relatively slow |
| 589 |
|
function on some processors, so all of the radial functions are |
| 590 |
|
precomputed on a fine grid and are spline-interpolated to provide |
| 791 |
|
|
| 792 |
|
\begin{figure} |
| 793 |
|
\centering |
| 794 |
< |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 794 |
> |
\includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
| 795 |
|
\caption{Statistical analysis of the quality of configurational |
| 796 |
|
energy differences for the real-space electrostatic methods |
| 797 |
|
compared with the reference Ewald sum. Results with a value equal |
| 864 |
|
|
| 865 |
|
\begin{figure} |
| 866 |
|
\centering |
| 867 |
< |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 867 |
> |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} |
| 868 |
|
\caption{Statistical analysis of the quality of the force vector |
| 869 |
|
magnitudes for the real-space electrostatic methods compared with |
| 870 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 878 |
|
|
| 879 |
|
\begin{figure} |
| 880 |
|
\centering |
| 881 |
< |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 881 |
> |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps} |
| 882 |
|
\caption{Statistical analysis of the quality of the torque vector |
| 883 |
|
magnitudes for the real-space electrostatic methods compared with |
| 884 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 936 |
|
|
| 937 |
|
\begin{figure} |
| 938 |
|
\centering |
| 939 |
< |
\includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} |
| 939 |
> |
\includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps} |
| 940 |
|
\caption{The circular variance of the direction of the force and |
| 941 |
|
torque vectors obtained from the real-space methods around the |
| 942 |
|
reference Ewald vectors. A variance equal to 0 (dashed line) |
| 956 |
|
in this series and provides the most comprehensive test of the new |
| 957 |
|
methods. A liquid-phase system was created with 2000 water molecules |
| 958 |
|
and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
| 959 |
< |
temperature of 300K. After equilibration, this liquid-phase system |
| 960 |
< |
was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with |
| 961 |
< |
a cutoff radius of 12\AA. The value of the damping coefficient was |
| 962 |
< |
also varied from the undamped case ($\alpha = 0$) to a heavily damped |
| 963 |
< |
case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods. A |
| 964 |
< |
sample was also run using the multipolar Ewald sum with the same |
| 965 |
< |
real-space cutoff. |
| 959 |
> |
temperature of 300K. After equilibration in the canonical (NVT) |
| 960 |
> |
ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system |
| 961 |
> |
was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald, |
| 962 |
> |
Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA. The |
| 963 |
> |
value of the damping coefficient was also varied from the undamped |
| 964 |
> |
case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$ |
| 965 |
> |
\AA$^{-1}$) for all of the real space methods. A sample was also run |
| 966 |
> |
using the multipolar Ewald sum with the same real-space cutoff. |
| 967 |
|
|
| 968 |
|
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
| 969 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
| 984 |
|
|
| 985 |
|
\begin{figure} |
| 986 |
|
\centering |
| 987 |
< |
\includegraphics[width=\textwidth]{newDrift_12.pdf} |
| 987 |
> |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
| 988 |
|
\label{fig:energyDrift} |
| 989 |
|
\caption{Analysis of the energy conservation of the real-space |
| 990 |
< |
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
| 991 |
< |
energy over time (in kcal / mol / particle / ns) and $\delta |
| 992 |
< |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 993 |
< |
around this drift (in kcal / mol / particle). All simulations were |
| 994 |
< |
of a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 990 |
> |
methods. $\delta \mathrm{E}_1$ is the linear drift in energy over |
| 991 |
> |
time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is |
| 992 |
> |
the standard deviation of energy fluctuations around this drift (in |
| 993 |
> |
kcal / mol / particle). Points that appear below the dashed grey |
| 994 |
> |
(Ewald) lines exhibit better energy conservation than commonly-used |
| 995 |
> |
parameters for Ewald-based electrostatics. All simulations were of |
| 996 |
> |
a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 997 |
|
300 K starting from the same initial configuration. All runs |
| 998 |
|
utilized the same real-space cutoff, $r_c = 12$\AA.} |
| 999 |
|
\end{figure} |
| 1000 |
|
|
| 1001 |
+ |
\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}} |
| 1002 |
+ |
The most important test of the modified interaction potentials is the |
| 1003 |
+ |
fidelity with which they can reproduce structural features and |
| 1004 |
+ |
dynamical properties in a liquid. One commonly-utilized measure of |
| 1005 |
+ |
structural ordering is the pair distribution function, $g(r)$, which |
| 1006 |
+ |
measures local density deviations in relation to the bulk density. In |
| 1007 |
+ |
the electrostatic approaches studied here, the short-range repulsion |
| 1008 |
+ |
from the Lennard-Jones potential is identical for the various |
| 1009 |
+ |
electrostatic methods, and since short range repulsion determines much |
| 1010 |
+ |
of the local liquid ordering, one would not expect to see many |
| 1011 |
+ |
differences in $g(r)$. Indeed, the pair distributions are essentially |
| 1012 |
+ |
identical for all of the electrostatic methods studied (for each of |
| 1013 |
+ |
the different systems under investigation). An example of this |
| 1014 |
+ |
agreement for the SSDQ water/ion system is shown in |
| 1015 |
+ |
Fig. \ref{fig:gofr}. |
| 1016 |
|
|
| 1017 |
+ |
\begin{figure} |
| 1018 |
+ |
\centering |
| 1019 |
+ |
\includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
| 1020 |
+ |
\label{fig:gofr} |
| 1021 |
+ |
\caption{The pair distribution functions, $g(r)$, for the SSDQ |
| 1022 |
+ |
water/ion system obtained using the different real-space methods are |
| 1023 |
+ |
essentially identical with the result from the Ewald |
| 1024 |
+ |
treatment.} |
| 1025 |
+ |
\end{figure} |
| 1026 |
+ |
|
| 1027 |
+ |
There is a very slight overstructuring of the first solvation shell |
| 1028 |
+ |
when using when using TSF at lower values of the damping coefficient |
| 1029 |
+ |
($\alpha \le 0.1$) or when using undamped GSF. With moderate damping, |
| 1030 |
+ |
GSF and SP produce pair distributions that are identical (within |
| 1031 |
+ |
numerical noise) to their Ewald counterparts. |
| 1032 |
+ |
|
| 1033 |
+ |
A structural property that is a more demanding test of modified |
| 1034 |
+ |
electrostatics is the mean value of the electrostatic energy $\langle |
| 1035 |
+ |
U_\mathrm{elect} \rangle / N$ which is obtained by sampling the |
| 1036 |
+ |
liquid-state configurations experienced by a liquid evolving entirely |
| 1037 |
+ |
under the influence of each of the methods. In table \ref{tab:Props} |
| 1038 |
+ |
we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with |
| 1039 |
+ |
the damping parameter, $\alpha$, for each of the methods. |
| 1040 |
+ |
|
| 1041 |
+ |
As in the crystals studied in the first paper, damping is important |
| 1042 |
+ |
for converging the mean electrostatic energy values, particularly for |
| 1043 |
+ |
the two shifted force methods (GSF and TSF). A value of $\alpha |
| 1044 |
+ |
\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF |
| 1045 |
+ |
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
| 1046 |
+ |
dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA). |
| 1047 |
+ |
Overdamping the real-space electrostatic methods occurs with $\alpha > |
| 1048 |
+ |
0.4$, causing the estimate of the energy to drop below the Ewald |
| 1049 |
+ |
results. |
| 1050 |
+ |
|
| 1051 |
+ |
These ``optimal'' values of the damping coefficient are slightly |
| 1052 |
+ |
larger than what were observed for DSF electrostatics for purely |
| 1053 |
+ |
point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for |
| 1054 |
+ |
$r_c = 12$\AA appears to be an excellent compromise for mixed charge |
| 1055 |
+ |
multipole systems. |
| 1056 |
+ |
|
| 1057 |
+ |
To test the fidelity of the electrostatic methods at reproducing |
| 1058 |
+ |
dynamics in a multipolar liquid, it is also useful to look at |
| 1059 |
+ |
transport properties, particularly the diffusion constant, |
| 1060 |
+ |
\begin{equation} |
| 1061 |
+ |
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
| 1062 |
+ |
\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle |
| 1063 |
+ |
\label{eq:diff} |
| 1064 |
+ |
\end{equation} |
| 1065 |
+ |
which measures long-time behavior and is sensitive to the forces on |
| 1066 |
+ |
the multipoles. For the soft dipolar fluid and the SSDQ liquid |
| 1067 |
+ |
systems, the self-diffusion constants (D) were calculated from linear |
| 1068 |
+ |
fits to the long-time portion of the mean square displacement, |
| 1069 |
+ |
$\langle r^{2}(t) \rangle$.\cite{Allen87} |
| 1070 |
+ |
|
| 1071 |
+ |
In addition to translational diffusion, orientational relaxation times |
| 1072 |
+ |
were calculated for comparisons with the Ewald simulations and with |
| 1073 |
+ |
experiments. These values were determined from the same 1~ns |
| 1074 |
+ |
microcanonical trajectories used for translational diffusion by |
| 1075 |
+ |
calculating the orientational time correlation function, |
| 1076 |
+ |
\begin{equation} |
| 1077 |
+ |
C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t) |
| 1078 |
+ |
\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle, |
| 1079 |
+ |
\label{eq:OrientCorr} |
| 1080 |
+ |
\end{equation} |
| 1081 |
+ |
where $P_l$ is the Legendre polynomial of order $l$ and |
| 1082 |
+ |
$\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis |
| 1083 |
+ |
$\gamma$ on a molecule.. Th body-fixed reference frame used for our |
| 1084 |
+ |
models has the $z$-axis running along the dipoles, and for the SSDQ |
| 1085 |
+ |
water model, the $y$-axis connects the two implied hydrogen atom |
| 1086 |
+ |
positions. From the orientation autocorrelation functions, we can |
| 1087 |
+ |
obtain time constants for rotational relaxation either by fitting an |
| 1088 |
+ |
exponential function or by integrating the entire correlation |
| 1089 |
+ |
function. In a good water model, these decay times would be |
| 1090 |
+ |
comparable to water orientational relaxation times from nuclear |
| 1091 |
+ |
magnetic resonance (NMR). The relaxation constant obtained from |
| 1092 |
+ |
$C_2^y(t)$ is normally of experimental interest because it describes |
| 1093 |
+ |
the relaxation of the principle axis connecting the hydrogen |
| 1094 |
+ |
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
| 1095 |
+ |
of the dipole-dipole relaxation from a proton NMR signal and should |
| 1096 |
+ |
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
| 1097 |
+ |
|
| 1098 |
+ |
Results for the diffusion constants and orientational relaxation times |
| 1099 |
+ |
are shown in figure \ref{tab:Props}. From this data, it is apparent |
| 1100 |
+ |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
| 1101 |
+ |
reproduced with reasonable fidelity by the GSF method. |
| 1102 |
+ |
|
| 1103 |
+ |
The $\tau_2$ results in \ref{tab:Props} show a much greater difference |
| 1104 |
+ |
between the real-space and the Ewald results. |
| 1105 |
+ |
|
| 1106 |
+ |
\begin{table} |
| 1107 |
+ |
\label{tab:Props} |
| 1108 |
+ |
\caption{Comparison of the structural and dynamic properties for the |
| 1109 |
+ |
soft dipolar liquid test for all of the real-space methods.} |
| 1110 |
+ |
\begin{tabular}{l|c|cccc|cccc|cccc} |
| 1111 |
+ |
& Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\ |
| 1112 |
+ |
$\alpha$ (\AA$^{-1}$) & & |
| 1113 |
+ |
0.0 & 0.1 & 0.2 & 0.3 & |
| 1114 |
+ |
0.0 & 0.1 & 0.2 & 0.3 & |
| 1115 |
+ |
0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14} |
| 1116 |
+ |
|
| 1117 |
+ |
$\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\ |
| 1118 |
+ |
D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)& |
| 1119 |
+ |
470.2(6) & |
| 1120 |
+ |
416.6(5) & |
| 1121 |
+ |
379.6(5) & |
| 1122 |
+ |
438.6(5) & |
| 1123 |
+ |
476.0(6) & |
| 1124 |
+ |
412.8(5) & |
| 1125 |
+ |
421.1(5) & |
| 1126 |
+ |
400.5(5) & |
| 1127 |
+ |
437.5(6) & |
| 1128 |
+ |
434.6(5) & |
| 1129 |
+ |
411.4(5) & |
| 1130 |
+ |
545.3(7) & |
| 1131 |
+ |
459.6(6) \\ |
| 1132 |
+ |
$\tau_2$ (fs) & |
| 1133 |
+ |
1.136 & |
| 1134 |
+ |
1.041 & |
| 1135 |
+ |
1.064 & |
| 1136 |
+ |
1.109 & |
| 1137 |
+ |
1.211 & |
| 1138 |
+ |
1.119 & |
| 1139 |
+ |
1.039 & |
| 1140 |
+ |
1.058 & |
| 1141 |
+ |
1.21 & |
| 1142 |
+ |
1.15 & |
| 1143 |
+ |
1.172 & |
| 1144 |
+ |
1.153 & |
| 1145 |
+ |
1.125 \\ |
| 1146 |
+ |
\end{tabular} |
| 1147 |
+ |
\end{table} |
| 1148 |
+ |
|
| 1149 |
+ |
|
| 1150 |
|
\section{CONCLUSION} |
| 1151 |
|
In the first paper in this series, we generalized the |
| 1152 |
|
charge-neutralized electrostatic energy originally developed by Wolf |
