multipole/multipole_2/multipole2.tex

% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,jcp,
 amsmath,amssymb,
preprint,
%reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
%\usepackage{bm}% bold math
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines
\usepackage{amsmath}
\usepackage{times}
\usepackage{mathptm}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}
\usepackage[english]{babel}


\begin{document}

\preprint{AIP/123-QED}

\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald
Sum. II. Comparison of Simulation Methodologies} % Force line breaks with \\

\author{Madan Lamichhane}
 \affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\

\author{Kathie E. Newman}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}

\author{J. Daniel Gezelter}%
 \email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash
}%

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
We have tested our recently developed shifted potential, gradient-shifted force, and Taylor-shifted force methods for the higher-order multipoles against Ewald’s method in different types of liquid and crystalline system. In this paper, we have also investigated the conservation of total energy in the molecular dynamic simulation using all of these methods. The shifted potential method shows better agreement with the Ewald in the energy differences between different configurations as compared to the direct truncation. Both the gradient shifted force and Taylor-shifted force methods reproduce very good energy conservation. But the absolute energy, force and torque evaluated from the gradient shifted force method shows better result as compared to taylor-shifted force method. Hence the gradient-shifted force method suitably mimics the electrostatic interaction in the molecular dynamic simulation.
\end{abstract}

\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                             % Classification Scheme.
\keywords{Electrostatics, Multipoles, Real-space}

\maketitle


\section{\label{sec:intro}Introduction}
Computing the interactions between electrostatic sites is one of the
most expensive aspects of molecular simulations, which is why there
have been significant efforts to develop practical, efficient and
convergent methods for handling these interactions. Ewald's method is
perhaps the best known and most accurate method for evaluating
energies, forces, and torques in explicitly-periodic simulation
cells. In this approach, the conditionally convergent electrostatic
energy is converted into two absolutely convergent contributions, one
which is carried out in real space with a cutoff radius, and one in
reciprocal space.\cite{Clarke:1986eu,Woodcock75}

When carried out as originally formulated, the reciprocal-space
portion of the Ewald sum exhibits relatively poor computational
scaling, making it prohibitive for large systems. By utilizing
particle meshes and three dimensional fast Fourier transforms (FFT),
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease
the computational cost from $O(N^2)$ down to $O(N \log
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}.

Because of the artificial periodicity required for the Ewald sum, the
method may require modification to compute interactions for
interfacial molecular systems such as membranes and liquid-vapor
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl}
To simulate interfacial systems, Parry’s extension of the 3D Ewald sum
is appropriate for slab geometries.\cite{Parry:1975if} The inherent
periodicity in the Ewald’s method can also be problematic for
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald
methods that were developed to handle two-dimensional (2D)
electrostatic interactions in interfacial systems have not had similar
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77,
  Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq}

\subsection{Real-space methods}
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$
method for calculating electrostatic interactions between point
charges. They argued that the effective Coulomb interaction in
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}.  For
an ordered lattice (e.g. when computing the Madelung constant of an
ionic solid), the material can be considered as a set of ions
interacting with neutral dipolar or quadrupolar ``molecules'' giving
an effective distance dependence for the electrostatic interactions of
$r^{-5}$ (see figure \ref{fig:NaCl}.  For this reason, careful
applications of Wolf's method are able to obtain accurate estimates of
Madelung constants using relatively short cutoff radii.  Recently,
Fukuda used neutralization of the higher order moments for the
calculation of the electrostatic interaction of the point charges
system.\cite{Fukuda:2013sf}

\begin{figure}[h!]
  \centering
  \includegraphics[width=0.50 \textwidth]{chargesystem.pdf}
  \caption{Top: NaCl crystal showing how spherical truncation can
    breaking effective charge ordering, and how complete \ce{(NaCl)4}
    molecules interact with the central ion.  Bottom: A dipolar
    crystal exhibiting similar behavior and illustrating how the
    effective dipole-octupole interactions can be disrupted by
    spherical truncation.}
  \label{fig:NaCl}
\end{figure}

The direct truncation of interactions at a cutoff radius creates
truncation defects. Wolf \textit{et al.} further argued that
truncation errors are due to net charge remaining inside the cutoff
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed
placing an image charge on the surface of the cutoff sphere for every
real charge inside the cutoff.  These charges are present for the
evaluation of both the pair interaction energy and the force, although
the force expression maintained a discontinuity at the cutoff sphere.
In the original Wolf formulation, the total energy for the charge and
image were not equal to the integral of their force expression, and as
a result, the total energy would not be conserved in molecular
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and
Fennel and Gezelter later proposed shifted force variants of the Wolf
method with commensurate force and energy expressions that do not
exhibit this problem.\cite{Fennell:2006lq}   Related real-space
methods were also proposed by Chen \textit{et
  al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw}
and by Wu and Brooks.\cite{Wu:044107}

Considering the interaction of one central ion in an ionic crystal
with a portion of the crystal at some distance, the effective Columbic
potential is found to be decreasing as $r^{-5}$. If one views the
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more
rapidly to the Madelung energy than the dipolar
approximation.\cite{Wolf92} To find the correct Madelung constant,
Lacman suggested that the NaCl structure could be constructed in a way
that the finite crystal terminates with complete \ce{(NaCl)4}
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded
by opposite charges. Similarly for each pair of charges, there is an
opposite pair of charge adjacent to it.  The central ion sees what is
effectively a set of octupoles at large distances. These facts suggest
that the Madelung constants are relatively short ranged for perfect
ionic crystals.\cite{Wolf:1999dn}

One can make a similar argument for crystals of point multipoles. The
Luttinger and Tisza treatment of energy constants for dipolar lattices
utilizes 24 basis vectors that contain dipoles at the eight corners of
a unit cube.  Only three of these basis vectors, $X_1, Y_1,
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have
zero net dipole and retain contributions only from higher order
multipoles.  The effective interaction between a dipole at the center
of a crystal and a group of eight dipoles farther away is
significantly shorter ranged than the $r^{-3}$ that one would expect
for raw dipole-dipole interactions.  Only in crystals which retain a
bulk dipole moment (e.g. ferroelectrics) does the analogy with the
ionic crystal break down -- ferroelectric dipolar crystals can exist,
while ionic crystals with net charge in each unit cell would be
unstable.

In ionic crystals, real-space truncation can break the effective
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant
swings in the electrostatic energy as the cutoff radius is increased
(or as individual ions move back and forth across the boundary).  This
is why the image charges were necessary for the Wolf sum to exhibit
rapid convergence.  Similarly, the real-space truncation of point
multipole interactions breaks higher order multipole arrangements, and
image multipoles are required for real-space treatments of
electrostatic energies.

% Because of this reason, although the nature of electrostatic
% interaction short ranged, the hard cutoff sphere creates very large
% fluctuation in the electrostatic energy for the perfect crystal. In
% addition, the charge neutralized potential proposed by Wolf et
% al. converged to correct Madelung constant but still holds oscillation
% in the energy about correct Madelung energy.\cite{Wolf:1999dn}.  This
% oscillation in the energy around its fully converged value can be due
% to the non-neutralized value of the higher order moments within the
% cutoff sphere.

The forces and torques acting on atomic sites are the fundamental
factors driving dynamics in molecular simulations. Fennell and
Gezelter proposed the damped shifted force (DSF) energy kernel to
obtain consistent energies and forces on the atoms within the cutoff
sphere. Both the energy and the force go smoothly to zero as an atom
aproaches the cutoff radius. The comparisons of the accuracy these
quantities between the DSF kernel and SPME was surprisingly
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for
calculating electrostatic interactions in molecular systems with
relatively uniform charge
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13}

\subsection{The damping function}
The damping function used in our research has been discussed in detail
in the first paper of this series.\cite{PaperI} The radial kernel
$1/r$ for the interactions between point charges can be replaced by
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to
accelerate the rate of convergence, where $\alpha$ is a damping
parameter with units of inverse distance.  Altering the value of
$\alpha$ is equivalent to changing the width of Gaussian charge
distributions that replace each point charge -- Gaussian overlap
integrals yield complementary error functions when truncated at a
finite distance.

By using suitable value of damping alpha ($\alpha \sim 0.2$) for a
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good
agreement with SPME for the interaction energies, forces and torques
for charge-charge interactions.\cite{Fennell:2006lq}

\subsection{Point multipoles in molecular modeling}
Coarse-graining approaches which treat entire molecular subsystems as
a single rigid body are now widely used. A common feature of many
coarse-graining approaches is simplification of the electrostatic
interactions between bodies so that fewer site-site interactions are
required to compute configurational energies.  Many coarse-grained
molecular structures would normally consist of equal positive and
negative charges, and rather than use multiple site-site interactions,
the interaction between higher order multipoles can also be used to
evaluate a single molecule-molecule
interaction.\cite{Ren06,Essex10,Essex11}

Because electrons in a molecule are not localized at specific points,
the assignment of partial charges to atomic centers is a relatively
rough approximation.  Atomic sites can also be assigned point
multipoles and polarizabilities to increase the accuracy of the
molecular model.  Recently, water has been modeled with point
multipoles up to octupolar
order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point
multipoles up to quadrupolar order have also been coupled with point
polarizabilities in the high-quality AMOEBA and iAMOEBA water
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}.  But
using point multipole with the real space truncation without
accounting for multipolar neutrality will create energy conservation
issues in molecular dynamics (MD) simulations.

In this paper we test a set of real-space methods that were developed
for point multipolar interactions.  These methods extend the damped
shifted force (DSF) and Wolf methods originally developed for
charge-charge interactions and generalize them for higher order
multipoles. The detailed mathematical development of these methods has
been presented in the first paper in this series, while this work
covers the testing the energies, forces, torques, and energy
conservation properties of the methods in realistic simulation
environments.  In all cases, the methods are compared with the
reference method, a full multipolar Ewald treatment.


%\subsection{Conservation of total energy } 
%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.

\section{\label{sec:method}Review of Methods}
Any real-space electrostatic method that is suitable for MD
simulations should have the electrostatic energy, forces and torques
between two sites go smoothly to zero as the distance between the
sites, $r_{\bf ab}$ approaches the cutoff radius, $r_c$.  Requiring
this continuity at the cutoff is essential for energy conservation in
MD simulations.  The mathematical details of the shifted potential
(SP), gradient-shifted-force (GSF) and Taylor shifted-force (TSF)
methods have been discussed in detail in the previous paper in this
series.\cite{PaperI} Here we briefly review the new methods and
describe their essential features.

\subsection{Taylor-shifted force (TSF)}

The electrostatic potential energy between point multipoles can be
expressed as the product of two multipole operators and a Coulombic
kernel,
\begin{equation}
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
\end{equation}
where the multipole operator for site $\bf a$,
\begin{equation}
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} 
+  Q_{{\bf a}\alpha\beta}
 \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
\end{equation}
is expressed in terms of the point charge, $C_{\bf a}$, dipole,
$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for
object $\bf a$.  Note that in this work, we use the primitive
quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k
r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site.

Interactions between multipoles can be expressed as higher derivatives
of the bare Coulomb potential, so one way of ensuring that the forces
and torques vanish at the cutoff distance is to include a larger
number of terms in the truncated Taylor expansion, e.g.,
%
\begin{equation}
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-R_c)^m}{m!} f^{(m)} \Big \lvert  _{R_c}  .
\end{equation}
%
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
Thus, for $f(r)=1/r$, we find
%
\begin{equation}
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} .
\end{equation}
This function is an approximate electrostatic potential that has
vanishing second derivatives at the cutoff radius, making it suitable
for shifting the forces and torques of charge-dipole interactions.

In general, the TSF potential for any multipole-multipole interaction
can be written
\begin{equation}
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r)
\label{generic}
\end{equation}
with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for
charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and
$n=4$ for quadrupole-quadrupole.  To ensure smooth convergence of the
energy, force, and torques, the required number of terms from Taylor
series expansion in $f_n(r)$ must be performed for different
multipole-multipole interactions. 

To carry out the same procedure for a damped electrostatic kernel, we
replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$.
Many of the derivatives of the damped kernel are well known from
Smith's early work on multipoles for the Ewald
summation.\cite{Smith82,Smith98} 

Note that increasing the value of $n$ will add additional terms to the
electrostatic potential, e.g., $f_2(r)$ includes orders up to
$(r-R_c)^3/R_c^4$, and so on.  Successive derivatives of the $f_n(r)$
functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) =
f^{\prime\prime}_2(r)$, etc.  These higher derivatives are required
for computing multipole energies, forces, and torques, and smooth
cutoffs of these quantities can be guaranteed as long as the number of
terms in the Taylor series exceeds the derivative order required.

For multipole-multipole interactions, following this procedure results
in separate radial functions for each distinct orientational
contribution to the potential, and ensures that the forces and torques
from {\it each} of these contributions will vanish at the cutoff
radius.  For example, the direct dipole dot product ($\mathbf{D}_{i}
\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance
dot products:
\begin{equation}
U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot
\mathbf{D}_{j} \right) \frac{g_2(r)}{r} 
-\frac{1}{4\pi \epsilon_0} 
\left( \mathbf{D}_{i} \cdot \hat{r} \right)
\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) -
  \frac{g_2(r)}{r} \right)
\end{equation}

The electrostatic forces and torques acting on the central multipole
site due to another site within cutoff sphere are derived from
Eq.~\ref{generic}, accounting for the appropriate number of
derivatives. Complete energy, force, and torque expressions are
presented in the first paper in this series (Reference
\citep{PaperI}).

\subsection{Gradient-shifted force (GSF)}

A second (and significantly simpler) method involves shifting the
gradient of the raw coulomb potential for each particular multipole
order.  For example, the raw dipole-dipole potential energy may be
shifted smoothly by finding the gradient for two interacting dipoles
which have been projected onto the surface of the cutoff sphere
without changing their relative orientation,
\begin{displaymath}
U_{D_{i}D_{j}}(r_{ij})  = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(R_c)
   - (r_{ij}-R_c) \hat{r}_{ij} \cdot
  \vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{R_c}
\end{displaymath}
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$
and $\mathbf{D}_{j}$, are retained at the cutoff distance (although
the signs are reversed for the dipole that has been projected onto the
cutoff sphere).  In many ways, this simpler approach is closer in
spirit to the original shifted force method, in that it projects a
neutralizing multipole (and the resulting forces from this multipole)
onto a cutoff sphere. The resulting functional forms for the
potentials, forces, and torques turn out to be quite similar in form
to the Taylor-shifted approach, although the radial contributions are
significantly less perturbed by the Gradient-shifted approach than
they are in the Taylor-shifted method.

In general, the gradient shifted potential between a central multipole
and any multipolar site inside the cutoff radius is given by,
\begin{equation}
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r}
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert  _{r_c} \right]
\label{generic2}
\end{equation}
where the sum describes a separate force-shifting that is applied to
each orientational contribution to the energy. 

The third term converges more rapidly than the first two terms as a
function of radius, hence the contribution of the third term is very
small for large cutoff radii.  The force and torque derived from
equation \ref{generic2} are consistent with the energy expression and
approach zero as $r \rightarrow R_c$.  Both the GSF and TSF methods
can be considered generalizations of the original DSF method for
higher order multipole interactions. GSF and TSF are also identical up
to the charge-dipole interaction but generate different expressions in
the energy, force and torque for higher order multipole-multipole
interactions. Complete energy, force, and torque expressions for the
GSF potential are presented in the first paper in this series
(Reference \citep{PaperI})


\subsection{Shifted potential (SP) }
A discontinuous truncation of the electrostatic potential at the
cutoff sphere introduces a severe artifact (oscillation in the
electrostatic energy) even for molecules with the higher-order
multipoles.\cite{PaperI} We have also formulated an extension of the
Wolf approach for point multipoles by simply projecting the image
multipole onto the surface of the cutoff sphere, and including the
interactions with the central multipole and the image. This
effectively shifts the total potential to zero at the cutoff radius,
\begin{equation}
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c)
\label{eq:SP}
\end{equation}          
where the sum describes separate potential shifting that is done for
each orientational contribution to the energy (e.g. the direct dipole
product contribution is shifted {\it separately} from the
dipole-distance terms in dipole-dipole interactions).  Note that this
is not a simple shifting of the total potential at $R_c$. Each radial
contribution is shifted separately.  One consequence of this is that
multipoles that reorient after leaving the cutoff sphere can re-enter
the cutoff sphere without perturbing the total energy.

The potential energy between a central multipole and other multipolar
sites then goes smoothly to zero as $r \rightarrow R_c$. However, the
force and torque obtained from the shifted potential (SP) are
discontinuous at $R_c$. Therefore, MD simulations will still
experience energy drift while operating under the SP potential, but it
may be suitable for Monte Carlo approaches where the configurational
energy differences are the primary quantity of interest.

\subsection{The Self term}
In the TSF, GSF, and SP methods, a self-interaction is retained for
the central multipole interacting with its own image on the surface of
the cutoff sphere.  This self interaction is nearly identical with the
self-terms that arise in the Ewald sum for multipoles.  Complete
expressions for the self terms are presented in the first paper in
this series (Reference \citep{PaperI})  


\section{\label{sec:methodology}Methodology}

To understand how the real-space multipole methods behave in computer
simulations, it is vital to test against established methods for
computing electrostatic interactions in periodic systems, and to
evaluate the size and sources of any errors that arise from the
real-space cutoffs.  In the first paper of this series, we compared
the dipolar and quadrupolar energy expressions against analytic
expressions for ordered dipolar and quadrupolar
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} This work uses the
multipolar Ewald sum as a reference method for comparing energies,
forces, and torques for molecular models that mimic disordered and
ordered condensed-phase systems.  These test-cases include:

\begin{itemize}
\item Soft Dipolar fluids ($\sigma = , \epsilon = , |D| = $)
\item Soft Dipolar solids ($\sigma = , \epsilon = , |D| = $)
\item Soft Quadrupolar fluids ($\sigma = , \epsilon = , Q_{xx} = ...$)
\item Soft Quadrupolar solids  ($\sigma = , \epsilon = , Q_{xx} = ...$)
\item A mixed multipole model for water 
\item A mixed multipole models for water with dissolved ions
\end{itemize}
This last test case exercises all levels of the multipole-multipole
interactions we have derived so far and represents the most complete
test of the new methods.

In the following section, we present results for the total
electrostatic energy, as well as the electrostatic contributions to
the force and torque on each molecule.  These quantities have been
computed using the SP, TSF, and GSF methods, as well as a hard cutoff,
and have been compared with the values obtaine from the multipolar
Ewald sum.  In Mote Carlo (MC) simulations, the energy differences
between two configurations is the primary quantity that governs how
the simulation proceeds. These differences are the most imporant
indicators of the reliability of a method even if the absolute
energies are not exact.  For each of the multipolar systems listed
above, we have compared the change in electrostatic potential energy
($\Delta E$) between 250 statistically-independent configurations.  In
molecular dynamics (MD) simulations, the forces and torques govern the
behavior of the simulation, so we also compute the electrostatic
contributions to the forces and torques.

