941 |
|
n_C = 4\pi\rho \int_{0}^{a}r^2\text{g}(r)dr, |
942 |
|
\end{equation} |
943 |
|
where $\rho$ is the number density of the site-site pair interactions, |
944 |
< |
$a$ and is the radial location of the minima following the first peak |
945 |
< |
in $g(r)$ ($a = 4.2$ \AA for the SSDQ water/ion system). The |
944 |
> |
and $a$ is the radial location of the minima following the first peak |
945 |
> |
in $g(r)$ ($a = 4.2$ \AA\ for the SSDQ water/ion system). The |
946 |
|
coordination number is shown as a function of the damping coefficient |
947 |
< |
for all of the real space methods in Fig. \ref{fig:Props}. |
947 |
> |
for all of the real space methods in Fig. \ref{fig:Props}. |
948 |
|
|
949 |
|
A more demanding test of modified electrostatics is the average value |
950 |
|
of the electrostatic energy $\langle U_\mathrm{elect} \rangle / N$ |
952 |
|
experienced by a liquid evolving entirely under the influence of each |
953 |
|
of the methods. In fig \ref{fig:Props} we demonstrate how $\langle |
954 |
|
U_\mathrm{elect} \rangle / N$ varies with the damping parameter, |
955 |
< |
$\alpha$, for each of the methods. |
955 |
> |
$\alpha$, for each of the methods. |
956 |
|
|
957 |
|
As in the crystals studied in the first paper, damping is important |
958 |
|
for converging the mean electrostatic energy values, particularly for |
964 |
|
0.3$, causing the estimate of the electrostatic energy to drop below |
965 |
|
the Ewald results. |
966 |
|
|
967 |
< |
These ``optimal'' values of the damping coefficient are slightly |
968 |
< |
larger than what were observed for DSF electrostatics for purely |
969 |
< |
point-charge systems, although the range $\alpha= 0.175 \rightarrow |
967 |
> |
These ``optimal'' values of the damping coefficient for structural |
968 |
> |
features are similar to those observed for DSF electrostatics for |
969 |
> |
purely point-charge systems, and the range $\alpha= 0.175 \rightarrow |
970 |
|
0.225$ \AA$^{-1}$ for $r_c = 12$\AA\ appears to be an excellent |
971 |
|
compromise for mixed charge/multipolar systems. |
972 |
|
|
979 |
|
\label{eq:diff} |
980 |
|
\end{equation} |
981 |
|
which measures long-time behavior and is sensitive to the forces on |
982 |
< |
the multipoles. The self-diffusion constants (D) were calculated from |
982 |
> |
the multipoles. The self-diffusion constants (D) were calculated from |
983 |
|
linear fits to the long-time portion of the mean square displacement, |
984 |
|
$\langle r^{2}(t) \rangle$.\cite{Allen87} In fig. \ref{fig:Props} we |
985 |
|
demonstrate how the diffusion constant depends on the choice of |
989 |
|
|
990 |
|
In addition to translational diffusion, orientational relaxation times |
991 |
|
were calculated for comparisons with the Ewald simulations and with |
992 |
< |
experiments. These values were determined from the same 1~ns |
993 |
< |
microcanonical trajectories used for translational diffusion by |
994 |
< |
calculating the orientational time correlation function, |
992 |
> |
experiments. These values were determined by calculating the |
993 |
> |
orientational time correlation function, |
994 |
|
\begin{equation} |
995 |
|
C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t) |
996 |
|
\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle, |
997 |
|
\label{eq:OrientCorr} |
998 |
|
\end{equation} |
999 |
< |
where $P_l$ is the Legendre polynomial of order $l$ and |
1000 |
< |
$\hat{\mathbf{A}}_\gamma$ is the unit vector for body axis $\gamma$. |
1001 |
< |
The reference frame used for our sample dipolar systems has the |
1002 |
< |
$z$-axis running along the dipoles, and for the SSDQ water model, the |
1003 |
< |
$y$-axis connects the two implied hydrogen atom positions. From the |
1004 |
< |
orientation autocorrelation functions, we can obtain time constants |
1005 |
< |
for rotational relaxation either by fitting an exponential function or |
1006 |
< |
by integrating the entire correlation function. In a good water |
1007 |
< |
model, these decay times would be comparable to water orientational |
1008 |
< |
relaxation times from nuclear magnetic resonance (NMR). The relaxation |
1010 |
< |
constant obtained from $C_2^y(t)$ is normally of experimental interest |
1011 |
< |
because it describes the relaxation of the principle axis connecting |
1012 |
< |
the hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the |
1013 |
< |
intermolecular portion of the dipole-dipole relaxation from a proton |
1014 |
< |
NMR signal and should provide an estimate of the NMR relaxation time |
1015 |
< |
constant.\cite{Impey82} |
999 |
> |
from the same 350 ps microcanonical trajectories that were used for |
1000 |
> |
translational diffusion. Here, $P_l$ is the Legendre polynomial of |
1001 |
> |
order $l$ and $\hat{\mathbf{A}}_\gamma$ is the unit vector for body |
1002 |
> |
axis $\gamma$. The reference frame used for our sample dipolar |
1003 |
> |
systems has the $z$-axis running along the dipoles, and for the SSDQ |
1004 |
> |
water model, the $y$-axis connects the two implied hydrogen atom |
1005 |
> |
positions. From the orientation autocorrelation functions, we can |
1006 |
> |
obtain time constants for rotational relaxation either by fitting to a |
1007 |
> |
multi-exponential model for the orientational relaxation, or by |
1008 |
> |
integrating the correlation functions. |
1009 |
|
|
1010 |
< |
Results for the diffusion constants and orientational relaxation times |
1011 |
< |
are shown in figure \ref{fig:Props}. From this data, it is apparent |
1012 |
< |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
1013 |
< |
reproduced with reasonable fidelity by the GSF method. |
1010 |
> |
In a good water model, the orientational decay times would be |
1011 |
> |
comparable to water orientational relaxation times from nuclear |
1012 |
> |
magnetic resonance (NMR). The relaxation constant obtained from |
1013 |
> |
$C_2^y(t)$ is normally of experimental interest because it describes |
1014 |
> |
the relaxation of the principle axis connecting the hydrogen |
1015 |
> |
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
1016 |
> |
of the dipole-dipole relaxation from a proton NMR signal and can |
1017 |
> |
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
1018 |
> |
In fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$ |
1019 |
> |
values for the various real-space methods over a range of different |
1020 |
> |
damping coefficients. The rotational relaxation for the $z$ axis |
1021 |
> |
primarily probes the torques on the dipoles, while the relaxation for |
1022 |
> |
the $y$ axis is sensitive primarily to the quadrupolar torques. |
1023 |
|
|
1024 |
|
\begin{figure} |
1025 |
|
\caption{Comparison of the structural and dynamic properties for the |
1026 |
|
combined multipolar liquid (SSDQ water + ions) for all of the |
1027 |
|
real-space methods with $r_c = 12$\AA. Electrostatic energies, |
1028 |
|
$\langle U_\mathrm{elect} \rangle / N$ (in kcal / mol), |
1029 |
< |
coordination numbers, $n_C$, diffusion constants (in cm$^2$ |
1030 |
< |
s$^{-1}$), and rotational correlation times (in fs) all show |
1031 |
< |
excellent agreement with Ewald results for damping coefficients in |
1032 |
< |
the range $\alpha= 0.175 \rightarrow 0.225$ |
1033 |
< |
\AA$^{-1}$. \label{fig:Props}} |
1029 |
> |
coordination numbers, $n_C$, diffusion constants (in $10^{-5} |
1030 |
> |
\mathrm{cm}^2\mathrm{s}^{-1}$), and rotational correlation times |
1031 |
> |
(in ps) all show excellent agreement with Ewald results for |
1032 |
> |
damping coefficients in the range $\alpha= 0.175 \rightarrow |
1033 |
> |
0.225$ \AA$^{-1}$. \label{fig:Props}} |
1034 |
|
\includegraphics[width=\textwidth]{properties.eps} |
1035 |
|
\end{figure} |
1036 |
|
|
1037 |
+ |
In Fig. \ref{fig:Props} it appears that values for $D$, $\tau_2^y$, |
1038 |
+ |
and $\tau_2^z$ using the Ewald sum are reproduced with excellent |
1039 |
+ |
fidelity by the GSF and SP methods. All of the real space methods can |
1040 |
+ |
be \textit{overdamped}, which reduces the effective range of multipole |
1041 |
+ |
interactions, causing structural and dynamical changes from the |
1042 |
+ |
correct behavior. Because overdamping weakens orientational |
1043 |
+ |
preferences between adjacent molecules, it manifests as too-rapid |
1044 |
+ |
orientational decay coupled with faster diffusion and |
1045 |
+ |
over-coordination of the liquid. Underdamping is less problematic for |
1046 |
+ |
the SP and GSF methods, as their structural and dynamical properties |
1047 |
+ |
still reproduce the Ewald results even in the completely undamped |
1048 |
+ |
($\alpha = 0$) case. An optimal range for the electrostatic damping |
1049 |
+ |
parameter appears to be $\alpha= 0.175 \rightarrow 0.225$ \AA$^{-1}$ |
1050 |
+ |
for $r_c = 12$\AA, which similar to the optimal range found for the |
1051 |
+ |
damped shifted force potential for point charges.\cite{Fennell:2006lq} |
1052 |
|
|
1053 |
|
\section{CONCLUSION} |
1054 |
|
In the first paper in this series, we generalized the |
1073 |
|
sphere to derive shifted force and torque expressions, and is a |
1074 |
|
significantly more gentle approach. |
1075 |
|
|
1076 |
< |
The GSF method produced quantitative agreement with Ewald energy, |
1077 |
< |
force, and torques. It also performs well in conserving energy in MD |
1076 |
> |
The GSF method produces quantitative agreement with Ewald energies, |
1077 |
> |
forcse, and torques. It also performs well in conserving energy in MD |
1078 |
|
simulations. The Taylor-shifted (TSF) model provides smooth dynamics, |
1079 |
|
but these take place on a potential energy surface that is |
1080 |
|
significantly perturbed from Ewald-based electrostatics. Because it |
1086 |
|
method also has the unique property that a large number of derivatives |
1087 |
|
can be made to vanish at the cutoff radius. This property has proven |
1088 |
|
useful in past treatments of the corrections to the fluctuation |
1089 |
< |
formula for dielectric constants.\cite{Izvekov:2008wo} |
1089 |
> |
formula for dielectric constants.\cite{Izvekov:2008wo} |
1090 |
|
|
1091 |
|
Reproduction of both structural and dynamical features in the liquid |
1092 |
|
systems is remarkably good for both the SP and GSF models. Pair |
1117 |
|
important. Both the SP and GSF methods retain excellent fidelity to |
1118 |
|
the Ewald energies, forces and torques. Additionally, the energy |
1119 |
|
drift and fluctuations from the GSF electrostatics are significantly |
1120 |
< |
better than a multipolar Ewald sum for finite-sized reciprocal spaces. |
1120 |
> |
better than a multipolar Ewald sum for finite-sized reciprocal spaces, |
1121 |
> |
and physical properties are reproduced accurately. |
1122 |
|
|
1123 |
|
As in all purely pairwise cutoff methods, the SP, GSF and TSF methods |
1124 |
|
are expected to scale approximately {\it linearly} with system size, |