| 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{times} |
| 38 |
+ |
\usepackage{mathptm} |
| 39 |
+ |
\usepackage{tabularx} |
| 40 |
+ |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 41 |
+ |
\usepackage{url} |
| 42 |
+ |
\usepackage[english]{babel} |
| 43 |
|
|
| 44 |
+ |
\newcolumntype{Y}{>{\centering\arraybackslash}X} |
| 45 |
+ |
|
| 46 |
|
\begin{document} |
| 47 |
|
|
| 48 |
< |
\preprint{AIP/123-QED} |
| 48 |
> |
%\preprint{AIP/123-QED} |
| 49 |
|
|
| 50 |
< |
\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald |
| 51 |
< |
Sum. II. performance in condensed phase simulations}% Force line breaks with \\ |
| 50 |
> |
\title{Real space alternatives to the Ewald |
| 51 |
> |
Sum. II. Comparison of Methods} % Force line breaks with \\ |
| 52 |
|
|
| 53 |
|
\author{Madan Lamichhane} |
| 54 |
|
\affiliation{Department of Physics, University |
| 67 |
|
% but any date may be explicitly specified |
| 68 |
|
|
| 69 |
|
\begin{abstract} |
| 70 |
< |
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. |
| 70 |
> |
We have tested the real-space shifted potential (SP), |
| 71 |
> |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
| 72 |
> |
for multipoles that were developed in the first paper in this series |
| 73 |
> |
against a reference method. The tests were carried out in a variety |
| 74 |
> |
of condensed-phase environments which were designed to test all |
| 75 |
> |
levels of the multipole-multipole interactions. Comparisons of the |
| 76 |
> |
energy differences between configurations, molecular forces, and |
| 77 |
> |
torques were used to analyze how well the real-space models perform |
| 78 |
> |
relative to the more computationally expensive Ewald sum. We have |
| 79 |
> |
also investigated the energy conservation properties of the new |
| 80 |
> |
methods in molecular dynamics simulations using all of these |
| 81 |
> |
methods. The SP method shows excellent agreement with |
| 82 |
> |
configurational energy differences, forces, and torques, and would |
| 83 |
> |
be suitable for use in Monte Carlo calculations. Of the two new |
| 84 |
> |
shifted-force methods, the GSF approach shows the best agreement |
| 85 |
> |
with Ewald-derived energies, forces, and torques and exhibits energy |
| 86 |
> |
conservation properties that make it an excellent choice for |
| 87 |
> |
efficiently computing electrostatic interactions in molecular |
| 88 |
> |
dynamics simulations. |
| 89 |
|
\end{abstract} |
| 90 |
|
|
| 91 |
< |
\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 91 |
> |
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 92 |
|
% Classification Scheme. |
| 93 |
< |
\keywords{Suggested keywords}%Use showkeys class option if keyword |
| 94 |
< |
%display desired |
| 93 |
> |
\keywords{Electrostatics, Multipoles, Real-space} |
| 94 |
> |
|
| 95 |
|
\maketitle |
| 96 |
|
|
| 97 |
|
|
| 98 |
< |
\section{\label{sec:level1}Introduction} |
| 99 |
< |
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} |
| 98 |
> |
\section{\label{sec:intro}Introduction} |
| 99 |
> |
Computing the interactions between electrostatic sites is one of the |
| 100 |
> |
most expensive aspects of molecular simulations, which is why there |
| 101 |
> |
have been significant efforts to develop practical, efficient and |
| 102 |
> |
convergent methods for handling these interactions. Ewald's method is |
| 103 |
> |
perhaps the best known and most accurate method for evaluating |
| 104 |
> |
energies, forces, and torques in explicitly-periodic simulation |
| 105 |
> |
cells. In this approach, the conditionally convergent electrostatic |
| 106 |
> |
energy is converted into two absolutely convergent contributions, one |
| 107 |
> |
which is carried out in real space with a cutoff radius, and one in |
| 108 |
> |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
| 109 |
|
|
| 110 |
< |
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} |
| 110 |
> |
When carried out as originally formulated, the reciprocal-space |
| 111 |
> |
portion of the Ewald sum exhibits relatively poor computational |
| 112 |
> |
scaling, making it prohibitive for large systems. By utilizing |
| 113 |
> |
particle meshes and three dimensional fast Fourier transforms (FFT), |
| 114 |
> |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
| 115 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
| 116 |
> |
the computational cost from $O(N^2)$ down to $O(N \log |
| 117 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
| 118 |
|
|
| 119 |
< |
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} |
| 120 |
< |
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). |
| 119 |
> |
Because of the artificial periodicity required for the Ewald sum, the |
| 120 |
> |
method may require modification to compute interactions for |
| 121 |
> |
interfacial molecular systems such as membranes and liquid-vapor |
| 122 |
> |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 123 |
> |
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
| 124 |
> |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
| 125 |
> |
periodicity in the Ewald’s method can also be problematic for |
| 126 |
> |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
| 127 |
> |
methods that were developed to handle two-dimensional (2D) |
| 128 |
> |
electrostatic interactions in interfacial systems have not had similar |
| 129 |
> |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 130 |
> |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
| 131 |
|
|
| 132 |
< |
\section{\label{sec:level2}Background} |
| 132 |
> |
\subsection{Real-space methods} |
| 133 |
> |
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
| 134 |
> |
method for calculating electrostatic interactions between point |
| 135 |
> |
charges. They argued that the effective Coulomb interaction in |
| 136 |
> |
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
| 137 |
> |
an ordered lattice (e.g., when computing the Madelung constant of an |
| 138 |
> |
ionic solid), the material can be considered as a set of ions |
| 139 |
> |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
| 140 |
> |
an effective distance dependence for the electrostatic interactions of |
| 141 |
> |
$r^{-5}$ (see figure \ref{fig:NaCl}). For this reason, careful |
| 142 |
> |
applications of Wolf's method are able to obtain accurate estimates of |
| 143 |
> |
Madelung constants using relatively short cutoff radii. Recently, |
| 144 |
> |
Fukuda used neutralization of the higher order moments for the |
| 145 |
> |
calculation of the electrostatic interaction of the point charges |
| 146 |
> |
system.\cite{Fukuda:2013sf} |
| 147 |
|
|
| 83 |
– |
\subsection{Short range nature of electrostatic interaction} |
| 84 |
– |
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. |
| 148 |
|
\begin{figure}[h!] |
| 149 |
< |
\centering |
| 150 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 151 |
< |
\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. } |
| 152 |
< |
\label{fig:NaCl} |
| 153 |
< |
\end{figure} |
| 154 |
< |
|
| 155 |
< |
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 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} |
| 149 |
> |
\centering |
| 150 |
> |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 151 |
> |
\caption{Top: NaCl crystal showing how spherical truncation can |
| 152 |
> |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
| 153 |
> |
molecules interact with the central ion. Bottom: A dipolar |
| 154 |
> |
crystal exhibiting similar behavior and illustrating how the |
| 155 |
> |
effective dipole-octupole interactions can be disrupted by |
| 156 |
> |
spherical truncation.} |
| 157 |
> |
\label{fig:NaCl} |
| 158 |
> |
\end{figure} |
| 159 |
|
|
| 160 |
< |
\subsection{Disordered system} |
| 161 |
< |
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 in Figure ~\ref{fig:rcutConvergence}. For quadrupolar liquid oscillation damped pretty quickly as seen in Figure ~\ref{fig:rcutConvergence_hardQuadrupole} because of short range nature of the quadrupole-quadrupole interaction. |
| 162 |
< |
\begin{figure} |
| 163 |
< |
\centering |
| 164 |
< |
\includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_dipolar.pdf} |
| 165 |
< |
\end{figure} |
| 166 |
< |
\begin{figure} |
| 167 |
< |
\centering |
| 168 |
< |
\includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_quadrupole.pdf} |
| 160 |
> |
The direct truncation of interactions at a cutoff radius creates |
| 161 |
> |
truncation defects. Wolf \textit{et al.} further argued that |
| 162 |
> |
truncation errors are due to net charge remaining inside the cutoff |
| 163 |
> |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
| 164 |
> |
placing an image charge on the surface of the cutoff sphere for every |
| 165 |
> |
real charge inside the cutoff. These charges are present for the |
| 166 |
> |
evaluation of both the pair interaction energy and the force, although |
| 167 |
> |
the force expression maintained a discontinuity at the cutoff sphere. |
| 168 |
> |
In the original Wolf formulation, the total energy for the charge and |
| 169 |
> |
image were not equal to the integral of their force expression, and as |
| 170 |
> |
a result, the total energy would not be conserved in molecular |
| 171 |
> |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
| 172 |
> |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
| 173 |
> |
method with commensurate force and energy expressions that do not |
| 174 |
> |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
| 175 |
> |
methods were also proposed by Chen \textit{et |
| 176 |
> |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
| 177 |
> |
and by Wu and Brooks.\cite{Wu:044107} |
| 178 |
|
|
| 179 |
< |
\caption{Total energy per molecule against cutoff radius, $r_c$ 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. } |
| 180 |
< |
\label{fig:rcutConvergence} |
| 181 |
< |
\end{figure} |
| 182 |
< |
\subsection{Oscillation in the electrostatic energy} |
| 183 |
< |
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. |
| 184 |
< |
\subsection{Conservation of total energy } |
| 185 |
< |
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. |
| 186 |
< |
\subsection{Damping function} |
| 187 |
< |
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} |
| 188 |
< |
\section{METHOD REVIEW} |
| 189 |
< |
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} |
| 179 |
> |
Considering the interaction of one central ion in an ionic crystal |
| 180 |
> |
with a portion of the crystal at some distance, the effective Columbic |
| 181 |
> |
potential is found to be decreasing as $r^{-5}$. If one views the |
| 182 |
> |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
| 183 |
> |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
| 184 |
> |
rapidly to the Madelung energy than the dipolar |
| 185 |
> |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
| 186 |
> |
Lacman suggested that the NaCl structure could be constructed in a way |
| 187 |
> |
that the finite crystal terminates with complete \ce{(NaCl)4} |
| 188 |
> |
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded |
| 189 |
> |
by opposite charges. Similarly for each pair of charges, there is an |
| 190 |
> |
opposite pair of charge adjacent to it. The central ion sees what is |
| 191 |
> |
effectively a set of octupoles at large distances. These facts suggest |
| 192 |
> |
that the Madelung constants are relatively short ranged for perfect |
| 193 |
> |
ionic crystals.\cite{Wolf:1999dn} |
| 194 |
|
|
| 195 |
< |
\subsection{Shifted potential (SP) } |
| 196 |
< |
As it shown in Figure ~\ref{fig:rcutConvergence}, 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, |
| 197 |
< |
\begin{equation} |
| 198 |
< |
U_{SP}(\vec r)=U(\vec r) - U(\vec r_c) |
| 199 |
< |
\label{eq:SP} |
| 200 |
< |
\end{equation} |
| 201 |
< |
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. |
| 202 |
< |
\subsection{Taylor-shifted force(TSF)} |
| 203 |
< |
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, |
| 195 |
> |
One can make a similar argument for crystals of point multipoles. The |
| 196 |
> |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
| 197 |
> |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
| 198 |
> |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
| 199 |
> |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
| 200 |
> |
zero net dipole and retain contributions only from higher order |
| 201 |
> |
multipoles. The effective interaction between a dipole at the center |
| 202 |
> |
of a crystal and a group of eight dipoles farther away is |
| 203 |
> |
significantly shorter ranged than the $r^{-3}$ that one would expect |
| 204 |
> |
for raw dipole-dipole interactions. Only in crystals which retain a |
| 205 |
> |
bulk dipole moment (e.g. ferroelectrics) does the analogy with the |
| 206 |
> |
ionic crystal break down -- ferroelectric dipolar crystals can exist, |
| 207 |
> |
while ionic crystals with net charge in each unit cell would be |
| 208 |
> |
unstable. |
| 209 |
|
|
| 210 |
< |
\begin{equation} |
| 211 |
< |
\begin{split} |
| 212 |
< |
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}})\\ |
| 213 |
< |
(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}})\\ |
| 214 |
< |
[(\frac{1}{r}-[\frac{1}{r_c}-(r-r_c)\frac{1}{r_c^2}+(r-r_c)^2\frac{1}{r_c^3}+...)] |
| 215 |
< |
\end{split} |
| 216 |
< |
\label{eq:TSF} |
| 217 |
< |
\end{equation} |
| 134 |
< |
|
| 135 |
< |
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 ~\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. |
| 136 |
< |
\subsection{Gradient-shifted force (GSF)} |
| 137 |
< |
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, |
| 138 |
< |
\begin{equation} |
| 139 |
< |
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} |
| 140 |
< |
\label{eq:GSF} |
| 141 |
< |
\end{equation} |
| 142 |
< |
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$. |
| 143 |
< |
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. |
| 144 |
< |
\section{Test} |
| 145 |
< |
\subsection{Test with Ewald} |
| 146 |
< |
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. |
| 147 |
< |
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. |
| 210 |
> |
In ionic crystals, real-space truncation can break the effective |
| 211 |
> |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
| 212 |
> |
swings in the electrostatic energy as individual ions move back and |
| 213 |
> |
forth across the boundary. This is why the image charges are |
| 214 |
> |
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
| 215 |
> |
the real-space truncation of point multipole interactions breaks |
| 216 |
> |
higher order multipole arrangements, and image multipoles are required |
| 217 |
> |
for real-space treatments of electrostatic energies. |
| 218 |
|
|
| 219 |
< |
R.A. Fisher has developed a probablity density function to analyse directional data sets is expressed as below, |
| 220 |
< |
\begin{equation} |
| 221 |
< |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta \exp(\kappa \cos\theta) |
| 222 |
< |
\label{eq:pdf} |
| 223 |
< |
\end{equation} |
| 224 |
< |
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{vector_statistics} 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, |
| 225 |
< |
\begin{equation} |
| 226 |
< |
R = \sqrt{(\sum\limits_{i=1}^N \sin\theta_i)^2 + (\sum\limits_{i=1}^N \sin\theta_i)^2} |
| 227 |
< |
\label{eq:displacement} |
| 158 |
< |
\end{equation} |
| 159 |
< |
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. |
| 219 |
> |
% Because of this reason, although the nature of electrostatic |
| 220 |
> |
% interaction short ranged, the hard cutoff sphere creates very large |
| 221 |
> |
% fluctuation in the electrostatic energy for the perfect crystal. In |
| 222 |
> |
% addition, the charge neutralized potential proposed by Wolf et |
| 223 |
> |
% al. converged to correct Madelung constant but still holds oscillation |
| 224 |
> |
% in the energy about correct Madelung energy.\cite{Wolf:1999dn}. This |
| 225 |
> |
% oscillation in the energy around its fully converged value can be due |
| 226 |
> |
% to the non-neutralized value of the higher order moments within the |
| 227 |
> |
% cutoff sphere. |
| 228 |
|
|
| 229 |
< |
\subsection{Modeled systems} |
| 230 |
< |
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. |
| 231 |
< |
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. |
| 232 |
< |
\subsection{Summation methods} |
| 233 |
< |
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. |
| 234 |
< |
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. |
| 235 |
< |
\subsection{Energy conservation and efficiency} |
| 236 |
< |
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. |
| 229 |
> |
The forces and torques acting on atomic sites are the fundamental |
| 230 |
> |
factors driving dynamics in molecular simulations. Fennell and |
| 231 |
> |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
| 232 |
> |
obtain consistent energies and forces on the atoms within the cutoff |
| 233 |
> |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 234 |
> |
aproaches the cutoff radius. The comparisons of the accuracy these |
| 235 |
> |
quantities between the DSF kernel and SPME was surprisingly |
| 236 |
> |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
| 237 |
> |
calculating electrostatic interactions in molecular systems with |
| 238 |
> |
relatively uniform charge |
| 239 |
> |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
| 240 |
|
|
| 241 |
< |
\section{RESULTS} |
| 242 |
< |
\subsection{Electrostatic energy and configurational energy differences} |
| 243 |
< |
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. |
| 244 |
< |
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. |
| 245 |
< |
\begin{figure}[h!] |
| 246 |
< |
\centering |
| 247 |
< |
\includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf} |
| 248 |
< |
\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. } |
| 249 |
< |
\label{fig:rcutConvergence_dipolarLiquid} |
| 250 |
< |
\end{figure} |
| 251 |
< |
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. |
| 181 |
< |
\begin{figure} |
| 182 |
< |
\centering |
| 183 |
< |
\includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 184 |
< |
\label{fig:barGraph1} |
| 185 |
< |
\end{figure} |
| 186 |
< |
\begin{figure} |
| 187 |
< |
\centering |
| 188 |
< |
\includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 189 |
< |
\caption{} |
| 190 |
< |
|
| 191 |
< |
\label{fig:barGraph2} |
| 192 |
< |
\end{figure} |
| 193 |
< |
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. |
| 194 |
< |
\begin{figure} |
| 195 |
< |
\centering |
| 196 |
< |
\includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 197 |
< |
\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)} |
| 198 |
< |
\label{fig:slopeCorr_energy} |
| 199 |
< |
\end{figure} |
| 200 |
< |
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. |
| 201 |
< |
\subsection{Magnitude of the force and torque vectors} |
| 202 |
< |
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 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. |
| 203 |
< |
\begin{figure} |
| 204 |
< |
\centering |
| 205 |
< |
\includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 206 |
< |
\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). } |
| 207 |
< |
\label{fig:slopeCorr_force} |
| 208 |
< |
\end{figure} |
| 209 |
< |
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. |
| 210 |
< |
\begin{figure} |
| 211 |
< |
\centering |
| 212 |
< |
\includegraphics[width=0.5 \textwidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 213 |
< |
\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).} |
| 214 |
< |
\label{fig:slopeCorr_torque} |
| 215 |
< |
\end{figure} |
| 216 |
< |
\subsection{Directionality of the force and torque vectors} |
| 217 |
< |
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$. |
| 241 |
> |
\subsection{The damping function} |
| 242 |
> |
The damping function used in our research has been discussed in detail |
| 243 |
> |
in the first paper of this series.\cite{PaperI} The radial kernel |
| 244 |
> |
$1/r$ for the interactions between point charges can be replaced by |
| 245 |
> |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
| 246 |
> |
accelerate the rate of convergence, where $\alpha$ is a damping |
| 247 |
> |
parameter with units of inverse distance. Altering the value of |
| 248 |
> |
$\alpha$ is equivalent to changing the width of Gaussian charge |
| 249 |
> |
distributions that replace each point charge -- Gaussian overlap |
| 250 |
> |
integrals yield complementary error functions when truncated at a |
| 251 |
> |
finite distance. |
| 252 |
|
|
| 253 |
< |
\begin{figure} |
| 254 |
< |
\centering |
| 255 |
< |
\includegraphics[width=0.5 \textwidth]{Variance_forceNtorque_modified-crop.pdf} |
| 256 |
< |
\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)} |
| 257 |
< |
\label{fig:slopeCorr_circularVariance} |
| 258 |
< |
\end{figure} |
| 259 |
< |
\subsection{Total energy conservation and efficiency } |
| 260 |
< |
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. |
| 261 |
< |
\begin{figure} |
| 262 |
< |
\centering |
| 263 |
< |
\includegraphics[width=0.45 \textwidth]{log(energyDrift).pdf} |
| 264 |
< |
\label{fig:energyDrift} |
| 265 |
< |
\end{figure} |
| 253 |
> |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
| 254 |
> |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
| 255 |
> |
agreement with SPME for the interaction energies, forces and torques |
| 256 |
> |
for charge-charge interactions.\cite{Fennell:2006lq} |
| 257 |
> |
|
| 258 |
> |
\subsection{Point multipoles in molecular modeling} |
| 259 |
> |
Coarse-graining approaches which treat entire molecular subsystems as |
| 260 |
> |
a single rigid body are now widely used. A common feature of many |
| 261 |
> |
coarse-graining approaches is simplification of the electrostatic |
| 262 |
> |
interactions between bodies so that fewer site-site interactions are |
| 263 |
> |
required to compute configurational energies. Many coarse-grained |
| 264 |
> |
molecular structures would normally consist of equal positive and |
| 265 |
> |
negative charges, and rather than use multiple site-site interactions, |
| 266 |
> |
the interaction between higher order multipoles can also be used to |
| 267 |
> |
evaluate a single molecule-molecule |
| 268 |
> |
interaction.\cite{Ren06,Essex10,Essex11} |
| 269 |
> |
|
| 270 |
> |
Because electrons in a molecule are not localized at specific points, |
| 271 |
> |
the assignment of partial charges to atomic centers is a relatively |
| 272 |
> |
rough approximation. Atomic sites can also be assigned point |
| 273 |
> |
multipoles and polarizabilities to increase the accuracy of the |
| 274 |
> |
molecular model. Recently, water has been modeled with point |
| 275 |
> |
multipoles up to octupolar order using the soft sticky |
| 276 |
> |
dipole-quadrupole-octupole (SSDQO) |
| 277 |
> |
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 278 |
> |
multipoles up to quadrupolar order have also been coupled with point |
| 279 |
> |
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
| 280 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But |
| 281 |
> |
using point multipole with the real space truncation without |
| 282 |
> |
accounting for multipolar neutrality will create energy conservation |
| 283 |
> |
issues in molecular dynamics (MD) simulations. |
| 284 |
> |
|
| 285 |
> |
In this paper we test a set of real-space methods that were developed |
| 286 |
> |
for point multipolar interactions. These methods extend the damped |
| 287 |
> |
shifted force (DSF) and Wolf methods originally developed for |
| 288 |
> |
charge-charge interactions and generalize them for higher order |
| 289 |
> |
multipoles. The detailed mathematical development of these methods has |
| 290 |
> |
been presented in the first paper in this series, while this work |
| 291 |
> |
covers the testing the energies, forces, torques, and energy |
| 292 |
> |
conservation properties of the methods in realistic simulation |
| 293 |
> |
environments. In all cases, the methods are compared with the |
| 294 |
> |
reference method, a full multipolar Ewald treatment. |
| 295 |
> |
|
| 296 |
> |
|
| 297 |
> |
%\subsection{Conservation of total energy } |
| 298 |
> |
%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. |
| 299 |
> |
|
| 300 |
> |
\section{\label{sec:method}Review of Methods} |
| 301 |
> |
Any real-space electrostatic method that is suitable for MD |
| 302 |
> |
simulations should have the electrostatic energy, forces and torques |
| 303 |
> |
between two sites go smoothly to zero as the distance between the |
| 304 |
> |
sites, $r_{\bf ab}$ approaches the cutoff radius, $r_c$. Requiring |
| 305 |
> |
this continuity at the cutoff is essential for energy conservation in |
| 306 |
> |
MD simulations. The mathematical details of the shifted potential |
| 307 |
> |
(SP), gradient-shifted-force (GSF) and Taylor shifted-force (TSF) |
| 308 |
> |
methods have been discussed in detail in the previous paper in this |
| 309 |
> |
series.\cite{PaperI} Here we briefly review the new methods and |
| 310 |
> |
describe their essential features. |
| 311 |
> |
|
| 312 |
> |
\subsection{Taylor-shifted force (TSF)} |
| 313 |
> |
|
| 314 |
> |
The electrostatic potential energy between point multipoles can be |
| 315 |
> |
expressed as the product of two multipole operators and a Coulombic |
| 316 |
> |
kernel, |
| 317 |
> |
\begin{equation} |
| 318 |
> |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
| 319 |
> |
\end{equation} |
| 320 |
> |
where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is |
| 321 |
> |
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
| 322 |
> |
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
| 323 |
> |
$\bf a$. |
| 324 |
> |
|
| 325 |
> |
% Interactions between multipoles can be expressed as higher derivatives |
| 326 |
> |
% of the bare Coulomb potential, so one way of ensuring that the forces |
| 327 |
> |
% and torques vanish at the cutoff distance is to include a larger |
| 328 |
> |
% number of terms in the truncated Taylor expansion, e.g., |
| 329 |
> |
% % |
| 330 |
> |
% \begin{equation} |
| 331 |
> |
% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
| 332 |
> |
% \end{equation} |
| 333 |
> |
% % |
| 334 |
> |
% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 335 |
> |
% Thus, for $f(r)=1/r$, we find |
| 336 |
> |
% % |
| 337 |
> |
% \begin{equation} |
| 338 |
> |
% f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . |
| 339 |
> |
% \end{equation} |
| 340 |
> |
% This function is an approximate electrostatic potential that has |
| 341 |
> |
% vanishing second derivatives at the cutoff radius, making it suitable |
| 342 |
> |
% for shifting the forces and torques of charge-dipole interactions. |
| 343 |
> |
|
| 344 |
> |
The TSF potential for any multipole-multipole interaction can be |
| 345 |
> |
written |
| 346 |
> |
\begin{equation} |
| 347 |
> |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
| 348 |
> |
\label{generic} |
| 349 |
> |
\end{equation} |
| 350 |
> |
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
| 351 |
> |
of the interaction, with $n=0$ for charge-charge, $n=1$ for |
| 352 |
> |
charge-dipole, $n=2$ for charge-quadrupole and dipole-dipole, $n=3$ |
| 353 |
> |
for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole. To ensure |
| 354 |
> |
smooth convergence of the energy, force, and torques, a Taylor |
| 355 |
> |
expansion with $n$ terms must be performed at cutoff radius ($r_c$) to |
| 356 |
> |
obtain $f_n(r)$. |
| 357 |
> |
|
| 358 |
> |
% To carry out the same procedure for a damped electrostatic kernel, we |
| 359 |
> |
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
| 360 |
> |
% Many of the derivatives of the damped kernel are well known from |
| 361 |
> |
% Smith's early work on multipoles for the Ewald |
| 362 |
> |
% summation.\cite{Smith82,Smith98} |
| 363 |
> |
|
| 364 |
> |
% Note that increasing the value of $n$ will add additional terms to the |
| 365 |
> |
% electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 366 |
> |
% $(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 367 |
> |
% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 368 |
> |
% f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 369 |
> |
% for computing multipole energies, forces, and torques, and smooth |
| 370 |
> |
% cutoffs of these quantities can be guaranteed as long as the number of |
| 371 |
> |
% terms in the Taylor series exceeds the derivative order required. |
| 372 |
> |
|
| 373 |
> |
For multipole-multipole interactions, following this procedure results |
| 374 |
> |
in separate radial functions for each of the distinct orientational |
| 375 |
> |
contributions to the potential, and ensures that the forces and |
| 376 |
> |
torques from each of these contributions will vanish at the cutoff |
| 377 |
> |
radius. For example, the direct dipole dot product |
| 378 |
> |
($\mathbf{D}_{\bf a} |
| 379 |
> |
\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance |
| 380 |
> |
dot products: |
| 381 |
> |
\begin{equation} |
| 382 |
> |
U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 383 |
> |
\mathbf{D}_{\bf a} \cdot |
| 384 |
> |
\mathbf{D}_{\bf b} \right) v_{21}(r) + |
| 385 |
> |
\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right) |
| 386 |
> |
\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right] |
| 387 |
> |
\end{equation} |
| 388 |
> |
|
| 389 |
> |
For the Taylor shifted (TSF) method with the undamped kernel, |
| 390 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} + |
| 391 |
> |
\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4} |
| 392 |
> |
- \frac{6}{r r_c^2}$. In these functions, one can easily see the |
| 393 |
> |
connection to unmodified electrostatics as well as the smooth |
| 394 |
> |
transition to zero in both these functions as $r\rightarrow r_c$. The |
| 395 |
> |
electrostatic forces and torques acting on the central multipole due |
| 396 |
> |
to another site within cutoff sphere are derived from |
| 397 |
> |
Eq.~\ref{generic}, accounting for the appropriate number of |
| 398 |
> |
derivatives. Complete energy, force, and torque expressions are |
| 399 |
> |
presented in the first paper in this series (Reference |
| 400 |
> |
\onlinecite{PaperI}). |
| 401 |
> |
|
| 402 |
> |
\subsection{Gradient-shifted force (GSF)} |
| 403 |
> |
|
| 404 |
> |
A second (and conceptually simpler) method involves shifting the |
| 405 |
> |
gradient of the raw Coulomb potential for each particular multipole |
| 406 |
> |
order. For example, the raw dipole-dipole potential energy may be |
| 407 |
> |
shifted smoothly by finding the gradient for two interacting dipoles |
| 408 |
> |
which have been projected onto the surface of the cutoff sphere |
| 409 |
> |
without changing their relative orientation, |
| 410 |
> |
\begin{displaymath} |
| 411 |
> |
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
| 412 |
> |
U_{D_{\bf a} D_{\bf b}}(r_c) |
| 413 |
> |
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
| 414 |
> |
\vec{\nabla} U_{D_{\bf a}D_{\bf b}}(r) \Big \lvert _{r_c} |
| 415 |
> |
\end{displaymath} |
| 416 |
> |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf |
| 417 |
> |
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
| 418 |
> |
(although the signs are reversed for the dipole that has been |
| 419 |
> |
projected onto the cutoff sphere). In many ways, this simpler |
| 420 |
> |
approach is closer in spirit to the original shifted force method, in |
| 421 |
> |
that it projects a neutralizing multipole (and the resulting forces |
| 422 |
> |
from this multipole) onto a cutoff sphere. The resulting functional |
| 423 |
> |
forms for the potentials, forces, and torques turn out to be quite |
| 424 |
> |
similar in form to the Taylor-shifted approach, although the radial |
| 425 |
> |
contributions are significantly less perturbed by the gradient-shifted |
| 426 |
> |
approach than they are in the Taylor-shifted method. |
| 427 |
> |
|
| 428 |
> |
For the gradient shifted (GSF) method with the undamped kernel, |
| 429 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
| 430 |
> |
$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. |
| 431 |
> |
Again, these functions go smoothly to zero as $r\rightarrow r_c$, and |
| 432 |
> |
because the Taylor expansion retains only one term, they are |
| 433 |
> |
significantly less perturbed than the TSF functions. |
| 434 |
> |
|
| 435 |
> |
In general, the gradient shifted potential between a central multipole |
| 436 |
> |
and any multipolar site inside the cutoff radius is given by, |
| 437 |
> |
\begin{equation} |
| 438 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 439 |
> |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
| 440 |
> |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
| 441 |
> |
\label{generic2} |
| 442 |
> |
\end{equation} |
| 443 |
> |
where the sum describes a separate force-shifting that is applied to |
| 444 |
> |
each orientational contribution to the energy. |
| 445 |
> |
|
| 446 |
> |
The third term converges more rapidly than the first two terms as a |
| 447 |
> |
function of radius, hence the contribution of the third term is very |
| 448 |
> |
small for large cutoff radii. The force and torque derived from |
| 449 |
> |
equation \ref{generic2} are consistent with the energy expression and |
| 450 |
> |
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
| 451 |
> |
can be considered generalizations of the original DSF method for |
| 452 |
> |
higher order multipole interactions. GSF and TSF are also identical up |
| 453 |
> |
to the charge-dipole interaction but generate different expressions in |
| 454 |
> |
the energy, force and torque for higher order multipole-multipole |
| 455 |
> |
interactions. Complete energy, force, and torque expressions for the |
| 456 |
> |
GSF potential are presented in the first paper in this series |
| 457 |
> |
(Reference~\onlinecite{PaperI}) |
| 458 |
> |
|
| 459 |
> |
|
| 460 |
> |
\subsection{Shifted potential (SP) } |
| 461 |
> |
A discontinuous truncation of the electrostatic potential at the |
| 462 |
> |
cutoff sphere introduces a severe artifact (oscillation in the |
| 463 |
> |
electrostatic energy) even for molecules with the higher-order |
| 464 |
> |
multipoles.\cite{PaperI} We have also formulated an extension of the |
| 465 |
> |
Wolf approach for point multipoles by simply projecting the image |
| 466 |
> |
multipole onto the surface of the cutoff sphere, and including the |
| 467 |
> |
interactions with the central multipole and the image. This |
| 468 |
> |
effectively shifts the total potential to zero at the cutoff radius, |
| 469 |
> |
\begin{equation} |
| 470 |
> |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 471 |
> |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 472 |
> |
\label{eq:SP} |
| 473 |
> |
\end{equation} |
| 474 |
> |
where the sum describes separate potential shifting that is done for |
| 475 |
> |
each orientational contribution to the energy (e.g. the direct dipole |
| 476 |
> |
product contribution is shifted {\it separately} from the |
| 477 |
> |
dipole-distance terms in dipole-dipole interactions). Note that this |
| 478 |
> |
is not a simple shifting of the total potential at $r_c$. Each radial |
| 479 |
> |
contribution is shifted separately. One consequence of this is that |
| 480 |
> |
multipoles that reorient after leaving the cutoff sphere can re-enter |
| 481 |
> |
the cutoff sphere without perturbing the total energy. |
| 482 |
> |
|
| 483 |
> |
For the shifted potential (SP) method with the undamped kernel, |
| 484 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) = |
| 485 |
> |
\frac{3}{r^3} - \frac{3}{r_c^3}$. The potential energy between a |
| 486 |
> |
central multipole and other multipolar sites goes smoothly to zero as |
| 487 |
> |
$r \rightarrow r_c$. However, the force and torque obtained from the |
| 488 |
> |
shifted potential (SP) are discontinuous at $r_c$. MD simulations |
| 489 |
> |
will still experience energy drift while operating under the SP |
| 490 |
> |
potential, but it may be suitable for Monte Carlo approaches where the |
| 491 |
> |
configurational energy differences are the primary quantity of |
| 492 |
> |
interest. |
| 493 |
> |
|
| 494 |
> |
\subsection{The Self Term} |
| 495 |
> |
In the TSF, GSF, and SP methods, a self-interaction is retained for |
| 496 |
> |
the central multipole interacting with its own image on the surface of |
| 497 |
> |
the cutoff sphere. This self interaction is nearly identical with the |
| 498 |
> |
self-terms that arise in the Ewald sum for multipoles. Complete |
| 499 |
> |
expressions for the self terms are presented in the first paper in |
| 500 |
> |
this series (Reference \onlinecite{PaperI}). |
| 501 |
> |
|
| 502 |
> |
|
| 503 |
> |
\section{\label{sec:methodology}Methodology} |
| 504 |
> |
|
| 505 |
> |
To understand how the real-space multipole methods behave in computer |
| 506 |
> |
simulations, it is vital to test against established methods for |
| 507 |
> |
computing electrostatic interactions in periodic systems, and to |
| 508 |
> |
evaluate the size and sources of any errors that arise from the |
| 509 |
> |
real-space cutoffs. In the first paper of this series, we compared |
| 510 |
> |
the dipolar and quadrupolar energy expressions against analytic |
| 511 |
> |
expressions for ordered dipolar and quadrupolar |
| 512 |
> |
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} In this work, we |
| 513 |
> |
used the multipolar Ewald sum as a reference method for comparing |
| 514 |
> |
energies, forces, and torques for molecular models that mimic |
| 515 |
> |
disordered and ordered condensed-phase systems. The parameters used |
| 516 |
> |
in the test cases are given in table~\ref{tab:pars}. |
| 517 |
> |
|
| 518 |
> |
\begin{table} |
| 519 |
> |
\label{tab:pars} |
| 520 |
> |
\caption{The parameters used in the systems used to evaluate the new |
| 521 |
> |
real-space methods. The most comprehensive test was a liquid |
| 522 |
> |
composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-} |
| 523 |
> |
ions). This test excercises all orders of the multipolar |
| 524 |
> |
interactions developed in the first paper.} |
| 525 |
> |
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
| 526 |
> |
& \multicolumn{2}{c|}{LJ parameters} & |
| 527 |
> |
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
| 528 |
> |
Test system & $\sigma$& $\epsilon$ & $C$ & $D$ & |
| 529 |
> |
$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass & $I_{xx}$ & $I_{yy}$ & |
| 530 |
> |
$I_{zz}$ \\ \cline{6-8}\cline{10-12} |
| 531 |
> |
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
| 532 |
> |
\AA\textsuperscript{2})} \\ \hline |
| 533 |
> |
Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ |
| 534 |
> |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ |
| 535 |
> |
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 536 |
> |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ |
| 537 |
> |
SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 538 |
> |
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
| 539 |
> |
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
| 540 |
> |
\end{tabularx} |
| 541 |
> |
\end{table} |
| 542 |
> |
The systems consist of pure multipolar solids (both dipole and |
| 543 |
> |
quadrupole), pure multipolar liquids (both dipole and quadrupole), a |
| 544 |
> |
fluid composed of sites containing both dipoles and quadrupoles |
| 545 |
> |
simultaneously, and a final test case that includes ions with point |
| 546 |
> |
charges in addition to the multipolar fluid. The solid-phase |
| 547 |
> |
parameters were chosen so that the systems can explore some |
| 548 |
> |
orientational freedom for the multipolar sites, while maintaining |
| 549 |
> |
relatively strict translational order. The SSDQ model used here is |
| 550 |
> |
not a particularly accurate water model, but it does test |
| 551 |
> |
dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole |
| 552 |
> |
interactions at roughly the same magnitudes. The last test case, SSDQ |
| 553 |
> |
water with dissolved ions, exercises \textit{all} levels of the |
| 554 |
> |
multipole-multipole interactions we have derived so far and represents |
| 555 |
> |
the most complete test of the new methods. |
| 556 |
> |
|
| 557 |
> |
In the following section, we present results for the total |
| 558 |
> |
electrostatic energy, as well as the electrostatic contributions to |
| 559 |
> |
the force and torque on each molecule. These quantities have been |
| 560 |
> |
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
| 561 |
> |
and have been compared with the values obtained from the multipolar |
| 562 |
> |
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
| 563 |
> |
between two configurations is the primary quantity that governs how |
| 564 |
> |
the simulation proceeds. These differences are the most imporant |
| 565 |
> |
indicators of the reliability of a method even if the absolute |
| 566 |
> |
energies are not exact. For each of the multipolar systems listed |
| 567 |
> |
above, we have compared the change in electrostatic potential energy |
| 568 |
> |
($\Delta E$) between 250 statistically-independent configurations. In |
| 569 |
> |
molecular dynamics (MD) simulations, the forces and torques govern the |
| 570 |
> |
behavior of the simulation, so we also compute the electrostatic |
| 571 |
> |
contributions to the forces and torques. |
| 572 |
> |
|
| 573 |
> |
\subsection{Implementation} |
| 574 |
> |
The real-space methods developed in the first paper in this series |
| 575 |
> |
have been implemented in our group's open source molecular simulation |
| 576 |
> |
program, OpenMD,\cite{openmd} which was used for all calculations in |
| 577 |
> |
this work. The complementary error function can be a relatively slow |
| 578 |
> |
function on some processors, so all of the radial functions are |
| 579 |
> |
precomputed on a fine grid and are spline-interpolated to provide |
| 580 |
> |
values when required. |
| 581 |
> |
|
| 582 |
> |
Using the same simulation code, we compare to a multipolar Ewald sum |
| 583 |
> |
with a reciprocal space cutoff, $k_\mathrm{max} = 7$. Our version of |
| 584 |
> |
the Ewald sum is a re-implementation of the algorithm originally |
| 585 |
> |
proposed by Smith that does not use the particle mesh or smoothing |
| 586 |
> |
approximations.\cite{Smith82,Smith98} In all cases, the quantities |
| 587 |
> |
being compared are the electrostatic contributions to energies, force, |
| 588 |
> |
and torques. All other contributions to these quantities (i.e. from |
| 589 |
> |
Lennard-Jones interactions) are removed prior to the comparisons. |
| 590 |
> |
|
| 591 |
> |
The convergence parameter ($\alpha$) also plays a role in the balance |
| 592 |
> |
of the real-space and reciprocal-space portions of the Ewald |
| 593 |
> |
calculation. Typical molecular mechanics packages set this to a value |
| 594 |
> |
that depends on the cutoff radius and a tolerance (typically less than |
| 595 |
> |
$1 \times 10^{-4}$ kcal/mol). Smaller tolerances are typically |
| 596 |
> |
associated with increasing accuracy at the expense of computational |
| 597 |
> |
time spent on the reciprocal-space portion of the |
| 598 |
> |
summation.\cite{Perram88,Essmann:1995pb} A default tolerance of $1 \times |
| 599 |
> |
10^{-8}$ kcal/mol was used in all Ewald calculations, resulting in |
| 600 |
> |
Ewald coefficient 0.3119 \AA$^{-1}$ for a cutoff radius of 12 \AA. |
| 601 |
> |
|
| 602 |
> |
The real-space models have self-interactions that provide |
| 603 |
> |
contributions to the energies only. Although the self interaction is |
| 604 |
> |
a rapid calculation, we note that in systems with fluctuating charges |
| 605 |
> |
or point polarizabilities, the self-term is not static and must be |
| 606 |
> |
recomputed at each time step. |
| 607 |
> |
|
| 608 |
> |
\subsection{Model systems} |
| 609 |
> |
To sample independent configurations of the multipolar crystals, body |
| 610 |
> |
centered cubic (bcc) crystals, which exhibit the minimum energy |
| 611 |
> |
structures for point dipoles, were generated using 3,456 molecules. |
| 612 |
> |
The multipoles were translationally locked in their respective crystal |
| 613 |
> |
sites for equilibration at a relatively low temperature (50K) so that |
| 614 |
> |
dipoles or quadrupoles could freely explore all accessible |
| 615 |
> |
orientations. The translational constraints were then removed, the |
| 616 |
> |
systems were re-equilibrated, and the crystals were simulated for an |
| 617 |
> |
additional 10 ps in the microcanonical (NVE) ensemble with an average |
| 618 |
> |
temperature of 50 K. The balance between moments of inertia and |
| 619 |
> |
particle mass were chosen to allow orientational sampling without |
| 620 |
> |
significant translational motion. Configurations were sampled at |
| 621 |
> |
equal time intervals in order to compare configurational energy |
| 622 |
> |
differences. The crystals were simulated far from the melting point |
| 623 |
> |
in order to avoid translational deformation away of the ideal lattice |
| 624 |
> |
geometry. |
| 625 |
> |
|
| 626 |
> |
For dipolar, quadrupolar, and mixed-multipole \textit{liquid} |
| 627 |
> |
simulations, each system was created with 2,048 randomly-oriented |
| 628 |
> |
molecules. These were equilibrated at a temperature of 300K for 1 ns. |
| 629 |
> |
Each system was then simulated for 1 ns in the microcanonical (NVE) |
| 630 |
> |
ensemble. We collected 250 different configurations at equal time |
| 631 |
> |
intervals. For the liquid system that included ionic species, we |
| 632 |
> |
converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24 |
| 633 |
> |
\ce{Cl-} ions and re-equilibrated. After equilibration, the system was |
| 634 |
> |
run under the same conditions for 1 ns. A total of 250 configurations |
| 635 |
> |
were collected. In the following comparisons of energies, forces, and |
| 636 |
> |
torques, the Lennard-Jones potentials were turned off and only the |
| 637 |
> |
purely electrostatic quantities were compared with the same values |
| 638 |
> |
obtained via the Ewald sum. |
| 639 |
> |
|
| 640 |
> |
\subsection{Accuracy of Energy Differences, Forces and Torques} |
| 641 |
> |
The pairwise summation techniques (outlined above) were evaluated for |
| 642 |
> |
use in MC simulations by studying the energy differences between |
| 643 |
> |
different configurations. We took the Ewald-computed energy |
| 644 |
> |
difference between two conformations to be the correct behavior. An |
| 645 |
> |
ideal performance by one of the new methods would reproduce these |
| 646 |
> |
energy differences exactly. The configurational energies being used |
| 647 |
> |
here contain only contributions from electrostatic interactions. |
| 648 |
> |
Lennard-Jones interactions were omitted from the comparison as they |
| 649 |
> |
should be identical for all methods. |
| 650 |
> |
|
| 651 |
> |
Since none of the real-space methods provide exact energy differences, |
| 652 |
> |
we used least square regressions analysis for the six different |
| 653 |
> |
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
| 654 |
> |
with the multipolar Ewald reference method. Unitary results for both |
| 655 |
> |
the correlation (slope) and correlation coefficient for these |
| 656 |
> |
regressions indicate perfect agreement between the real-space method |
| 657 |
> |
and the multipolar Ewald sum. |
| 658 |
> |
|
| 659 |
> |
Molecular systems were run long enough to explore independent |
| 660 |
> |
configurations and 250 configurations were recorded for comparison. |
| 661 |
> |
Each system provided 31,125 energy differences for a total of 186,750 |
| 662 |
> |
data points. Similarly, the magnitudes of the forces and torques have |
| 663 |
> |
also been compared using least squares regression analysis. In the |
| 664 |
> |
forces and torques comparison, the magnitudes of the forces acting in |
| 665 |
> |
each molecule for each configuration were evaluated. For example, our |
| 666 |
> |
dipolar liquid simulation contains 2048 molecules and there are 250 |
| 667 |
> |
different configurations for each system resulting in 3,072,000 data |
| 668 |
> |
points for comparison of forces and torques. |
| 669 |
> |
|
| 670 |
> |
\subsection{Analysis of vector quantities} |
| 671 |
> |
Getting the magnitudes of the force and torque vectors correct is only |
| 672 |
> |
part of the issue for carrying out accurate molecular dynamics |
| 673 |
> |
simulations. Because the real space methods reweight the different |
| 674 |
> |
orientational contributions to the energies, it is also important to |
| 675 |
> |
understand how the methods impact the \textit{directionality} of the |
| 676 |
> |
force and torque vectors. Fisher developed a probablity density |
| 677 |
> |
function to analyse directional data sets, |
| 678 |
> |
\begin{equation} |
| 679 |
> |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta e^{\kappa \cos\theta} |
| 680 |
> |
\label{eq:pdf} |
| 681 |
> |
\end{equation} |
| 682 |
> |
where $\kappa$ measures directional dispersion of the data around the |
| 683 |
> |
mean direction.\cite{fisher53} This quantity $(\kappa)$ can be |
| 684 |
> |
estimated as a reciprocal of the circular variance.\cite{Allen91} To |
| 685 |
> |
quantify the directional error, forces obtained from the Ewald sum |
| 686 |
> |
were taken as the mean (or correct) direction and the angle between |
| 687 |
> |
the forces obtained via the Ewald sum and the real-space methods were |
| 688 |
> |
evaluated, |
| 689 |
> |
\begin{equation} |
| 690 |
> |
\cos\theta_i = \frac{\vec{f}_i^\mathrm{~Ewald} \cdot |
| 691 |
> |
\vec{f}_i^\mathrm{~GSF}}{\left|\vec{f}_i^\mathrm{~Ewald}\right| \left|\vec{f}_i^\mathrm{~GSF}\right|} |
| 692 |
> |
\end{equation} |
| 693 |
> |
The total angular displacement of the vectors was calculated as, |
| 694 |
> |
\begin{equation} |
| 695 |
> |
R = \sqrt{\left(\sum\limits_{i=1}^N \cos\theta_i\right)^2 + \left(\sum\limits_{i=1}^N \sin\theta_i\right)^2} |
| 696 |
> |
\label{eq:displacement} |
| 697 |
> |
\end{equation} |
| 698 |
> |
where $N$ is number of force vectors. The circular variance is |
| 699 |
> |
defined as |
| 700 |
> |
\begin{equation} |
| 701 |
> |
\mathrm{Var}(\theta) \approx 1/\kappa = 1 - R/N |
| 702 |
> |
\end{equation} |
| 703 |
> |
The circular variance takes on values between from 0 to 1, with 0 |
| 704 |
> |
indicating a perfect directional match between the Ewald force vectors |
| 705 |
> |
and the real-space forces. Lower values of $\mathrm{Var}(\theta)$ |
| 706 |
> |
correspond to higher values of $\kappa$, which indicates tighter |
| 707 |
> |
clustering of the real-space force vectors around the Ewald forces. |
| 708 |
> |
|
| 709 |
> |
A similar analysis was carried out for the electrostatic contribution |
| 710 |
> |
to the molecular torques as well as forces. |
| 711 |
> |
|
| 712 |
> |
\subsection{Energy conservation} |
| 713 |
> |
To test conservation the energy for the methods, the mixed molecular |
| 714 |
> |
system of 2000 SSDQ water molecules with 24 \ce{Na+} and 24 \ce{Cl-} |
| 715 |
> |
ions was run for 1 ns in the microcanonical ensemble at an average |
| 716 |
> |
temperature of 300K. Each of the different electrostatic methods |
| 717 |
> |
(Ewald, Hard, SP, GSF, and TSF) was tested for a range of different |
| 718 |
> |
damping values. The molecular system was started with same initial |
| 719 |
> |
positions and velocities for all cutoff methods. The energy drift |
| 720 |
> |
($\delta E_1$) and standard deviation of the energy about the slope |
| 721 |
> |
($\delta E_0$) were evaluated from the total energy of the system as a |
| 722 |
> |
function of time. Although both measures are valuable at |
| 723 |
> |
investigating new methods for molecular dynamics, a useful interaction |
| 724 |
> |
model must allow for long simulation times with minimal energy drift. |
| 725 |
> |
|
| 726 |
> |
\section{\label{sec:result}RESULTS} |
| 727 |
> |
\subsection{Configurational energy differences} |
| 728 |
> |
%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. |
| 729 |
> |
%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. |
| 730 |
> |
%\begin{figure}[h!] |
| 731 |
> |
% \centering |
| 732 |
> |
% \includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf} |
| 733 |
> |
% \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. } |
| 734 |
> |
% \label{fig:rcutConvergence_dipolarLiquid} |
| 735 |
> |
% \end{figure} |
| 736 |
> |
%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. |
| 737 |
> |
%\begin{figure} |
| 738 |
> |
% \centering |
| 739 |
> |
% \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 740 |
> |
% \label{fig:barGraph1} |
| 741 |
> |
% \end{figure} |
| 742 |
> |
% \begin{figure} |
| 743 |
> |
% \centering |
| 744 |
> |
% \includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 745 |
> |
% \caption{} |
| 746 |
> |
|
| 747 |
> |
% \label{fig:barGraph2} |
| 748 |
> |
% \end{figure} |
| 749 |
> |
%The correlation coefficient ($R^2$) and slope of the linear |
| 750 |
> |
%regression plots for the energy differences for all six different |
| 751 |
> |
%molecular systems is shown in figure 4a and 4b.The plot shows that |
| 752 |
> |
%the correlation coefficient improves for the SP cutoff method as |
| 753 |
> |
%compared to the undamped hard cutoff method in the case of SSDQC, |
| 754 |
> |
%SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar |
| 755 |
> |
%crystal and liquid, the correlation coefficient is almost unchanged |
| 756 |
> |
%and close to 1. The correlation coefficient is smallest (0.696276 |
| 757 |
> |
%for $r_c$ = 9 $A^\circ$) for the SSDQC liquid because of the presence of |
| 758 |
> |
%charge-charge and charge-multipole interactions. Since the |
| 759 |
> |
%charge-charge and charge-multipole interaction is long ranged, there |
| 760 |
> |
%is huge deviation of correlation coefficient from 1. Similarly, the |
| 761 |
> |
%quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with |
| 762 |
> |
%compared to interactions in the other multipolar systems, thus the |
| 763 |
> |
%correlation coefficient very close to 1 even for hard cutoff |
| 764 |
> |
%method. The idea of placing image multipole on the surface of the |
| 765 |
> |
%cutoff sphere improves the correlation coefficient and makes it close |
| 766 |
> |
%to 1 for all types of multipolar systems. Similarly the slope is |
| 767 |
> |
%hugely deviated from the correct value for the lower order |
| 768 |
> |
%multipole-multipole interaction and slightly deviated for higher |
| 769 |
> |
%order multipole – multipole interaction. The SP method improves both |
| 770 |
> |
%correlation coefficient ($R^2$) and slope significantly in SSDQC and |
| 771 |
> |
%dipolar systems. The Slope is found to be deviated more in dipolar |
| 772 |
> |
%crystal as compared to liquid which is associated with the large |
| 773 |
> |
%fluctuation in the electrostatic energy in crystal. The GSF also |
| 774 |
> |
%produced better values of correlation coefficient and slope with the |
| 775 |
> |
%proper selection of the damping alpha (Interested reader can consult |
| 776 |
> |
%accompanying supporting material). The TSF method gives good value of |
| 777 |
> |
%correlation coefficient for the dipolar crystal, dipolar liquid, |
| 778 |
> |
%SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the |
| 779 |
> |
%regression slopes are significantly deviated. |
| 780 |
> |
|
| 781 |
|
\begin{figure} |
| 782 |
< |
\centering |
| 783 |
< |
\includegraphics[width=0.45 \textwidth]{logSD.pdf} |
| 784 |
< |
\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. } |
| 785 |
< |
\label{fig:fluctuation} |
| 786 |
< |
\end{figure} |
| 782 |
> |
\centering |
| 783 |
> |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 784 |
> |
\caption{Statistical analysis of the quality of configurational |
| 785 |
> |
energy differences for the real-space electrostatic methods |
| 786 |
> |
compared with the reference Ewald sum. Results with a value equal |
| 787 |
> |
to 1 (dashed line) indicate $\Delta E$ values indistinguishable |
| 788 |
> |
from those obtained using the multipolar Ewald sum. Different |
| 789 |
> |
values of the cutoff radius are indicated with different symbols |
| 790 |
> |
(9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted |
| 791 |
> |
triangles).} |
| 792 |
> |
\label{fig:slopeCorr_energy} |
| 793 |
> |
\end{figure} |
| 794 |
> |
|
| 795 |
> |
The combined correlation coefficient and slope for all six systems is |
| 796 |
> |
shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the methods |
| 797 |
> |
reproduce the Ewald configurational energy differences with remarkable |
| 798 |
> |
fidelity. Undamped hard cutoffs introduce a significant amount of |
| 799 |
> |
random scatter in the energy differences which is apparent in the |
| 800 |
> |
reduced value of the correlation coefficient for this method. This |
| 801 |
> |
can be easily understood as configurations which exhibit small |
| 802 |
> |
traversals of a few dipoles or quadrupoles out of the cutoff sphere |
| 803 |
> |
will see large energy jumps when hard cutoffs are used. The |
| 804 |
> |
orientations of the multipoles (particularly in the ordered crystals) |
| 805 |
> |
mean that these energy jumps can go in either direction, producing a |
| 806 |
> |
significant amount of random scatter, but no systematic error. |
| 807 |
> |
|
| 808 |
> |
The TSF method produces energy differences that are highly correlated |
| 809 |
> |
with the Ewald results, but it also introduces a significant |
| 810 |
> |
systematic bias in the values of the energies, particularly for |
| 811 |
> |
smaller cutoff values. The TSF method alters the distance dependence |
| 812 |
> |
of different orientational contributions to the energy in a |
| 813 |
> |
non-uniform way, so the size of the cutoff sphere can have a large |
| 814 |
> |
effect, particularly for the crystalline systems. |
| 815 |
> |
|
| 816 |
> |
Both the SP and GSF methods appear to reproduce the Ewald results with |
| 817 |
> |
excellent fidelity, particularly for moderate damping ($\alpha = |
| 818 |
> |
0.1-0.2$\AA$^{-1}$) and with a commonly-used cutoff value ($r_c = |
| 819 |
> |
12$\AA). With the exception of the undamped hard cutoff, and the TSF |
| 820 |
> |
method with short cutoffs, all of the methods would be appropriate for |
| 821 |
> |
use in Monte Carlo simulations. |
| 822 |
> |
|
| 823 |
> |
\subsection{Magnitude of the force and torque vectors} |
| 824 |
> |
|
| 825 |
> |
The comparisons of the magnitudes of the forces and torques for the |
| 826 |
> |
data accumulated from all six systems are shown in Figures |
| 827 |
> |
~\ref{fig:slopeCorr_force} and \ref{fig:slopeCorr_torque}. The |
| 828 |
> |
correlation and slope for the forces agree well with the Ewald sum |
| 829 |
> |
even for the hard cutoffs. |
| 830 |
> |
|
| 831 |
> |
For systems of molecules with only multipolar interactions, the pair |
| 832 |
> |
energy contributions are quite short ranged. Moreover, the force |
| 833 |
> |
decays more rapidly than the electrostatic energy, hence the hard |
| 834 |
> |
cutoff method can also produce reasonable agreement for this quantity. |
| 835 |
> |
Although the pure cutoff gives reasonably good electrostatic forces |
| 836 |
> |
for pairs of molecules included within each other's cutoff spheres, |
| 837 |
> |
the discontinuity in the force at the cutoff radius can potentially |
| 838 |
> |
cause energy conservation problems as molecules enter and leave the |
| 839 |
> |
cutoff spheres. This is discussed in detail in section |
| 840 |
> |
\ref{sec:conservation}. |
| 841 |
> |
|
| 842 |
> |
The two shifted-force methods (GSF and TSF) exhibit a small amount of |
| 843 |
> |
systematic variation and scatter compared with the Ewald forces. The |
| 844 |
> |
shifted-force models intentionally perturb the forces between pairs of |
| 845 |
> |
molecules inside each other's cutoff spheres in order to correct the |
| 846 |
> |
energy conservation issues, and this perturbation is evident in the |
| 847 |
> |
statistics accumulated for the molecular forces. The GSF |
| 848 |
> |
perturbations are minimal, particularly for moderate damping and |
| 849 |
> |
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
| 850 |
> |
reasonable agreement in the correlation coefficient but again the |
| 851 |
> |
systematic error in the forces is concerning if replication of Ewald |
| 852 |
> |
forces is desired. |
| 853 |
> |
|
| 854 |
> |
\begin{figure} |
| 855 |
> |
\centering |
| 856 |
> |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 857 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 858 |
> |
magnitudes for the real-space electrostatic methods compared with |
| 859 |
> |
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 860 |
> |
line) indicate force magnitude values indistinguishable from those |
| 861 |
> |
obtained using the multipolar Ewald sum. Different values of the |
| 862 |
> |
cutoff radius are indicated with different symbols (9\AA\ = |
| 863 |
> |
circles, 12\AA\ = squares, and 15\AA\ = inverted triangles). } |
| 864 |
> |
\label{fig:slopeCorr_force} |
| 865 |
> |
\end{figure} |
| 866 |
> |
|
| 867 |
> |
|
| 868 |
> |
\begin{figure} |
| 869 |
> |
\centering |
| 870 |
> |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 871 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 872 |
> |
magnitudes for the real-space electrostatic methods compared with |
| 873 |
> |
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 874 |
> |
line) indicate force magnitude values indistinguishable from those |
| 875 |
> |
obtained using the multipolar Ewald sum. Different values of the |
| 876 |
> |
cutoff radius are indicated with different symbols (9\AA\ = |
| 877 |
> |
circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 878 |
> |
\label{fig:slopeCorr_torque} |
| 879 |
> |
\end{figure} |
| 880 |
> |
|
| 881 |
> |
The torques (Fig. \ref{fig:slopeCorr_torque}) appear to be |
| 882 |
> |
significantly influenced by the choice of real-space method. The |
| 883 |
> |
torque expressions have the same distance dependence as the energies, |
| 884 |
> |
which are naturally longer-ranged expressions than the inter-site |
| 885 |
> |
forces. Torques are also quite sensitive to orientations of |
| 886 |
> |
neighboring molecules, even those that are near the cutoff distance. |
| 887 |
> |
|
| 888 |
> |
The results shows that the torque from the hard cutoff method |
| 889 |
> |
reproduces the torques in quite good agreement with the Ewald sum. |
| 890 |
> |
The other real-space methods can cause some deviations, but excellent |
| 891 |
> |
agreement with the Ewald sum torques is recovered at moderate values |
| 892 |
> |
of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1}$) and cutoff |
| 893 |
> |
radius ($r_c \ge 12$\AA). The TSF method exhibits only fair agreement |
| 894 |
> |
in the slope when compared with the Ewald torques even for larger |
| 895 |
> |
cutoff radii. It appears that the severity of the perturbations in |
| 896 |
> |
the TSF method are most in evidence for the torques. |
| 897 |
> |
|
| 898 |
> |
\subsection{Directionality of the force and torque vectors} |
| 899 |
> |
|
| 900 |
> |
The accurate evaluation of force and torque directions is just as |
| 901 |
> |
important for molecular dynamics simulations as the magnitudes of |
| 902 |
> |
these quantities. Force and torque vectors for all six systems were |
| 903 |
> |
analyzed using Fisher statistics, and the quality of the vector |
| 904 |
> |
directionality is shown in terms of circular variance |
| 905 |
> |
($\mathrm{Var}(\theta)$) in figure |
| 906 |
> |
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
| 907 |
> |
from the new real-space methods exhibit nearly-ideal Fisher probability |
| 908 |
> |
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
| 909 |
> |
exhibit the best vectorial agreement with the Ewald sum. The force and |
| 910 |
> |
torque vectors from GSF method also show good agreement with the Ewald |
| 911 |
> |
method, which can also be systematically improved by using moderate |
| 912 |
> |
damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c = |
| 913 |
> |
12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds |
| 914 |
> |
to a distribution with 95\% of force vectors within $6.37^\circ$ of |
| 915 |
> |
the corresponding Ewald forces. The TSF method produces the poorest |
| 916 |
> |
agreement with the Ewald force directions. |
| 917 |
> |
|
| 918 |
> |
Torques are again more perturbed than the forces by the new real-space |
| 919 |
> |
methods, but even here the variance is reasonably small. For the same |
| 920 |
> |
method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA), |
| 921 |
> |
the circular variance was 0.01415, corresponds to a distribution which |
| 922 |
> |
has 95\% of torque vectors are within $16.75^\circ$ of the Ewald |
| 923 |
> |
results. Again, the direction of the force and torque vectors can be |
| 924 |
> |
systematically improved by varying $\alpha$ and $r_c$. |
| 925 |
> |
|
| 926 |
> |
\begin{figure} |
| 927 |
> |
\centering |
| 928 |
> |
\includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} |
| 929 |
> |
\caption{The circular variance of the direction of the force and |
| 930 |
> |
torque vectors obtained from the real-space methods around the |
| 931 |
> |
reference Ewald vectors. A variance equal to 0 (dashed line) |
| 932 |
> |
indicates direction of the force or torque vectors are |
| 933 |
> |
indistinguishable from those obtained from the Ewald sum. Here |
| 934 |
> |
different symbols represent different values of the cutoff radius |
| 935 |
> |
(9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)} |
| 936 |
> |
\label{fig:slopeCorr_circularVariance} |
| 937 |
> |
\end{figure} |
| 938 |
> |
|
| 939 |
> |
\subsection{Energy conservation\label{sec:conservation}} |
| 940 |
> |
|
| 941 |
> |
We have tested the conservation of energy one can expect to see with |
| 942 |
> |
the new real-space methods using the SSDQ water model with a small |
| 943 |
> |
fraction of solvated ions. This is a test system which exercises all |
| 944 |
> |
orders of multipole-multipole interactions derived in the first paper |
| 945 |
> |
in this series and provides the most comprehensive test of the new |
| 946 |
> |
methods. A liquid-phase system was created with 2000 water molecules |
| 947 |
> |
and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
| 948 |
> |
temperature of 300K. After equilibration, this liquid-phase system |
| 949 |
> |
was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with |
| 950 |
> |
a cutoff radius of 12\AA. The value of the damping coefficient was |
| 951 |
> |
also varied from the undamped case ($\alpha = 0$) to a heavily damped |
| 952 |
> |
case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods. A |
| 953 |
> |
sample was also run using the multipolar Ewald sum with the same |
| 954 |
> |
real-space cutoff. |
| 955 |
> |
|
| 956 |
> |
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
| 957 |
> |
energy over time, $\delta E_1$, and the standard deviation of energy |
| 958 |
> |
fluctuations around this drift $\delta E_0$. Both of the |
| 959 |
> |
shifted-force methods (GSF and TSF) provide excellent energy |
| 960 |
> |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
| 961 |
> |
while the hard cutoff is essentially unusable for molecular dynamics. |
| 962 |
> |
SP provides some benefit over the hard cutoff because the energetic |
| 963 |
> |
jumps that happen as particles leave and enter the cutoff sphere are |
| 964 |
> |
somewhat reduced, but like the Wolf method for charges, the SP method |
| 965 |
> |
would not be as useful for molecular dynamics as either of the |
| 966 |
> |
shifted-force methods. |
| 967 |
> |
|
| 968 |
> |
We note that for all tested values of the cutoff radius, the new |
| 969 |
> |
real-space methods can provide better energy conservation behavior |
| 970 |
> |
than the multipolar Ewald sum, even when utilizing a relatively large |
| 971 |
> |
$k$-space cutoff values. |
| 972 |
> |
|
| 973 |
> |
\begin{figure} |
| 974 |
> |
\centering |
| 975 |
> |
\includegraphics[width=\textwidth]{newDrift_12.pdf} |
| 976 |
> |
\label{fig:energyDrift} |
| 977 |
> |
\caption{Analysis of the energy conservation of the real-space |
| 978 |
> |
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
| 979 |
> |
energy over time (in kcal / mol / particle / ns) and $\delta |
| 980 |
> |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 981 |
> |
around this drift (in kcal / mol / particle). All simulations were |
| 982 |
> |
of a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 983 |
> |
300 K starting from the same initial configuration. All runs |
| 984 |
> |
utilized the same real-space cutoff, $r_c = 12$\AA.} |
| 985 |
> |
\end{figure} |
| 986 |
> |
|
| 987 |
> |
|
| 988 |
|
\section{CONCLUSION} |
| 989 |
< |
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. |
| 990 |
< |
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. |
| 991 |
< |
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. |
| 992 |
< |
\bibliographystyle{rev4-1} |
| 989 |
> |
In the first paper in this series, we generalized the |
| 990 |
> |
charge-neutralized electrostatic energy originally developed by Wolf |
| 991 |
> |
\textit{et al.}\cite{Wolf:1999dn} to multipole-multipole interactions |
| 992 |
> |
up to quadrupolar order. The SP method is essentially a |
| 993 |
> |
multipole-capable version of the Wolf model. The SP method for |
| 994 |
> |
multipoles provides excellent agreement with Ewald-derived energies, |
| 995 |
> |
forces and torques, and is suitable for Monte Carlo simulations, |
| 996 |
> |
although the forces and torques retain discontinuities at the cutoff |
| 997 |
> |
distance that prevents its use in molecular dynamics. |
| 998 |
> |
|
| 999 |
> |
We also developed two natural extensions of the damped shifted-force |
| 1000 |
> |
(DSF) model originally proposed by Fennel and |
| 1001 |
> |
Gezelter.\cite{Fennell:2006lq} The GSF and TSF approaches provide |
| 1002 |
> |
smooth truncation of energies, forces, and torques at the real-space |
| 1003 |
> |
cutoff, and both converge to DSF electrostatics for point-charge |
| 1004 |
> |
interactions. The TSF model is based on a high-order truncated Taylor |
| 1005 |
> |
expansion which can be relatively perturbative inside the cutoff |
| 1006 |
> |
sphere. The GSF model takes the gradient from an images of the |
| 1007 |
> |
interacting multipole that has been projected onto the cutoff sphere |
| 1008 |
> |
to derive shifted force and torque expressions, and is a significantly |
| 1009 |
> |
more gentle approach. |
| 1010 |
> |
|
| 1011 |
> |
Of the two newly-developed shifted force models, the GSF method |
| 1012 |
> |
produced quantitative agreement with Ewald energy, force, and torques. |
| 1013 |
> |
It also performs well in conserving energy in MD simulations. The |
| 1014 |
> |
Taylor-shifted (TSF) model provides smooth dynamics, but these take |
| 1015 |
> |
place on a potential energy surface that is significantly perturbed |
| 1016 |
> |
from Ewald-based electrostatics. |
| 1017 |
> |
|
| 1018 |
> |
% The direct truncation of any electrostatic potential energy without |
| 1019 |
> |
% multipole neutralization creates large fluctuations in molecular |
| 1020 |
> |
% simulations. This fluctuation in the energy is very large for the case |
| 1021 |
> |
% of crystal because of long range of multipole ordering (Refer paper |
| 1022 |
> |
% I).\cite{PaperI} This is also significant in the case of the liquid |
| 1023 |
> |
% because of the local multipole ordering in the molecules. If the net |
| 1024 |
> |
% multipole within cutoff radius neutralized within cutoff sphere by |
| 1025 |
> |
% placing image multiples on the surface of the sphere, this fluctuation |
| 1026 |
> |
% in the energy reduced significantly. Also, the multipole |
| 1027 |
> |
% neutralization in the generalized SP method showed very good agreement |
| 1028 |
> |
% with the Ewald as compared to direct truncation for the evaluation of |
| 1029 |
> |
% the $\triangle E$ between the configurations. In MD simulations, the |
| 1030 |
> |
% energy conservation is very important. The conservation of the total |
| 1031 |
> |
% energy can be ensured by i) enforcing the smooth truncation of the |
| 1032 |
> |
% energy, force and torque in the cutoff radius and ii) making the |
| 1033 |
> |
% energy, force and torque consistent with each other. The GSF and TSF |
| 1034 |
> |
% methods ensure the consistency and smooth truncation of the energy, |
| 1035 |
> |
% force and torque at the cutoff radius, as a result show very good |
| 1036 |
> |
% total energy conservation. But the TSF method does not show good |
| 1037 |
> |
% agreement in the absolute value of the electrostatic energy, force and |
| 1038 |
> |
% torque with the Ewald. The GSF method has mimicked Ewald’s force, |
| 1039 |
> |
% energy and torque accurately and also conserved energy. |
| 1040 |
> |
|
| 1041 |
> |
The only cases we have found where the new GSF and SP real-space |
| 1042 |
> |
methods can be problematic are those which retain a bulk dipole moment |
| 1043 |
> |
at large distances (e.g. the $Z_1$ dipolar lattice). In ferroelectric |
| 1044 |
> |
materials, uniform weighting of the orientational contributions can be |
| 1045 |
> |
important for converging the total energy. In these cases, the |
| 1046 |
> |
damping function which causes the non-uniform weighting can be |
| 1047 |
> |
replaced by the bare electrostatic kernel, and the energies return to |
| 1048 |
> |
the expected converged values. |
| 1049 |
> |
|
| 1050 |
> |
Based on the results of this work, the GSF method is a suitable and |
| 1051 |
> |
efficient replacement for the Ewald sum for evaluating electrostatic |
| 1052 |
> |
interactions in MD simulations. Both methods retain excellent |
| 1053 |
> |
fidelity to the Ewald energies, forces and torques. Additionally, the |
| 1054 |
> |
energy drift and fluctuations from the GSF electrostatics are better |
| 1055 |
> |
than a multipolar Ewald sum for finite-sized reciprocal spaces. |
| 1056 |
> |
Because they use real-space cutoffs with moderate cutoff radii, the |
| 1057 |
> |
GSF and SP models approach $\mathcal{O}(N)$ scaling as the system size |
| 1058 |
> |
increases. Additionally, they can be made extremely efficient using |
| 1059 |
> |
spline interpolations of the radial functions. They require no |
| 1060 |
> |
Fourier transforms or $k$-space sums, and guarantee the smooth |
| 1061 |
> |
handling of energies, forces, and torques as multipoles cross the |
| 1062 |
> |
real-space cutoff boundary. |
| 1063 |
> |
|
| 1064 |
> |
\begin{acknowledgments} |
| 1065 |
> |
JDG acknowledges helpful discussions with Christopher |
| 1066 |
> |
Fennell. Support for this project was provided by the National |
| 1067 |
> |
Science Foundation under grant CHE-1362211. Computational time was |
| 1068 |
> |
provided by the Center for Research Computing (CRC) at the |
| 1069 |
> |
University of Notre Dame. |
| 1070 |
> |
\end{acknowledgments} |
| 1071 |
> |
|
| 1072 |
> |
%\bibliographystyle{aip} |
| 1073 |
> |
\newpage |
| 1074 |
|
\bibliography{references} |
| 1075 |
|
\end{document} |
| 1076 |
|
|