\subsection{Model systems}
To sample independent configurations of multipolar crystals, a body
centered cubic (BCC) crystal which is a minimum energy structure for
point dipoles was generated using 3,456 molecules.  The multipoles
were translationally locked in their respective crystal sites for
equilibration at a relatively low temperature (50K), so that dipoles
or quadrupoles could freely explore all accessible orientations.  The
translational constraints were removed, and the crystals were
simulated for 10 ps in the microcanonical (NVE) ensemble with an
average temperature of 50 K.  Configurations were sampled at equal
time intervals for the comparison of the configurational energy
differences.  The crystals were not simulated close to the melting
points in order to avoid translational deformation away of the ideal
lattice geometry.

For dipolar, quadrupolar, and mixed-multipole liquid simulations, each
system was created with 2048 molecules oriented randomly.  These were 

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 \ce{Na+} and 24
\ce{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, Lennard-Jones potentials were turned off and purely
electrostatic interaction had been compared.

\subsection{Accuracy of Energy Differences, Forces and Torques}
The pairwise summation techniques (outlined above) were evaluated for
use in MC simulations by studying the energy differences between
different configurations.  We took the Ewald-computed energy
difference between two conformations to be the correct behavior. An
ideal performance by one of the new methods would reproduce these
energy differences exactly. The configurational energies being used
here contain only contributions from electrostatic interactions.
Lennard-Jones interactions were omitted from the comparison as they
should be identical for all methods.

Since none of the real-space methods provide exact energy differences,
we used least square regressions analysiss for the six different
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF
with the multipolar Ewald reference method.  Unitary results for both
the correlation (slope) and correlation coefficient for these
regressions indicate perfect agreement between the real-space method
and the multipolar Ewald sum.

Molecular systems were run long enough to explore independent
configurations and 250 configurations were recorded for comparison.
Each system provided 31,125 energy differences for a total of 186,750
data points.  Similarly, the magnitudes of the forces and torques have
also been compared by using least squares 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 resulting in 3,072,000 data
points for comparison of forces and torques.

\subsection{Analysis of vector quantities}
Getting the magnitudes of the force and torque vectors correct is only
part of the issue for carrying out accurate molecular dynamics
simulations.  Because the real space methods reweight the different
orientational contributions to the energies, it is also important to
understand how the methods impact the \textit{directionality} of the
force and torque vectors. Fisher developed a probablity density
function to analyse directional data sets,
\begin{equation}
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta e^{\kappa \cos\theta}
\label{eq:pdf}
\end{equation}
where $\kappa$ measures directional dispersion of the data around the
mean direction.\cite{fisher53} This quantity $(\kappa)$ can be
estimated as a reciprocal of the circular variance.\cite{Allen91} To
quantify the directional error, forces obtained from the Ewald sum
were taken as the mean (or correct) direction and the angle between
the forces obtained via the Ewald sum and the real-space methods were
evaluated,
\begin{equation}
\cos\theta_i =  \frac{\vec{f}_i^\mathrm{~Ewald} \cdot
  \vec{f}_i^\mathrm{~GSF}}{\left|\vec{f}_i^\mathrm{~Ewald}\right| \left|\vec{f}_i^\mathrm{~GSF}\right|}
\end{equation}
The total angular displacement of the vectors was calculated as,
\begin{equation}
R = \sqrt{\left(\sum\limits_{i=1}^N \cos\theta_i\right)^2 + \left(\sum\limits_{i=1}^N \sin\theta_i\right)^2}
\label{eq:displacement}
\end{equation} 
where $N$ is number of force vectors.  The circular variance is
defined as
\begin{equation}
\mathrm{Var}(\theta) \approx 1/\kappa = 1 - R/N
\end{equation}
The circular variance takes on values between from 0 to 1, with 0
indicating a perfect directional match between the Ewald force vectors
and the real-space forces. Lower values of $\mathrm{Var}(\theta)$
correspond to higher values of $\kappa$, which indicates tighter
clustering of the real-space force vectors around the Ewald forces.

A similar analysis was carried out for the electrostatic contribution
to the molecular torques as well as forces.  

\subsection{Energy conservation}
To test conservation the energy for the methods, the mixed molecular
system of 2000 SSDQ water molecules with 24 \ce{Na+} and 24 \ce{Cl-}
ions was run for 1 ns in the microcanonical ensemble at an average
temperature of 300K.  Each of the different electrostatic methods
(Ewald, Hard, SP, GSF, and TSF) was tested for a range of different
damping values. The molecular system was started with same initial
positions and velocities for all cutoff methods. The energy drift
($\delta E_1$) and standard deviation of the energy about the slope
($\delta E_0$) were evaluated from the total energy of the system as a
function of time.  Although both measures are valuable at
investigating new methods for molecular dynamics, a useful interaction
model must allow for long simulation times with minimal energy drift.

\section{\label{sec:result}RESULTS}
\subsection{Configurational energy differences}
%The magnitude of the fluctuation in the total electrostatic energy per molecule for a dipolar crystal is very high as shown in (Fig … paper I).\cite{PaperI}As soon as, the net dipole moment within a cutoff radius is neutralized in the SP method, the magnitude of the fluctuation in the total electrostatic energy per molecule reduced significantly and rapidly converged to the correct energy constant (Refer figure … Paper I).\cite{PaperI} The GSF potential energy also converged to the correct energy constant for the cutoff radius rc = 6a for the undamped case. The potential energy from the TSF method converges towards the correct value for a very large cutoff radius. The speed of convergence for the all the cutoff methods can be increased by using damping function as shown in figure … Paper I\cite{PaperI}. For the quadrupolar crystals, the fluctuation in the total electrostatic energy for the hard cutoff method is small and short ranged as compared to the dipolar crystals (figure … in the paper I).\cite{PaperI}  Similar to the dipolar crystals, the net quadrupolar neutralization of the cutoff sphere in the SP method reduces oscillation rapidly and converge electrostatic energy to the correct energy constant.
%The oscillation in the the electrostatic energy for the hard cutoff method is even true for the dipolar liquids as shown in Figure ~\ref{fig:rcutConvergence_dipolarLiquid}. As we placed image on the surface of the cutoff sphere in SP method, the oscillation in the energy is reduced. The fluctuation in the energy in liquid is much smaller as compared to the crystal (This result is similar to the results observed by \textit{Wolf et al.} in the case of ionic crystal and Mgo melt). The large magnitude in the fluctuation of the electrostatic energy in the crystal is because of large range of multipole ordering in the crystal. When the energy is evaluated by the direct truncation, it breaks up large number of multipolar ordering leaving behind net multipole moment. But in the case of liquid, there is only local ordering of the multipoles and their ordering disappears in the long range. Therefore, the direct truncation results a small oscillation in the electrostatic energy (which is smaller than deviation SP energy from the Ewald Figure ~\ref{fig:rcutConvergence_dipolarLiquid}) in the case of liquid. Although, the oscillation in the energy is very small for the case of liquid, this affects the change in potential energy ($\triangle E$), which is observed when $\triangle E$ evaluated from the SP method compared with Ewald as shown in figure 4a and 4b.
%\begin{figure}[h!]
%        \centering
%        \includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf}
%        \caption{The energy per molecule plotted against cutoff radius, rc for i) Hard ii) SP iii) GSF, and iv) TSF method. The hard cutoff method shows fluctuation in the electorstatic energy and it disappers in all other methods.  }
%        \label{fig:rcutConvergence_dipolarLiquid}
%    \end{figure}
%In MC simulations, the electrostatic differences between the molecules are important parameter for sampling. We have compared $\triangle E$ from the different methods (Hard, SP, GSF, and TSF) with the Ewald using linear regression analysis. The correlation coefficient ($R^2$) of the regression line measures the deviation of the evaluated quantities from the mean slope. We know that Ewald method evaluates accurate value of the electrostatic energy. Hence, if the proposed methods can quantify the electrostatic energy as good as Ewald then the correlation coefficient is 1.The correlation coefficient is 0 for the completely random result for any physical quantity measured by the proposed method. The slope is a measure of the accuracy of the average of a physical quantity obtained from the proposed methods. If the slope is 1 then we can conclude that the average of the physical quantity measured by the method is as good as Ewald. The deviation of the slope from 1 state that the method used in quantifying physical quantities is statistically biased as compared to Ewald.
%\begin{figure}
%        \centering
%        \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf}
%        \label{fig:barGraph1}
%        \end{figure}
%        \begin{figure}
%        \centering
 %       \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf}
%        \caption{}
       
%        \label{fig:barGraph2}
%      \end{figure} 
%The correlation coefficient ($R^2$) and slope of the linear regression plots for the energy differences for all six different molecular systems is shown in figure 4a and 4b.The plot shows that the correlation coefficient improves for the SP cutoff method as compared to the undamped hard cutoff method in the case of SSDQC, SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar crystal and liquid, the correlation coefficient is almost unchanged and close to 1.  The correlation coefficient is smallest (0.696276 for $r_c$ = 9 $A^o$) for the SSDQC liquid because of the presence of charge-charge and charge-multipole interactions. Since the charge-charge and charge-multipole interaction is long ranged, there is huge deviation of correlation coefficient from 1. Similarly, the quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with compared to interactions in the other multipolar systems, thus the correlation coefficient very close to 1 even for hard cutoff method. The idea of placing image multipole on the surface of the cutoff sphere improves the correlation coefficient and makes it close to 1 for all types of multipolar systems. Similarly the slope is hugely deviated from the correct value for the lower order multipole-multipole interaction and slightly deviated for higher order multipole – multipole interaction. The SP method improves both correlation coefficient ($R^2$) and slope significantly in SSDQC and dipolar systems.  The Slope is found to be deviated more in dipolar crystal as compared to liquid which is associated with the large fluctuation in the electrostatic energy in crystal. The GSF also produced better values of correlation coefficient and slope with the proper selection of the damping alpha (Interested reader can consult accompanying supporting material). The TSF method gives good value of correlation coefficient for the dipolar crystal, dipolar liquid, SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the regression slopes are significantly deviated. 
\begin{figure}
        \centering
        \includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf}
        \caption{The correlation coefficient and regression slope of configurational energy differences for a given method with compared with the reference Ewald method. The value of result equal to 1(dashed line) indicates energy difference is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\  = circle, 12 \AA\  = square 15 \AA\  = inverted triangle)}
        \label{fig:slopeCorr_energy}
    \end{figure} 
The combined correlation coefficient and slope for all six systems is shown in Figure ~\ref{fig:slopeCorr_energy}. The correlation coefficient for the undamped hard cutoff method is does not have good agreement with the Ewald because of the fluctuation of the electrostatic energy in the direct truncation method. This deviation in correlation coefficient is improved by using SP, GSF, and TSF method. But the TSF method worsens the regression slope stating that this method produces statistically more biased result as compared to Ewald. Also the GSF method slightly deviate slope but it can be alleviated by using proper value of damping alpha and cutoff radius. The SP method shows good agreement with Ewald method for all values of damping alpha and radii. 
\subsection{Magnitude of the force and torque vectors}
The comparison of the magnitude of the combined forces and torques for the data accumulated from all system types are shown in Figure ~\ref{fig:slopeCorr_force}. The correlation and slope for the forces agree with the Ewald even for the hard cutoff method. For the system of molecules with higher order multipoles, the interaction is short ranged. Moreover, the force decays more rapidly than the electrostatic energy hence the hard cutoff method also produces good results. Although the pure cutoff gives the good match of the electrostatic force, the discontinuity in the force at the cutoff radius causes problem in the total energy conservation in MD simulations, which will be discussed in detail in subsection D. The correlation coefficient for GSF method also perfectly matches with Ewald but the slope is slightly deviated (due to extra term obtained from the angular differentiation). This deviation in the slope can be alleviated with proper selection of the damping alpha and radii ($\alpha = 0.2$ and $r_c = 12 A^o$ are good choice). The TSF method shows good agreement in the correlation coefficient but the slope is not good as compared to the Ewald.
\begin{figure}
        \centering
        \includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf}
        \caption{The correlation coefficient and regression slope of the magnitude of the force for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\  = circle, 12 \AA\  = square 15 \AA\  = inverted triangle). }
        \label{fig:slopeCorr_force}
    \end{figure} 
The torques appears to be very influenced because of extra term generated when the potential energy is modified to get consistent force and torque.  The result shows that the torque from the hard cutoff method has good agreement with Ewald. As the potential is modified to make it consistent with the force and torque, the correlation and slope is deviated as shown in Figure~\ref{fig:slopeCorr_torque} for SP, GSF and TSF cutoff methods.  But the proper value of the damping alpha and radius can improve the agreement of the GSF with the Ewald method. The TSF method shows worst agreement in the slope as compared to Ewald even for larger cutoff radii.
\begin{figure}
        \centering
        \includegraphics[width=0.5 \textwidth]{torquePlot_slopeCorrelation_combined-crop.pdf}
        \caption{The correlation coefficient and regression slope of the magnitude of the torque for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle).}
        \label{fig:slopeCorr_torque}
    \end{figure}
\subsection{Directionality of the force and torque vectors}   
The accurate evaluation of the direction of the force and torques are also important for the dynamic simulation.In our research, the direction data sets were computed from the purposed method and compared with Ewald using Fisher statistics and results are expressed in terms of circular variance ($Var(\theta$).The force and torque vectors from the purposed method followed Fisher probability distribution function expressed in equation~\ref{eq:pdf}. The circular variance for the force and torque vectors of each molecule in the 250 configurations for all system types is shown in Figure~\ref{fig:slopeCorr_circularVariance}. The direction of the force and torque vectors from hard and SP cutoff methods showed best directional agreement with the Ewald. The force and torque vectors from GSF method also showed good agreement with the Ewald method, which can also be improved by varying damping alpha and cutoff radius.For $\alpha = 0.2$ and $r_c = 12 A^o$, $ Var(\theta) $ for direction of the force was found to be 0.002061 and corresponding value of $\kappa $ was 485.20. Integration of equation ~\ref{eq:pdf} for that corresponding value of $\kappa$ showed that 95\% of force vectors are with in $6.37^o$. The TSF method is the poorest in evaluating accurate direction with compared to Hard, SP, and GSF methods. The circular variance for the direction of the torques is larger as compared to force. For same $\alpha = 0.2, r_c = 12 A^o$ and GSF method, the circular variance was 0.01415, which showed 95\% of torque vectors are within $16.75^o$.The direction of the force and torque vectors can be improved by varying $\alpha$ and $r_c$.

\begin{figure}
        \centering
        \includegraphics[width=0.5 \textwidth]{Variance_forceNtorque_modified-crop.pdf}
        \caption{The circular variance of the data sets of the
          direction of the  force and torque vectors obtained from a
          given method about reference Ewald method. The result equal
          to 0 (dashed line) indicates direction of the vectors are
          indistinguishable from the Ewald method. Here different
          symbols represent different value of the cutoff radius (9
          \AA\ = circle, 12 \AA\ = square, 15 \AA\  = inverted triangle)}
        \label{fig:slopeCorr_circularVariance}
    \end{figure}
\subsection{Total energy conservation}
We have tested the conservation of energy in the SSDQC liquid system
by running system for 1ns in the Hard, SP, GSF and TSF method. The
Hard cutoff method shows very high energy drifts 433.53
KCal/Mol/ns/particle and energy fluctuation 32574.53 Kcal/Mol
(measured by the SD from the slope) for the undamped case, which makes
it completely unusable in MD simulations. The SP method also shows
large value of energy drift 1.289 Kcal/Mol/ns/particle and energy
fluctuation 0.01953 Kcal/Mol. The energy fluctuation in the SP method
is due to the non-vanishing nature of the torque and force at the
cutoff radius. We can improve the energy conservation in some extent
by the proper selection of the damping alpha but the improvement is
not good enough, which can be observed in Figure 9a and 9b .The GSF
and TSF shows very low value of energy drift 0.09016, 0.07371
KCal/Mol/ns/particle and fluctuation 3.016E-6, 1.217E-6 Kcal/Mol
respectively for the undamped case. Since the absolute value of the
evaluated electrostatic energy, force and torque from TSF method are
deviated from the Ewald, it does not mimic MD simulations
appropriately. The electrostatic energy, force and torque from the GSF
method have very good agreement with the Ewald. In addition, the
energy drift and energy fluctuation from the GSF method is much better
than Ewald’s method for reciprocal space vector value ($k_f$) equal to
7 as shown in Figure~\ref{fig:energyDrift} and
~\ref{fig:fluctuation}. We can improve the total energy fluctuation
and drift for the Ewald’s method by increasing size of the reciprocal
space, which extremely increseses the simulation time. In our current
simulation, the simulation time for the Hard, SP, and GSF methods are
about 5.5 times faster than the Ewald method. 

In Fig.~\ref{fig:energyDrift}, $\delta \mbox{E}_1$ is a measure of the
linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle
over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the
standard deviation of the energy fluctuations in units of $\mbox{kcal
  mol}^{-1}$ per particle. In the bottom plot, it is apparent that the
energy drift is reduced by a significant amount (2 to 3 orders of
magnitude improvement at all values of the damping coefficient) by
chosing either of the shifted-force methods over the hard or SP
methods.  We note that the two shifted-force method can give
significantly better energy conservation than the multipolar Ewald sum
with the same choice of real-space cutoffs.

\begin{figure}
  \centering
  \includegraphics[width=\textwidth]{newDrift.pdf}
\label{fig:energyDrift}        
\caption{Analysis of the energy conservation of the real space
  electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in
  energy over time and $\delta \mathrm{E}_0$ is the standard deviation
  of energy fluctuations around this drift.  All simulations were of a
  2000-molecule simulation of SSDQ water with 48 ionic charges at 300
  K starting from the same initial configuration.}
\end{figure} 

\section{CONCLUSION}
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.
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. 
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. 

Note that the TSF, GSF, and SP models are $\mathcal{O}(N)$ methods
that can be made extremely efficient using spline interpolations of
the radial functions.  They require no Fourier transforms or $k$-space
sums, and guarantee the smooth handling of energies, forces, and
torques as multipoles cross the real-space cutoff boundary.  

%\bibliographystyle{aip}
\newpage
\bibliography{references}
\end{document}

%
% ****** End of file aipsamp.tex ******
Revision:	4171
Committed:	Wed Jun 4 19:31:06 2014 UTC (10 years ago) by gezelter
Content type:	application/x-tex
File size:	53105 byte(s)
Log Message:	New figure