| 20 |
|
% Use this file as a source of example code for your aip document. |
| 21 |
|
% Use the file aiptemplate.tex as a template for your document. |
| 22 |
|
\documentclass[% |
| 23 |
< |
aip, |
| 24 |
< |
jmp, |
| 23 |
> |
aip,jcp, |
| 24 |
|
amsmath,amssymb, |
| 25 |
< |
%preprint,% |
| 26 |
< |
reprint,% |
| 25 |
> |
preprint, |
| 26 |
> |
%reprint,% |
| 27 |
|
%author-year,% |
| 28 |
|
%author-numerical,% |
| 29 |
|
]{revtex4-1} |
| 30 |
|
|
| 31 |
|
\usepackage{graphicx}% Include figure files |
| 32 |
|
\usepackage{dcolumn}% Align table columns on decimal point |
| 33 |
< |
\usepackage{bm}% bold math |
| 33 |
> |
%\usepackage{bm}% bold math |
| 34 |
|
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math |
| 35 |
|
%\linenumbers\relax % Commence numbering lines |
| 36 |
|
\usepackage{amsmath} |
| 37 |
+ |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 38 |
+ |
\usepackage{url} |
| 39 |
+ |
\usepackage[english]{babel} |
| 40 |
|
|
| 41 |
+ |
|
| 42 |
|
\begin{document} |
| 43 |
|
|
| 44 |
|
\preprint{AIP/123-QED} |
| 45 |
|
|
| 46 |
|
\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald |
| 47 |
< |
Sum. II. performance in condensed phase simulations}% Force line breaks with \\ |
| 47 |
> |
Sum. II. Comparison of Simulation Methodologies} % Force line breaks with \\ |
| 48 |
|
|
| 49 |
|
\author{Madan Lamichhane} |
| 50 |
|
\affiliation{Department of Physics, University |
| 68 |
|
|
| 69 |
|
\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 70 |
|
% Classification Scheme. |
| 71 |
< |
\keywords{Suggested keywords}%Use showkeys class option if keyword |
| 72 |
< |
%display desired |
| 71 |
> |
\keywords{Electrostatics, Multipoles, Real-space} |
| 72 |
> |
|
| 73 |
|
\maketitle |
| 74 |
|
|
| 75 |
|
|
| 76 |
|
\section{\label{sec:intro}Introduction} |
| 77 |
< |
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 molecular 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} |
| 77 |
> |
Computing the interactions between electrostatic sites is one of the |
| 78 |
> |
most expensive aspects of molecular simulations, which is why there |
| 79 |
> |
have been significant efforts to develop practical, efficient and |
| 80 |
> |
convergent methods for handling these interactions. Ewald's method is |
| 81 |
> |
perhaps the best known and most accurate method for evaluating |
| 82 |
> |
energies, forces, and torques in explicitly-periodic simulation |
| 83 |
> |
cells. In this approach, the conditionally convergent electrostatic |
| 84 |
> |
energy is converted into two absolutely convergent contributions, one |
| 85 |
> |
which is carried out in real space with a cutoff radius, and one in |
| 86 |
> |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
| 87 |
|
|
| 88 |
+ |
When carried out as originally formulated, the reciprocal-space |
| 89 |
+ |
portion of the Ewald sum exhibits relatively poor computational |
| 90 |
+ |
scaling, making it prohibitive for large systems. By utilizing |
| 91 |
+ |
particle meshes and three dimensional fast Fourier transforms (FFT), |
| 92 |
+ |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 93 |
+ |
($P^3ME$), and smooth particle mesh Ewald (SPME) methods can decrease |
| 94 |
+ |
the computational cost from $O(N^2)$ down to $O(N \log |
| 95 |
+ |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
| 96 |
+ |
|
| 97 |
+ |
Because of the artificial periodicity required for the Ewald sum, the |
| 98 |
+ |
method may require modification to compute interactions for |
| 99 |
+ |
interfacial molecular systems such as membranes and liquid-vapor |
| 100 |
+ |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 101 |
+ |
To simulate interfacial systems, Parry’s extension of the 3D Ewald sum |
| 102 |
+ |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
| 103 |
+ |
periodicity in the Ewald’s method can also be problematic for |
| 104 |
+ |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
| 105 |
+ |
methods that were developed to handle two-dimensional (2D) |
| 106 |
+ |
electrostatic interactions in interfacial systems have not had similar |
| 107 |
+ |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 108 |
+ |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
| 109 |
+ |
|
| 110 |
|
\subsection{Real-space methods} |
| 111 |
< |
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} |
| 111 |
> |
Recently, Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space |
| 112 |
> |
$O(N)$ method for calculating electrostatic interactions between point |
| 113 |
> |
charges. They argued that the effective Coulomb interaction in |
| 114 |
> |
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
| 115 |
> |
an ordered lattice (e.g. when computing the Madelung constant of an |
| 116 |
> |
ionic solid), the material can be considered as a set of ions |
| 117 |
> |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
| 118 |
> |
an effective distance dependence for the electrostatic interactions of |
| 119 |
> |
$r^{-5}$ (see figure \ref{fig:NaCl}. For this reason, careful |
| 120 |
> |
applications of Wolf's method are able to obtain accurate estimates of |
| 121 |
> |
Madelung constants using relatively short cutoff radii. Recently, |
| 122 |
> |
Fukuda used neutralization of the higher order moments for the |
| 123 |
> |
calculation of the electrostatic interaction of the point charges |
| 124 |
> |
system.\cite{Fukuda:2013sf} |
| 125 |
|
|
| 79 |
– |
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. |
| 126 |
|
\begin{figure}[h!] |
| 127 |
< |
\centering |
| 128 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 129 |
< |
\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. } |
| 130 |
< |
\label{fig:NaCl} |
| 131 |
< |
\end{figure} |
| 132 |
< |
|
| 133 |
< |
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 Figure ~\ref{fig:NaCl}. 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 can be due to the non-neutralized value of the 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} |
| 127 |
> |
\centering |
| 128 |
> |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 129 |
> |
\caption{Top: NaCl crystal showing how spherical truncation can |
| 130 |
> |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
| 131 |
> |
molecules interact with the central ion. Bottom: A dipolar |
| 132 |
> |
crystal exhibiting similar behavior and illustrating how the |
| 133 |
> |
effective dipole-octupole interactions can be disrupted by |
| 134 |
> |
spherical truncation.} |
| 135 |
> |
\label{fig:NaCl} |
| 136 |
> |
\end{figure} |
| 137 |
|
|
| 138 |
< |
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} |
| 138 |
> |
|
| 139 |
> |
The direct truncation of interactions at a cutoff radius creates |
| 140 |
> |
truncation defects. Wolf \textit{et al.} further argued that |
| 141 |
> |
truncation errors are due to net charge remaining inside the cutoff |
| 142 |
> |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
| 143 |
> |
placing an image charge on the surface of the cutoff sphere for every |
| 144 |
> |
real charge inside the cutoff. These charges are present for the |
| 145 |
> |
evaluation of both the pair interaction energy and the force, although |
| 146 |
> |
the force expression maintained a discontinuity at the cutoff sphere. |
| 147 |
> |
In the original Wolf formulation, the total energy for the charge and |
| 148 |
> |
image were not equal to the integral of their force expression, and as |
| 149 |
> |
a result, the total energy would not be conserved in molecular |
| 150 |
> |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
| 151 |
> |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
| 152 |
> |
method with commensurate force and energy expressions that do not |
| 153 |
> |
exhibit this problem.\cite{Fennell:2006lq} |
| 154 |
|
|
| 155 |
+ |
Considering the interaction of one central ion in an ionic crystal |
| 156 |
+ |
with a portion of the crystal at some distance, the effective Columbic |
| 157 |
+ |
potential is found to be decreasing as $r^{-5}$. If one views the |
| 158 |
+ |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
| 159 |
+ |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
| 160 |
+ |
rapidly to the Madelung energy than the dipolar |
| 161 |
+ |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
| 162 |
+ |
Lacman suggested that the NaCl structure could be constructed in a way |
| 163 |
+ |
that the finite crystal terminates with complete \ce{(NaCl)4} |
| 164 |
+ |
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded |
| 165 |
+ |
by opposite charges. Similarly for each pair of charges, there is an |
| 166 |
+ |
opposite pair of charge adjacent to it. The central ion sees what is |
| 167 |
+ |
effectively a set of octupoles at large distances. These facts suggest |
| 168 |
+ |
that the Madelung constants are relatively short ranged for perfect |
| 169 |
+ |
ionic crystals.\cite{Wolf:1999dn} |
| 170 |
+ |
|
| 171 |
+ |
One can make a similar argument for crystals of point multipoles. The |
| 172 |
+ |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
| 173 |
+ |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
| 174 |
+ |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
| 175 |
+ |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
| 176 |
+ |
zero net dipole and retain contributions only from higher order |
| 177 |
+ |
multipoles. The effective interaction between a dipole at the center |
| 178 |
+ |
of a crystal and a group of eight dipoles farther away is |
| 179 |
+ |
significantly shorter ranged than the $r^{-3}$ that one would expect |
| 180 |
+ |
for raw dipole-dipole interactions. Only in crystals which retain a |
| 181 |
+ |
bulk dipole moment (e.g. ferroelectrics) does the analogy with the |
| 182 |
+ |
ionic crystal break down -- ferroelectric dipolar crystals can exist, |
| 183 |
+ |
while ionic crystals with net charge in each unit cell would be |
| 184 |
+ |
unstable. |
| 185 |
+ |
|
| 186 |
+ |
In ionic crystals, real-space truncation can break the effective |
| 187 |
+ |
octupolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
| 188 |
+ |
swings in the electrostatic energy as the cutoff radius is increased |
| 189 |
+ |
(or as individual ions move back and forth across the boundary). This |
| 190 |
+ |
is why the image charges were necessary for the Wolf sum to exhibit |
| 191 |
+ |
rapid convergence. Similarly, the real-space truncation of point |
| 192 |
+ |
multipole interactions breaks dipole-octupole arrangements, and image |
| 193 |
+ |
multipoles are required for real-space treatments of electrostatic |
| 194 |
+ |
energies. |
| 195 |
+ |
|
| 196 |
+ |
% Because of this reason, although the nature of electrostatic |
| 197 |
+ |
% interaction short ranged, the hard cutoff sphere creates very large |
| 198 |
+ |
% fluctuation in the electrostatic energy for the perfect crystal. In |
| 199 |
+ |
% addition, the charge neutralized potential proposed by Wolf et |
| 200 |
+ |
% al. converged to correct Madelung constant but still holds oscillation |
| 201 |
+ |
% in the energy about correct Madelung energy.\cite{Wolf:1999dn}. This |
| 202 |
+ |
% oscillation in the energy around its fully converged value can be due |
| 203 |
+ |
% to the non-neutralized value of the higher order moments within the |
| 204 |
+ |
% cutoff sphere. |
| 205 |
+ |
|
| 206 |
+ |
The forces and torques acting on atomic sites are the fundamental |
| 207 |
+ |
factors driving dynamics in molecular simulations. Fennell and |
| 208 |
+ |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 209 |
+ |
obtain consistent energies and forces on the atoms within the cutoff |
| 210 |
+ |
sphere which both smoothly go to zero as an atom aproaches the cutoff |
| 211 |
+ |
radius. Also, the comparisons of the accuracy of the potential energy |
| 212 |
+ |
and force between the DSF method and SPME was surprisingly |
| 213 |
+ |
good.\cite{Fennell:2006lq} The DSF method has seen wide use in to |
| 214 |
+ |
calculated electrostatic interactions in molecular systems with |
| 215 |
+ |
relatively uniform charge densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 216 |
+ |
|
| 217 |
|
\subsection{Damping function} |
| 218 |
< |
The damping function used in our research has 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} |
| 218 |
> |
The damping function used in our research has been discussed in detail |
| 219 |
> |
in the first paper of this series.\cite{PaperI} The radial kernel |
| 220 |
> |
$1/r$ for the interactions between point charges is replaced by the |
| 221 |
> |
complementary error function $\erfc(\alpha r)/r$ to accelerate the |
| 222 |
> |
rate of convergence, where $\alpha$ is damping parameter with units of |
| 223 |
> |
inverse distance. Altering the value of $\alpha$ is equivalent to |
| 224 |
> |
changing the width of the small Gaussian charge distributions that are |
| 225 |
> |
replacing each point charge -- Gaussian overlap integrals yield |
| 226 |
> |
complementary error functions when truncated at a finite distance. |
| 227 |
|
|
| 228 |
+ |
e can perform necessary mathematical manipulation |
| 229 |
+ |
by varying $\alpha$ in the damping function for the calculation of the |
| 230 |
+ |
electrostatic energy, force and torque\cite{Wolf:1999dn}. By using |
| 231 |
+ |
suitable value of damping alpha ($\alpha = 0.2$) for a cutoff radius |
| 232 |
+ |
($r_{c}=9 A$), \textit{Fennel and Gezelter}\cite{Fennell:2006lq} |
| 233 |
+ |
produced very good agreement of the interaction energies, forces and |
| 234 |
+ |
torques for charge-charge interactions.\cite{Fennell:2006lq} |
| 235 |
+ |
|
| 236 |
|
\subsection{Point multipoles for CG modeling} |
| 237 |
< |
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). |
| 237 |
> |
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). |
| 238 |
|
|
| 239 |
|
|
| 240 |
|
%\subsection{Conservation of total energy } |
| 241 |
< |
%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. |
| 241 |
> |
%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. |
| 242 |
|
|
| 243 |
|
\section{\label{sec:method}REVIEW OF METHODS} |
| 244 |
|
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} |
| 322 |
|
% \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 323 |
|
% \caption{} |
| 324 |
|
|
| 325 |
< |
\label{fig:barGraph2} |
| 326 |
< |
\end{figure} |
| 325 |
> |
% \label{fig:barGraph2} |
| 326 |
> |
% \end{figure} |
| 327 |
|
%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. |
| 328 |
|
\begin{figure} |
| 329 |
|
\centering |
| 330 |
|
\includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 331 |
< |
\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)} |
| 331 |
> |
\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)} |
| 332 |
|
\label{fig:slopeCorr_energy} |
| 333 |
|
\end{figure} |
| 334 |
|
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. |
| 337 |
|
\begin{figure} |
| 338 |
|
\centering |
| 339 |
|
\includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 340 |
< |
\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). } |
| 340 |
> |
\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). } |
| 341 |
|
\label{fig:slopeCorr_force} |
| 342 |
|
\end{figure} |
| 343 |
|
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. |
| 370 |
|
\label{fig:fluctuation} |
| 371 |
|
\end{figure} |
| 372 |
|
\section{CONCLUSION} |
| 373 |
< |
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. |
| 373 |
> |
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. |
| 374 |
|
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. |
| 375 |
|
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. |
| 376 |
< |
\bibliographystyle{rev4-1} |
| 376 |
> |
%\bibliographystyle{aip} |
| 377 |
|
\bibliography{references} |
| 378 |
|
\end{document} |
| 379 |
|
|
