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 2643 by gezelter, Mon Mar 20 17:32:33 2006 UTC vs.
Revision 2645 by gezelter, Mon Mar 20 21:40:49 2006 UTC

# Line 492 | Line 492 | In MC, the potential energy difference between two sub
492   techniques utilize pairwise summations of interactions between
493   particle sites, but they use these summations in different ways.
494  
495 < In MC, the potential energy difference between two subsequent
496 < configurations dictates the progression of MC sampling.  Going back to
497 < the origins of this method, the acceptance criterion for the canonical
498 < ensemble laid out by Metropolis \textit{et al.} states that a
499 < subsequent configuration is accepted if $\Delta E < 0$ or if $\xi <
500 < \exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and
501 < 1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an
502 < alternate method for handling the long-range electrostatics will
503 < ensure proper sampling from the ensemble.
495 > In MC, the potential energy difference between configurations dictates
496 > the progression of MC sampling.  Going back to the origins of this
497 > method, the acceptance criterion for the canonical ensemble laid out
498 > by Metropolis \textit{et al.} states that a subsequent configuration
499 > is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where
500 > $\xi$ is a random number between 0 and 1.\cite{Metropolis53}
501 > Maintaining the correct $\Delta E$ when using an alternate method for
502 > handling the long-range electrostatics will ensure proper sampling
503 > from the ensemble.
504  
505   In MD, the derivative of the potential governs how the system will
506   progress in time.  Consequently, the force and torque vectors on each
# Line 513 | Line 513 | The pairwise summation techniques (outlined in section
513   vectors will diverge from each other more rapidly.
514  
515   \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
516 +
517   The pairwise summation techniques (outlined in section
518   \ref{sec:ESMethods}) were evaluated for use in MC simulations by
519   studying the energy differences between conformations.  We took the
520   SPME-computed energy difference between two conformations to be the
521   correct behavior. An ideal performance by an alternative method would
522 < reproduce these energy differences exactly.  Since none of the methods
523 < provide exact energy differences, we used linear least squares
524 < regressions of the $\Delta E$ values between configurations using SPME
525 < against $\Delta E$ values using tested methods provides a quantitative
526 < comparison of this agreement.  Unitary results for both the
527 < correlation and correlation coefficient for these regressions indicate
528 < equivalent energetic results between the method under consideration
529 < and electrostatics handled using SPME.  Sample correlation plots for
530 < two alternate methods are shown in Fig. \ref{fig:linearFit}.
522 > reproduce these energy differences exactly (even if the absolute
523 > energies calculated by the methods are different).  Since none of the
524 > methods provide exact energy differences, we used linear least squares
525 > regressions of energy gap data to evaluate how closely the methods
526 > mimicked the Ewald energy gaps.  Unitary results for both the
527 > correlation (slope) and correlation coefficient for these regressions
528 > indicate perfect agreement between the alternative method and SPME.
529 > Sample correlation plots for two alternate methods are shown in
530 > Fig. \ref{fig:linearFit}.
531  
532   \begin{figure}
533   \centering
534   \includegraphics[width = \linewidth]{./dualLinear.pdf}
535 < \caption{Example least squares regressions of the configuration energy differences for SPC/E water systems. The upper plot shows a data set with a poor correlation coefficient ($R^2$), while the lower plot shows a data set with a good correlation coefficient.}
536 < \label{fig:linearFit}
535 > \caption{Example least squares regressions of the configuration energy
536 > differences for SPC/E water systems. The upper plot shows a data set
537 > with a poor correlation coefficient ($R^2$), while the lower plot
538 > shows a data set with a good correlation coefficient.}
539 > \label{fig:linearFit}
540   \end{figure}
541  
542   Each system type (detailed in section \ref{sec:RepSims}) was
543   represented using 500 independent configurations.  Additionally, we
544 < used seven different system types, so each of the alternate
544 > used seven different system types, so each of the alternative
545   (non-Ewald) electrostatic summation methods was evaluated using
546   873,250 configurational energy differences.
547  
# Line 567 | Line 571 | between those computed from the particular method and
571   investigated through measurement of the angle ($\theta$) formed
572   between those computed from the particular method and those from SPME,
573   \begin{equation}
574 < \theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right),
574 > \theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right),
575   \end{equation}
576 < where $\hat{f}_\textrm{M}$ is the unit vector pointing along the
577 < force vector computed using method $M$.  
576 > where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force
577 > vector computed using method M.
578  
579   Each of these $\theta$ values was accumulated in a distribution
580 < function, weighted by the area on the unit sphere.  Non-linear
580 > function and weighted by the area on the unit sphere.  Non-linear
581   Gaussian fits were used to measure the width of the resulting
582   distributions.
583  
584   \begin{figure}
585   \centering
586   \includegraphics[width = \linewidth]{./gaussFit.pdf}
587 < \caption{Sample fit of the angular distribution of the force vectors over all of the studied systems.  Gaussian fits were used to obtain values for the variance in force and torque vectors used in the following figure.}
587 > \caption{Sample fit of the angular distribution of the force vectors
588 > accumulated using all of the studied systems.  Gaussian fits were used
589 > to obtain values for the variance in force and torque vectors.}
590   \label{fig:gaussian}
591   \end{figure}
592  
# Line 597 | Line 603 | Evaluation of the short-time dynamics of charged syste
603   when using the reference method (SPME).
604  
605   \subsection{Short-time Dynamics}
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.  All systems were started
611 < with the same initial positions and velocities.  Simulations were
612 < performed under the microcanonical ensemble, and velocity
606 >
607 > The effects of the alternative electrostatic summation methods on the
608 > short-time dynamics of charged systems were evaluated by considering a
609 > NaCl crystal at a temperature of 1000 K.  A subset of the best
610 > performing pairwise methods was used in this comparison.  The NaCl
611 > crystal was chosen to avoid possible complications from the treatment
612 > of orientational motion in molecular systems.  All systems were
613 > started with the same initial positions and velocities.  Simulations
614 > were performed under the microcanonical ensemble, and velocity
615   autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each
616   of the trajectories,
617   \begin{equation}
# Line 617 | Line 625 | Evaluation of the long-time dynamics of charged system
625   functions was used for comparisons.
626  
627   \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
628 < Evaluation of the long-time dynamics of charged systems was performed
629 < by considering the NaCl crystal system, again while using a subset of
630 < the best performing pairwise methods.  To enhance the atomic motion,
631 < these crystals were equilibrated at 1000 K, near the experimental
632 < $T_m$ for NaCl.  Simulations were performed under the microcanonical
633 < ensemble, and velocity information was saved every 5 fs over 100 ps
634 < trajectories.  The power spectrum ($I(\omega)$) was obtained via
627 < Fourier transform of the velocity autocorrelation function
628 < \begin{equation}
629 < I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
628 >
629 > The effects of the same subset of alternative electrostatic methods on
630 > the {\it long-time} dynamics of charged systems were evaluated using
631 > the same model system (NaCl crystals at 1000K).  The power spectrum
632 > ($I(\omega)$) was obtained via Fourier transform of the velocity
633 > autocorrelation function, \begin{equation} I(\omega) =
634 > \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
635   \label{eq:powerSpec}
636   \end{equation}
637   where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
638   NaCl crystal is composed of two different atom types, the average of
639 < the two resulting power spectra was used for comparisons.
639 > the two resulting power spectra was used for comparisons. Simulations
640 > were performed under the microcanonical ensemble, and velocity
641 > information was saved every 5 fs over 100 ps trajectories.
642  
643   \subsection{Representative Simulations}\label{sec:RepSims}
644 < A variety of common and representative simulations were analyzed to
645 < determine the relative effectiveness of the pairwise summation
646 < techniques in reproducing the energetics and dynamics exhibited by
647 < SPME.  The studied systems were as follows:
644 > A variety of representative simulations were analyzed to determine the
645 > relative effectiveness of the pairwise summation techniques in
646 > reproducing the energetics and dynamics exhibited by SPME.  We wanted
647 > to span the space of modern simulations (i.e. from liquids of neutral
648 > molecules to ionic crystals), so the systems studied were:
649   \begin{enumerate}
650 < \item Liquid Water
651 < \item Crystalline Water (Ice I$_\textrm{c}$)
652 < \item NaCl Crystal
653 < \item NaCl Melt
654 < \item Low Ionic Strength Solution of NaCl in Water
655 < \item High Ionic Strength Solution of NaCl in Water
656 < \item 6 \AA\  Radius Sphere of Argon in Water
650 > \item liquid water (SPC/E),\cite{Berendsen87}
651 > \item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E),
652 > \item NaCl crystals,
653 > \item NaCl melts,
654 > \item a low ionic strength solution of NaCl in water (0.11 M),
655 > \item a high ionic strength solution of NaCl in water (1.1 M), and
656 > \item a 6 \AA\  radius sphere of Argon in water.
657   \end{enumerate}
658   By utilizing the pairwise techniques (outlined in section
659   \ref{sec:ESMethods}) in systems composed entirely of neutral groups,
660 < charged particles, and mixtures of the two, we can comment on possible
661 < system dependence and/or universal applicability of the techniques.
660 > charged particles, and mixtures of the two, we hope to discern under
661 > which conditions it will be possible to use one of the alternative
662 > summation methodologies instead of the Ewald sum.
663  
664 < Generation of the system configurations was dependent on the system
665 < type.  For the solid and liquid water configurations, configuration
666 < snapshots were taken at regular intervals from higher temperature 1000
667 < SPC/E water molecule trajectories and each equilibrated
668 < individually.\cite{Berendsen87} The solid and liquid NaCl systems
669 < consisted of 500 Na+ and 500 Cl- ions and were selected and
670 < equilibrated in the same fashion as the water systems.  For the low
671 < and high ionic strength NaCl solutions, 4 and 40 ions were first
672 < solvated in a 1000 water molecule boxes respectively.  Ion and water
673 < positions were then randomly swapped, and the resulting configurations
674 < were again equilibrated individually.  Finally, for the Argon/Water
675 < "charge void" systems, the identities of all the SPC/E waters within 6
676 < \AA\ of the center of the equilibrated water configurations were
677 < converted to argon (Fig. \ref{fig:argonSlice}).
664 > For the solid and liquid water configurations, configurations were
665 > taken at regular intervals from high temperature trajectories of 1000
666 > SPC/E water molecules.  Each configuration was equilibrated
667 > independently at a lower temperature (300~K for the liquid, 200~K for
668 > the crystal).  The solid and liquid NaCl systems consisted of 500
669 > $\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions.  Configurations for
670 > these systems were selected and equilibrated in the same manner as the
671 > water systems.  The equilibrated temperatures were 1000~K for the NaCl
672 > crystal and 7000~K for the liquid. The ionic solutions were made by
673 > solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water
674 > molecules.  Ion and water positions were then randomly swapped, and
675 > the resulting configurations were again equilibrated individually.
676 > Finally, for the Argon / Water ``charge void'' systems, the identities
677 > of all the SPC/E waters within 6 \AA\ of the center of the
678 > equilibrated water configurations were converted to argon
679 > (Fig. \ref{fig:argonSlice}).
680  
681 + These procedures guaranteed us a set of representative configurations
682 + from chemically-relevant systems sampled from an appropriate
683 + ensemble. Force field parameters for the ions and Argon were taken
684 + from the force field utilized by {\sc oopse}.\cite{Meineke05}
685 +
686   \begin{figure}
687   \centering
688   \includegraphics[width = \linewidth]{./slice.pdf}
689 < \caption{A slice from the center of a water box used in a charge void simulation.  The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.}
689 > \caption{A slice from the center of a water box used in a charge void
690 > simulation.  The darkened region represents the boundary sphere within
691 > which the water molecules were converted to argon atoms.}
692   \label{fig:argonSlice}
693   \end{figure}
694  
695 < \subsection{Electrostatic Summation Methods}\label{sec:ESMethods}
696 < Electrostatic summation method comparisons were performed using SPME,
697 < the {\sc sp} and {\sc sf} methods - both with damping
698 < parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak,
699 < moderate, and strong damping respectively), reaction field with an
700 < infinite dielectric constant, and an unmodified cutoff.  Group-based
701 < cutoffs with a fifth-order polynomial switching function were
702 < necessary for the reaction field simulations and were utilized in the
703 < SP, SF, and pure cutoff methods for comparison to the standard lack of
704 < group-based cutoffs with a hard truncation.  The SPME calculations
705 < were performed using the TINKER implementation of SPME,\cite{Ponder87}
706 < while all other method calculations were performed using the OOPSE
707 < molecular mechanics package.\cite{Meineke05}
695 > \subsection{Comparison of Summation Methods}\label{sec:ESMethods}
696 > We compared the following alternative summation methods with results
697 > from the reference method (SPME):
698 > \begin{itemize}
699 > \item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
700 > and 0.3 \AA$^{-1}$,
701 > \item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
702 > and 0.3 \AA$^{-1}$,
703 > \item reaction field with an infinite dielectric constant, and
704 > \item an unmodified cutoff.
705 > \end{itemize}
706 > Group-based cutoffs with a fifth-order polynomial switching function
707 > were utilized for the reaction field simulations.  Additionally, we
708 > investigated the use of these cutoffs with the SP, SF, and pure
709 > cutoff.  The SPME electrostatics were performed using the TINKER
710 > implementation of SPME,\cite{Ponder87} while all other method
711 > calculations were performed using the OOPSE molecular mechanics
712 > package.\cite{Meineke05} All other portions of the energy calculation
713 > (i.e. Lennard-Jones interactions) were handled in exactly the same
714 > manner across all systems and configurations.
715  
716 < These methods were additionally evaluated with three different cutoff
717 < radii (9, 12, and 15 \AA) to investigate possible cutoff radius
718 < dependence.  It should be noted that the damping parameter chosen in
719 < SPME, or so called ``Ewald Coefficient", has a significant effect on
720 < the energies and forces calculated.  Typical molecular mechanics
721 < packages default this to a value dependent on the cutoff radius and a
722 < tolerance (typically less than $1 \times 10^{-4}$ kcal/mol).  Smaller
723 < tolerances are typically associated with increased accuracy, but this
724 < usually means more time spent calculating the reciprocal-space portion
725 < of the summation.\cite{Perram88,Essmann95} The default TINKER
726 < tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
727 < calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and
728 < 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.
716 > The althernative methods were also evaluated with three different
717 > cutoff radii (9, 12, and 15 \AA).  It should be noted that the damping
718 > parameter chosen in SPME, or so called ``Ewald Coefficient'', has a
719 > significant effect on the energies and forces calculated.  Typical
720 > molecular mechanics packages set this to a value dependent on the
721 > cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$
722 > kcal/mol).  Smaller tolerances are typically associated with increased
723 > accuracy at the expense of increased time spent calculating the
724 > reciprocal-space portion of the summation.\cite{Perram88,Essmann95}
725 > The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used
726 > in all SPME calculations, resulting in Ewald Coefficients of 0.4200,
727 > 0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\
728 > respectively.
729  
730   \section{Results and Discussion}
731  
# Line 714 | Line 739 | figure \ref{fig:delE}.
739   \begin{figure}
740   \centering
741   \includegraphics[width=5.5in]{./delEplot.pdf}
742 < \caption{Statistical analysis of the quality of configurational energy differences for a given electrostatic method compared with the reference Ewald sum.  Results with a value equal to 1 (dashed line) indicate $\Delta E$ values indistinguishable from those obtained using SPME.  Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
742 > \caption{Statistical analysis of the quality of configurational energy
743 > differences for a given electrostatic method compared with the
744 > reference Ewald sum.  Results with a value equal to 1 (dashed line)
745 > indicate $\Delta E$ values indistinguishable from those obtained using
746 > SPME.  Different values of the cutoff radius are indicated with
747 > different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
748 > inverted triangles).}
749   \label{fig:delE}
750   \end{figure}
751  
752 < In this figure, it is apparent that it is unreasonable to expect
753 < realistic results using an unmodified cutoff.  This is not all that
754 < surprising since this results in large energy fluctuations as atoms or
755 < molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79}
756 < These fluctuations can be alleviated to some degree by using group
757 < based cutoffs with a switching
727 < function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff
728 < row doesn't show a significant improvement in this plot because the
729 < salt and salt solution systems contain non-neutral groups, see the
730 < accompanying supporting information for a comparison where all groups
731 < are neutral.
752 > The most striking feature of this plot is how well the Shifted Force
753 > ({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
754 > differences.  For the undamped {\sc sf} method, and the
755 > moderately-damped {\sc sp} methods, the results are nearly
756 > indistinguishable from the Ewald results.  The other common methods do
757 > significantly less well.  
758  
759 < Correcting the resulting charged cutoff sphere is one of the purposes
760 < of the damped Coulomb summation proposed by Wolf \textit{et
761 < al.},\cite{Wolf99} and this correction indeed improves the results as
762 < seen in the {\sc sp} rows.  While the undamped case of this
763 < method is a significant improvement over the pure cutoff, it still
764 < doesn't correlate that well with SPME.  Inclusion of potential damping
765 < improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
759 > The unmodified cutoff method is essentially unusable.  This is not
760 > surprising since hard cutoffs give large energy fluctuations as atoms
761 > or molecules move in and out of the cutoff
762 > radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to
763 > some degree by using group based cutoffs with a switching
764 > function.\cite{Adams79,Steinbach94,Leach01} However, we do not see
765 > significant improvement using the group-switched cutoff because the
766 > salt and salt solution systems contain non-neutral groups.  Interested
767 > readers can consult the accompanying supporting information for a
768 > comparison where all groups are neutral.
769 >
770 > For the {\sc sp} method, inclusion of potential damping improves the
771 > agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
772   an excellent correlation and quality of fit with the SPME results,
773 < particularly with a cutoff radius greater than 12 \AA .  Use of a
774 < larger damping parameter is more helpful for the shortest cutoff
775 < shown, but it has a detrimental effect on simulations with larger
776 < cutoffs.  In the {\sc sf} sets, increasing damping results in
777 < progressively poorer correlation.  Overall, the undamped case is the
778 < best performing set, as the correlation and quality of fits are
779 < consistently superior regardless of the cutoff distance.  This result
780 < is beneficial in that the undamped case is less computationally
781 < prohibitive do to the lack of complimentary error function calculation
782 < when performing the electrostatic pair interaction.  The reaction
783 < field results illustrates some of that method's limitations, primarily
752 < that it was developed for use in homogenous systems; although it does
753 < provide results that are an improvement over those from an unmodified
754 < cutoff.
773 > particularly with a cutoff radius greater than 12
774 > \AA .  Use of a larger damping parameter is more helpful for the
775 > shortest cutoff shown, but it has a detrimental effect on simulations
776 > with larger cutoffs.  
777 >
778 > In the {\sc sf} sets, increasing damping results in progressively
779 > worse correlation with Ewald.  Overall, the undamped case is the best
780 > performing set, as the correlation and quality of fits are
781 > consistently superior regardless of the cutoff distance.  The undamped
782 > case is also less computationally demanding (because no evaluation of
783 > the complementary error function is required).
784  
785 + The reaction field results illustrates some of that method's
786 + limitations, primarily that it was developed for use in homogenous
787 + systems; although it does provide results that are an improvement over
788 + those from an unmodified cutoff.
789 +
790   \subsection{Magnitudes of the Force and Torque Vectors}
791  
792   Evaluation of pairwise methods for use in Molecular Dynamics

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines