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[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^3ME$), 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}
Recently, 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}

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
octupolar 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 dipole-octupole 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 which both smoothly go to zero as an atom aproaches the cutoff
radius. Also, the comparisons of the accuracy of the potential energy
and force between the DSF method and SPME was surprisingly
good.\cite{Fennell:2006lq} The DSF method has seen wide use in to
calculated electrostatic interactions in molecular systems with
relatively uniform charge densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13}

\subsection{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 is replaced by the
complementary error function $\erfc(\alpha r)/r$ to accelerate the
rate of convergence, where $\alpha$ is damping parameter with units of
inverse distance.  Altering the value of $\alpha$ is equivalent to
changing the width of the small Gaussian charge distributions that are
replacing each point charge -- Gaussian overlap integrals yield
complementary error functions when truncated at a finite distance.

e can perform necessary mathematical manipulation
by varying $\alpha$ in the damping function for the calculation of the
electrostatic energy, force and torque\cite{Wolf:1999dn}. By using
suitable value of damping alpha ($\alpha = 0.2$) for a cutoff radius
($r_{c}=9 A$), \textit{Fennel and Gezelter}\cite{Fennell:2006lq}
produced very good agreement of the interaction energies, forces and
torques for charge-charge interactions.\cite{Fennell:2006lq}

\subsection{Point multipoles for CG modeling}
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). 


%\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 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} 

\subsection{Taylor-shifted force(TSF)}
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,

\begin{equation}
\begin{split}
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}})\\
(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}})\\
[(\frac{1}{r}-[\frac{1}{r_c}-(r-r_c)\frac{1}{r_c^2}+(r-r_c)^2\frac{1}{r_c^3}+...)]
\end{split}
\label{eq:TSF}
\end{equation}
   
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.

\subsection{Shifted potential (SP) }
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, 
\begin{equation}
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c)
\label{eq:SP}
\end{equation}          
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.

\subsection{Gradient-shifted force (GSF)}
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,
        \begin{equation}
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}
\label{eq:GSF}
\end{equation}    
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$.
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. 
\subsection{Self term}

\section{\label{sec:test}Test systems}
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.

\subsection{Modeled systems}
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.
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. 

\subsection{Statistical analysis} 
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. 

\subsection{Analysis of vector quantities}
R.A. Fisher has developed a probablity density function to analyse directional data sets is expressed as below,\cite{fisher53}
\begin{equation}
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta \exp(\kappa \cos\theta)
\label{eq:pdf}
\end{equation}
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,
\begin{equation}
R = \sqrt{(\sum\limits_{i=1}^N \sin\theta_i)^2 + (\sum\limits_{i=1}^N \sin\theta_i)^2}
\label{eq:displacement}
\end{equation} 
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. 

\subsection{Energy conservation}
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.

\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 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = 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. 
\begin{figure}
        \centering
        \includegraphics[width=0.5 \textwidth]{log(energyDrift)-crop.pdf}
\label{fig:energyDrift}        
        \end{figure}
\begin{figure}
        \centering
        \includegraphics[width=0.5 \textwidth]{logSD-crop.pdf}       
        \caption{The plot showing (a) standard deviation, and (b) total energy drift in the total energy conservation plot for different values of the damping alpha for different cut off methods. }
        \label{fig:fluctuation}
    \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. 
%\bibliographystyle{aip}
\bibliography{references}
\end{document}

%
% ****** End of file aipsamp.tex ******
Revision:	4167
Committed:	Tue Jun 3 20:31:58 2014 UTC (10 years ago) by gezelter
Content type:	application/x-tex
File size:	44639 byte(s)
Log Message:	Edits