| 34 |
|
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math |
| 35 |
|
%\linenumbers\relax % Commence numbering lines |
| 36 |
|
\usepackage{amsmath} |
| 37 |
+ |
\usepackage{times} |
| 38 |
+ |
\usepackage{mathptm} |
| 39 |
|
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 40 |
|
\usepackage{url} |
| 41 |
|
\usepackage[english]{babel} |
| 92 |
|
scaling, making it prohibitive for large systems. By utilizing |
| 93 |
|
particle meshes and three dimensional fast Fourier transforms (FFT), |
| 94 |
|
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 95 |
< |
($P^3ME$), and smooth particle mesh Ewald (SPME) methods can decrease |
| 95 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
| 96 |
|
the computational cost from $O(N^2)$ down to $O(N \log |
| 97 |
|
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
| 98 |
|
|
| 110 |
|
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
| 111 |
|
|
| 112 |
|
\subsection{Real-space methods} |
| 113 |
< |
Recently, Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space |
| 114 |
< |
$O(N)$ method for calculating electrostatic interactions between point |
| 113 |
> |
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
| 114 |
> |
method for calculating electrostatic interactions between point |
| 115 |
|
charges. They argued that the effective Coulomb interaction in |
| 116 |
|
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
| 117 |
|
an ordered lattice (e.g. when computing the Madelung constant of an |
| 137 |
|
\label{fig:NaCl} |
| 138 |
|
\end{figure} |
| 139 |
|
|
| 138 |
– |
|
| 140 |
|
The direct truncation of interactions at a cutoff radius creates |
| 141 |
|
truncation defects. Wolf \textit{et al.} further argued that |
| 142 |
|
truncation errors are due to net charge remaining inside the cutoff |
| 151 |
|
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
| 152 |
|
Fennel and Gezelter later proposed shifted force variants of the Wolf |
| 153 |
|
method with commensurate force and energy expressions that do not |
| 154 |
< |
exhibit this problem.\cite{Fennell:2006lq} |
| 154 |
> |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
| 155 |
> |
methods were also proposed by Chen \textit{et |
| 156 |
> |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
| 157 |
> |
and by Wu and Brooks.\cite{Wu:044107} |
| 158 |
|
|
| 159 |
|
Considering the interaction of one central ion in an ionic crystal |
| 160 |
|
with a portion of the crystal at some distance, the effective Columbic |
| 188 |
|
unstable. |
| 189 |
|
|
| 190 |
|
In ionic crystals, real-space truncation can break the effective |
| 191 |
< |
octupolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
| 191 |
> |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
| 192 |
|
swings in the electrostatic energy as the cutoff radius is increased |
| 193 |
|
(or as individual ions move back and forth across the boundary). This |
| 194 |
|
is why the image charges were necessary for the Wolf sum to exhibit |
| 195 |
|
rapid convergence. Similarly, the real-space truncation of point |
| 196 |
< |
multipole interactions breaks dipole-octupole arrangements, and image |
| 197 |
< |
multipoles are required for real-space treatments of electrostatic |
| 198 |
< |
energies. |
| 196 |
> |
multipole interactions breaks higher order multipole arrangements, and |
| 197 |
> |
image multipoles are required for real-space treatments of |
| 198 |
> |
electrostatic energies. |
| 199 |
|
|
| 200 |
|
% Because of this reason, although the nature of electrostatic |
| 201 |
|
% interaction short ranged, the hard cutoff sphere creates very large |
| 211 |
|
factors driving dynamics in molecular simulations. Fennell and |
| 212 |
|
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 213 |
|
obtain consistent energies and forces on the atoms within the cutoff |
| 214 |
< |
sphere which both smoothly go to zero as an atom aproaches the cutoff |
| 215 |
< |
radius. Also, the comparisons of the accuracy of the potential energy |
| 216 |
< |
and force between the DSF method and SPME was surprisingly |
| 217 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen wide use in to |
| 218 |
< |
calculated electrostatic interactions in molecular systems with |
| 219 |
< |
relatively uniform charge densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 214 |
> |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 215 |
> |
aproaches the cutoff radius. The comparisons of the accuracy these |
| 216 |
> |
quantities between the DSF kernel and SPME was surprisingly |
| 217 |
> |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
| 218 |
> |
calculating electrostatic interactions in molecular systems with |
| 219 |
> |
relatively uniform charge |
| 220 |
> |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 221 |
|
|
| 222 |
< |
\subsection{Damping function} |
| 222 |
> |
\subsection{The damping function} |
| 223 |
|
The damping function used in our research has been discussed in detail |
| 224 |
|
in the first paper of this series.\cite{PaperI} The radial kernel |
| 225 |
< |
$1/r$ for the interactions between point charges is replaced by the |
| 226 |
< |
complementary error function $\erfc(\alpha r)/r$ to accelerate the |
| 227 |
< |
rate of convergence, where $\alpha$ is damping parameter with units of |
| 228 |
< |
inverse distance. Altering the value of $\alpha$ is equivalent to |
| 229 |
< |
changing the width of the small Gaussian charge distributions that are |
| 230 |
< |
replacing each point charge -- Gaussian overlap integrals yield |
| 231 |
< |
complementary error functions when truncated at a finite distance. |
| 225 |
> |
$1/r$ for the interactions between point charges can be replaced by |
| 226 |
> |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
| 227 |
> |
accelerate the rate of convergence, where $\alpha$ is a damping |
| 228 |
> |
parameter with units of inverse distance. Altering the value of |
| 229 |
> |
$\alpha$ is equivalent to changing the width of Gaussian charge |
| 230 |
> |
distributions that replace each point charge -- Gaussian overlap |
| 231 |
> |
integrals yield complementary error functions when truncated at a |
| 232 |
> |
finite distance. |
| 233 |
|
|
| 234 |
< |
e can perform necessary mathematical manipulation |
| 235 |
< |
by varying $\alpha$ in the damping function for the calculation of the |
| 236 |
< |
electrostatic energy, force and torque\cite{Wolf:1999dn}. By using |
| 237 |
< |
suitable value of damping alpha ($\alpha = 0.2$) for a cutoff radius |
| 232 |
< |
($r_{c}=9 A$), \textit{Fennel and Gezelter}\cite{Fennell:2006lq} |
| 233 |
< |
produced very good agreement of the interaction energies, forces and |
| 234 |
< |
torques for charge-charge interactions.\cite{Fennell:2006lq} |
| 234 |
> |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
| 235 |
> |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
| 236 |
> |
agreement with SPME for the interaction energies, forces and torques |
| 237 |
> |
for charge-charge interactions.\cite{Fennell:2006lq} |
| 238 |
|
|
| 239 |
< |
\subsection{Point multipoles for CG modeling} |
| 240 |
< |
Since a molecule consists of equal positive and negative charges, instead taking of the most common case of atomic site-site interaction, the interaction between higher order multipoles can also be used to evaluate molecule-molecule interactions. The short-ranged interaction between the molecules is dominated by Lennard-Jones repulsion. Also, electrons in a molecule is not localized at a specific point, thus a molecule can be coarse-grained to approximate as point multipole.\cite{Ren06, Essex10, Essex11}Recently, water has been modeled with point multipoles up to octupolar order.\cite{Ichiye10_1, Ichiye10_2, Ichiye10_3}. The point multipoles method has also been used in the AMOEBA water model.\cite{Ponder:2010vl, Gordon07,Smith80}. But using point multipole in the real space cutoff method without account of multipolar neutrality creates problem in the total energy conservation in MD simulations. In this paper we extended the original idea of the charge neutrality by Wolf’s into point dipoles and quadrupoles. Also, we used the previously developed idea of the damped shifted potential (DSF) for the charge-charge interaction\cite{Fennell:2006lq}and generalized it into higher order multipoles to conserve the total energy in the molecular dynamic simulation (The detail mathematical development of the purposed methods have been discussed in paper I). |
| 239 |
> |
\subsection{Point multipoles in molecular modeling} |
| 240 |
> |
Coarse-graining approaches which treat entire molecular subsystems as |
| 241 |
> |
a single rigid body are now widely used. A common feature of many |
| 242 |
> |
coarse-graining approaches is simplification of the electrostatic |
| 243 |
> |
interactions between bodies so that fewer site-site interactions are |
| 244 |
> |
required to compute configurational energies. Many coarse-grained |
| 245 |
> |
molecular structures would normally consist of equal positive and |
| 246 |
> |
negative charges, and rather than use multiple site-site interactions, |
| 247 |
> |
the interaction between higher order multipoles can also be used to |
| 248 |
> |
evaluate a single molecule-molecule |
| 249 |
> |
interaction.\cite{Ren06,Essex10,Essex11} |
| 250 |
|
|
| 251 |
+ |
Because electrons in a molecule are not localized at specific points, |
| 252 |
+ |
the assignment of partial charges to atomic centers is a relatively |
| 253 |
+ |
rough approximation. Atomic sites can also be assigned point |
| 254 |
+ |
multipoles and polarizabilities to increase the accuracy of the |
| 255 |
+ |
molecular model. Recently, water has been modeled with point |
| 256 |
+ |
multipoles up to octupolar |
| 257 |
+ |
order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 258 |
+ |
multipoles up to quadrupolar order have also been coupled with point |
| 259 |
+ |
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
| 260 |
+ |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}. But |
| 261 |
+ |
using point multipole with the real space truncation without |
| 262 |
+ |
accounting for multipolar neutrality will create energy conservation |
| 263 |
+ |
issues in molecular dynamics (MD) simulations. |
| 264 |
|
|
| 265 |
+ |
In this paper we test a set of real-space methods that were developed |
| 266 |
+ |
for point multipolar interactions. These methods extend the damped |
| 267 |
+ |
shifted force (DSF) and Wolf methods originally developed for |
| 268 |
+ |
charge-charge interactions and generalize them for higher order |
| 269 |
+ |
multipoles. The detailed mathematical development of these methods has |
| 270 |
+ |
been presented in the first paper in this series, while this work |
| 271 |
+ |
covers the testing the energies, forces, torques, and energy |
| 272 |
+ |
conservation properties of the methods in realistic simulation |
| 273 |
+ |
environments. In all cases, the methods are compared with the |
| 274 |
+ |
reference method, a full multipolar Ewald treatment. |
| 275 |
+ |
|
| 276 |
+ |
|
| 277 |
|
%\subsection{Conservation of total energy } |
| 278 |
|
%To conserve the total energy in MD simulations, the energy, force, and torque on a central molecule due to another molecule should smoothly tends to zero as second molecule approaches to cutoff radius. In addition, the force should be derivable from the energy and vice versa. If only the first condition holds but not the second, the total energy does not conserve.\cite{Fennell:2006lq}. The hard cutoff method does not ensure the smooth transition of the energy, force, and torque at the cutoff radius.\cite{Wolf:1999dn} By placing image charge on the surface of the cutoff sphere, the smooth transition of the energy can be ensured but the force and torque remains discontinuous. Therefore the purposed methods should have smooth transition of the energy, force and torque to ensure the total energy conservation and the expression should be close to idea of placing image multipole on the surface of the cutoff sphere. |
| 279 |
|
|
| 280 |
< |
\section{\label{sec:method}REVIEW OF METHODS} |
| 281 |
< |
Any force field associated with MD simulation should have the electrostatic energy, force and the torque between central molecule and any other molecule within cutoff radius should smoothly approach to zero as $r$ tends to $r_c$. This issue of continuous nature of the electrostatic interaction at the cutoff radius is associated with the conservation of total energy in the MD simulation. The mathematical detail for the SP, GSF and TSF has already been discussed in detail in previous paper I.\cite{PaperI} |
| 280 |
> |
\section{\label{sec:method}Review of Methods} |
| 281 |
> |
Any real-space electrostatic method that is suitable for MD |
| 282 |
> |
simulations should have the electrostatic energy, forces and torques |
| 283 |
> |
between two sites go smoothly to zero as the distance between the |
| 284 |
> |
sites, $r_{\bf ab}$ approaches the cutoff radius, $r_c$. Requiring |
| 285 |
> |
this continuity at the cutoff is essential for energy conservation in |
| 286 |
> |
MD simulations. The mathematical details of the shifted potential |
| 287 |
> |
(SP), gradient-shifted-force (GSF) and Taylor shifted-force (TSF) |
| 288 |
> |
methods have been discussed in detail in the previous paper in this |
| 289 |
> |
series.\cite{PaperI} Here we briefly review the new methods and |
| 290 |
> |
describe their essential features. |
| 291 |
|
|
| 292 |
< |
\subsection{Taylor-shifted force(TSF)} |
| 247 |
< |
The detail mathematical expression for the multipole-multipole interaction by the TSF method has been described in paper I.\cite{PaperI}. The electrostatic potential energy between groups of charges or multipoles is expressed as the product of operator and potential due to point charge as shown in \textit{equation 4 in Paper I}.\cite{PaperI} In the Taylor Shifted Force (TSF) method, we shifted kernel $1/r$ (the potential due to a point charge) by $1/r_c$ and performed Taylor Series expansion of the shifted part about the cutoff radius before operating with the operators. To ensure smooth convergence of the energy, force, and torque to zero at the cut off radius, the required number of terms from Taylor Series expansion are performed for different multipole-multipole interactions. Also, the mathematical consistency between the energy, force and the torque has been established. The potential energy for the multipole-multipole interaction is given by, |
| 292 |
> |
\subsection{Taylor-shifted force (TSF)} |
| 293 |
|
|
| 294 |
+ |
The electrostatic potential energy between point multipoles can be |
| 295 |
+ |
expressed as the product of two multipole operators and a Coulombic |
| 296 |
+ |
kernel, |
| 297 |
|
\begin{equation} |
| 298 |
< |
\begin{split} |
| 251 |
< |
U_{TSF}(\vec r)=\sum_{\alpha=1}^3\sum_{\beta=1}^3(C_a - D_{a \alpha }\frac{\partial}{\partial r_{a \alpha}}+Q_{a \alpha \beta }\frac{\partial}{\partial r_{a \alpha}\partial r_{a \beta}})\\ |
| 252 |
< |
(C_b - D_{b \alpha }\frac{\partial}{\partial r_{b \alpha}}+Q_{b \alpha \beta }\frac{\partial}{\partial r_{b \alpha}\partial r_{b \beta}})\\ |
| 253 |
< |
[(\frac{1}{r}-[\frac{1}{r_c}-(r-r_c)\frac{1}{r_c^2}+(r-r_c)^2\frac{1}{r_c^3}+...)] |
| 254 |
< |
\end{split} |
| 255 |
< |
\label{eq:TSF} |
| 298 |
> |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
| 299 |
|
\end{equation} |
| 300 |
< |
|
| 258 |
< |
where $C_a = \sum_{k\;in\; a}q_k$ , $D_{a\alpha}=\sum_{k \;in\;a}q_k r_k\alpha$, and $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k r_{k\alpha}r_{k\beta}$ stand for charge, dipole and quadrupole moment respectively (detail in paperI\cite{PaperI}). The electrostatic force and torque acting on the central molecule due to a molecule within cutoff sphere are derived from the equation ~\ref{eq:TSF} with the account of appropriate number of terms. This method is developed on the basis of using kernel potential due to the point charge ($1/r$) and their image charge potential ($1/r_c$) with its Taylor series expansion and considering that the expression for multipole-multipole interaction can be obtained operating the modified kernel by their corresponding operators. |
| 259 |
< |
|
| 260 |
< |
\subsection{Shifted potential (SP) } |
| 261 |
< |
A discontinuous truncation of the electrostatic potential at the cutoff sphere introduces severe artifact(Oscillation in the electrostatic energy) even for molecules with the higher-order multipoles.\cite{Paper I} This artifact is due to the existence of multipole moments within the cutoff spheres contributed by the breaking of the multipole ordering at the the surface of the cutoff sphere. The multipole moments of the cutoff sphere can be neutralized by placing image multipole for every multipole within the cutoff sphere. The electrostatic potential between multipoles for the SP method is given by, |
| 300 |
> |
where the multipole operator for site $\bf a$, |
| 301 |
|
\begin{equation} |
| 302 |
< |
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) |
| 303 |
< |
\label{eq:SP} |
| 304 |
< |
\end{equation} |
| 305 |
< |
The SP method compensates the artifact created by truncation of the multipole ordering by placing image on the cutoff surface. Also, the potential energy between central multipole and other multipole within sphere approaches smoothly to zero as $r$ tends to $r_c$. But the force and torque obtained from the shifted potential are discontinuous at $r_c$. Therefore, the MD simulation will still have the total energy drift for a longer simulation. If we derive the force and torque from the direct shifting about $r_c$ like in shifted potential then inconsistency between the force, torque, and potential fails the energy conservation in the dynamic simulation. |
| 302 |
> |
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} |
| 303 |
> |
+ Q_{{\bf a}\alpha\beta} |
| 304 |
> |
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
| 305 |
> |
\end{equation} |
| 306 |
> |
is expressed in terms of the point charge, $C_{\bf a}$, dipole, |
| 307 |
> |
$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for |
| 308 |
> |
object $\bf a$. Note that in this work, we use the primitive |
| 309 |
> |
quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k |
| 310 |
> |
r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site. |
| 311 |
|
|
| 312 |
< |
\subsection{Gradient-shifted force (GSF)} |
| 313 |
< |
As we mentioned earlier, in the MD simulation the electrostatic energy, force and torque should approach to zero as r tends to $r_c$. Also, the energy, force and torque should be consistent with each other for the total energy conservation. The GSF method is developed to address both the issues of consistency and convergence of the energy, force and the torque. Furthermore, the compensating of charge or multipole ordering breakage in the SP method due to direct spherical truncation will remain intact for large $r_c$. The electrostatic potential energy between central molecule and any molecule inside cutoff radius is given by, |
| 314 |
< |
\begin{equation} |
| 315 |
< |
U_{SF}(\vec r)=\sum U(\vec r) - U(\vec r_c)-(\vec r-\vec r_c)\cdot\vec \nabla U(\vec r)|_{r=r_c} |
| 316 |
< |
\label{eq:GSF} |
| 317 |
< |
\end{equation} |
| 318 |
< |
where the third term converges more rapidly as compared to first two terms hence the contribution of the third term is very small for large $r_c$ value. Hence the GSF method similar to SP method for large $r_c$. Moreover, the force and torque derived from equation 3 are consistent with the energy and approaches to zero as $r$ tends to $r_c$. |
| 319 |
< |
Both GSF and TSF methods are the generalization of the original DSF method to higher order multipole-multipole interactions. These two methods are same up to charge-dipole interaction level but generate different expressions in the energy, force and torque for the higher order multipole-multipole interactions. |
| 320 |
< |
\subsection{Self term} |
| 312 |
> |
Interactions between multipoles can be expressed as higher derivatives |
| 313 |
> |
of the bare Coulomb potential, so one way of ensuring that the forces |
| 314 |
> |
and torques vanish at the cutoff distance is to include a larger |
| 315 |
> |
number of terms in the truncated Taylor expansion, e.g., |
| 316 |
> |
% |
| 317 |
> |
\begin{equation} |
| 318 |
> |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-R_c)^m}{m!} f^{(m)} \Big \lvert _{R_c} . |
| 319 |
> |
\end{equation} |
| 320 |
> |
% |
| 321 |
> |
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 322 |
> |
Thus, for $f(r)=1/r$, we find |
| 323 |
> |
% |
| 324 |
> |
\begin{equation} |
| 325 |
> |
f_1(r)=\frac{1}{r}- \frac{1}{R_c} + (r - R_c) \frac{1}{R_c^2} - \frac{(r-R_c)^2}{R_c^3} . |
| 326 |
> |
\end{equation} |
| 327 |
> |
This function is an approximate electrostatic potential that has |
| 328 |
> |
vanishing second derivatives at the cutoff radius, making it suitable |
| 329 |
> |
for shifting the forces and torques of charge-dipole interactions. |
| 330 |
|
|
| 331 |
< |
\section{\label{sec:test}Test systems} |
| 332 |
< |
We have compared the electrostatic force and torque of each molecule from SP, TSF and GSF method with the multipolar-Ewald method. Furthermore, total electrostatic energies of a molecular system from the different methods have also been compared with total energy from the Ewald. In Mote Carlo (MC) simulation, the energy difference between different configurations of the molecular system is important, even though absolute energies are not accurate. We have compared the change in electrostatic potential energy ($\triangle E$) of 250 different configurations of the various multipolar molecular systems (Section IV B) calculated from the Hard, SP, GSF, and TSF methods with the well-known Ewald method. In MD simulations, the force and torque acting on the molecules drives the whole dynamics of the molecules in a system. The magnitudes of the electrostatic force, torque and their direction for each molecule of the all 250 configurations have also been compared against the Ewald’s method. |
| 331 |
> |
In general, the TSF potential for any multipole-multipole interaction |
| 332 |
> |
can be written |
| 333 |
> |
\begin{equation} |
| 334 |
> |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
| 335 |
> |
\label{generic} |
| 336 |
> |
\end{equation} |
| 337 |
> |
with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for |
| 338 |
> |
charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and |
| 339 |
> |
$n=4$ for quadrupole-quadrupole. To ensure smooth convergence of the |
| 340 |
> |
energy, force, and torques, the required number of terms from Taylor |
| 341 |
> |
series expansion in $f_n(r)$ must be performed for different |
| 342 |
> |
multipole-multipole interactions. |
| 343 |
|
|
| 344 |
< |
\subsection{Modeled systems} |
| 345 |
< |
We studied the comparison of the energy differences, forces and torques for six different systems; i) dipolar liquid, ii) quadrupolar liquid, iii) dipolar crystal, iv) quadrupolar crystal v) dipolar-quadrupolar liquid(SSDQ), and vi) ions in dipolar-qudrupolar liquid(SSDQC). To simulate different configurations of the crystals, the body centered cubic (BCC) minimum energy crystal with 3,456 molecules was taken and translationally locked in their respective crystal sites. The thermal energy was supplied to the rotational motion so that dipoles or quadrupoles can freely explore all possible orientation. The crystals were simulated for 10,000 fs in NVE ensemble at 50 K and 250 different configurations was taken in equal time interval for the comparative study. The crystals were not simulated at high temperature and for a long run time to avoid possible translational deformation of the crystal sites. |
| 346 |
< |
For dipolar, quadrupolar, and dipolar-quadrupolar liquids simulation, each molecular system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K temperature after equilibration. We collected 250 different configurations in equal interval of time. For the ions mixed liquid system, we converted 48 different molecules into 24 $Na^+$ and $24 Cl^-$ ions and equilibrated. After equilibration, the system was run at the same environment for 1ns and 250 configurations were collected. While comparing energies, forces, and torques with Ewald method, Lennad Jone’s potentials were turned off and purely electrostatic interaction had been compared. |
| 344 |
> |
To carry out the same procedure for a damped electrostatic kernel, we |
| 345 |
> |
replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
| 346 |
> |
Many of the derivatives of the damped kernel are well known from |
| 347 |
> |
Smith's early work on multipoles for the Ewald |
| 348 |
> |
summation.\cite{Smith82,Smith98} |
| 349 |
|
|
| 350 |
< |
\subsection{Statistical analysis} |
| 351 |
< |
We have used least square regression analyses for six different molecular systems to compare $\triangle E$ from Hard, SP, GSF, and TSF with the reference method. Molecular systems were run longer enough to explore various configurations and 250 independent configurations were recorded for comparison. The total numbers of 31,125 energy differences from the proposed methods have been compared with the Ewald. Similarly, the magnitudes of the forces and torques have also been compared by using least square regression analyses. In the forces and torques comparison, the magnitudes of the forces acting in each molecule for each configuration were evaluated. For example, our dipolar liquid simulation contains 2048 molecules and there are 250 different configurations for each system thus there are 512,000 force and torque comparisons. The correlation coefficient and correlation slope varies from 0 to 1, where 1 is the best agreement between the two methods. |
| 350 |
> |
Note that increasing the value of $n$ will add additional terms to the |
| 351 |
> |
electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 352 |
> |
$(r-R_c)^3/R_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 353 |
> |
functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 354 |
> |
f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 355 |
> |
for computing multipole energies, forces, and torques, and smooth |
| 356 |
> |
cutoffs of these quantities can be guaranteed as long as the number of |
| 357 |
> |
terms in the Taylor series exceeds the derivative order required. |
| 358 |
|
|
| 359 |
< |
\subsection{Analysis of vector quantities} |
| 360 |
< |
R.A. Fisher has developed a probablity density function to analyse directional data sets is expressed as below,\cite{fisher53} |
| 359 |
> |
For multipole-multipole interactions, following this procedure results |
| 360 |
> |
in separate radial functions for each distinct orientational |
| 361 |
> |
contribution to the potential, and ensures that the forces and torques |
| 362 |
> |
from {\it each} of these contributions will vanish at the cutoff |
| 363 |
> |
radius. For example, the direct dipole dot product ($\mathbf{D}_{i} |
| 364 |
> |
\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance |
| 365 |
> |
dot products: |
| 366 |
|
\begin{equation} |
| 367 |
< |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta \exp(\kappa \cos\theta) |
| 368 |
< |
\label{eq:pdf} |
| 367 |
> |
U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot |
| 368 |
> |
\mathbf{D}_{j} \right) \frac{g_2(r)}{r} |
| 369 |
> |
-\frac{1}{4\pi \epsilon_0} |
| 370 |
> |
\left( \mathbf{D}_{i} \cdot \hat{r} \right) |
| 371 |
> |
\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) - |
| 372 |
> |
\frac{g_2(r)}{r} \right) |
| 373 |
|
\end{equation} |
| 374 |
< |
where $\kappa$ measures directional dispersion of the data about mean direction can be estimated as a reciprocal of the circular variance for large number of directional data sets.\cite{Allen91} In our calculation, the unit vector from the Ewald method was considered as mean direction and the angle between the vectors from Ewald and the purposed method were evaluated.The total displacement of the unit vectors from the purposed method was calculated as, |
| 374 |
> |
|
| 375 |
> |
The electrostatic forces and torques acting on the central multipole |
| 376 |
> |
site due to another site within cutoff sphere are derived from |
| 377 |
> |
Eq.~\ref{generic}, accounting for the appropriate number of |
| 378 |
> |
derivatives. Complete energy, force, and torque expressions are |
| 379 |
> |
presented in the first paper in this series (Reference |
| 380 |
> |
\citep{PaperI}). |
| 381 |
> |
|
| 382 |
> |
\subsection{Gradient-shifted force (GSF)} |
| 383 |
> |
|
| 384 |
> |
A second (and significantly simpler) method involves shifting the |
| 385 |
> |
gradient of the raw coulomb potential for each particular multipole |
| 386 |
> |
order. For example, the raw dipole-dipole potential energy may be |
| 387 |
> |
shifted smoothly by finding the gradient for two interacting dipoles |
| 388 |
> |
which have been projected onto the surface of the cutoff sphere |
| 389 |
> |
without changing their relative orientation, |
| 390 |
> |
\begin{displaymath} |
| 391 |
> |
U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(R_c) |
| 392 |
> |
- (r_{ij}-R_c) \hat{r}_{ij} \cdot |
| 393 |
> |
\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{R_c} |
| 394 |
> |
\end{displaymath} |
| 395 |
> |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$ |
| 396 |
> |
and $\mathbf{D}_{j}$, are retained at the cutoff distance (although |
| 397 |
> |
the signs are reversed for the dipole that has been projected onto the |
| 398 |
> |
cutoff sphere). In many ways, this simpler approach is closer in |
| 399 |
> |
spirit to the original shifted force method, in that it projects a |
| 400 |
> |
neutralizing multipole (and the resulting forces from this multipole) |
| 401 |
> |
onto a cutoff sphere. The resulting functional forms for the |
| 402 |
> |
potentials, forces, and torques turn out to be quite similar in form |
| 403 |
> |
to the Taylor-shifted approach, although the radial contributions are |
| 404 |
> |
significantly less perturbed by the Gradient-shifted approach than |
| 405 |
> |
they are in the Taylor-shifted method. |
| 406 |
> |
|
| 407 |
> |
In general, the gradient shifted potential between a central multipole |
| 408 |
> |
and any multipolar site inside the cutoff radius is given by, |
| 409 |
|
\begin{equation} |
| 410 |
< |
R = \sqrt{(\sum\limits_{i=1}^N \sin\theta_i)^2 + (\sum\limits_{i=1}^N \sin\theta_i)^2} |
| 410 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 411 |
> |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
| 412 |
> |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
| 413 |
> |
\label{generic2} |
| 414 |
> |
\end{equation} |
| 415 |
> |
where the sum describes a separate force-shifting that is applied to |
| 416 |
> |
each orientational contribution to the energy. |
| 417 |
> |
|
| 418 |
> |
The third term converges more rapidly than the first two terms as a |
| 419 |
> |
function of radius, hence the contribution of the third term is very |
| 420 |
> |
small for large cutoff radii. The force and torque derived from |
| 421 |
> |
equation \ref{generic2} are consistent with the energy expression and |
| 422 |
> |
approach zero as $r \rightarrow R_c$. Both the GSF and TSF methods |
| 423 |
> |
can be considered generalizations of the original DSF method for |
| 424 |
> |
higher order multipole interactions. GSF and TSF are also identical up |
| 425 |
> |
to the charge-dipole interaction but generate different expressions in |
| 426 |
> |
the energy, force and torque for higher order multipole-multipole |
| 427 |
> |
interactions. Complete energy, force, and torque expressions for the |
| 428 |
> |
GSF potential are presented in the first paper in this series |
| 429 |
> |
(Reference \citep{PaperI}) |
| 430 |
> |
|
| 431 |
> |
|
| 432 |
> |
\subsection{Shifted potential (SP) } |
| 433 |
> |
A discontinuous truncation of the electrostatic potential at the |
| 434 |
> |
cutoff sphere introduces a severe artifact (oscillation in the |
| 435 |
> |
electrostatic energy) even for molecules with the higher-order |
| 436 |
> |
multipoles.\cite{PaperI} We have also formulated an extension of the |
| 437 |
> |
Wolf approach for point multipoles by simply projecting the image |
| 438 |
> |
multipole onto the surface of the cutoff sphere, and including the |
| 439 |
> |
interactions with the central multipole and the image. This |
| 440 |
> |
effectively shifts the total potential to zero at the cutoff radius, |
| 441 |
> |
\begin{equation} |
| 442 |
> |
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) |
| 443 |
> |
\label{eq:SP} |
| 444 |
> |
\end{equation} |
| 445 |
> |
where the sum describes separate potential shifting that is done for |
| 446 |
> |
each orientational contribution to the energy (e.g. the direct dipole |
| 447 |
> |
product contribution is shifted {\it separately} from the |
| 448 |
> |
dipole-distance terms in dipole-dipole interactions). Note that this |
| 449 |
> |
is not a simple shifting of the total potential at $R_c$. Each radial |
| 450 |
> |
contribution is shifted separately. One consequence of this is that |
| 451 |
> |
multipoles that reorient after leaving the cutoff sphere can re-enter |
| 452 |
> |
the cutoff sphere without perturbing the total energy. |
| 453 |
> |
|
| 454 |
> |
The potential energy between a central multipole and other multipolar |
| 455 |
> |
sites then goes smoothly to zero as $r \rightarrow R_c$. However, the |
| 456 |
> |
force and torque obtained from the shifted potential (SP) are |
| 457 |
> |
discontinuous at $R_c$. Therefore, MD simulations will still |
| 458 |
> |
experience energy drift while operating under the SP potential, but it |
| 459 |
> |
may be suitable for Monte Carlo approaches where the configurational |
| 460 |
> |
energy differences are the primary quantity of interest. |
| 461 |
> |
|
| 462 |
> |
\subsection{The Self term} |
| 463 |
> |
In the TSF, GSF, and SP methods, a self-interaction is retained for |
| 464 |
> |
the central multipole interacting with its own image on the surface of |
| 465 |
> |
the cutoff sphere. This self interaction is nearly identical with the |
| 466 |
> |
self-terms that arise in the Ewald sum for multipoles. Complete |
| 467 |
> |
expressions for the self terms are presented in the first paper in |
| 468 |
> |
this series (Reference \citep{PaperI}) |
| 469 |
> |
|
| 470 |
> |
|
| 471 |
> |
\section{\label{sec:methodology}Methodology} |
| 472 |
> |
|
| 473 |
> |
To understand how the real-space multipole methods behave in computer |
| 474 |
> |
simulations, it is vital to test against established methods for |
| 475 |
> |
computing electrostatic interactions in periodic systems, and to |
| 476 |
> |
evaluate the size and sources of any errors that arise from the |
| 477 |
> |
real-space cutoffs. In the first paper of this series, we compared |
| 478 |
> |
the dipolar and quadrupolar energy expressions against analytic |
| 479 |
> |
expressions for ordered dipolar and quadrupolar |
| 480 |
> |
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} This work uses the |
| 481 |
> |
multipolar Ewald sum as a reference method for comparing energies, |
| 482 |
> |
forces, and torques for molecular models that mimic disordered and |
| 483 |
> |
ordered condensed-phase systems. These test-cases include: |
| 484 |
> |
|
| 485 |
> |
\begin{itemize} |
| 486 |
> |
\item Soft Dipolar fluids ($\sigma = , \epsilon = , |D| = $) |
| 487 |
> |
\item Soft Dipolar solids ($\sigma = , \epsilon = , |D| = $) |
| 488 |
> |
\item Soft Quadrupolar fluids ($\sigma = , \epsilon = , Q_{xx} = ...$) |
| 489 |
> |
\item Soft Quadrupolar solids ($\sigma = , \epsilon = , Q_{xx} = ...$) |
| 490 |
> |
\item A mixed multipole model for water |
| 491 |
> |
\item A mixed multipole models for water with dissolved ions |
| 492 |
> |
\end{itemize} |
| 493 |
> |
This last test case exercises all levels of the multipole-multipole |
| 494 |
> |
interactions we have derived so far and represents the most complete |
| 495 |
> |
test of the new methods. |
| 496 |
> |
|
| 497 |
> |
In the following section, we present results for the total |
| 498 |
> |
electrostatic energy, as well as the electrostatic contributions to |
| 499 |
> |
the force and torque on each molecule. These quantities have been |
| 500 |
> |
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
| 501 |
> |
and have been compared with the values obtaine from the multipolar |
| 502 |
> |
Ewald sum. In Mote Carlo (MC) simulations, the energy differences |
| 503 |
> |
between two configurations is the primary quantity that governs how |
| 504 |
> |
the simulation proceeds. These differences are the most imporant |
| 505 |
> |
indicators of the reliability of a method even if the absolute |
| 506 |
> |
energies are not exact. For each of the multipolar systems listed |
| 507 |
> |
above, we have compared the change in electrostatic potential energy |
| 508 |
> |
($\Delta E$) between 250 statistically-independent configurations. In |
| 509 |
> |
molecular dynamics (MD) simulations, the forces and torques govern the |
| 510 |
> |
behavior of the simulation, so we also compute the electrostatic |
| 511 |
> |
contributions to the forces and torques. |
| 512 |
> |
|
| 513 |
> |
\subsection{Model systems} |
| 514 |
> |
To sample independent configurations of multipolar crystals, a body |
| 515 |
> |
centered cubic (BCC) crystal which is a minimum energy structure for |
| 516 |
> |
point dipoles was generated using 3,456 molecules. The multipoles |
| 517 |
> |
were translationally locked in their respective crystal sites for |
| 518 |
> |
equilibration at a relatively low temperature (50K), so that dipoles |
| 519 |
> |
or quadrupoles could freely explore all accessible orientations. The |
| 520 |
> |
translational constraints were removed, and the crystals were |
| 521 |
> |
simulated for 10 ps in the microcanonical (NVE) ensemble with an |
| 522 |
> |
average temperature of 50 K. Configurations were sampled at equal |
| 523 |
> |
time intervals for the comparison of the configurational energy |
| 524 |
> |
differences. The crystals were not simulated close to the melting |
| 525 |
> |
points in order to avoid translational deformation away of the ideal |
| 526 |
> |
lattice geometry. |
| 527 |
> |
|
| 528 |
> |
For dipolar, quadrupolar, and mixed-multipole liquid simulations, each |
| 529 |
> |
system was created with 2048 molecules oriented randomly. These were |
| 530 |
> |
|
| 531 |
> |
system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K |
| 532 |
> |
temperature after equilibration. We collected 250 different |
| 533 |
> |
configurations in equal interval of time. For the ions mixed liquid |
| 534 |
> |
system, we converted 48 different molecules into 24 \ce{Na+} and 24 |
| 535 |
> |
\ce{Cl-} ions and equilibrated. After equilibration, the system was run |
| 536 |
> |
at the same environment for 1ns and 250 configurations were |
| 537 |
> |
collected. While comparing energies, forces, and torques with Ewald |
| 538 |
> |
method, Lennard-Jones potentials were turned off and purely |
| 539 |
> |
electrostatic interaction had been compared. |
| 540 |
> |
|
| 541 |
> |
\subsection{Accuracy of Energy Differences, Forces and Torques} |
| 542 |
> |
The pairwise summation techniques (outlined above) were evaluated for |
| 543 |
> |
use in MC simulations by studying the energy differences between |
| 544 |
> |
different configurations. We took the Ewald-computed energy |
| 545 |
> |
difference between two conformations to be the correct behavior. An |
| 546 |
> |
ideal performance by one of the new methods would reproduce these |
| 547 |
> |
energy differences exactly. The configurational energies being used |
| 548 |
> |
here contain only contributions from electrostatic interactions. |
| 549 |
> |
Lennard-Jones interactions were omitted from the comparison as they |
| 550 |
> |
should be identical for all methods. |
| 551 |
> |
|
| 552 |
> |
Since none of the real-space methods provide exact energy differences, |
| 553 |
> |
we used least square regressions analysiss for the six different |
| 554 |
> |
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
| 555 |
> |
with the multipolar Ewald reference method. Unitary results for both |
| 556 |
> |
the correlation (slope) and correlation coefficient for these |
| 557 |
> |
regressions indicate perfect agreement between the real-space method |
| 558 |
> |
and the multipolar Ewald sum. |
| 559 |
> |
|
| 560 |
> |
Molecular systems were run long enough to explore independent |
| 561 |
> |
configurations and 250 configurations were recorded for comparison. |
| 562 |
> |
Each system provided 31,125 energy differences for a total of 186,750 |
| 563 |
> |
data points. Similarly, the magnitudes of the forces and torques have |
| 564 |
> |
also been compared by using least squares regression analyses. In the |
| 565 |
> |
forces and torques comparison, the magnitudes of the forces acting in |
| 566 |
> |
each molecule for each configuration were evaluated. For example, our |
| 567 |
> |
dipolar liquid simulation contains 2048 molecules and there are 250 |
| 568 |
> |
different configurations for each system resulting in 3,072,000 data |
| 569 |
> |
points for comparison of forces and torques. |
| 570 |
> |
|
| 571 |
> |
\subsection{Analysis of vector quantities} |
| 572 |
> |
Getting the magnitudes of the force and torque vectors correct is only |
| 573 |
> |
part of the issue for carrying out accurate molecular dynamics |
| 574 |
> |
simulations. Because the real space methods reweight the different |
| 575 |
> |
orientational contributions to the energies, it is also important to |
| 576 |
> |
understand how the methods impact the \textit{directionality} of the |
| 577 |
> |
force and torque vectors. Fisher developed a probablity density |
| 578 |
> |
function to analyse directional data sets, |
| 579 |
> |
\begin{equation} |
| 580 |
> |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta e^{\kappa \cos\theta} |
| 581 |
> |
\label{eq:pdf} |
| 582 |
> |
\end{equation} |
| 583 |
> |
where $\kappa$ measures directional dispersion of the data around the |
| 584 |
> |
mean direction.\cite{fisher53} This quantity $(\kappa)$ can be |
| 585 |
> |
estimated as a reciprocal of the circular variance.\cite{Allen91} To |
| 586 |
> |
quantify the directional error, forces obtained from the Ewald sum |
| 587 |
> |
were taken as the mean (or correct) direction and the angle between |
| 588 |
> |
the forces obtained via the Ewald sum and the real-space methods were |
| 589 |
> |
evaluated, |
| 590 |
> |
\begin{equation} |
| 591 |
> |
\cos\theta_i = \frac{\vec{f}_i^\mathrm{~Ewald} \cdot |
| 592 |
> |
\vec{f}_i^\mathrm{~GSF}}{\left|\vec{f}_i^\mathrm{~Ewald}\right| \left|\vec{f}_i^\mathrm{~GSF}\right|} |
| 593 |
> |
\end{equation} |
| 594 |
> |
The total angular displacement of the vectors was calculated as, |
| 595 |
> |
\begin{equation} |
| 596 |
> |
R = \sqrt{\left(\sum\limits_{i=1}^N \cos\theta_i\right)^2 + \left(\sum\limits_{i=1}^N \sin\theta_i\right)^2} |
| 597 |
|
\label{eq:displacement} |
| 598 |
|
\end{equation} |
| 599 |
< |
where N is number of directional data sets and $theta_i$ are the angles between unit vectors evaluated from the Ewald and the purposed methods. The circular variance is defined as $ Var(\theta) = 1 -R/N$. The value of circular variance varies from 0 to 1. The lower the value of $Var{\theta}$ is higher the value of $\kappa$, which expresses tighter clustering of the direction sets around Ewald direction. |
| 599 |
> |
where $N$ is number of force vectors. The circular variance is |
| 600 |
> |
defined as |
| 601 |
> |
\begin{equation} |
| 602 |
> |
\mathrm{Var}(\theta) \approx 1/\kappa = 1 - R/N |
| 603 |
> |
\end{equation} |
| 604 |
> |
The circular variance takes on values between from 0 to 1, with 0 |
| 605 |
> |
indicating a perfect directional match between the Ewald force vectors |
| 606 |
> |
and the real-space forces. Lower values of $\mathrm{Var}(\theta)$ |
| 607 |
> |
correspond to higher values of $\kappa$, which indicates tighter |
| 608 |
> |
clustering of the real-space force vectors around the Ewald forces. |
| 609 |
|
|
| 610 |
+ |
A similar analysis was carried out for the electrostatic contribution |
| 611 |
+ |
to the molecular torques as well as forces. |
| 612 |
+ |
|
| 613 |
|
\subsection{Energy conservation} |
| 614 |
< |
To test conservation of the energy, the mixed molecular system of 2000 dipolar-quadrupolar molecules with 24 $Na^+$, and 24 $Cl^-$ was run for 1ns in the microcanonical ensemble at 300 K temperature for different cutoff methods (Ewald, Hard, SP, GSF, and TSF). The molecular system was run in 12 parallel computers and started with same initial positions and velocities for all cutoff methods. The slope and Standard Deviation of the energy about the slope (SD) were evaluated in the total energy versus time plot, where the slope evaluates the total energy drift and SD calculates the energy fluctuation in MD simulations. Also, the time duration for the simulation was recorded to compare efficiency of the purposed methods with the Ewald. |
| 614 |
> |
To test conservation the energy for the methods, the mixed molecular |
| 615 |
> |
system of 2000 SSDQ water molecules with 24 \ce{Na+} and 24 \ce{Cl-} |
| 616 |
> |
ions was run for 1 ns in the microcanonical ensemble at an average |
| 617 |
> |
temperature of 300K. Each of the different electrostatic methods |
| 618 |
> |
(Ewald, Hard, SP, GSF, and TSF) was tested for a range of different |
| 619 |
> |
damping values. The molecular system was started with same initial |
| 620 |
> |
positions and velocities for all cutoff methods. The energy drift |
| 621 |
> |
($\delta E_1$) and standard deviation of the energy about the slope |
| 622 |
> |
($\delta E_0$) were evaluated from the total energy of the system as a |
| 623 |
> |
function of time. Although both measures are valuable at |
| 624 |
> |
investigating new methods for molecular dynamics, a useful interaction |
| 625 |
> |
model must allow for long simulation times with minimal energy drift. |
| 626 |
|
|
| 627 |
|
\section{\label{sec:result}RESULTS} |
| 628 |
|
\subsection{Configurational energy differences} |
| 695 |
|
\section{CONCLUSION} |
| 696 |
|
We have generalized the charged neutralized potential energy originally developed by the Wolf et al.\cite{Wolf:1999dn} for the charge-charge interaction to the charge-multipole and multipole-multipole interaction in the SP method for higher order multipoles. Also, we have developed GSF and TSF methods by implementing the modification purposed by Fennel and Gezelter\cite{Fennell:2006lq} for the charge-charge interaction to the higher order multipoles to ensure consistency and smooth truncation of the electrostatic energy, force, and torque for the spherical truncation. The SP methods for multipoles proved its suitability in MC simulations. On the other hand, the results from the GSF method produced good agreement with the Ewald's energy, force, and torque. Also, it shows very good energy conservation in MD simulations. |
| 697 |
|
The direct truncation of any molecular system without multipole neutralization creates the fluctuation in the electrostatic energy. This fluctuation in the energy is very large for the case of crystal because of long range of multipole ordering (Refer paper I).\cite{PaperI} This is also significant in the case of the liquid because of the local multipole ordering in the molecules. If the net multipole within cutoff radius neutralized within cutoff sphere by placing image multiples on the surface of the sphere, this fluctuation in the energy reduced significantly. Also, the multipole neutralization in the generalized SP method showed very good agreement with the Ewald as compared to direct truncation for the evaluation of the $\triangle E$ between the configurations. |
| 698 |
< |
In MD simulations, the energy conservation is very important. The conservation of the total energy can be ensured by i) enforcing the smooth truncation of the energy, force and torque in the cutoff radius and ii) making the energy, force and torque consistent with each other. The GSF and TSF methods ensure the consistency and smooth truncation of the energy, force and torque at the cutoff radius, as a result show very good total energy conservation. But the TSF method does not show good agreement in the absolute value of the electrostatic energy, force and torque with the Ewald. The GSF method has mimicked Ewald’s force, energy and torque accurately and also conserved energy. Therefore, the GSF method is the suitable method for evaluating required force field in MD simulations. In addition, the energy drift and fluctuation from the GSF method is much better than Ewald’s method for finite-sized reciprocal space. |
| 698 |
> |
In MD simulations, the energy conservation is very important. The |
| 699 |
> |
conservation of the total energy can be ensured by i) enforcing the |
| 700 |
> |
smooth truncation of the energy, force and torque in the cutoff radius |
| 701 |
> |
and ii) making the energy, force and torque consistent with each |
| 702 |
> |
other. The GSF and TSF methods ensure the consistency and smooth |
| 703 |
> |
truncation of the energy, force and torque at the cutoff radius, as a |
| 704 |
> |
result show very good total energy conservation. But the TSF method |
| 705 |
> |
does not show good agreement in the absolute value of the |
| 706 |
> |
electrostatic energy, force and torque with the Ewald. The GSF method |
| 707 |
> |
has mimicked Ewald’s force, energy and torque accurately and also |
| 708 |
> |
conserved energy. Therefore, the GSF method is the suitable method for |
| 709 |
> |
evaluating required force field in MD simulations. In addition, the |
| 710 |
> |
energy drift and fluctuation from the GSF method is much better than |
| 711 |
> |
Ewald’s method for finite-sized reciprocal space. |
| 712 |
> |
|
| 713 |
> |
Note that the TSF, GSF, and SP models are $\mathcal{O}(N)$ methods |
| 714 |
> |
that can be made extremely efficient using spline interpolations of |
| 715 |
> |
the radial functions. They require no Fourier transforms or $k$-space |
| 716 |
> |
sums, and guarantee the smooth handling of energies, forces, and |
| 717 |
> |
torques as multipoles cross the real-space cutoff boundary. |
| 718 |
> |
|
| 719 |
|
%\bibliographystyle{aip} |
| 720 |
+ |
\newpage |
| 721 |
|
\bibliography{references} |
| 722 |
|
\end{document} |
| 723 |
|
|
