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

Comparing trunk/electrostaticMethodsPaper/electrostaticMethods.tex (file contents):
Revision 2637 by chrisfen, Sun Mar 19 02:48:19 2006 UTC vs.
Revision 2641 by chrisfen, Mon Mar 20 15:43:13 2006 UTC

# Line 77 | Line 77 | accurately incorporate their effect, and since the com
77   leading to an effect excluded from the pair interactions within a unit
78   box.  In large systems, excessively large cutoffs need to be used to
79   accurately incorporate their effect, and since the computational cost
80 < increases proportionally with the cutoff sphere, it quickly becomes an
81 < impractical task to perform these calculations.
80 > increases proportionally with the cutoff sphere, it quickly becomes
81 > very time-consuming to perform these calculations.
82  
83 + There have been many efforts to address this issue of both proper and
84 + practical handling of electrostatic interactions, and these have
85 + resulted in the availability of a variety of
86 + techniques.\cite{Roux99,Sagui99,Tobias01} These are typically
87 + classified as implicit methods (i.e., continuum dielectrics, static
88 + dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e.,
89 + Ewald summations, interaction shifting or
90 + truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
91 + reaction field type methods, fast multipole
92 + methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
93 + often preferred because they incorporate dynamic solvent molecules in
94 + the system of interest, but these methods are sometimes difficult to
95 + utilize because of their high computational cost.\cite{Roux99} In
96 + addition to this cost, there has been some question of the inherent
97 + periodicity of the explicit Ewald summation artificially influencing
98 + systems dynamics.\cite{Tobias01}
99 +
100 + In this paper, we focus on the common mixed and explicit methods of
101 + reaction filed and smooth particle mesh
102 + Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods
103 + devised by Wolf {\it et al.} which we further extend.\cite{Wolf99}
104 + These new methods for handling electrostatics are quite
105 + computationally efficient, since they involve only a simple
106 + modification to the direct pairwise sum, and they lack the added
107 + periodicity of the Ewald sum. Below, these methods are evaluated using
108 + a variety of model systems and comparison methodologies to establish
109 + their usability in molecular simulations.
110 +
111   \subsection{The Ewald Sum}
112 < The complete accumulation electrostatic interactions in a system with periodic boundary conditions (PBC) requires the consideration of the effect of all charges within a simulation box, as well as those in the periodic replicas,
112 > The complete accumulation electrostatic interactions in a system with
113 > periodic boundary conditions (PBC) requires the consideration of the
114 > effect of all charges within a simulation box, as well as those in the
115 > periodic replicas,
116   \begin{equation}
117   V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime \left[ \sum_{i=1}^N\sum_{j=1}^N \phi\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) \right],
118   \label{eq:PBCSum}
119   \end{equation}
120 < where the sum over $\mathbf{n}$ is a sum over all periodic box replicas
121 < with integer coordinates $\mathbf{n} = (l,m,n)$, and the prime indicates
122 < $i = j$ are neglected for $\mathbf{n} = 0$.\cite{deLeeuw80} Within the
123 < sum, $N$ is the number of electrostatic particles, $\mathbf{r}_{ij}$ is
124 < $\mathbf{r}_j - \mathbf{r}_i$, $L$ is the cell length, $\bm{\Omega}_{i,j}$ are
125 < the Euler angles for $i$ and $j$, and $\phi$ is Poisson's equation
126 < ($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge
127 < interactions). In the case of monopole electrostatics,
128 < eq. (\ref{eq:PBCSum}) is conditionally convergent and is discontiuous
129 < for non-neutral systems.
120 > where the sum over $\mathbf{n}$ is a sum over all periodic box
121 > replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the
122 > prime indicates $i = j$ are neglected for $\mathbf{n} =
123 > 0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
124 > particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
125 > the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
126 > $j$, and $\phi$ is Poisson's equation ($\phi(\mathbf{r}_{ij}) = q_i
127 > q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge interactions). In the
128 > case of monopole electrostatics, eq. (\ref{eq:PBCSum}) is
129 > conditionally convergent and is discontinuous for non-neutral systems.
130  
131   This electrostatic summation problem was originally studied by Ewald
132   for the case of an infinite crystal.\cite{Ewald21}. The approach he
# Line 145 | Line 176 | real-space or reciprocal space portion of the summatio
176   direct and reciprocal-space portions of the summation.  The choice of
177   the magnitude of this value allows one to select whether the
178   real-space or reciprocal space portion of the summation is an
179 < $\mathscr{O}(N^2)$ calcualtion (with the other being
179 > $\mathscr{O}(N^2)$ calculation (with the other being
180   $\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$
181   and thoughtful algorithm development, this cost can be brought down to
182   $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to
# Line 181 | Line 212 | considering the use of the Ewald summation where the i
212   artificially stabilized by the periodic replicas introduced by the
213   Ewald summation.\cite{Weber00} Thus, care ought to be taken when
214   considering the use of the Ewald summation where the intrinsic
215 < perodicity may negatively affect the system dynamics.
215 > periodicity may negatively affect the system dynamics.
216  
217  
218   \subsection{The Wolf and Zahn Methods}
# Line 199 | Line 230 | the real-space portion of the Ewald sum) to aid conver
230   and a distance-dependent damping function (identical to that seen in
231   the real-space portion of the Ewald sum) to aid convergence
232   \begin{equation}
233 < V_{\textrm{Wolf}}(r_{ij})= \frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}.
233 > V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}.
234   \label{eq:WolfPot}
235   \end{equation}
236   Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
# Line 542 | Line 573 | between those computed from the particular method and
573   investigated through measurement of the angle ($\theta$) formed
574   between those computed from the particular method and those from SPME,
575   \begin{equation}
576 < \theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method},
576 > \theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right),
577   \end{equation}
578   where $\hat{f}_\textrm{M}$ is the unit vector pointing along the
579   force vector computed using method $M$.  
# Line 572 | Line 603 | when using the reference method (SPME).
603   when using the reference method (SPME).
604  
605   \subsection{Short-time Dynamics}
606 <
607 < \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
577 < Evaluation of the long-time dynamics of charged systems was performed
578 < by considering the NaCl crystal system while using a subset of the
606 > Evaluation of the short-time dynamics of charged systems was performed
607 > by considering the 1000 K NaCl crystal system while using a subset of the
608   best performing pairwise methods.  The NaCl crystal was chosen to
609   avoid possible complications involving the propagation techniques of
610 < orientational motion in molecular systems.  To enhance the atomic
611 < motion, these crystals were equilibrated at 1000 K, near the
612 < experimental $T_m$ for NaCl.  Simulations were performed under the
613 < microcanonical ensemble, and velocity autocorrelation functions
614 < (Eq. \ref{eq:vCorr}) were computed for each of the trajectories,
610 > orientational motion in molecular systems.  All systems were started
611 > with the same initial positions and velocities.  Simulations were
612 > performed under the microcanonical ensemble, and velocity
613 > autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each
614 > of the trajectories,
615   \begin{equation}
616 < C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle.
616 > C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}.
617   \label{eq:vCorr}
618   \end{equation}
619 < Velocity autocorrelation functions require detailed short time data
620 < and long trajectories for good statistics, thus velocity information
621 < was saved every 5 fs over 100 ps trajectories.  The power spectrum
622 < ($I(\omega)$) is obtained via Fourier transform of the autocorrelation
623 < function
619 > Velocity autocorrelation functions require detailed short time data,
620 > thus velocity information was saved every 2 fs over 10 ps
621 > trajectories. Because the NaCl crystal is composed of two different
622 > atom types, the average of the two resulting velocity autocorrelation
623 > functions was used for comparisons.
624 >
625 > \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
626 > Evaluation of the long-time dynamics of charged systems was performed
627 > by considering the NaCl crystal system, again while using a subset of
628 > the best performing pairwise methods.  To enhance the atomic motion,
629 > these crystals were equilibrated at 1000 K, near the experimental
630 > $T_m$ for NaCl.  Simulations were performed under the microcanonical
631 > ensemble, and velocity information was saved every 5 fs over 100 ps
632 > trajectories.  The power spectrum ($I(\omega)$) was obtained via
633 > Fourier transform of the velocity autocorrelation function
634   \begin{equation}
635   I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
636   \label{eq:powerSpec}
637   \end{equation}
638 < where the frequency, $\omega=0,\ 1,\ ...,\ N-1$.
638 > where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
639 > NaCl crystal is composed of two different atom types, the average of
640 > the two resulting power spectra was used for comparisons.
641  
642   \subsection{Representative Simulations}\label{sec:RepSims}
643   A variety of common and representative simulations were analyzed to
# Line 620 | Line 661 | snapshots were taken at regular intervals from higher
661   Generation of the system configurations was dependent on the system
662   type.  For the solid and liquid water configurations, configuration
663   snapshots were taken at regular intervals from higher temperature 1000
664 < SPC/E water molecule trajectories and each equilibrated individually.
665 < The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl-
666 < ions and were selected and equilibrated in the same fashion as the
667 < water systems.  For the low and high ionic strength NaCl solutions, 4
668 < and 40 ions were first solvated in a 1000 water molecule boxes
669 < respectively.  Ion and water positions were then randomly swapped, and
670 < the resulting configurations were again equilibrated individually.
671 < Finally, for the Argon/Water "charge void" systems, the identities of
672 < all the SPC/E waters within 6 \AA\ of the center of the equilibrated
673 < water configurations were converted to argon
674 < (Fig. \ref{fig:argonSlice}).
664 > SPC/E water molecule trajectories and each equilibrated
665 > individually.\cite{Berendsen87} The solid and liquid NaCl systems
666 > consisted of 500 Na+ and 500 Cl- ions and were selected and
667 > equilibrated in the same fashion as the water systems.  For the low
668 > and high ionic strength NaCl solutions, 4 and 40 ions were first
669 > solvated in a 1000 water molecule boxes respectively.  Ion and water
670 > positions were then randomly swapped, and the resulting configurations
671 > were again equilibrated individually.  Finally, for the Argon/Water
672 > "charge void" systems, the identities of all the SPC/E waters within 6
673 > \AA\ of the center of the equilibrated water configurations were
674 > converted to argon (Fig. \ref{fig:argonSlice}).
675  
676   \begin{figure}
677   \centering
# Line 685 | Line 726 | realistic results using an unmodified cutoff.  This is
726  
727   In this figure, it is apparent that it is unreasonable to expect
728   realistic results using an unmodified cutoff.  This is not all that
729 < surprising since this results in large energy fluctuations as atoms
730 < move in and out of the cutoff radius.  These fluctuations can be
731 < alleviated to some degree by using group based cutoffs with a
732 < switching function.\cite{Steinbach94} The Group Switch Cutoff row
733 < doesn't show a significant improvement in this plot because the salt
734 < and salt solution systems contain non-neutral groups, see the
729 > surprising since this results in large energy fluctuations as atoms or
730 > molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79}
731 > These fluctuations can be alleviated to some degree by using group
732 > based cutoffs with a switching
733 > function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff
734 > row doesn't show a significant improvement in this plot because the
735 > salt and salt solution systems contain non-neutral groups, see the
736   accompanying supporting information for a comparison where all groups
737   are neutral.
738  
739   Correcting the resulting charged cutoff sphere is one of the purposes
740   of the damped Coulomb summation proposed by Wolf \textit{et
741   al.},\cite{Wolf99} and this correction indeed improves the results as
742 < seen in the Shifted-Potental rows.  While the undamped case of this
742 > seen in the {\sc sp} rows.  While the undamped case of this
743   method is a significant improvement over the pure cutoff, it still
744   doesn't correlate that well with SPME.  Inclusion of potential damping
745   improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
# Line 882 | Line 924 | unnecessary when using the {\sc sf} method.
924   up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably
925   unnecessary when using the {\sc sf} method.
926  
927 < \subsection{Collective Motion: Power Spectra of NaCl Crystals}
927 > \subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}
928  
929   In the previous studies using a {\sc sf} variant of the damped
930   Wolf coulomb potential, the structure and dynamics of water were
# Line 897 | Line 939 | summation methods from the above results.
939  
940   \begin{figure}
941   \centering
942 + \includegraphics[width = \linewidth]{./vCorrPlot.pdf}
943 + \caption{Velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first trough. The times to first collision are nearly identical, but the differences can be seen in the peaks and troughs, where the undamped to weakly damped methods are stiffer than the moderately damped and SPME methods.}
944 + \label{fig:vCorrPlot}
945 + \end{figure}
946 +
947 + The short-time decays through the first collision are nearly identical
948 + in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the
949 + functions show how the methods differ.  The undamped {\sc sf} method
950 + has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher
951 + peaks than any of the other methods.  As the damping function is
952 + increased, these peaks are smoothed out, and approach the SPME
953 + curve. The damping acts as a distance dependent Gaussian screening of
954 + the point charges for the pairwise summation methods; thus, the
955 + collisions are more elastic in the undamped {\sc sf} potential, and the
956 + stiffness of the potential is diminished as the electrostatic
957 + interactions are softened by the damping function.  With $\alpha$
958 + values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are
959 + nearly identical and track the SPME features quite well.  This is not
960 + too surprising in that the differences between the {\sc sf} and {\sc
961 + sp} potentials are mitigated with increased damping.  However, this
962 + appears to indicate that once damping is utilized, the form of the
963 + potential seems to play a lesser role in the crystal dynamics.
964 +
965 + \subsection{Collective Motion: Power Spectra of NaCl Crystals}
966 +
967 + The short time dynamics were extended to evaluate how the differences
968 + between the methods affect the collective long-time motion.  The same
969 + electrostatic summation methods were used as in the short time
970 + velocity autocorrelation function evaluation, but the trajectories
971 + were sampled over a much longer time. The power spectra of the
972 + resulting velocity autocorrelation functions were calculated and are
973 + displayed in figure \ref{fig:methodPS}.
974 +
975 + \begin{figure}
976 + \centering
977   \includegraphics[width = \linewidth]{./spectraSquare.pdf}
978   \caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2).  Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude.  The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.}
979   \label{fig:methodPS}
980   \end{figure}
981  
982 < Figure \ref{fig:methodPS} shows the power spectra for the NaCl
983 < crystals (from averaged Na and Cl ion velocity autocorrelation
984 < functions) using the stated electrostatic summation methods.  While
985 < high frequency peaks of all the spectra overlap, showing the same
986 < general features, the low frequency region shows how the summation
987 < methods differ.  Considering the low-frequency inset (expanded in the
988 < upper frame of figure \ref{fig:dampInc}), at frequencies below 100
989 < cm$^{-1}$, the correlated motions are blue-shifted when using undamped
990 < or weakly damped {\sc sf}.  When using moderate damping ($\alpha
991 < = 0.2$ \AA$^{-1}$) both the {\sc sf} and {\sc sp}
992 < methods give near identical correlated motion behavior as the Ewald
993 < method (which has a damping value of 0.3119).  The damping acts as a
994 < distance dependent Gaussian screening of the point charges for the
995 < pairwise summation methods.  This weakening of the electrostatic
996 < interaction with distance explains why the long-ranged correlated
920 < motions are at lower frequencies for the moderately damped methods
921 < than for undamped or weakly damped methods.  To see this effect more
922 < clearly, we show how damping strength affects a simple real-space
923 < electrostatic potential,
982 > While high frequency peaks of the spectra in this figure overlap,
983 > showing the same general features, the low frequency region shows how
984 > the summation methods differ.  Considering the low-frequency inset
985 > (expanded in the upper frame of figure \ref{fig:dampInc}), at
986 > frequencies below 100 cm$^{-1}$, the correlated motions are
987 > blue-shifted when using undamped or weakly damped {\sc sf}.  When
988 > using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf}
989 > and {\sc sp} methods give near identical correlated motion behavior as
990 > the Ewald method (which has a damping value of 0.3119).  This
991 > weakening of the electrostatic interaction with increased damping
992 > explains why the long-ranged correlated motions are at lower
993 > frequencies for the moderately damped methods than for undamped or
994 > weakly damped methods.  To see this effect more clearly, we show how
995 > damping strength alone affects a simple real-space electrostatic
996 > potential,
997   \begin{equation}
998   V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r),
999   \end{equation}
# Line 935 | Line 1008 | blue-shifted such that the lowest frequency peak resid
1008   shift to higher frequency in exponential fashion.  Though not shown,
1009   the spectrum for the simple undamped electrostatic potential is
1010   blue-shifted such that the lowest frequency peak resides near 325
1011 < cm$^{-1}$.  In light of these results, the undamped {\sc sf}
1012 < method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is
1013 < quite respectable; however, it appears as though moderate damping is
1014 < required for accurate reproduction of crystal dynamics.
1011 > cm$^{-1}$.  In light of these results, the undamped {\sc sf} method
1012 > producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite
1013 > respectable and shows that the shifted force procedure accounts for
1014 > most of the effect afforded through use of the Ewald summation.
1015 > However, it appears as though moderate damping is required for
1016 > accurate reproduction of crystal dynamics.
1017   \begin{figure}
1018   \centering
1019   \includegraphics[width = \linewidth]{./comboSquare.pdf}
# Line 982 | Line 1057 | today, the Ewald summation may no longer be required t
1057   standard by which these simple pairwise sums are judged.  However,
1058   these results do suggest that in the typical simulations performed
1059   today, the Ewald summation may no longer be required to obtain the
1060 < level of accuracy most researcher have come to expect
1060 > level of accuracy most researchers have come to expect
1061  
1062   \section{Acknowledgments}
1063   \newpage

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines