| 140 |
|
\end{split} |
| 141 |
|
\label{eq:EwaldSum} |
| 142 |
|
\end{equation} |
| 143 |
< |
where $\alpha$ is a damping parameter, or separation constant, with |
| 144 |
< |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are |
| 145 |
< |
equal to $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the |
| 146 |
< |
dielectric constant of the surrounding medium. The final two terms of |
| 143 |
> |
where $\alpha$ is the damping or convergence parameter with units of |
| 144 |
> |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 145 |
> |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 146 |
> |
constant of the surrounding medium. The final two terms of |
| 147 |
|
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 148 |
|
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 149 |
|
dipolar term was neglected in early applications in molecular |
| 154 |
|
system is said to be using conducting (or ``tin-foil'') boundary |
| 155 |
|
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 156 |
|
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 157 |
< |
time. Initially, due to the small sizes of the systems that could be |
| 158 |
< |
feasibly simulated, the entire simulation box was replicated to |
| 159 |
< |
convergence. In more modern simulations, the simulation boxes have |
| 160 |
< |
grown large enough that a real-space cutoff could potentially give |
| 161 |
< |
convergent behavior. Indeed, it has often been observed that the |
| 162 |
< |
reciprocal-space portion of the Ewald sum can be vanishingly |
| 163 |
< |
small compared to the real-space portion.\cite{XXX} |
| 157 |
> |
time. Initially, due to the small system sizes that could be |
| 158 |
> |
simulated feasibly, the entire simulation box was replicated to |
| 159 |
> |
convergence. In more modern simulations, the systems have grown large |
| 160 |
> |
enough that a real-space cutoff could potentially give convergent |
| 161 |
> |
behavior. Indeed, it has been observed that with the choice of a |
| 162 |
> |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 163 |
> |
rapidly convergent and small relative to the real-space |
| 164 |
> |
portion.\cite{Karasawa89,Kolafa92} |
| 165 |
|
|
| 166 |
|
\begin{figure} |
| 167 |
|
\centering |
| 168 |
< |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 169 |
< |
\caption{How the application of the Ewald summation has changed with |
| 170 |
< |
the increase in computer power. Initially, only small numbers of |
| 171 |
< |
particles could be studied, and the Ewald sum acted to replicate the |
| 172 |
< |
unit cell charge distribution out to convergence. Now, much larger |
| 173 |
< |
systems of charges are investigated with fixed distance cutoffs. The |
| 173 |
< |
calculated structure factor is used to sum out to great distance, and |
| 174 |
< |
a surrounding dielectric term is included.} |
| 168 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression2.pdf} |
| 169 |
> |
\caption{The change in the application of the Ewald sum with |
| 170 |
> |
increasing computational power. Initially, only small systems could |
| 171 |
> |
be studied, and the Ewald sum replicated the simulation box to |
| 172 |
> |
convergence. Now, much larger systems of charges are investigated |
| 173 |
> |
with fixed-distance cutoffs.} |
| 174 |
|
\label{fig:ewaldTime} |
| 175 |
|
\end{figure} |
| 176 |
|
|
| 177 |
|
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 178 |
< |
separation constant $(\alpha)$ plays an important role in balancing |
| 178 |
> |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 179 |
|
the computational cost between the direct and reciprocal-space |
| 180 |
|
portions of the summation. The choice of this value allows one to |
| 181 |
|
select whether the real-space or reciprocal space portion of the |
| 573 |
|
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 574 |
|
\end{equation} |
| 575 |
|
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 576 |
< |
vector computed using method M. |
| 577 |
< |
|
| 578 |
< |
Each of these $\theta$ values was accumulated in a distribution |
| 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 |
| 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 |
< |
|
| 593 |
< |
Figure \ref{fig:gaussian} shows an example distribution with applied |
| 594 |
< |
non-linear fits. The solid line is a Gaussian profile, while the |
| 595 |
< |
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 596 |
< |
Lorentzian. Since this distribution is a measure of angular error |
| 576 |
> |
vector computed using method M. Each of these $\theta$ values was |
| 577 |
> |
accumulated in a distribution function and weighted by the area on the |
| 578 |
> |
unit sphere. Since this distribution is a measure of angular error |
| 579 |
|
between two different electrostatic summation methods, there is no |
| 580 |
< |
{\it a priori} reason for the profile to adhere to any specific shape. |
| 581 |
< |
Gaussian fits was used to compare all the tested methods. The |
| 582 |
< |
variance ($\sigma^2$) was extracted from each of these fits and was |
| 583 |
< |
used to compare distribution widths. Values of $\sigma^2$ near zero |
| 584 |
< |
indicate vector directions indistinguishable from those calculated |
| 585 |
< |
when using the reference method (SPME). |
| 580 |
> |
{\it a priori} reason for the profile to adhere to any specific |
| 581 |
> |
shape. Thus, gaussian fits were used to measure the width of the |
| 582 |
> |
resulting 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 |
| 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 |
> |
% |
| 593 |
> |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 594 |
> |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 595 |
> |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 596 |
> |
%Lorentzian. |
| 597 |
> |
%Since this distribution is a measure of angular error between two |
| 598 |
> |
%different electrostatic summation methods, there is no {\it a priori} |
| 599 |
> |
%reason for the profile to adhere to any specific shape. |
| 600 |
> |
%Gaussian fits was used to compare all the tested methods. |
| 601 |
> |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 602 |
> |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 603 |
> |
zero indicate vector directions indistinguishable from those |
| 604 |
> |
calculated when using the reference method (SPME). |
| 605 |
|
|
| 606 |
|
\subsection{Short-time Dynamics} |
| 607 |
|
|
| 676 |
|
the resulting configurations were again equilibrated individually. |
| 677 |
|
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 678 |
|
of all the SPC/E waters within 6 \AA\ of the center of the |
| 679 |
< |
equilibrated water configurations were converted to argon |
| 680 |
< |
(Fig. \ref{fig:argonSlice}). |
| 679 |
> |
equilibrated water configurations were converted to argon. |
| 680 |
> |
%(Fig. \ref{fig:argonSlice}). |
| 681 |
|
|
| 682 |
|
These procedures guaranteed us a set of representative configurations |
| 683 |
< |
from chemically-relevant systems sampled from an appropriate |
| 684 |
< |
ensemble. Force field parameters for the ions and Argon were taken |
| 683 |
> |
from chemically-relevant systems sampled from appropriate |
| 684 |
> |
ensembles. Force field parameters for the ions and Argon were taken |
| 685 |
|
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 686 |
|
|
| 687 |
< |
\begin{figure} |
| 688 |
< |
\centering |
| 689 |
< |
\includegraphics[width = \linewidth]{./slice.pdf} |
| 690 |
< |
\caption{A slice from the center of a water box used in a charge void |
| 691 |
< |
simulation. The darkened region represents the boundary sphere within |
| 692 |
< |
which the water molecules were converted to argon atoms.} |
| 693 |
< |
\label{fig:argonSlice} |
| 694 |
< |
\end{figure} |
| 687 |
> |
%\begin{figure} |
| 688 |
> |
%\centering |
| 689 |
> |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 690 |
> |
%\caption{A slice from the center of a water box used in a charge void |
| 691 |
> |
%simulation. The darkened region represents the boundary sphere within |
| 692 |
> |
%which the water molecules were converted to argon atoms.} |
| 693 |
> |
%\label{fig:argonSlice} |
| 694 |
> |
%\end{figure} |
| 695 |
|
|
| 696 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 697 |
|
We compared the following alternative summation methods with results |
| 708 |
|
were utilized for the reaction field simulations. Additionally, we |
| 709 |
|
investigated the use of these cutoffs with the SP, SF, and pure |
| 710 |
|
cutoff. The SPME electrostatics were performed using the TINKER |
| 711 |
< |
implementation of SPME,\cite{Ponder87} while all other method |
| 712 |
< |
calculations were performed using the OOPSE molecular mechanics |
| 711 |
> |
implementation of SPME,\cite{Ponder87} while all other calculations |
| 712 |
> |
were performed using the {\sc oopse} molecular mechanics |
| 713 |
|
package.\cite{Meineke05} All other portions of the energy calculation |
| 714 |
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 715 |
|
manner across all systems and configurations. |
| 716 |
|
|
| 717 |
|
The althernative methods were also evaluated with three different |
| 718 |
< |
cutoff radii (9, 12, and 15 \AA). It should be noted that the damping |
| 719 |
< |
parameter chosen in SPME, or so called ``Ewald Coefficient'', has a |
| 720 |
< |
significant effect on the energies and forces calculated. Typical |
| 721 |
< |
molecular mechanics packages set this to a value dependent on the |
| 722 |
< |
cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$ |
| 723 |
< |
kcal/mol). Smaller tolerances are typically associated with increased |
| 724 |
< |
accuracy at the expense of increased time spent calculating the |
| 718 |
> |
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 719 |
> |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 720 |
> |
real-space and reciprocal-space portions of the Ewald calculation. |
| 721 |
> |
Typical molecular mechanics packages set this to a value dependent on |
| 722 |
> |
the cutoff radius and a tolerance (typically less than $1 \times |
| 723 |
> |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 724 |
> |
increasing accuracy at the expense of computational time spent on the |
| 725 |
|
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 726 |
|
The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 727 |
< |
in all SPME calculations, resulting in Ewald Coefficients of 0.4200, |
| 727 |
> |
in all SPME calculations, resulting in Ewald coefficients of 0.4200, |
| 728 |
|
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 729 |
|
respectively. |
| 730 |
|
|
| 768 |
|
readers can consult the accompanying supporting information for a |
| 769 |
|
comparison where all groups are neutral. |
| 770 |
|
|
| 771 |
< |
For the {\sc sp} method, inclusion of potential damping improves the |
| 772 |
< |
agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 773 |
< |
an excellent correlation and quality of fit with the SPME results, |
| 774 |
< |
particularly with a cutoff radius greater than 12 |
| 771 |
> |
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 772 |
> |
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 773 |
> |
shows an excellent correlation and quality of fit with the SPME |
| 774 |
> |
results, particularly with a cutoff radius greater than 12 |
| 775 |
|
\AA . Use of a larger damping parameter is more helpful for the |
| 776 |
|
shortest cutoff shown, but it has a detrimental effect on simulations |
| 777 |
|
with larger cutoffs. |
| 778 |
|
|
| 779 |
< |
In the {\sc sf} sets, increasing damping results in progressively |
| 780 |
< |
worse correlation with Ewald. Overall, the undamped case is the best |
| 779 |
> |
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 780 |
> |
worse} correlation with Ewald. Overall, the undamped case is the best |
| 781 |
|
performing set, as the correlation and quality of fits are |
| 782 |
|
consistently superior regardless of the cutoff distance. The undamped |
| 783 |
|
case is also less computationally demanding (because no evaluation of |
| 792 |
|
|
| 793 |
|
Evaluation of pairwise methods for use in Molecular Dynamics |
| 794 |
|
simulations requires consideration of effects on the forces and |
| 795 |
< |
torques. Investigation of the force and torque vector magnitudes |
| 796 |
< |
provides a measure of the strength of these values relative to SPME. |
| 797 |
< |
Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the |
| 798 |
< |
force and torque vector magnitude regression results for the |
| 798 |
< |
accumulated analysis over all the system types. |
| 795 |
> |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
| 796 |
> |
regression results for the force and torque vector magnitudes, |
| 797 |
> |
respectively. The data in these figures was generated from an |
| 798 |
> |
accumulation of the statistics from all of the system types. |
| 799 |
|
|
| 800 |
|
\begin{figure} |
| 801 |
|
\centering |
| 802 |
|
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 803 |
< |
\caption{Statistical analysis of the quality of the force vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate force magnitude 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).} |
| 803 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 804 |
> |
magnitudes for a given electrostatic method compared with the |
| 805 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 806 |
> |
indicate force magnitude values indistinguishable from those obtained |
| 807 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 808 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 809 |
> |
inverted triangles).} |
| 810 |
|
\label{fig:frcMag} |
| 811 |
|
\end{figure} |
| 812 |
|
|
| 813 |
+ |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 814 |
+ |
methods are doing at reproducing the SPME forces. The undamped and |
| 815 |
+ |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 816 |
+ |
This is perhaps expected because this method explicitly incorporates a |
| 817 |
+ |
smooth transition in the forces at the cutoff radius as well as the |
| 818 |
+ |
neutralizing image charges. |
| 819 |
+ |
|
| 820 |
|
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 821 |
|
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 822 |
|
poor, but using group based cutoffs and a switching function provides |
| 823 |
< |
a improvement much more significant than what was seen with $\Delta |
| 824 |
< |
E$. Looking at the {\sc sp} sets, the slope and $R^2$ |
| 825 |
< |
improve with the use of damping to an optimal result of 0.2 \AA |
| 826 |
< |
$^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, |
| 823 |
> |
an improvement much more significant than what was seen with $\Delta |
| 824 |
> |
E$. |
| 825 |
> |
|
| 826 |
> |
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 827 |
> |
method is generating usable forces. Further increases in damping, |
| 828 |
|
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 829 |
< |
detrimental to simulations with larger cutoff radii. The undamped |
| 830 |
< |
{\sc sf} method gives forces in line with those obtained using |
| 831 |
< |
SPME, and use of a damping function results in minor improvement. The |
| 818 |
< |
reaction field results are surprisingly good, considering the poor |
| 829 |
> |
detrimental to simulations with larger cutoff radii. |
| 830 |
> |
|
| 831 |
> |
The reaction field results are surprisingly good, considering the poor |
| 832 |
|
quality of the fits for the $\Delta E$ results. There is still a |
| 833 |
< |
considerable degree of scatter in the data, but it correlates well in |
| 834 |
< |
general. To be fair, we again note that the reaction field |
| 835 |
< |
calculations do not encompass NaCl crystal and melt systems, so these |
| 833 |
> |
considerable degree of scatter in the data, but the forces correlate |
| 834 |
> |
well with the Ewald forces in general. We note that the reaction |
| 835 |
> |
field calculations do not include the pure NaCl systems, so these |
| 836 |
|
results are partly biased towards conditions in which the method |
| 837 |
|
performs more favorably. |
| 838 |
|
|
| 839 |
|
\begin{figure} |
| 840 |
|
\centering |
| 841 |
|
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 842 |
< |
\caption{Statistical analysis of the quality of the torque vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate torque magnitude 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).} |
| 842 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 843 |
> |
magnitudes for a given electrostatic method compared with the |
| 844 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 845 |
> |
indicate torque magnitude values indistinguishable from those obtained |
| 846 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 847 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 848 |
> |
inverted triangles).} |
| 849 |
|
\label{fig:trqMag} |
| 850 |
|
\end{figure} |
| 851 |
|
|
| 852 |
< |
To evaluate the torque vector magnitudes, the data set from which |
| 853 |
< |
values are drawn is limited to rigid molecules in the systems |
| 854 |
< |
(i.e. water molecules). In spite of this smaller sampling pool, the |
| 836 |
< |
torque vector magnitude results in figure \ref{fig:trqMag} are still |
| 837 |
< |
similar to those seen for the forces; however, they more clearly show |
| 838 |
< |
the improved behavior that comes with increasing the cutoff radius. |
| 839 |
< |
Moderate damping is beneficial to the {\sc sp} and helpful |
| 840 |
< |
yet possibly unnecessary with the {\sc sf} method, and they also |
| 841 |
< |
show that over-damping adversely effects all cutoff radii rather than |
| 842 |
< |
showing an improvement for systems with short cutoffs. The reaction |
| 843 |
< |
field method performs well when calculating the torques, better than |
| 844 |
< |
the Shifted Force method over this limited data set. |
| 852 |
> |
Molecular torques were only available from the systems which contained |
| 853 |
> |
rigid molecules (i.e. the systems containing water). The data in |
| 854 |
> |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 855 |
|
|
| 856 |
+ |
Torques appear to be much more sensitive to charges at a longer |
| 857 |
+ |
distance. The striking feature in comparing the new electrostatic |
| 858 |
+ |
methods with SPME is how much the agreement improves with increasing |
| 859 |
+ |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 860 |
+ |
appears to be reproducing the SPME torques most accurately. |
| 861 |
+ |
|
| 862 |
+ |
Water molecules are dipolar, and the reaction field method reproduces |
| 863 |
+ |
the effect of the surrounding polarized medium on each of the |
| 864 |
+ |
molecular bodies. Therefore it is not surprising that reaction field |
| 865 |
+ |
performs best of all of the methods on molecular torques. |
| 866 |
+ |
|
| 867 |
|
\subsection{Directionality of the Force and Torque Vectors} |
| 868 |
|
|
| 869 |
< |
Having force and torque vectors with magnitudes that are well |
| 870 |
< |
correlated to SPME is good, but if they are not pointing in the proper |
| 871 |
< |
direction the results will be incorrect. These vector directions were |
| 872 |
< |
investigated through measurement of the angle formed between them and |
| 873 |
< |
those from SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared |
| 874 |
< |
through the variance ($\sigma^2$) of the Gaussian fits of the angle |
| 875 |
< |
error distributions of the combined set over all system types. |
| 869 |
> |
It is clearly important that a new electrostatic method can reproduce |
| 870 |
> |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 871 |
> |
sum. However, the {\it directionality} of these vectors will also be |
| 872 |
> |
vital in calculating dynamical quantities accurately. Force and |
| 873 |
> |
torque directionalities were investigated by measuring the angles |
| 874 |
> |
formed between these vectors and the same vectors calculated using |
| 875 |
> |
SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 876 |
> |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 877 |
> |
distributions of the combined set over all system types. |
| 878 |
|
|
| 879 |
|
\begin{figure} |
| 880 |
|
\centering |
| 881 |
|
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 882 |
< |
\caption{Statistical analysis of the quality of the Gaussian fit of the force and torque vector angular distributions for a given electrostatic method compared with the reference Ewald sum. Results with a variance ($\sigma^2$) equal to zero (dashed line) indicate force and torque directions 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).} |
| 882 |
> |
\caption{Statistical analysis of the width of the angular distribution |
| 883 |
> |
that the force and torque vectors from a given electrostatic method |
| 884 |
> |
make with their counterparts obtained using the reference Ewald sum. |
| 885 |
> |
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 886 |
> |
indicate force and torque directions indistinguishable from those |
| 887 |
> |
obtained using SPME. Different values of the cutoff radius are |
| 888 |
> |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 889 |
> |
and 15\AA\ = inverted triangles).} |
| 890 |
|
\label{fig:frcTrqAng} |
| 891 |
|
\end{figure} |
| 892 |
|
|
| 893 |
|
Both the force and torque $\sigma^2$ results from the analysis of the |
| 894 |
|
total accumulated system data are tabulated in figure |
| 895 |
< |
\ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case |
| 896 |
< |
show the improvement afforded by choosing a longer simulation cutoff. |
| 897 |
< |
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving |
| 898 |
< |
of the distribution widths, with a similar improvement going from 12 |
| 899 |
< |
to 15 \AA . The undamped {\sc sf}, Group Based Cutoff, and |
| 870 |
< |
Reaction Field methods all do equivalently well at capturing the |
| 871 |
< |
direction of both the force and torque vectors. Using damping |
| 872 |
< |
improves the angular behavior significantly for the {\sc sp} |
| 873 |
< |
and moderately for the {\sc sf} methods. Increasing the damping |
| 874 |
< |
too far is destructive for both methods, particularly to the torque |
| 875 |
< |
vectors. Again it is important to recognize that the force vectors |
| 876 |
< |
cover all particles in the systems, while torque vectors are only |
| 877 |
< |
available for neutral molecular groups. Damping appears to have a |
| 878 |
< |
more beneficial effect on non-neutral bodies, and this observation is |
| 879 |
< |
investigated further in the accompanying supporting information. |
| 895 |
> |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 896 |
> |
sp}) method would be essentially unusable for molecular dynamics until |
| 897 |
> |
the damping function is added. The Shifted Force ({\sc sf}) method, |
| 898 |
> |
however, is generating force and torque vectors which are within a few |
| 899 |
> |
degrees of the Ewald results even with weak (or no) damping. |
| 900 |
|
|
| 901 |
+ |
All of the sets (aside from the over-damped case) show the improvement |
| 902 |
+ |
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 903 |
+ |
from 9 to 12 \AA\ typically results in a halving of the width of the |
| 904 |
+ |
distribution, with a similar improvement going from 12 to 15 |
| 905 |
+ |
\AA . |
| 906 |
+ |
|
| 907 |
+ |
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 908 |
+ |
all do equivalently well at capturing the direction of both the force |
| 909 |
+ |
and torque vectors. Using damping improves the angular behavior |
| 910 |
+ |
significantly for the {\sc sp} and moderately for the {\sc sf} |
| 911 |
+ |
methods. Overdamping is detrimental to both methods. Again it is |
| 912 |
+ |
important to recognize that the force vectors cover all particles in |
| 913 |
+ |
the systems, while torque vectors are only available for neutral |
| 914 |
+ |
molecular groups. Damping appears to have a more beneficial effect on |
| 915 |
+ |
charged bodies, and this observation is investigated further in the |
| 916 |
+ |
accompanying supporting information. |
| 917 |
+ |
|
| 918 |
+ |
Although not discussed previously, group based cutoffs can be applied |
| 919 |
+ |
to both the {\sc sp} and {\sc sf} methods. Use of a switching |
| 920 |
+ |
function corrects for the discontinuities that arise when atoms of the |
| 921 |
+ |
two groups exit the cutoff radius before the group centers leave each |
| 922 |
+ |
other's cutoff. Though there are no significant benefits or drawbacks |
| 923 |
+ |
observed in $\Delta E$ and vector magnitude results when doing this, |
| 924 |
+ |
there is a measurable improvement in the vector angle results. Table |
| 925 |
+ |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 926 |
+ |
group based cutoffs and a switching function alongside the results |
| 927 |
+ |
seen in figure \ref{fig:frcTrqAng}. The {\sc sp} shows much narrower |
| 928 |
+ |
angular distributions for both the force and torque vectors when using |
| 929 |
+ |
an $\alpha$ of 0.2 \AA$^{-1}$ or less, while {\sc sf} shows |
| 930 |
+ |
improvements in the undamped and lightly damped cases. Thus, by |
| 931 |
+ |
calculating the electrostatic interactions in terms of molecular pairs |
| 932 |
+ |
rather than atomic pairs, the direction of the force and torque |
| 933 |
+ |
vectors can be determined more accurately. |
| 934 |
+ |
|
| 935 |
|
\begin{table}[htbp] |
| 936 |
|
\centering |
| 937 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 937 |
> |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 938 |
> |
(bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 939 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 940 |
|
\\ |
| 941 |
|
\toprule |
| 965 |
|
\end{tabular} |
| 966 |
|
\label{tab:groupAngle} |
| 967 |
|
\end{table} |
| 913 |
– |
|
| 914 |
– |
Although not discussed previously, group based cutoffs can be applied |
| 915 |
– |
to both the {\sc sp} and {\sc sf} methods. Use off a |
| 916 |
– |
switching function corrects for the discontinuities that arise when |
| 917 |
– |
atoms of a group exit the cutoff before the group's center of mass. |
| 918 |
– |
Though there are no significant benefit or drawbacks observed in |
| 919 |
– |
$\Delta E$ and vector magnitude results when doing this, there is a |
| 920 |
– |
measurable improvement in the vector angle results. Table |
| 921 |
– |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 922 |
– |
group based cutoffs and a switching function alongside the standard |
| 923 |
– |
results seen in figure \ref{fig:frcTrqAng} for comparison purposes. |
| 924 |
– |
The {\sc sp} shows much narrower angular distributions for |
| 925 |
– |
both the force and torque vectors when using an $\alpha$ of 0.2 |
| 926 |
– |
\AA$^{-1}$ or less, while {\sc sf} shows improvements in the |
| 927 |
– |
undamped and lightly damped cases. Thus, by calculating the |
| 928 |
– |
electrostatic interactions in terms of molecular pairs rather than |
| 929 |
– |
atomic pairs, the direction of the force and torque vectors are |
| 930 |
– |
determined more accurately. |
| 968 |
|
|
| 969 |
|
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 970 |
< |
that the $\sigma^2$ values for both {\sc sp} and |
| 971 |
< |
{\sc sf} converge as $\alpha$ increases, something that is easier |
| 972 |
< |
to see when using group based cutoffs. Looking back on figures |
| 973 |
< |
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this |
| 974 |
< |
behavior clearly at large $\alpha$ and cutoff values. The reason for |
| 975 |
< |
this is that the complimentary error function inserted into the |
| 976 |
< |
potential weakens the electrostatic interaction as $\alpha$ increases. |
| 977 |
< |
Thus, at larger values of $\alpha$, both the summation method types |
| 978 |
< |
progress toward non-interacting functions, so care is required in |
| 979 |
< |
choosing large damping functions lest one generate an undesirable loss |
| 980 |
< |
in the pair interaction. Kast \textit{et al.} developed a method for |
| 981 |
< |
choosing appropriate $\alpha$ values for these types of electrostatic |
| 982 |
< |
summation methods by fitting to $g(r)$ data, and their methods |
| 983 |
< |
indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff |
| 984 |
< |
values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear |
| 985 |
< |
to be reasonable choices to obtain proper MC behavior |
| 986 |
< |
(Fig. \ref{fig:delE}); however, based on these findings, choices this |
| 987 |
< |
high would introduce error in the molecular torques, particularly for |
| 988 |
< |
the shorter cutoffs. Based on the above findings, empirical damping |
| 952 |
< |
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 953 |
< |
unnecessary when using the {\sc sf} method. |
| 970 |
> |
that the $\sigma^2$ values for both {\sc sp} and {\sc sf} converge as |
| 971 |
> |
$\alpha$ increases, something that is easier to see when using group |
| 972 |
> |
based cutoffs. The reason for this is that the complimentary error |
| 973 |
> |
function inserted into the potential weakens the electrostatic |
| 974 |
> |
interaction as $\alpha$ increases. Thus, at larger values of |
| 975 |
> |
$\alpha$, both summation methods progress toward non-interacting |
| 976 |
> |
functions, so care is required in choosing large damping functions |
| 977 |
> |
lest one generate an undesirable loss in the pair interaction. Kast |
| 978 |
> |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 979 |
> |
values for these types of electrostatic summation methods by fitting |
| 980 |
> |
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 981 |
> |
0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ |
| 982 |
> |
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 983 |
> |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 984 |
> |
these findings, choices this high would introduce error in the |
| 985 |
> |
molecular torques, particularly for the shorter cutoffs. Based on the |
| 986 |
> |
above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be |
| 987 |
> |
beneficial, but damping may be unnecessary when using the {\sc sf} |
| 988 |
> |
method. |
| 989 |
|
|
| 990 |
|
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 991 |
|
|
| 992 |
< |
In the previous studies using a {\sc sf} variant of the damped |
| 993 |
< |
Wolf coulomb potential, the structure and dynamics of water were |
| 994 |
< |
investigated rather extensively.\cite{Zahn02,Kast03} Their results |
| 995 |
< |
indicated that the damped {\sc sf} method results in properties |
| 996 |
< |
very similar to those obtained when using the Ewald summation. |
| 997 |
< |
Considering the statistical results shown above, the good performance |
| 998 |
< |
of this method is not that surprising. Rather than consider the same |
| 999 |
< |
systems and simply recapitulate their results, we decided to look at |
| 1000 |
< |
the solid state dynamical behavior obtained using the best performing |
| 1001 |
< |
summation methods from the above results. |
| 992 |
> |
Zahn {\it et al.} investigated the structure and dynamics of water |
| 993 |
> |
using eqs. (\ref{eq:ZahnPot}) and |
| 994 |
> |
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 995 |
> |
that a method similar (but not identical with) the damped {\sc sf} |
| 996 |
> |
method resulted in properties very similar to those obtained when |
| 997 |
> |
using the Ewald summation. The properties they studied (pair |
| 998 |
> |
distribution functions, diffusion constants, and velocity and |
| 999 |
> |
orientational correlation functions) may not be particularly sensitive |
| 1000 |
> |
to the long-range and collective behavior that governs the |
| 1001 |
> |
low-frequency behavior in crystalline systems. |
| 1002 |
|
|
| 1003 |
+ |
We are using two separate measures to probe the effects of these |
| 1004 |
+ |
alternative electrostatic methods on the dynamics in crystalline |
| 1005 |
+ |
materials. For short- and intermediate-time dynamics, we are |
| 1006 |
+ |
computing the velocity autocorrelation function, and for long-time |
| 1007 |
+ |
and large length-scale collective motions, we are looking at the |
| 1008 |
+ |
low-frequency portion of the power spectrum. |
| 1009 |
+ |
|
| 1010 |
|
\begin{figure} |
| 1011 |
|
\centering |
| 1012 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 1013 |
< |
\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.} |
| 1013 |
> |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 1014 |
> |
1000 K using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1015 |
> |
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 1016 |
> |
the first minimum. The times to first collision are nearly identical, |
| 1017 |
> |
but differences can be seen in the peaks and troughs, where the |
| 1018 |
> |
undamped and weakly damped methods are stiffer than the moderately |
| 1019 |
> |
damped and SPME methods.} |
| 1020 |
|
\label{fig:vCorrPlot} |
| 1021 |
|
\end{figure} |
| 1022 |
|
|
| 1051 |
|
\begin{figure} |
| 1052 |
|
\centering |
| 1053 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1054 |
< |
\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.} |
| 1054 |
> |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1055 |
> |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1056 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1057 |
> |
Apodization of the correlation functions via a cubic switching |
| 1058 |
> |
function between 40 and 50 ps was used to clear up the spectral noise |
| 1059 |
> |
resulting from data truncation, and had no noticeable effect on peak |
| 1060 |
> |
location or magnitude. The inset shows the frequency region below 100 |
| 1061 |
> |
cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1062 |
|
\label{fig:methodPS} |
| 1063 |
|
\end{figure} |
| 1064 |
|
|
| 1100 |
|
\begin{figure} |
| 1101 |
|
\centering |
| 1102 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1103 |
< |
\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods. The upper plot is a zoomed inset from figure \ref{fig:methodPS}. As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift. The lower plot is of spectra when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 1103 |
> |
\caption{Regions of spectra showing the low-frequency correlated |
| 1104 |
> |
motions for NaCl crystals at 1000 K using various electrostatic |
| 1105 |
> |
summation methods. The upper plot is a zoomed inset from figure |
| 1106 |
> |
\ref{fig:methodPS}. As the damping value for the {\sc sf} potential |
| 1107 |
> |
increases, the low-frequency peaks red-shift. The lower plot is of |
| 1108 |
> |
spectra when using SPME and a simple damped Coulombic sum with damping |
| 1109 |
> |
coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As |
| 1110 |
> |
$\alpha$ increases, the peaks are red-shifted toward and eventually |
| 1111 |
> |
beyond the values given by SPME. The larger $\alpha$ values weaken |
| 1112 |
> |
the real-space electrostatics, explaining this shift towards less |
| 1113 |
> |
strongly correlated motions in the crystal.} |
| 1114 |
|
\label{fig:dampInc} |
| 1115 |
|
\end{figure} |
| 1116 |
|
|
