--- trunk/multipole/multipole_2/multipole2.tex 2014/06/06 16:01:36 4175 +++ trunk/multipole/multipole_2/multipole2.tex 2014/08/06 19:10:04 4203 @@ -35,7 +35,7 @@ preprint, %\linenumbers\relax % Commence numbering lines \usepackage{amsmath} \usepackage{times} -\usepackage{mathptm} +\usepackage{mathptmx} \usepackage{tabularx} \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions \usepackage{url} @@ -47,158 +47,158 @@ preprint, %\preprint{AIP/123-QED} -\title{Real space alternatives to the Ewald -Sum. II. Comparison of Methods} % Force line breaks with \\ +\title{Real space electrostatics for multipoles. II. Comparisons with + the Ewald Sum} \author{Madan Lamichhane} - \affiliation{Department of Physics, University -of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\ + \affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} \author{Kathie E. Newman} -\affiliation{Department of Physics, University -of Notre Dame, Notre Dame, IN 46556} +\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} \author{J. Daniel Gezelter}% \email{gezelter@nd.edu.} -\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash -}% +\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556 +} -\date{\today}% It is always \today, today, - % but any date may be explicitly specified +\date{\today} \begin{abstract} - We have tested the real-space shifted potential (SP), - gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods - for multipoles that were developed in the first paper in this series - against a reference method. The tests were carried out in a variety - of condensed-phase environments which were designed to test all - levels of the multipole-multipole interactions. Comparisons of the + We report on tests of the shifted potential (SP), gradient shifted + force (GSF), and Taylor shifted force (TSF) real-space methods for + multipole interactions developed in the first paper in this series, + using the multipolar Ewald sum as a reference method. The tests were + carried out in a variety of condensed-phase environments designed to + test up to quadrupole-quadrupole interactions. Comparisons of the energy differences between configurations, molecular forces, and torques were used to analyze how well the real-space models perform - relative to the more computationally expensive Ewald sum. We have - also investigated the energy conservation properties of the new - methods in molecular dynamics simulations using all of these - methods. The SP method shows excellent agreement with - configurational energy differences, forces, and torques, and would - be suitable for use in Monte Carlo calculations. Of the two new - shifted-force methods, the GSF approach shows the best agreement - with Ewald-derived energies, forces, and torques and exhibits energy - conservation properties that make it an excellent choice for - efficiently computing electrostatic interactions in molecular - dynamics simulations. + relative to the more computationally expensive Ewald treatment. We + have also investigated the energy conservation properties of the new + methods in molecular dynamics simulations. The SP method shows + excellent agreement with configurational energy differences, forces, + and torques, and would be suitable for use in Monte Carlo + calculations. Of the two new shifted-force methods, the GSF + approach shows the best agreement with Ewald-derived energies, + forces, and torques and also exhibits energy conservation properties + that make it an excellent choice for efficient computation of + electrostatic interactions in molecular dynamics simulations. \end{abstract} %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy % Classification Scheme. -\keywords{Electrostatics, Multipoles, Real-space} +%\keywords{Electrostatics, Multipoles, Real-space} \maketitle - \section{\label{sec:intro}Introduction} Computing the interactions between electrostatic sites is one of the -most expensive aspects of molecular simulations, which is why there -have been significant efforts to develop practical, efficient and -convergent methods for handling these interactions. Ewald's method is -perhaps the best known and most accurate method for evaluating -energies, forces, and torques in explicitly-periodic simulation -cells. In this approach, the conditionally convergent electrostatic -energy is converted into two absolutely convergent contributions, one -which is carried out in real space with a cutoff radius, and one in -reciprocal space.\cite{Clarke:1986eu,Woodcock75} +most expensive aspects of molecular simulations. There have been +significant efforts to develop practical, efficient and convergent +methods for handling these interactions. Ewald's method is perhaps the +best known and most accurate method for evaluating energies, forces, +and torques in explicitly-periodic simulation cells. In this approach, +the conditionally convergent electrostatic energy is converted into +two absolutely convergent contributions, one which is carried out in +real space with a cutoff radius, and one in reciprocal +space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} When carried out as originally formulated, the reciprocal-space portion of the Ewald sum exhibits relatively poor computational -scaling, making it prohibitive for large systems. By utilizing -particle meshes and three dimensional fast Fourier transforms (FFT), -the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald -(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease -the computational cost from $O(N^2)$ down to $O(N \log -N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. +scaling, making it prohibitive for large systems. By utilizing a +particle mesh and three dimensional fast Fourier transforms (FFT), the +particle-mesh Ewald (PME), particle-particle particle-mesh Ewald +(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) +methods can decrease the computational cost from $O(N^2)$ down to $O(N +\log +N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} -Because of the artificial periodicity required for the Ewald sum, the -method may require modification to compute interactions for +Because of the artificial periodicity required for the Ewald sum, interfacial molecular systems such as membranes and liquid-vapor -interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} -To simulate interfacial systems, Parry’s extension of the 3D Ewald sum -is appropriate for slab geometries.\cite{Parry:1975if} The inherent -periodicity in the Ewald’s method can also be problematic for -interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald -methods that were developed to handle two-dimensional (2D) -electrostatic interactions in interfacial systems have not had similar -particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, - Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} +interfaces require modifications to the method. Parry's extension of +the three dimensional Ewald sum is appropriate for slab +geometries.\cite{Parry:1975if} Modified Ewald methods that were +developed to handle two-dimensional (2-D) electrostatic +interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} +These methods were originally quite computationally +expensive.\cite{Spohr97,Yeh99} There have been several successful +efforts that reduced the computational cost of 2-D lattice summations, +bringing them more in line with the scaling for the full 3-D +treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The +inherent periodicity required by the Ewald method can also be +problematic in a number of protein/solvent and ionic solution +environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} \subsection{Real-space methods} Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ method for calculating electrostatic interactions between point -charges. They argued that the effective Coulomb interaction in -condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For -an ordered lattice (e.g. when computing the Madelung constant of an -ionic solid), the material can be considered as a set of ions -interacting with neutral dipolar or quadrupolar ``molecules'' giving -an effective distance dependence for the electrostatic interactions of -$r^{-5}$ (see figure \ref{fig:NaCl}. For this reason, careful -applications of Wolf's method are able to obtain accurate estimates of -Madelung constants using relatively short cutoff radii. Recently, -Fukuda used neutralization of the higher order moments for the -calculation of the electrostatic interaction of the point charges -system.\cite{Fukuda:2013sf} +charges. They argued that the effective Coulomb interaction in most +condensed phase systems is effectively short +ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when +computing the Madelung constant of an ionic solid), the material can +be considered as a set of ions interacting with neutral dipolar or +quadrupolar ``molecules'' giving an effective distance dependence for +the electrostatic interactions of $r^{-5}$ (see figure +\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple +cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the +electrostatic energy per ion converges more rapidly to the Madelung +energy than the dipolar approximation.\cite{Wolf92} To find the +correct Madelung constant, Lacman suggested that the NaCl structure +could be constructed in a way that the finite crystal terminates with +complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees +what is effectively a set of octupoles at large distances. These facts +suggest that the Madelung constants are relatively short ranged for +perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful +application of Wolf's method can provide accurate estimates of +Madelung constants using relatively short cutoff radii. -\begin{figure}[h!] +Direct truncation of interactions at a cutoff radius creates numerical +errors. Wolf \textit{et al.} suggest that truncation errors are due +to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To +neutralize this charge they proposed placing an image charge on the +surface of the cutoff sphere for every real charge inside the cutoff. +These charges are present for the evaluation of both the pair +interaction energy and the force, although the force expression +maintains a discontinuity at the cutoff sphere. In the original Wolf +formulation, the total energy for the charge and image were not equal +to the integral of the force expression, and as a result, the total +energy would not be conserved in molecular dynamics (MD) +simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and +Gezelter later proposed shifted force variants of the Wolf method with +commensurate force and energy expressions that do not exhibit this +problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods +were also proposed by Chen \textit{et + al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} +and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly +used additional neutralization of higher order moments for systems of +point charges.\cite{Fukuda:2013sf} + +\begin{figure} \centering - \includegraphics[width=0.50 \textwidth]{chargesystem.pdf} - \caption{Top: NaCl crystal showing how spherical truncation can - breaking effective charge ordering, and how complete \ce{(NaCl)4} - molecules interact with the central ion. Bottom: A dipolar - crystal exhibiting similar behavior and illustrating how the - effective dipole-octupole interactions can be disrupted by - spherical truncation.} - \label{fig:NaCl} + \includegraphics[width=\linewidth]{schematic.eps} + \caption{Top: Ionic systems exhibit local clustering of dissimilar + charges (in the smaller grey circle), so interactions are + effectively charge-multipole at longer distances. With hard + cutoffs, motion of individual charges in and out of the cutoff + sphere can break the effective multipolar ordering. Bottom: + dipolar crystals and fluids have a similar effective + \textit{quadrupolar} ordering (in the smaller grey circles), and + orientational averaging helps to reduce the effective range of the + interactions in the fluid. Placement of reversed image multipoles + on the surface of the cutoff sphere recovers the effective + higher-order multipole behavior.} + \label{fig:schematic} \end{figure} -The direct truncation of interactions at a cutoff radius creates -truncation defects. Wolf \textit{et al.} further argued that -truncation errors are due to net charge remaining inside the cutoff -sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed -placing an image charge on the surface of the cutoff sphere for every -real charge inside the cutoff. These charges are present for the -evaluation of both the pair interaction energy and the force, although -the force expression maintained a discontinuity at the cutoff sphere. -In the original Wolf formulation, the total energy for the charge and -image were not equal to the integral of their force expression, and as -a result, the total energy would not be conserved in molecular -dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and -Fennel and Gezelter later proposed shifted force variants of the Wolf -method with commensurate force and energy expressions that do not -exhibit this problem.\cite{Fennell:2006lq} Related real-space -methods were also proposed by Chen \textit{et - al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} -and by Wu and Brooks.\cite{Wu:044107} - -Considering the interaction of one central ion in an ionic crystal -with a portion of the crystal at some distance, the effective Columbic -potential is found to be decreasing as $r^{-5}$. If one views the -\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar -\ce{(NaCl)4} basis, the electrostatic energy per ion converges more -rapidly to the Madelung energy than the dipolar -approximation.\cite{Wolf92} To find the correct Madelung constant, -Lacman suggested that the NaCl structure could be constructed in a way -that the finite crystal terminates with complete \ce{(NaCl)4} -molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded -by opposite charges. Similarly for each pair of charges, there is an -opposite pair of charge adjacent to it. The central ion sees what is -effectively a set of octupoles at large distances. These facts suggest -that the Madelung constants are relatively short ranged for perfect -ionic crystals.\cite{Wolf:1999dn} - -One can make a similar argument for crystals of point multipoles. The -Luttinger and Tisza treatment of energy constants for dipolar lattices -utilizes 24 basis vectors that contain dipoles at the eight corners of -a unit cube. Only three of these basis vectors, $X_1, Y_1, -\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have -zero net dipole and retain contributions only from higher order -multipoles. The effective interaction between a dipole at the center +One can make a similar effective range argument for crystals of point +\textit{multipoles}. The Luttinger and Tisza treatment of energy +constants for dipolar lattices utilizes 24 basis vectors that contain +dipoles at the eight corners of a unit cube.\cite{LT} Only three of +these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole +moments, while the rest have zero net dipole and retain contributions +only from higher order multipoles. The lowest-energy crystalline +structures are built out of basis vectors that have only residual +quadrupolar moments (e.g. the $Z_5$ array). In these low energy +structures, the effective interaction between a dipole at the center of a crystal and a group of eight dipoles farther away is significantly shorter ranged than the $r^{-3}$ that one would expect for raw dipole-dipole interactions. Only in crystals which retain a @@ -208,15 +208,30 @@ multipolar arrangements (see Fig. \ref{fig:NaCl}), cau unstable. In ionic crystals, real-space truncation can break the effective -multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant -swings in the electrostatic energy as the cutoff radius is increased -(or as individual ions move back and forth across the boundary). This -is why the image charges were necessary for the Wolf sum to exhibit -rapid convergence. Similarly, the real-space truncation of point -multipole interactions breaks higher order multipole arrangements, and -image multipoles are required for real-space treatments of -electrostatic energies. +multipolar arrangements (see Fig. \ref{fig:schematic}), causing +significant swings in the electrostatic energy as individual ions move +back and forth across the boundary. This is why the image charges are +necessary for the Wolf sum to exhibit rapid convergence. Similarly, +the real-space truncation of point multipole interactions breaks +higher order multipole arrangements, and image multipoles are required +for real-space treatments of electrostatic energies. +The shorter effective range of electrostatic interactions is not +limited to perfect crystals, but can also apply in disordered fluids. +Even at elevated temperatures, there is local charge balance in an +ionic liquid, where each positive ion has surroundings dominated by +negaitve ions and vice versa. The reversed-charge images on the +cutoff sphere that are integral to the Wolf and DSF approaches retain +the effective multipolar interactions as the charges traverse the +cutoff boundary. + +In multipolar fluids (see Fig. \ref{fig:schematic}) there is +significant orientational averaging that additionally reduces the +effect of long-range multipolar interactions. The image multipoles +that are introduced in the TSF, GSF, and SP methods mimic this effect +and reduce the effective range of the multipolar interactions as +interacting molecules traverse each other's cutoff boundaries. + % 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 @@ -227,71 +242,66 @@ The forces and torques acting on atomic sites are the % to the non-neutralized value of the higher order moments within the % cutoff sphere. -The forces and torques acting on atomic sites are the fundamental -factors driving dynamics in molecular simulations. Fennell and -Gezelter proposed the damped shifted force (DSF) energy kernel to -obtain consistent energies and forces on the atoms within the cutoff -sphere. Both the energy and the force go smoothly to zero as an atom -aproaches the cutoff radius. The comparisons of the accuracy these -quantities between the DSF kernel and SPME was surprisingly -good.\cite{Fennell:2006lq} The DSF method has seen increasing use for -calculating electrostatic interactions in molecular systems with -relatively uniform charge -densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} +Forces and torques acting on atomic sites are fundamental in driving +dynamics in molecular simulations, and the damped shifted force (DSF) +energy kernel provides consistent energies and forces on charged atoms +within the cutoff sphere. Both the energy and the force go smoothly to +zero as an atom aproaches the cutoff radius. The comparisons of the +accuracy these quantities between the DSF kernel and SPME was +surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method +has seen increasing use in molecular systems with relatively uniform +charge +densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} \subsection{The damping function} -The damping function used in our research has been discussed in detail -in the first paper of this series.\cite{PaperI} The radial kernel -$1/r$ for the interactions between point charges can be replaced by -the complementary error function $\mathrm{erfc}(\alpha r)/r$ to -accelerate the rate of convergence, where $\alpha$ is a damping -parameter with units of inverse distance. Altering the value of -$\alpha$ is equivalent to changing the width of Gaussian charge -distributions that replace each point charge -- Gaussian overlap -integrals yield complementary error functions when truncated at a -finite distance. +The damping function has been discussed in detail in the first paper +of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the +interactions between point charges can be replaced by the +complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate +convergence, where $\alpha$ is a damping parameter with units of +inverse distance. Altering the value of $\alpha$ is equivalent to +changing the width of Gaussian charge distributions that replace each +point charge, as Coulomb integrals with Gaussian charge distributions +produce complementary error functions when truncated at a finite +distance. -By using suitable value of damping alpha ($\alpha \sim 0.2$) for a -cutoff radius ($r_{­c}=9 A$), Fennel and Gezelter produced very good -agreement with SPME for the interaction energies, forces and torques -for charge-charge interactions.\cite{Fennell:2006lq} +With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method +produced very good agreement with SPME for interaction energies, +forces and torques for charge-charge +interactions.\cite{Fennell:2006lq} \subsection{Point multipoles in molecular modeling} Coarse-graining approaches which treat entire molecular subsystems as a single rigid body are now widely used. A common feature of many coarse-graining approaches is simplification of the electrostatic interactions between bodies so that fewer site-site interactions are -required to compute configurational energies. Many coarse-grained -molecular structures would normally consist of equal positive and -negative charges, and rather than use multiple site-site interactions, -the interaction between higher order multipoles can also be used to -evaluate a single molecule-molecule -interaction.\cite{Ren06,Essex10,Essex11} +required to compute configurational +energies.\cite{Ren06,Essex10,Essex11} -Because electrons in a molecule are not localized at specific points, -the assignment of partial charges to atomic centers is a relatively -rough approximation. Atomic sites can also be assigned point -multipoles and polarizabilities to increase the accuracy of the -molecular model. Recently, water has been modeled with point -multipoles up to octupolar -order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point +Additionally, because electrons in a molecule are not localized at +specific points, the assignment of partial charges to atomic centers +is always an approximation. For increased accuracy, atomic sites can +also be assigned point multipoles and polarizabilities. Recently, +water has been modeled with point multipoles up to octupolar order +using the soft sticky dipole-quadrupole-octupole (SSDQO) +model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point multipoles up to quadrupolar order have also been coupled with point polarizabilities in the high-quality AMOEBA and iAMOEBA water -models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}. But -using point multipole with the real space truncation without -accounting for multipolar neutrality will create energy conservation -issues in molecular dynamics (MD) simulations. +models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, +truncating point multipoles without smoothing the forces and torques +can create energy conservation issues in molecular dynamics +simulations. In this paper we test a set of real-space methods that were developed for point multipolar interactions. These methods extend the damped shifted force (DSF) and Wolf methods originally developed for charge-charge interactions and generalize them for higher order -multipoles. The detailed mathematical development of these methods has -been presented in the first paper in this series, while this work -covers the testing the energies, forces, torques, and energy +multipoles. The detailed mathematical development of these methods +has been presented in the first paper in this series, while this work +covers the testing of energies, forces, torques, and energy conservation properties of the methods in realistic simulation environments. In all cases, the methods are compared with the -reference method, a full multipolar Ewald treatment. +reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} %\subsection{Conservation of total energy } @@ -317,83 +327,83 @@ where the multipole operator for site $\bf a$, \begin{equation} U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. \end{equation} -where the multipole operator for site $\bf a$, -\begin{equation} -\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} -+ Q_{{\bf a}\alpha\beta} - \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots -\end{equation} -is expressed in terms of the point charge, $C_{\bf a}$, dipole, -$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for -object $\bf a$. Note that in this work, we use the primitive -quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k -r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site. +where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is +expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf + a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object +$\bf a$, etc. -Interactions between multipoles can be expressed as higher derivatives -of the bare Coulomb potential, so one way of ensuring that the forces -and torques vanish at the cutoff distance is to include a larger -number of terms in the truncated Taylor expansion, e.g., -% -\begin{equation} -f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . -\end{equation} -% -The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. -Thus, for $f(r)=1/r$, we find -% -\begin{equation} -f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . -\end{equation} -This function is an approximate electrostatic potential that has -vanishing second derivatives at the cutoff radius, making it suitable -for shifting the forces and torques of charge-dipole interactions. +% Interactions between multipoles can be expressed as higher derivatives +% of the bare Coulomb potential, so one way of ensuring that the forces +% and torques vanish at the cutoff distance is to include a larger +% number of terms in the truncated Taylor expansion, e.g., +% % +% \begin{equation} +% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . +% \end{equation} +% % +% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. +% Thus, for $f(r)=1/r$, we find +% % +% \begin{equation} +% f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . +% \end{equation} +% This function is an approximate electrostatic potential that has +% vanishing second derivatives at the cutoff radius, making it suitable +% for shifting the forces and torques of charge-dipole interactions. -In general, the TSF potential for any multipole-multipole interaction -can be written +The TSF potential for any multipole-multipole interaction can be +written \begin{equation} U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) \label{generic} \end{equation} -with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for -charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and -$n=4$ for quadrupole-quadrupole. To ensure smooth convergence of the -energy, force, and torques, the required number of terms from Taylor -series expansion in $f_n(r)$ must be performed for different -multipole-multipole interactions. +where $f_n(r)$ is a shifted kernel that is appropriate for the order +of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for +charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole +and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for +quadrupole-quadrupole. To ensure smooth convergence of the energy, +force, and torques, a Taylor expansion with $n$ terms must be +performed at cutoff radius ($r_c$) to obtain $f_n(r)$. -To carry out the same procedure for a damped electrostatic kernel, we -replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. -Many of the derivatives of the damped kernel are well known from -Smith's early work on multipoles for the Ewald -summation.\cite{Smith82,Smith98} +% To carry out the same procedure for a damped electrostatic kernel, we +% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. +% Many of the derivatives of the damped kernel are well known from +% Smith's early work on multipoles for the Ewald +% summation.\cite{Smith82,Smith98} -Note that increasing the value of $n$ will add additional terms to the -electrostatic potential, e.g., $f_2(r)$ includes orders up to -$(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ -functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = -f^{\prime\prime}_2(r)$, etc. These higher derivatives are required -for computing multipole energies, forces, and torques, and smooth -cutoffs of these quantities can be guaranteed as long as the number of -terms in the Taylor series exceeds the derivative order required. +% Note that increasing the value of $n$ will add additional terms to the +% electrostatic potential, e.g., $f_2(r)$ includes orders up to +% $(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ +% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = +% f^{\prime\prime}_2(r)$, etc. These higher derivatives are required +% for computing multipole energies, forces, and torques, and smooth +% cutoffs of these quantities can be guaranteed as long as the number of +% terms in the Taylor series exceeds the derivative order required. For multipole-multipole interactions, following this procedure results -in separate radial functions for each distinct orientational -contribution to the potential, and ensures that the forces and torques -from {\it each} of these contributions will vanish at the cutoff -radius. For example, the direct dipole dot product ($\mathbf{D}_{i} -\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance +in separate radial functions for each of the distinct orientational +contributions to the potential, and ensures that the forces and +torques from each of these contributions will vanish at the cutoff +radius. For example, the direct dipole dot product +($\mathbf{D}_{\bf a} +\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance dot products: \begin{equation} -U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot -\mathbf{D}_{j} \right) \frac{g_2(r)}{r} --\frac{1}{4\pi \epsilon_0} -\left( \mathbf{D}_{i} \cdot \hat{r} \right) -\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) - - \frac{g_2(r)}{r} \right) +U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( + \mathbf{D}_{\bf a} \cdot +\mathbf{D}_{\bf b} \right) v_{21}(r) + +\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right) +\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right] \end{equation} -The electrostatic forces and torques acting on the central multipole -site due to another site within cutoff sphere are derived from +For the Taylor shifted (TSF) method with the undamped kernel, +$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} + +\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4} +- \frac{6}{r r_c^2}$. In these functions, one can easily see the +connection to unmodified electrostatics as well as the smooth +transition to zero in both these functions as $r\rightarrow r_c$. The +electrostatic forces and torques acting on the central multipole due +to another site within the cutoff sphere are derived from Eq.~\ref{generic}, accounting for the appropriate number of derivatives. Complete energy, force, and torque expressions are presented in the first paper in this series (Reference @@ -401,44 +411,55 @@ A second (and significantly simpler) method involves s \subsection{Gradient-shifted force (GSF)} -A second (and significantly simpler) method involves shifting the -gradient of the raw coulomb potential for each particular multipole +A second (and conceptually simpler) method involves shifting the +gradient of the raw Coulomb potential for each particular multipole order. For example, the raw dipole-dipole potential energy may be shifted smoothly by finding the gradient for two interacting dipoles which have been projected onto the surface of the cutoff sphere without changing their relative orientation, -\begin{displaymath} -U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(r_c) - - (r_{ij}-r_c) \hat{r}_{ij} \cdot - \vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{r_c} -\end{displaymath} -Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$ -and $\mathbf{D}_{j}$, are retained at the cutoff distance (although -the signs are reversed for the dipole that has been projected onto the -cutoff sphere). In many ways, this simpler approach is closer in -spirit to the original shifted force method, in that it projects a -neutralizing multipole (and the resulting forces from this multipole) -onto a cutoff sphere. The resulting functional forms for the -potentials, forces, and torques turn out to be quite similar in form -to the Taylor-shifted approach, although the radial contributions are -significantly less perturbed by the Gradient-shifted approach than -they are in the Taylor-shifted method. +\begin{equation} +U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - +U_{D_{\bf a} D_{\bf b}}(r_c) + - (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot + \nabla U_{D_{\bf a}D_{\bf b}}(r_c). +\end{equation} +Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf + a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance +(although the signs are reversed for the dipole that has been +projected onto the cutoff sphere). In many ways, this simpler +approach is closer in spirit to the original shifted force method, in +that it projects a neutralizing multipole (and the resulting forces +from this multipole) onto a cutoff sphere. The resulting functional +forms for the potentials, forces, and torques turn out to be quite +similar in form to the Taylor-shifted approach, although the radial +contributions are significantly less perturbed by the gradient-shifted +approach than they are in the Taylor-shifted method. +For the gradient shifted (GSF) method with the undamped kernel, +$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and +$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. +Again, these functions go smoothly to zero as $r\rightarrow r_c$, and +because the Taylor expansion retains only one term, they are +significantly less perturbed than the TSF functions. + In general, the gradient shifted potential between a central multipole and any multipolar site inside the cutoff radius is given by, \begin{equation} -U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - -U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} -\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] + U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - + U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} + \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] \label{generic2} \end{equation} where the sum describes a separate force-shifting that is applied to -each orientational contribution to the energy. +each orientational contribution to the energy. In this expression, +$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles +($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ +represent the orientations the multipoles. The third term converges more rapidly than the first two terms as a function of radius, hence the contribution of the third term is very small for large cutoff radii. The force and torque derived from -equation \ref{generic2} are consistent with the energy expression and +Eq. \ref{generic2} are consistent with the energy expression and approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods can be considered generalizations of the original DSF method for higher order multipole interactions. GSF and TSF are also identical up @@ -446,7 +467,7 @@ GSF potential are presented in the first paper in this the energy, force and torque for higher order multipole-multipole interactions. Complete energy, force, and torque expressions for the GSF potential are presented in the first paper in this series -(Reference~\onlinecite{PaperI}) +(Reference~\onlinecite{PaperI}). \subsection{Shifted potential (SP) } @@ -459,7 +480,8 @@ U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) interactions with the central multipole and the image. This effectively shifts the total potential to zero at the cutoff radius, \begin{equation} -U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) +U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - +U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] \label{eq:SP} \end{equation} where the sum describes separate potential shifting that is done for @@ -471,15 +493,18 @@ The potential energy between a central multipole and o multipoles that reorient after leaving the cutoff sphere can re-enter the cutoff sphere without perturbing the total energy. -The potential energy between a central multipole and other multipolar -sites then goes smoothly to zero as $r \rightarrow r_c$. However, the -force and torque obtained from the shifted potential (SP) are -discontinuous at $r_c$. Therefore, MD simulations will still -experience energy drift while operating under the SP potential, but it -may be suitable for Monte Carlo approaches where the configurational -energy differences are the primary quantity of interest. +For the shifted potential (SP) method with the undamped kernel, +$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) = +\frac{3}{r^3} - \frac{3}{r_c^3}$. The potential energy between a +central multipole and other multipolar sites goes smoothly to zero as +$r \rightarrow r_c$. However, the force and torque obtained from the +shifted potential (SP) are discontinuous at $r_c$. MD simulations +will still experience energy drift while operating under the SP +potential, but it may be suitable for Monte Carlo approaches where the +configurational energy differences are the primary quantity of +interest. -\subsection{The Self term} +\subsection{The Self Term} In the TSF, GSF, and SP methods, a self-interaction is retained for the central multipole interacting with its own image on the surface of the cutoff sphere. This self interaction is nearly identical with the @@ -501,7 +526,7 @@ in the test-cases are given in table~\ref{tab:pars}. used the multipolar Ewald sum as a reference method for comparing energies, forces, and torques for molecular models that mimic disordered and ordered condensed-phase systems. The parameters used -in the test-cases are given in table~\ref{tab:pars}. +in the test cases are given in table~\ref{tab:pars}. \begin{table} \label{tab:pars} @@ -519,9 +544,9 @@ in the test-cases are given in table~\ref{tab:pars}. & (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu \AA\textsuperscript{2})} \\ \hline Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ - Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & 10,000 & 17.6 &17.6 & 0 \\ + Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ -Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & 10,000 & 17.6&17.6&0 \\ +Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ \ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ \ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline @@ -546,10 +571,10 @@ and have been compared with the values obtaine from th electrostatic energy, as well as the electrostatic contributions to the force and torque on each molecule. These quantities have been computed using the SP, TSF, and GSF methods, as well as a hard cutoff, -and have been compared with the values obtaine from the multipolar -Ewald sum. In Mote Carlo (MC) simulations, the energy differences +and have been compared with the values obtained from the multipolar +Ewald sum. In Monte Carlo (MC) simulations, the energy differences between two configurations is the primary quantity that governs how -the simulation proceeds. These differences are the most imporant +the simulation proceeds. These differences are the most important indicators of the reliability of a method even if the absolute energies are not exact. For each of the multipolar systems listed above, we have compared the change in electrostatic potential energy @@ -561,7 +586,7 @@ program, OpenMD,\cite{openmd} which was used for all c \subsection{Implementation} The real-space methods developed in the first paper in this series have been implemented in our group's open source molecular simulation -program, OpenMD,\cite{openmd} which was used for all calculations in +program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in this work. The complementary error function can be a relatively slow function on some processors, so all of the radial functions are precomputed on a fine grid and are spline-interpolated to provide @@ -594,33 +619,37 @@ To sample independent configurations of multipolar cry recomputed at each time step. \subsection{Model systems} -To sample independent configurations of multipolar crystals, a body -centered cubic (bcc) crystal which is a minimum energy structure for -point dipoles was generated using 3,456 molecules. The multipoles -were translationally locked in their respective crystal sites for -equilibration at a relatively low temperature (50K), so that dipoles -or quadrupoles could freely explore all accessible orientations. The -translational constraints were removed, and the crystals were -simulated for 10 ps in the microcanonical (NVE) ensemble with an -average temperature of 50 K. Configurations were sampled at equal -time intervals for the comparison of the configurational energy -differences. The crystals were not simulated close to the melting -points in order to avoid translational deformation away of the ideal -lattice geometry. +To sample independent configurations of the multipolar crystals, body +centered cubic (bcc) crystals, which exhibit the minimum energy +structures for point dipoles, were generated using 3,456 molecules. +The multipoles were translationally locked in their respective crystal +sites for equilibration at a relatively low temperature (50K) so that +dipoles or quadrupoles could freely explore all accessible +orientations. The translational constraints were then removed, the +systems were re-equilibrated, and the crystals were simulated for an +additional 10 ps in the microcanonical (NVE) ensemble with an average +temperature of 50 K. The balance between moments of inertia and +particle mass were chosen to allow orientational sampling without +significant translational motion. Configurations were sampled at +equal time intervals in order to compare configurational energy +differences. The crystals were simulated far from the melting point +in order to avoid translational deformation away of the ideal lattice +geometry. -For dipolar, quadrupolar, and mixed-multipole liquid simulations, each -system was created with 2048 molecules oriented randomly. These were +For dipolar, quadrupolar, and mixed-multipole \textit{liquid} +simulations, each system was created with 2,048 randomly-oriented +molecules. These were equilibrated at a temperature of 300K for 1 ns. +Each system was then simulated for 1 ns in the microcanonical (NVE) +ensemble. We collected 250 different configurations at equal time +intervals. For the liquid system that included ionic species, we +converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24 +\ce{Cl-} ions and re-equilibrated. After equilibration, the system was +run under the same conditions for 1 ns. A total of 250 configurations +were collected. In the following comparisons of energies, forces, and +torques, the Lennard-Jones potentials were turned off and only the +purely electrostatic quantities were compared with the same values +obtained via the Ewald sum. -system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K -temperature after equilibration. We collected 250 different -configurations in equal interval of time. For the ions mixed liquid -system, we converted 48 different molecules into 24 \ce{Na+} and 24 -\ce{Cl-} ions and equilibrated. After equilibration, the system was run -at the same environment for 1ns and 250 configurations were -collected. While comparing energies, forces, and torques with Ewald -method, Lennard-Jones potentials were turned off and purely -electrostatic interaction had been compared. - \subsection{Accuracy of Energy Differences, Forces and Torques} The pairwise summation techniques (outlined above) were evaluated for use in MC simulations by studying the energy differences between @@ -633,7 +662,7 @@ we used least square regressions analysiss for the six should be identical for all methods. Since none of the real-space methods provide exact energy differences, -we used least square regressions analysiss for the six different +we used least square regressions analysis for the six different molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF with the multipolar Ewald reference method. Unitary results for both the correlation (slope) and correlation coefficient for these @@ -644,7 +673,7 @@ also been compared by using least squares regression a configurations and 250 configurations were recorded for comparison. Each system provided 31,125 energy differences for a total of 186,750 data points. Similarly, the magnitudes of the forces and torques have -also been compared by using least squares regression analyses. In the +also been compared using least squares regression analysis. 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 @@ -764,7 +793,7 @@ model must allow for long simulation times with minima \begin{figure} \centering - \includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} + \includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps} \caption{Statistical analysis of the quality of configurational energy differences for the real-space electrostatic methods compared with the reference Ewald sum. Results with a value equal @@ -829,7 +858,7 @@ perturbations are minimal, particularly for moderate d molecules inside each other's cutoff spheres in order to correct the energy conservation issues, and this perturbation is evident in the statistics accumulated for the molecular forces. The GSF -perturbations are minimal, particularly for moderate damping and and +perturbations are minimal, particularly for moderate damping and commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows reasonable agreement in the correlation coefficient but again the systematic error in the forces is concerning if replication of Ewald @@ -837,7 +866,7 @@ forces is desired. \begin{figure} \centering - \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} + \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} \caption{Statistical analysis of the quality of the force vector magnitudes for the real-space electrostatic methods compared with the reference Ewald sum. Results with a value equal to 1 (dashed @@ -851,7 +880,7 @@ forces is desired. \begin{figure} \centering - \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} + \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps} \caption{Statistical analysis of the quality of the torque vector magnitudes for the real-space electrostatic methods compared with the reference Ewald sum. Results with a value equal to 1 (dashed @@ -886,7 +915,7 @@ directionality is shown in terms of circular variance these quantities. Force and torque vectors for all six systems were analyzed using Fisher statistics, and the quality of the vector directionality is shown in terms of circular variance -($\mathrm{Var}(\theta$) in figure +($\mathrm{Var}(\theta)$) in figure \ref{fig:slopeCorr_circularVariance}. The force and torque vectors from the new real-space methods exhibit nearly-ideal Fisher probability distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods @@ -909,7 +938,7 @@ systematically improved by varying $\alpha$ and $r_c$. \begin{figure} \centering - \includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} + \includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps} \caption{The circular variance of the direction of the force and torque vectors obtained from the real-space methods around the reference Ewald vectors. A variance equal to 0 (dashed line) @@ -941,7 +970,7 @@ conservation (drift less than $10^{-6}$ kcal / mol / n energy over time, $\delta E_1$, and the standard deviation of energy fluctuations around this drift $\delta E_0$. Both of the shifted-force methods (GSF and TSF) provide excellent energy -conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), +conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), while the hard cutoff is essentially unusable for molecular dynamics. SP provides some benefit over the hard cutoff because the energetic jumps that happen as particles leave and enter the cutoff sphere are @@ -956,18 +985,124 @@ $k$-space cutoff values. \begin{figure} \centering - \includegraphics[width=\textwidth]{newDrift.pdf} + \includegraphics[width=\textwidth]{newDrift_12.eps} \label{fig:energyDrift} \caption{Analysis of the energy conservation of the real-space electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in - energy over time and $\delta \mathrm{E}_0$ is the standard deviation - of energy fluctuations around this drift. All simulations were of a - 2000-molecule simulation of SSDQ water with 48 ionic charges at 300 - K starting from the same initial configuration. All runs utilized - the same real-space cutoff, $r_c = 12$\AA.} + energy over time (in kcal / mol / particle / ns) and $\delta + \mathrm{E}_0$ is the standard deviation of energy fluctuations + around this drift (in kcal / mol / particle). All simulations were + of a 2000-molecule simulation of SSDQ water with 48 ionic charges at + 300 K starting from the same initial configuration. All runs + utilized the same real-space cutoff, $r_c = 12$\AA.} \end{figure} +\subsection{Reproduction of Structural Features\label{sec:structure}} +One of the best tests of modified interaction potentials is the +fidelity with which they can reproduce structural features in a +liquid. One commonly-utilized measure of structural ordering is the +pair distribution function, $g(r)$, which measures local density +deviations in relation to the bulk density. In the electrostatic +approaches studied here, the short-range repulsion from the +Lennard-Jones potential is identical for the various electrostatic +methods, and since short range repulsion determines much of the local +liquid ordering, one would not expect to see any differences in +$g(r)$. Indeed, the pair distributions are essentially identical for +all of the electrostatic methods studied (for each of the different +systems under investigation). Interested readers may consult the +supplementary information for plots of these pair distribution +functions. +A direct measure of the structural features that is a more +enlightening test of the modified electrostatic methods is the average +value of the electrostatic energy $\langle U_\mathrm{elect} \rangle$ +which is obtained by sampling the liquid-state configurations +experienced by a liquid evolving entirely under the influence of the +methods being investigated. In figure \ref{fig:Uelect} we show how +$\langle U_\mathrm{elect} \rangle$ for varies with the damping parameter, +$\alpha$, for each of the methods. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{averagePotentialEnergy_r9_12.eps} +\label{fig:Uelect} +\caption{The average electrostatic potential energy, + $\langle U_\mathrm{elect} \rangle$ for the SSDQ water with ions as a function + of the damping parameter, $\alpha$, for each of the real-space + electrostatic methods. Top panel: simulations run with a real-space + cutoff, $r_c = 9$\AA. Bottom panel: the same quantity, but with a + larger cutoff, $r_c = 12$\AA.} +\end{figure} + +It is clear that moderate damping is important for converging the mean +potential energy values, particularly for the two shifted force +methods (GSF and TSF). A value of $\alpha \approx 0.18$ \AA$^{-1}$ is +sufficient to converge the SP and GSF energies with a cutoff of 12 +\AA, while shorter cutoffs require more dramatic damping ($\alpha +\approx 0.36$ \AA$^{-1}$ for $r_c = 9$ \AA). It is also clear from +fig. \ref{fig:Uelect} that it is possible to overdamp the real-space +electrostatic methods, causing the estimate of the energy to drop +below the Ewald results. + +These ``optimal'' values of the damping coefficient are slightly +larger than what were observed for DSF electrostatics for purely +point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for +$r_c = 12$\AA appears to be an excellent compromise for mixed charge +multipole systems. + +\subsection{Reproduction of Dynamic Properties\label{sec:structure}} +To test the fidelity of the electrostatic methods at reproducing +dynamics in a multipolar liquid, it is also useful to look at +transport properties, particularly the diffusion constant, +\begin{equation} +D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| + \mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle +\label{eq:diff} +\end{equation} +which measures long-time behavior and is sensitive to the forces on +the multipoles. For the soft dipolar fluid, and the SSDQ liquid +systems, the self-diffusion constants (D) were calculated from linear +fits to the long-time portion of the mean square displacement +($\langle r^{2}(t) \rangle$).\cite{Allen87} + +In addition to translational diffusion, orientational relaxation times +were calculated for comparisons with the Ewald simulations and with +experiments. These values were determined from the same 1~ns $NVE$ +trajectories used for translational diffusion by calculating the +orientational time correlation function, +\begin{equation} +C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t) + \cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle, +\label{eq:OrientCorr} +\end{equation} +where $P_l$ is the Legendre polynomial of order $l$ and +$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along +axis $\gamma$. The body-fixed reference frame used for our +orientational correlation functions has the $z$-axis running along the +dipoles, and for the SSDQ water model, the $y$-axis connects the two +implied hydrogen atoms. + +From the orientation autocorrelation functions, we can obtain time +constants for rotational relaxation either by fitting an exponential +function or by integrating the entire correlation function. These +decay times are directly comparable to water orientational relaxation +times from nuclear magnetic resonance (NMR). The relaxation constant +obtained from $C_2^y(t)$ is normally of experimental interest because +it describes the relaxation of the principle axis connecting the +hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular +portion of the dipole-dipole relaxation from a proton NMR signal and +should provide an estimate of the NMR relaxation time +constant.\cite{Impey82} + +Results for the diffusion constants and orientational relaxation times +are shown in figure \ref{fig:dynamics}. From this data, it is apparent +that the values for both $D$ and $\tau_2$ using the Ewald sum are +reproduced with high fidelity by the GSF method. + +The $\tau_2$ results in \ref{fig:dynamics} show a much greater +difference between the real-space and the Ewald results. + + \section{CONCLUSION} In the first paper in this series, we generalized the charge-neutralized electrostatic energy originally developed by Wolf @@ -1044,6 +1179,14 @@ real-space cutoff boundary. handling of energies, forces, and torques as multipoles cross the real-space cutoff boundary. +\begin{acknowledgments} + JDG acknowledges helpful discussions with Christopher + Fennell. Support for this project was provided by the National + Science Foundation under grant CHE-1362211. Computational time was + provided by the Center for Research Computing (CRC) at the + University of Notre Dame. +\end{acknowledgments} + %\bibliographystyle{aip} \newpage \bibliography{references}