ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipole/multipole_2/multipole2.tex
(Generate patch)

Comparing trunk/multipole/multipole_2/multipole2.tex (file contents):
Revision 4188 by gezelter, Sun Jun 15 16:27:42 2014 UTC vs.
Revision 4206 by gezelter, Thu Aug 7 20:53:27 2014 UTC

# Line 47 | Line 47 | preprint,
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}
# Line 117 | Line 118 | interactions,\cite{Parry:1975if,Parry:1976fq,Clarke77,
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 < but these methods were originally quite computationally
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
124 < summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
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.  The inherent periodicity in the Ewald method can also
127 < be problematic for interfacial molecular
128 < systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq}
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)$
# Line 173 | Line 174 | point charges.\cite{Fukuda:2013sf}
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
# Line 324 | Line 325 | U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1
325   expressed as the product of two multipole operators and a Coulombic
326   kernel,
327   \begin{equation}
328 < U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
328 > U_{ab}(r)= M_{a} M_{b} \frac{1}{r}  \label{kernel}.
329   \end{equation}
330 < where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is
331 < expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf
332 <    a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object
332 < $\bf a$, etc.
330 > where the multipole operator for site $a$, $M_{a}$, is
331 > expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object
332 > $a$, etc.
333  
334   % Interactions between multipoles can be expressed as higher derivatives
335   % of the bare Coulomb potential, so one way of ensuring that the forces
# Line 384 | Line 384 | radius.  For example, the direct dipole dot product
384   contributions to the potential, and ensures that the forces and
385   torques from each of these contributions will vanish at the cutoff
386   radius.  For example, the direct dipole dot product
387 < ($\mathbf{D}_{\bf a}
388 < \cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance
387 > ($\mathbf{D}_{a}
388 > \cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance
389   dot products:
390   \begin{equation}
391 < U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
392 <  \mathbf{D}_{\bf a} \cdot
393 < \mathbf{D}_{\bf b} \right) v_{21}(r) +
394 < \left( \mathbf{D}_{\bf a} \cdot \hat{r} \right)
395 < \left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right]
391 > U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
392 >  \mathbf{D}_{a} \cdot
393 > \mathbf{D}_{b} \right) v_{21}(r) +
394 > \left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right)
395 > \left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right]
396   \end{equation}
397  
398   For the Taylor shifted (TSF) method with the undamped kernel,
# Line 417 | Line 417 | U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r)
417   which have been projected onto the surface of the cutoff sphere
418   without changing their relative orientation,
419   \begin{equation}
420 < U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r) -
421 < U_{D_{\bf a} D_{\bf b}}(r_c)
422 <   - (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot
423 <  \nabla U_{D_{\bf a}D_{\bf b}}(r_c).
420 > U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
421 > U_{D_{a}D_{b}}(r_c)
422 >   - (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot
423 >  \nabla U_{D_{a}D_{b}}(r_c).
424   \end{equation}
425 < Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf
426 <  a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance
425 > Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance
426   (although the signs are reversed for the dipole that has been
427   projected onto the cutoff sphere).  In many ways, this simpler
428   approach is closer in spirit to the original shifted force method, in
# Line 444 | Line 443 | and any multipolar site inside the cutoff radius is gi
443   In general, the gradient shifted potential between a central multipole
444   and any multipolar site inside the cutoff radius is given by,
445   \begin{equation}
446 <  U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
447 <    U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}}
448 <    \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
446 > U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
447 > U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c)
448 > \hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
449   \label{generic2}
450   \end{equation}
451   where the sum describes a separate force-shifting that is applied to
452   each orientational contribution to the energy.  In this expression,
453   $\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles
454 < ($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$
454 > ($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$
455   represent the orientations the multipoles.
456  
457   The third term converges more rapidly than the first two terms as a
# Line 479 | Line 478 | U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf
478   interactions with the central multipole and the image. This
479   effectively shifts the total potential to zero at the cutoff radius,
480   \begin{equation}
481 < U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
482 < U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
481 > U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
482 > U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
483   \label{eq:SP}
484   \end{equation}          
485   where the sum describes separate potential shifting that is done for
# Line 573 | Line 572 | the simulation proceeds. These differences are the mos
572   and have been compared with the values obtained from the multipolar
573   Ewald sum.  In Monte Carlo (MC) simulations, the energy differences
574   between two configurations is the primary quantity that governs how
575 < the simulation proceeds. These differences are the most imporant
575 > the simulation proceeds. These differences are the most important
576   indicators of the reliability of a method even if the absolute
577   energies are not exact.  For each of the multipolar systems listed
578   above, we have compared the change in electrostatic potential energy
# Line 792 | Line 791 | model must allow for long simulation times with minima
791  
792   \begin{figure}
793    \centering
794 <  \includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf}
794 >  \includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps}
795    \caption{Statistical analysis of the quality of configurational
796      energy differences for the real-space electrostatic methods
797      compared with the reference Ewald sum.  Results with a value equal
# Line 865 | Line 864 | forces is desired.
864  
865   \begin{figure}
866    \centering
867 <  \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf}
867 >  \includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps}
868    \caption{Statistical analysis of the quality of the force vector
869      magnitudes for the real-space electrostatic methods compared with
870      the reference Ewald sum. Results with a value equal to 1 (dashed
# Line 879 | Line 878 | forces is desired.
878  
879   \begin{figure}
880    \centering
881 <  \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf}
881 >  \includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps}
882    \caption{Statistical analysis of the quality of the torque vector
883      magnitudes for the real-space electrostatic methods compared with
884      the reference Ewald sum. Results with a value equal to 1 (dashed
# Line 937 | Line 936 | systematically improved by varying $\alpha$ and $r_c$.
936  
937   \begin{figure}
938    \centering
939 <  \includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf}
939 >  \includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps}
940    \caption{The circular variance of the direction of the force and
941      torque vectors obtained from the real-space methods around the
942      reference Ewald vectors. A variance equal to 0 (dashed line)
# Line 957 | Line 956 | temperature of 300K.  After equilibration, this liquid
956   in this series and provides the most comprehensive test of the new
957   methods.  A liquid-phase system was created with 2000 water molecules
958   and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a
959 < temperature of 300K.  After equilibration, this liquid-phase system
960 < was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with
961 < a cutoff radius of 12\AA.  The value of the damping coefficient was
962 < also varied from the undamped case ($\alpha = 0$) to a heavily damped
963 < case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods.  A
964 < sample was also run using the multipolar Ewald sum with the same
965 < real-space cutoff.
959 > temperature of 300K.  After equilibration in the canonical (NVT)
960 > ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system
961 > was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald,
962 > Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA.  The
963 > value of the damping coefficient was also varied from the undamped
964 > case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$
965 > \AA$^{-1}$) for all of the real space methods.  A sample was also run
966 > using the multipolar Ewald sum with the same real-space cutoff.
967  
968   In figure~\ref{fig:energyDrift} we show the both the linear drift in
969   energy over time, $\delta E_1$, and the standard deviation of energy
# Line 984 | Line 984 | $k$-space cutoff values.
984  
985   \begin{figure}
986    \centering
987 <  \includegraphics[width=\textwidth]{newDrift_12.pdf}
987 >  \includegraphics[width=\textwidth]{newDrift_12.eps}
988   \label{fig:energyDrift}        
989   \caption{Analysis of the energy conservation of the real-space
990 <  electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in
991 <  energy over time (in kcal / mol / particle / ns) and $\delta
992 <  \mathrm{E}_0$ is the standard deviation of energy fluctuations
993 <  around this drift (in kcal / mol / particle).  All simulations were
994 <  of a 2000-molecule simulation of SSDQ water with 48 ionic charges at
990 >  methods. $\delta \mathrm{E}_1$ is the linear drift in energy over
991 >  time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is
992 >  the standard deviation of energy fluctuations around this drift (in
993 >  kcal / mol / particle).  Points that appear below the dashed grey
994 >  (Ewald) lines exhibit better energy conservation than commonly-used
995 >  parameters for Ewald-based electrostatics.  All simulations were of
996 >  a 2000-molecule simulation of SSDQ water with 48 ionic charges at
997    300 K starting from the same initial configuration. All runs
998    utilized the same real-space cutoff, $r_c = 12$\AA.}
999 + \end{figure}
1000 +
1001 + \subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}}
1002 + The most important test of the modified interaction potentials is the
1003 + fidelity with which they can reproduce structural features and
1004 + dynamical properties in a liquid.  One commonly-utilized measure of
1005 + structural ordering is the pair distribution function, $g(r)$, which
1006 + measures local density deviations in relation to the bulk density.  In
1007 + the electrostatic approaches studied here, the short-range repulsion
1008 + from the Lennard-Jones potential is identical for the various
1009 + electrostatic methods, and since short range repulsion determines much
1010 + of the local liquid ordering, one would not expect to see many
1011 + differences in $g(r)$.  Indeed, the pair distributions are essentially
1012 + identical for all of the electrostatic methods studied (for each of
1013 + the different systems under investigation).  An example of this
1014 + agreement for the SSDQ water/ion system is shown in
1015 + Fig. \ref{fig:gofr}.
1016 +
1017 + \begin{figure}
1018 +  \centering
1019 +  \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
1020 + \label{fig:gofr}        
1021 + \caption{The pair distribution functions, $g(r)$, for the SSDQ
1022 +  water/ion system obtained using the different real-space methods are
1023 +  essentially identical with the result from the Ewald
1024 +  treatment.}
1025   \end{figure}
1026  
1027 + There is a very slight overstructuring of the first solvation shell
1028 + when using when using TSF at lower values of the damping coefficient
1029 + ($\alpha \le 0.1$) or when using undamped GSF.  With moderate damping,
1030 + GSF and SP produce pair distributions that are identical (within
1031 + numerical noise) to their Ewald counterparts.
1032  
1033 + A structural property that is a more demanding test of modified
1034 + electrostatics is the mean value of the electrostatic energy $\langle
1035 + U_\mathrm{elect} \rangle / N$ which is obtained by sampling the
1036 + liquid-state configurations experienced by a liquid evolving entirely
1037 + under the influence of each of the methods.  In table \ref{tab:Props}
1038 + we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with
1039 + the damping parameter, $\alpha$, for each of the methods.
1040 +
1041 + As in the crystals studied in the first paper, damping is important
1042 + for converging the mean electrostatic energy values, particularly for
1043 + the two shifted force methods (GSF and TSF).  A value of $\alpha
1044 + \approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF
1045 + energies with a cutoff of 12 \AA, while shorter cutoffs require more
1046 + dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA).
1047 + Overdamping the real-space electrostatic methods occurs with $\alpha >
1048 + 0.4$, causing the estimate of the energy to drop below the Ewald
1049 + results.
1050 +
1051 + These ``optimal'' values of the damping coefficient are slightly
1052 + larger than what were observed for DSF electrostatics for purely
1053 + point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for
1054 + $r_c = 12$\AA appears to be an excellent compromise for mixed charge
1055 + multipole systems.
1056 +
1057 + To test the fidelity of the electrostatic methods at reproducing
1058 + dynamics in a multipolar liquid, it is also useful to look at
1059 + transport properties, particularly the diffusion constant,
1060 + \begin{equation}
1061 + D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left|
1062 +  \mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle
1063 + \label{eq:diff}
1064 + \end{equation}
1065 + which measures long-time behavior and is sensitive to the forces on
1066 + the multipoles.  For the soft dipolar fluid and the SSDQ liquid
1067 + systems, the self-diffusion constants (D) were calculated from linear
1068 + fits to the long-time portion of the mean square displacement,
1069 + $\langle r^{2}(t) \rangle$.\cite{Allen87}
1070 +
1071 + In addition to translational diffusion, orientational relaxation times
1072 + were calculated for comparisons with the Ewald simulations and with
1073 + experiments. These values were determined from the same 1~ns
1074 + microcanonical trajectories used for translational diffusion by
1075 + calculating the orientational time correlation function,
1076 + \begin{equation}
1077 + C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t)
1078 +                \cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle,
1079 + \label{eq:OrientCorr}
1080 + \end{equation}
1081 + where $P_l$ is the Legendre polynomial of order $l$ and
1082 + $\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis
1083 + $\gamma$ on a molecule..  Th body-fixed reference frame used for our
1084 + models has the $z$-axis running along the dipoles, and for the SSDQ
1085 + water model, the $y$-axis connects the two implied hydrogen atom
1086 + positions.  From the orientation autocorrelation functions, we can
1087 + obtain time constants for rotational relaxation either by fitting an
1088 + exponential function or by integrating the entire correlation
1089 + function.  In a good water model, these decay times would be
1090 + comparable to water orientational relaxation times from nuclear
1091 + magnetic resonance (NMR). The relaxation constant obtained from
1092 + $C_2^y(t)$ is normally of experimental interest because it describes
1093 + the relaxation of the principle axis connecting the hydrogen
1094 + atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
1095 + of the dipole-dipole relaxation from a proton NMR signal and should
1096 + provide an estimate of the NMR relaxation time constant.\cite{Impey82}
1097 +
1098 + Results for the diffusion constants and orientational relaxation times
1099 + are shown in figure \ref{tab:Props}. From this data, it is apparent
1100 + that the values for both $D$ and $\tau_2$ using the Ewald sum are
1101 + reproduced with reasonable fidelity by the GSF method.
1102 +
1103 + The $\tau_2$ results in \ref{tab:Props} show a much greater difference
1104 + between the real-space and the Ewald results.
1105 +
1106 + \begin{table}
1107 + \label{tab:Props}
1108 + \caption{Comparison of the structural and dynamic properties for the
1109 +  soft dipolar liquid test for all of the real-space methods.}
1110 + \begin{tabular}{l|c|cccc|cccc|cccc}
1111 +         & Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\
1112 + $\alpha$ (\AA$^{-1}$) & &      
1113 + 0.0 & 0.1 & 0.2 & 0.3 &
1114 + 0.0 & 0.1 & 0.2 & 0.3 &
1115 + 0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14}
1116 +
1117 + $\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\
1118 + D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)&
1119 + 470.2(6) &
1120 + 416.6(5) &
1121 + 379.6(5) &
1122 + 438.6(5) &
1123 + 476.0(6) &
1124 + 412.8(5) &
1125 + 421.1(5) &
1126 + 400.5(5) &
1127 + 437.5(6) &
1128 + 434.6(5) &
1129 + 411.4(5) &
1130 + 545.3(7) &
1131 + 459.6(6) \\
1132 + $\tau_2$ (fs) &
1133 + 1.136 &
1134 + 1.041 &
1135 + 1.064 &
1136 + 1.109 &
1137 + 1.211 &
1138 + 1.119 &
1139 + 1.039 &
1140 + 1.058 &
1141 + 1.21  &
1142 + 1.15  &
1143 + 1.172 &
1144 + 1.153 &
1145 + 1.125 \\
1146 + \end{tabular}
1147 + \end{table}
1148 +
1149 +
1150   \section{CONCLUSION}
1151   In the first paper in this series, we generalized the
1152   charge-neutralized electrostatic energy originally developed by Wolf

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines