47 |
|
|
48 |
|
%\preprint{AIP/123-QED} |
49 |
|
|
50 |
< |
\title{Real space alternatives to the Ewald Sum. II. Comparison of Methods} |
50 |
> |
\title{Real space electrostatics for multipoles. II. Comparisons with |
51 |
> |
the Ewald Sum} |
52 |
|
|
53 |
|
\author{Madan Lamichhane} |
54 |
|
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
64 |
|
\date{\today} |
65 |
|
|
66 |
|
\begin{abstract} |
67 |
< |
We report on tests of the real-space shifted potential (SP), |
68 |
< |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
69 |
< |
for multipole interactions developed in the first paper in this |
70 |
< |
series, using the multipolar Ewald sum as a reference method. The |
71 |
< |
tests were carried out in a variety of condensed-phase environments |
72 |
< |
designed to test up to quadrupole-quadrupole interactions. |
73 |
< |
Comparisons of the energy differences between configurations, |
74 |
< |
molecular forces, and torques were used to analyze how well the |
75 |
< |
real-space models perform relative to the more computationally |
76 |
< |
expensive Ewald treatment. We have also investigated the energy |
77 |
< |
conservation properties of the new methods in molecular dynamics |
78 |
< |
simulations. The SP method shows excellent agreement with |
79 |
< |
configurational energy differences, forces, and torques, and would |
80 |
< |
be suitable for use in Monte Carlo calculations. Of the two new |
81 |
< |
shifted-force methods, the GSF approach shows the best agreement |
82 |
< |
with Ewald-derived energies, forces, and torques and also exhibits |
83 |
< |
energy conservation properties that make it an excellent choice for |
84 |
< |
efficient computation of electrostatic interactions in molecular |
84 |
< |
dynamics simulations. |
67 |
> |
We report on tests of the shifted potential (SP), gradient shifted |
68 |
> |
force (GSF), and Taylor shifted force (TSF) real-space methods for |
69 |
> |
multipole interactions developed in the first paper in this series, |
70 |
> |
using the multipolar Ewald sum as a reference method. The tests were |
71 |
> |
carried out in a variety of condensed-phase environments designed to |
72 |
> |
test up to quadrupole-quadrupole interactions. Comparisons of the |
73 |
> |
energy differences between configurations, molecular forces, and |
74 |
> |
torques were used to analyze how well the real-space models perform |
75 |
> |
relative to the more computationally expensive Ewald treatment. We |
76 |
> |
have also investigated the energy conservation properties of the new |
77 |
> |
methods in molecular dynamics simulations. The SP method shows |
78 |
> |
excellent agreement with configurational energy differences, forces, |
79 |
> |
and torques, and would be suitable for use in Monte Carlo |
80 |
> |
calculations. Of the two new shifted-force methods, the GSF |
81 |
> |
approach shows the best agreement with Ewald-derived energies, |
82 |
> |
forces, and torques and also exhibits energy conservation properties |
83 |
> |
that make it an excellent choice for efficient computation of |
84 |
> |
electrostatic interactions in molecular dynamics simulations. |
85 |
|
\end{abstract} |
86 |
|
|
87 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
100 |
|
the conditionally convergent electrostatic energy is converted into |
101 |
|
two absolutely convergent contributions, one which is carried out in |
102 |
|
real space with a cutoff radius, and one in reciprocal |
103 |
< |
space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
103 |
> |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
104 |
|
|
105 |
|
When carried out as originally formulated, the reciprocal-space |
106 |
|
portion of the Ewald sum exhibits relatively poor computational |
107 |
< |
scaling, making it prohibitive for large systems. By utilizing |
108 |
< |
particle meshes and three dimensional fast Fourier transforms (FFT), |
109 |
< |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
107 |
> |
scaling, making it prohibitive for large systems. By utilizing a |
108 |
> |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
109 |
> |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
110 |
|
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
111 |
|
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
112 |
|
\log |
113 |
< |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
113 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
114 |
|
|
115 |
|
Because of the artificial periodicity required for the Ewald sum, |
116 |
|
interfacial molecular systems such as membranes and liquid-vapor |
117 |
< |
interfaces require modifications to the |
118 |
< |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
119 |
< |
Parry's extension of the three dimensional Ewald sum is appropriate |
120 |
< |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
121 |
< |
were developed to handle two-dimensional (2D) electrostatic |
122 |
< |
interactions in interfacial systems have not seen similar |
123 |
< |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
124 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
125 |
< |
with system size. The inherent periodicity in the Ewald’s method can |
126 |
< |
also be problematic for interfacial molecular |
127 |
< |
systems.\cite{Fennell:2006lq} |
117 |
> |
interfaces require modifications to the method. Parry's extension of |
118 |
> |
the three dimensional Ewald sum is appropriate for slab |
119 |
> |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
120 |
> |
developed to handle two-dimensional (2-D) electrostatic |
121 |
> |
interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
122 |
> |
These methods were originally quite computationally |
123 |
> |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
124 |
> |
efforts that reduced the computational cost of 2-D lattice summations, |
125 |
> |
bringing them more in line with the scaling for the full 3-D |
126 |
> |
treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The |
127 |
> |
inherent periodicity required by the Ewald method can also be |
128 |
> |
problematic in a number of protein/solvent and ionic solution |
129 |
> |
environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
130 |
|
|
131 |
|
\subsection{Real-space methods} |
132 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
174 |
|
|
175 |
|
\begin{figure} |
176 |
|
\centering |
177 |
< |
\includegraphics[width=\linewidth]{schematic.pdf} |
177 |
> |
\includegraphics[width=\linewidth]{schematic.eps} |
178 |
|
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
179 |
|
charges (in the smaller grey circle), so interactions are |
180 |
|
effectively charge-multipole at longer distances. With hard |
574 |
|
and have been compared with the values obtained from the multipolar |
575 |
|
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
576 |
|
between two configurations is the primary quantity that governs how |
577 |
< |
the simulation proceeds. These differences are the most imporant |
577 |
> |
the simulation proceeds. These differences are the most important |
578 |
|
indicators of the reliability of a method even if the absolute |
579 |
|
energies are not exact. For each of the multipolar systems listed |
580 |
|
above, we have compared the change in electrostatic potential energy |
586 |
|
\subsection{Implementation} |
587 |
|
The real-space methods developed in the first paper in this series |
588 |
|
have been implemented in our group's open source molecular simulation |
589 |
< |
program, OpenMD,\cite{openmd} which was used for all calculations in |
589 |
> |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
590 |
|
this work. The complementary error function can be a relatively slow |
591 |
|
function on some processors, so all of the radial functions are |
592 |
|
precomputed on a fine grid and are spline-interpolated to provide |
793 |
|
|
794 |
|
\begin{figure} |
795 |
|
\centering |
796 |
< |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
796 |
> |
\includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
797 |
|
\caption{Statistical analysis of the quality of configurational |
798 |
|
energy differences for the real-space electrostatic methods |
799 |
|
compared with the reference Ewald sum. Results with a value equal |
866 |
|
|
867 |
|
\begin{figure} |
868 |
|
\centering |
869 |
< |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
869 |
> |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} |
870 |
|
\caption{Statistical analysis of the quality of the force vector |
871 |
|
magnitudes for the real-space electrostatic methods compared with |
872 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
880 |
|
|
881 |
|
\begin{figure} |
882 |
|
\centering |
883 |
< |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
883 |
> |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps} |
884 |
|
\caption{Statistical analysis of the quality of the torque vector |
885 |
|
magnitudes for the real-space electrostatic methods compared with |
886 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
938 |
|
|
939 |
|
\begin{figure} |
940 |
|
\centering |
941 |
< |
\includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} |
941 |
> |
\includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps} |
942 |
|
\caption{The circular variance of the direction of the force and |
943 |
|
torque vectors obtained from the real-space methods around the |
944 |
|
reference Ewald vectors. A variance equal to 0 (dashed line) |
985 |
|
|
986 |
|
\begin{figure} |
987 |
|
\centering |
988 |
< |
\includegraphics[width=\textwidth]{newDrift_12.pdf} |
988 |
> |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
989 |
|
\label{fig:energyDrift} |
990 |
|
\caption{Analysis of the energy conservation of the real-space |
991 |
|
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
997 |
|
utilized the same real-space cutoff, $r_c = 12$\AA.} |
998 |
|
\end{figure} |
999 |
|
|
1000 |
+ |
\subsection{Reproduction of Structural Features\label{sec:structure}} |
1001 |
+ |
One of the best tests of modified interaction potentials is the |
1002 |
+ |
fidelity with which they can reproduce structural features in a |
1003 |
+ |
liquid. One commonly-utilized measure of structural ordering is the |
1004 |
+ |
pair distribution function, $g(r)$, which measures local density |
1005 |
+ |
deviations in relation to the bulk density. In the electrostatic |
1006 |
+ |
approaches studied here, the short-range repulsion from the |
1007 |
+ |
Lennard-Jones potential is identical for the various electrostatic |
1008 |
+ |
methods, and since short range repulsion determines much of the local |
1009 |
+ |
liquid ordering, one would not expect to see any differences in |
1010 |
+ |
$g(r)$. Indeed, the pair distributions are essentially identical for |
1011 |
+ |
all of the electrostatic methods studied (for each of the different |
1012 |
+ |
systems under investigation). Interested readers may consult the |
1013 |
+ |
supplementary information for plots of these pair distribution |
1014 |
+ |
functions. |
1015 |
|
|
1016 |
+ |
A direct measure of the structural features that is a more |
1017 |
+ |
enlightening test of the modified electrostatic methods is the average |
1018 |
+ |
value of the electrostatic energy $\langle U_\mathrm{elect} \rangle$ |
1019 |
+ |
which is obtained by sampling the liquid-state configurations |
1020 |
+ |
experienced by a liquid evolving entirely under the influence of the |
1021 |
+ |
methods being investigated. In figure \ref{fig:Uelect} we show how |
1022 |
+ |
$\langle U_\mathrm{elect} \rangle$ for varies with the damping parameter, |
1023 |
+ |
$\alpha$, for each of the methods. |
1024 |
+ |
|
1025 |
+ |
\begin{figure} |
1026 |
+ |
\centering |
1027 |
+ |
\includegraphics[width=\textwidth]{averagePotentialEnergy_r9_12.eps} |
1028 |
+ |
\label{fig:Uelect} |
1029 |
+ |
\caption{The average electrostatic potential energy, |
1030 |
+ |
$\langle U_\mathrm{elect} \rangle$ for the SSDQ water with ions as a function |
1031 |
+ |
of the damping parameter, $\alpha$, for each of the real-space |
1032 |
+ |
electrostatic methods. Top panel: simulations run with a real-space |
1033 |
+ |
cutoff, $r_c = 9$\AA. Bottom panel: the same quantity, but with a |
1034 |
+ |
larger cutoff, $r_c = 12$\AA.} |
1035 |
+ |
\end{figure} |
1036 |
+ |
|
1037 |
+ |
It is clear that moderate damping is important for converging the mean |
1038 |
+ |
potential energy values, particularly for the two shifted force |
1039 |
+ |
methods (GSF and TSF). A value of $\alpha \approx 0.18$ \AA$^{-1}$ is |
1040 |
+ |
sufficient to converge the SP and GSF energies with a cutoff of 12 |
1041 |
+ |
\AA, while shorter cutoffs require more dramatic damping ($\alpha |
1042 |
+ |
\approx 0.36$ \AA$^{-1}$ for $r_c = 9$ \AA). It is also clear from |
1043 |
+ |
fig. \ref{fig:Uelect} that it is possible to overdamp the real-space |
1044 |
+ |
electrostatic methods, causing the estimate of the energy to drop |
1045 |
+ |
below the Ewald results. |
1046 |
+ |
|
1047 |
+ |
These ``optimal'' values of the damping coefficient are slightly |
1048 |
+ |
larger than what were observed for DSF electrostatics for purely |
1049 |
+ |
point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for |
1050 |
+ |
$r_c = 12$\AA appears to be an excellent compromise for mixed charge |
1051 |
+ |
multipole systems. |
1052 |
+ |
|
1053 |
+ |
\subsection{Reproduction of Dynamic Properties\label{sec:structure}} |
1054 |
+ |
To test the fidelity of the electrostatic methods at reproducing |
1055 |
+ |
dynamics in a multipolar liquid, it is also useful to look at |
1056 |
+ |
transport properties, particularly the diffusion constant, |
1057 |
+ |
\begin{equation} |
1058 |
+ |
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
1059 |
+ |
\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle |
1060 |
+ |
\label{eq:diff} |
1061 |
+ |
\end{equation} |
1062 |
+ |
which measures long-time behavior and is sensitive to the forces on |
1063 |
+ |
the multipoles. For the soft dipolar fluid, and the SSDQ liquid |
1064 |
+ |
systems, the self-diffusion constants (D) were calculated from linear |
1065 |
+ |
fits to the long-time portion of the mean square displacement |
1066 |
+ |
($\langle r^{2}(t) \rangle$).\cite{Allen87} |
1067 |
+ |
|
1068 |
+ |
In addition to translational diffusion, orientational relaxation times |
1069 |
+ |
were calculated for comparisons with the Ewald simulations and with |
1070 |
+ |
experiments. These values were determined from the same 1~ns $NVE$ |
1071 |
+ |
trajectories used for translational diffusion by calculating the |
1072 |
+ |
orientational time correlation function, |
1073 |
+ |
\begin{equation} |
1074 |
+ |
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t) |
1075 |
+ |
\cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle, |
1076 |
+ |
\label{eq:OrientCorr} |
1077 |
+ |
\end{equation} |
1078 |
+ |
where $P_l$ is the Legendre polynomial of order $l$ and |
1079 |
+ |
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along |
1080 |
+ |
axis $\gamma$. The body-fixed reference frame used for our |
1081 |
+ |
orientational correlation functions has the $z$-axis running along the |
1082 |
+ |
dipoles, and for the SSDQ water model, the $y$-axis connects the two |
1083 |
+ |
implied hydrogen atoms. |
1084 |
+ |
|
1085 |
+ |
From the orientation autocorrelation functions, we can obtain time |
1086 |
+ |
constants for rotational relaxation either by fitting an exponential |
1087 |
+ |
function or by integrating the entire correlation function. These |
1088 |
+ |
decay times are directly comparable to water orientational relaxation |
1089 |
+ |
times from nuclear magnetic resonance (NMR). The relaxation constant |
1090 |
+ |
obtained from $C_2^y(t)$ is normally of experimental interest because |
1091 |
+ |
it describes the relaxation of the principle axis connecting the |
1092 |
+ |
hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular |
1093 |
+ |
portion of the dipole-dipole relaxation from a proton NMR signal and |
1094 |
+ |
should provide an estimate of the NMR relaxation time |
1095 |
+ |
constant.\cite{Impey82} |
1096 |
+ |
|
1097 |
+ |
Results for the diffusion constants and orientational relaxation times |
1098 |
+ |
are shown in figure \ref{fig:dynamics}. From this data, it is apparent |
1099 |
+ |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
1100 |
+ |
reproduced with high fidelity by the GSF method. |
1101 |
+ |
|
1102 |
+ |
The $\tau_2$ results in \ref{fig:dynamics} show a much greater |
1103 |
+ |
difference between the real-space and the Ewald results. |
1104 |
+ |
|
1105 |
+ |
|
1106 |
|
\section{CONCLUSION} |
1107 |
|
In the first paper in this series, we generalized the |
1108 |
|
charge-neutralized electrostatic energy originally developed by Wolf |