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,
 jmp,
 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}

\begin{document}

\preprint{AIP/123-QED}

\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald
Sum. II. performance in condensed phase simulations}% 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{Suggested keywords}%Use showkeys class option if keyword
                              %display desired
\maketitle


\section{\label{sec:level1}Introduction}
The interaction between charges has always been the most expensive part in molecular simulations.  There have been many efforts to develop practical and efficient method for handling electrostatic interactions. The Ewald’s method has always been accepted as the most precise method for evaluating electrostatic energies, forces and torques. In this method, the conditionally convergent electrostatic energy is converted into the sum of the rapidly converging real and reciprocal space contribution of artificially made periodic system.\cite{Woodcock86, Woodcock75} Because of this artificially created periodicity, Ewald’s sum is not a suitable method to calculate electrostatic interaction in the interfacial molecular systems such as bicrystals, free surfaces, and liquid-vapor interfaces.\cite{Wolf99}To simulate such interfacial systems, the Parry’s extension of the 3D Ewald sum appropriate for the slab geometry is used,\cite{Parry75} which is computationally very expensive.  Also, the reciprocal part of the Ewald’s sum is computationally expensive which makes it inappropriate to use for the larger system. By using Fast Fourier Transform(FFT) in the  particle-mesh Ewald (PME) and particle-particle particle-mesh  Ewald ($P^3ME$) in the reciprocal space term, the computational cost has been decreased from $O(N^2)$ down to $O(Nlog N)$.\cite{Takada93, Gunsteren94, Gunsteren95, Pedersen93, Pedersen95}. Although the computational time has been reduced, the inherent periodicity in the Ewald’s method can be problematic for the interfacial molecular system.\cite{Gezelter06}  Furthermore, the modified Ewald’s methods developed to handle two-dimensional (2D) electrostatic interactions\cite{Parry75, Parry76, Clarke77, Perram79,Rahman89} in the interfacial systems are also computationally expensive.\cite{Spohr97,Berkowitz99}

Recently, \textit{Wolf et al.}\cite{Wolf99} proposed a real space $O(N)$ method for calculating electrostatic interaction between charges. They showed that the effective Coulomb interaction in the condensed system is actually short ranged.\cite{Wolf92, Wolf95}. Furthermore, the Madelung energy of an ion considering lattice summation over neutral dipolar molecules decreases as $r^{-5}$.\cite{Wolf92, Wolf95}. Thus, the careful application of the real-space method for a calculation of the electrostatic energy should be able to obtain correct Madelung energy for a significant size of the cutoff sphere. But the direct truncation of the cutoff sphere for the evaluation of the electrostatic energy always create truncation defect. This cutoff defect in the electrostatic energy is due to the existence of the net charge within the cutoff sphere.\cite{Wolf99} To neutralize net charge within the cutoff sphere, \textit{Wolf et al.}\cite{Wolf99} proposed a method of placing an image charge, for every charge within a cutoff sphere, on the surface to evaluate the electrostatic energy and force. Both the electrostatic energy and force for the central charge are evaluated separately from the interaction of the configuration of real charges within the cutoff sphere and image charges on the surface of the sphere. But the energy of an individual charge due to another charge within the cutoff sphere and its image on the surface is not an integral of their force, as a result the total energy does not conserve in molecular dynamic (MD) simulations.\cite{Zahn02}

The force and torque acting on molecules are the fundamental factors to drive the dynamics of the molecular simulation. \textit{Fennell and Gezelter} proposed the damped shifted force (DSF) potential energy to obtain consistent configurational force on the central charge by the charges within the cutoff sphere and their image charge on the surface. Since the force is consistent with the energy, MD simulations conserve the total energy. Also, the comparison of accuracy of the potential energy and force from DSF method with the Ewald shows surprisingly good results.\cite{Gezelter06}Now a days, the DSF method is being used in several molecular systems with uniform charge density to calculate electrostatic interaction.\cite{Luebke13, Daivis13, Acevedo13, Space12,English08, Lawrence13, Vergne13} 
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{Gordon10, 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{Gezelter06}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). 

\section{\label{sec:level2}Background}

\subsection{Short range nature of electrostatic interaction}
Considering the interaction of an ion with dipolar molecular shell, the effective Columbic potential for a perfect ionic crystal is found to be decreasing as $r^{-5}$.\cite{Wolf99} Furthermore, viewing the NaCl crystal as simple cubic (SC) structure with octupolar $(NaCl)_{4}$ basis, the electrostatic energy per ion converges more rapidly to Madelong than the dipolar approximation.\cite{Wolf92} Also, to find the correct Madelung constant, Lacman.\cite{Lacman65}suggested that the NaCl structure should be constructed in a such way that the finite crystal terminates with only complete $(NaCl)_4$ molecules.  These facts suggest that the Madelung energy is short ranged for a perfect ionic crystal.  
\begin{figure}[h!]
        \centering
        \includegraphics[width=0.50 \textwidth]{chargesystem.pdf}
        \caption{NaCl crystal showing (a) breaking of the charge ordering in the direct spherical truncation, and (b) complete $(NaCl)_{4}$ molecule interacting with the central ion. }
        \label{figure1}
    \end{figure}
 
Any charge in a NaCl crystal is surrounded by opposite charges. Similarly for each pair of charges, there is an opposite pair of charge to its adjacent as shown in Fig 1a.  Furthermore for each group of four charges, there should be an oppositely aligned group of four charges as shown in Fig 1b.  If we consider any group of charges, suppose $(NaCl)_4$, far from the central charge, they have little electrostatic interaction with  the central charge (acts like point octopole when it is far from the center ). But if the cutoff sphere passes through the $(NaCl)_4$ molecule leaving behind net positive or negative charge, the electrostatic contribution due to these broken charges going to be very large (for point charge  radial function $1/r_c$ and for point octupole $1/r_c$). 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{Wolf99}.  This oscillation in the energy around its fully converged value should be due to the non-neutralized value of the dipole and higher order moments within the cutoff sphere.  Recently, \textit{Ikuo Fukuda} used neutralization of the higher order moments for the calculation of the electrostatic interaction of the point charges system.\cite{Fukuda13}

\subsection{Disordered system}
The $r ^{-5}$ convergence behaviors is not only limited to the perfect crystals but also applied in the highly disordered crystal.\cite{Wolf99} At high temperature there should be local ordering of the charge and higher multipole moments in the liquids (To form the structure which is electrostaticaly neutral) but this ordering disappears at the long range. As in ionic crystal, even for liquid positive ion tends to be surrounded by the negative ion and vice versa, so the spherical truncation breaks the short range charge ordering present in the liquid system which results in oscillation (smaller amplitude in electrostatic energy of liquid as compared to crystal).\cite{Wolf99} This idea can also be generalized in molecule with multipole moments assuming local ordering is even true for multipoles, which is supported by the presence of the oscillation of the electrostatic energy as it plotted against the cutoff radius for dipolar liquid 2a. For quadrupolar liquid oscillation damped pretty quickly as seen in Fig 2c because of short range nature of the quadrupole-quadrupole interaction.
\begin{figure}
        \centering
        \includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_dipolar.pdf}
    \end{figure}
\begin{figure}
        \centering
        \includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_quadrupole.pdf}
        \caption{Total energy per molecule against cutoff radius, rc for (a) dipolar liquid (b) dipolar crystal, (c) quadrupolar liquid and (b) quadrupolar crystal to compare the oscillation in the electrostatic energy. The crystalline system shows larger oscillation as compared to liquid. Also the fluctuation in the dipolar system is very large as compared to quadupolar system. }
        \label{figure1}
    \end{figure}
\subsection{Oscillation in the electrostatic energy}
The oscillation of the electrostatic potential energy per molecule for a direct truncation method is associated with the charge neutrality.\cite{Wolf99} The electrostatic energy of a central molecule due to all other molecules within cutoff sphere is plotted against the cutoff radius for (i)dipolar liquid (i) perfect dipolar crystal, (iii) quadrupolar liquid, and (iv)quadrupolar cyrstal in Figure 2a, 2b, 2c and 2d. The larger amplitude in the oscillation of electrostatic energy for the perfect crystal is because of the long range of multipolar ordering.\cite{Wolf99} Moreover, the oscillation damped much faster for the system of higher order multipoles (Compare range of oscillation between the dipolar and quadrupolar system in figure 2a, 2b, 2c and 2d). As in the case of the charge system, this oscillating nature of the electrostatic energy of the central molecule should be due to the net charge-multipole within the cut off sphere. If the amplitude of the oscillation is very large such as in ionic\cite{Wolf99} and dipolar crystal (Figure 3) then it will introduce huge error in the calculation of the absolute energy. On the other hand, if the oscillation is small, it can contribute error in the evaluation of the energy differences between configurations.
\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{Gezelter06}. The hard cutoff method does not ensure the smooth transition of the energy, force, and torque at the cutoff radius.\cite{Wolf99} 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.
\subsection{Damping function}
The damping function used in our research have already been discussed in detail in the paper I.\cite{PaperI} The radial function $1/r$ of the interactions between the charges can be replaced by the complementary error function $erfc(\alpha r)/r$  to accelerate the rate of convergence, where $\alpha$ is damping parameter. We can perform necessary mathematical manipulation by varying $\alpha$ in the damping function for the calculation of the electrostatic energy, force and torque\cite{Wolf99}. By using suitable value of damping alpha ($\alpha = 0.2$) for a cutoff radius ($r_{c}=9 A$), \textit{Fennel and Gezelter}\cite{Gezelter06} produced very good agreement of the interaction energies, forces and torques for charge-charge interactions.\cite{Gezelter06}
\section{METHOD REVIEW}
Any force field associated with MD simulation should address two major issues in the electrostatic interaction. First, it should deal with the breaking of the charge or multipole ordering due to direct spherical truncation. Second, the electrostatic energy, force and the torque between central molecule and any other molecule should smoothly approaches to zero as $r$ tends to $r_c$. The first issue is associated with the oscillation of the total electrostatic potential energy of the central molecule due to all other molecules within cutoff sphere and second issue is related with the continuous nature of the electrostatic interaction at the cutoff radius, which eventually related 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{Shifted potential (SP) }
As it shown in Fig 2, 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. This artifact is due to the existence of higher order moments within the cutoff spheres.The net multipole moment within cutoff sphere is contributed by the breaking of the multipole ordering in direct truncation of the cutoff sphere. The multipole moments of the cutoff sphere can be neutralized by placing image multipole on the surface, for every each multipole within it. The electrostatic potential between multipoles for the SP method is given by, 
\begin{equation}
U_{SF}(\vec r)=U(\vec r) - U(\vec r_c)
\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{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_{i=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}
\end{equation}
   
where $C$ , $D_{\alpha,\beta}$, and $Q_{\alpha,\beta}$ stands for charge, component of the 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 (3) 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{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)=U(\vec r) - U(\vec r_c)-(\vec r-\vec r_c)\cdot\vec \nabla U(\vec r)|_{r=r_c}
\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. 
\section{Test}
\subsection{Test with Ewald}
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.  The Metropolis Monte Carlo algorithm states that the new configuration of the molecular system is accepted if the energy difference between the new and previous configuration $(\triangle E) < 0$ or if  any random number  $R< exp(-\triangle E/kT)$, where R is between 0 to 1, for the case $\triangle E >0$. 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. 
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. The directions of the forces and torques were compared by using formula, , where  and are unit vectors along direction of the force (or torque) evaluated by Ewald and new methods respectively and  is angle between them. The standard deviation about the zero for the angle between the forces and torques was computed for every molecule in each configuration. 
\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{Summation methods}
The Ewald summation for charge, dipole, and quadurpole was performed by using multipolar Ewald’s code in OpendMD/2.1. For different types of multipolar systems, damping alpha (α) for Ewald’s method was derived by plotting total electrostatic energy versus cutoff radius for the different values of $\alpha$ and compared the converged potential energy with the converged electrostatic energy from the pure cutoff method for very large cutoff radius (20 $A^o$). We found 0.3  $A^o^{-1}$  damping alpha for 12 $A^o$ cutoff radius is suitable for all kind of multipolar systems.
The energies, forces and torques for all methods i) Hard, ii) SP, iii) GSF, and iv) TSF   are evaluated for different cutoff radii i) 9, ii) 12, and iii) 15 $A^o$ with damping parameter ($\alpha$) 0.0, 0.1, 0.2, and 0.3 then compared with Ewald’s method. All the simulations for the various systems were conducted in the OpenMD/2.1.
\subsection{Energy conservation and efficiency}
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{RESULTS}
\subsection{Electrostatic energy and 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 3. 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 -refer figure 3) 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 ($\nabla E$), which is observed when $\nabla 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.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{figure1}
    \end{figure}
In MC simulations, the electrostatic differences between the molecules are important parameter for sampling. We have compared $\nabla 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}
        \end{figure}
        \begin{figure}
        \centering
        \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.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{figure1}
    \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]{energy_combined_slope_correlation.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 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle)}
        \label{figure1}
    \end{figure} 
The combined correlation coefficient and slope for all six systems is shown in Figure 5. 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 6. 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 section 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]{force_combined_slope_correlation.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 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle). }
        \label{figure1}
    \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 7 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.50 \textwidth]{torque_combined_slope_correlation.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 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle).}
        \label{figure1}
    \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 computed from the new methods was compared with the Ewald by measuring angle between them. And, the standard deviation of the angles between the force vectors from the proposed methods and Ewald evaluates the accuracy of the new methods with compared to Ewald.  Figure 8 shows that the standard deviation of the combined force and torque vectors for all molecules of the 250 configurations in all system types. 
The result shows that hard cutoff and SP methods show good directional agreement with Ewald in the case of force. The GSF method generate force within few degrees of Ewald method, for example alpha = 0.2 and $r = 12A^o$ the standard deviation is equal to 4.51$^o$.  The TSF method is the poorest in evaluating accurate direction as compared to Hard, SP and GSF methods. The direction of the torques has larger deviation of the angle as compared to the direction of the forces. For the same damping $\alpha = 0.2$ and $rc = 12 A^o$, the standard deviation of the torque is found to be 12.01$^o$ for the GSF method. The standard deviation of the both force and torque can also be minimized by varying damping alpha.
\begin{figure}
        \centering
        \includegraphics[width=0.50 \textwidth]{angle_forceNtorque_modified.pdf}
        \caption{The standard deviation of the angle between force and torque vectors from a given method and the reference Ewald method. The value of 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{figure1}
    \end{figure}
\subsection{Total energy conservation and efficiency }
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 9a and 9b. 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.45 \textwidth]{log(energyDrift).pdf}
        \end{figure}
\begin{figure}
        \centering
        \includegraphics[width=0.45 \textwidth]{logSD.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{figure1}
    \end{figure} 
\section{CONCLUSION}
We have generalized the charged neutralized potential energy originally developed by the Wolf et al.\cite{Wolf99} 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{Gezelter06} 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 $\nabla 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. 
\bibliography{references}
\end{document}

%
% ****** End of file aipsamp.tex ******
Revision:	4114
Committed:	Tue May 13 16:14:48 2014 UTC (10 years, 1 month ago) by mlamichh
Content type:	application/x-tex
File size:	46006 byte(s)
Log Message: