| 389 |
|
approach than they are in the Taylor-shifted method. |
| 390 |
|
|
| 391 |
|
For the gradient shifted (GSF) method with the undamped kernel, |
| 392 |
< |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
| 392 |
> |
$v_{21}(r) = -\frac{1}{r^3} - \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
| 393 |
|
$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. |
| 394 |
|
Again, these functions go smoothly to zero as $r\rightarrow r_c$, and |
| 395 |
|
because the Taylor expansion retains only one term, they are |
| 896 |
|
\begin{figure} |
| 897 |
|
\centering |
| 898 |
|
\includegraphics[width=\textwidth]{newDrift_12.eps} |
| 899 |
< |
\caption{Analysis of the energy conservation of the real-space methods |
| 900 |
< |
for the SSDQ water/ion system. $\delta \mathrm{E}_1$ is the linear |
| 901 |
< |
drift in energy over time (in kcal/mol/particle/ns) and $\delta |
| 902 |
< |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 903 |
< |
around this drift (in kcal/mol/particle). Points that appear in the |
| 904 |
< |
green region at the bottom exhibit better energy conservation than |
| 905 |
< |
would be obtained using common parameters for Ewald-based |
| 906 |
< |
electrostatics.\label{fig:energyDrift}} |
| 899 |
> |
\caption{Energy conservation of the real-space methods for the SSDQ |
| 900 |
> |
water/ion system. $\delta \mathrm{E}_1$ is the linear drift in |
| 901 |
> |
energy over time (in kcal/mol/particle/ns) and $\delta |
| 902 |
> |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 903 |
> |
around this drift (in kcal/mol/particle). Points that appear in |
| 904 |
> |
the green region at the bottom exhibit better energy conservation |
| 905 |
> |
than would be obtained using common parameters for Ewald-based |
| 906 |
> |
electrostatics.\label{fig:energyDrift}} |
| 907 |
|
\end{figure} |
| 908 |
|
|
| 909 |
|
\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}} |
| 931 |
|
treatment.\label{fig:gofr}} |
| 932 |
|
\end{figure} |
| 933 |
|
|
| 934 |
< |
There is a very slight overstructuring of the first solvation shell |
| 935 |
< |
when using when using TSF at lower values of the damping coefficient |
| 936 |
< |
($\alpha \le 0.1$) or when using undamped GSF. With moderate damping, |
| 937 |
< |
GSF and SP produce pair distributions that are identical (within |
| 938 |
< |
numerical noise) to their Ewald counterparts. |
| 934 |
> |
There is a minor overstructuring of the first solvation shell when |
| 935 |
> |
using TSF or when overdamping with any of the real-space methods. |
| 936 |
> |
With moderate damping, GSF and SP produce pair distributions that are |
| 937 |
> |
identical (within numerical noise) to their Ewald counterparts. The |
| 938 |
> |
degree of overstructuring can be measured most easily using the |
| 939 |
> |
coordination number, |
| 940 |
> |
\begin{equation} |
| 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 |
> |
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}. |
| 948 |
|
|
| 949 |
< |
A structural property that is a more demanding test of modified |
| 950 |
< |
electrostatics is the mean value of the electrostatic energy $\langle |
| 951 |
< |
U_\mathrm{elect} \rangle / N$ which is obtained by sampling the |
| 952 |
< |
liquid-state configurations experienced by a liquid evolving entirely |
| 953 |
< |
under the influence of each of the methods. In table \ref{tab:Props} |
| 954 |
< |
we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with |
| 955 |
< |
the damping parameter, $\alpha$, for each of the methods. |
| 949 |
> |
A more demanding test of modified electrostatics is the average value |
| 950 |
> |
of the electrostatic energy $\langle U_\mathrm{elect} \rangle / N$ |
| 951 |
> |
which is obtained by sampling the liquid-state configurations |
| 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. |
| 956 |
|
|
| 957 |
|
As in the crystals studied in the first paper, damping is important |
| 958 |
|
for converging the mean electrostatic energy values, particularly for |
| 959 |
|
the two shifted force methods (GSF and TSF). A value of $\alpha |
| 960 |
|
\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF |
| 961 |
|
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
| 962 |
< |
dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA). |
| 962 |
> |
dramatic damping ($\alpha \approx 0.28$ \AA$^{-1}$ for $r_c = 9$ \AA). |
| 963 |
|
Overdamping the real-space electrostatic methods occurs with $\alpha > |
| 964 |
< |
0.4$, causing the estimate of the energy to drop below the Ewald |
| 965 |
< |
results. |
| 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 a value of $\alpha=0.18$ \AA$^{-1}$ for |
| 970 |
< |
$r_c = 12$\AA\ appears to be an excellent compromise for mixed |
| 971 |
< |
charge/multipolar systems. |
| 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 |
|
|
| 973 |
|
To test the fidelity of the electrostatic methods at reproducing |
| 974 |
< |
dynamics in a multipolar liquid, it is also useful to look at |
| 974 |
> |
\textit{dynamics} in a multipolar liquid, it is also useful to look at |
| 975 |
|
transport properties, particularly the diffusion constant, |
| 976 |
|
\begin{equation} |
| 977 |
|
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
| 979 |
|
\label{eq:diff} |
| 980 |
|
\end{equation} |
| 981 |
|
which measures long-time behavior and is sensitive to the forces on |
| 982 |
< |
the multipoles. For the soft dipolar fluid and the SSDQ liquid |
| 983 |
< |
systems, the self-diffusion constants (D) were calculated from linear |
| 984 |
< |
fits to the long-time portion of the mean square displacement, |
| 985 |
< |
$\langle r^{2}(t) \rangle$.\cite{Allen87} |
| 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 |
| 986 |
> |
real-space methods and the damping coefficient. Both the SP and GSF |
| 987 |
> |
methods can obtain excellent agreement with Ewald again using moderate |
| 988 |
> |
damping. |
| 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 |
| 982 |
< |
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 |
| 998 |
< |
constant obtained from $C_2^y(t)$ is normally of experimental interest |
| 999 |
< |
because it describes the relaxation of the principle axis connecting |
| 1000 |
< |
the hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the |
| 1001 |
< |
intermolecular portion of the dipole-dipole relaxation from a proton |
| 1002 |
< |
NMR signal and should provide an estimate of the NMR relaxation time |
| 1003 |
< |
constant.\cite{Impey82} |
| 1004 |
< |
|
| 1005 |
< |
Results for the diffusion constants and orientational relaxation times |
| 1006 |
< |
are shown in figure \ref{tab:Props}. From this data, it is apparent |
| 1007 |
< |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
| 1008 |
< |
reproduced with reasonable fidelity by the GSF method. |
| 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 |
< |
The $\tau_2$ results in \ref{tab:Props} show a much greater difference |
| 1011 |
< |
between the real-space and the Ewald results. |
| 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{table} |
| 1025 |
< |
\caption{Comparison of the structural and dynamic properties for the |
| 1026 |
< |
soft dipolar liquid test for all of the real-space methods.\label{tab:Props}} |
| 1027 |
< |
\begin{tabular}{l|c|cccc|cccc|cccc} |
| 1028 |
< |
& Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\ |
| 1029 |
< |
$\alpha$ (\AA$^{-1}$) & & |
| 1030 |
< |
0.0 & 0.1 & 0.2 & 0.3 & |
| 1031 |
< |
0.0 & 0.1 & 0.2 & 0.3 & |
| 1032 |
< |
0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14} |
| 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 $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 |
< |
$\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\ |
| 1038 |
< |
D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)& |
| 1039 |
< |
470.2(6) & |
| 1040 |
< |
416.6(5) & |
| 1041 |
< |
379.6(5) & |
| 1042 |
< |
438.6(5) & |
| 1043 |
< |
476.0(6) & |
| 1044 |
< |
412.8(5) & |
| 1045 |
< |
421.1(5) & |
| 1046 |
< |
400.5(5) & |
| 1047 |
< |
437.5(6) & |
| 1048 |
< |
434.6(5) & |
| 1049 |
< |
411.4(5) & |
| 1050 |
< |
545.3(7) & |
| 1051 |
< |
459.6(6) \\ |
| 1038 |
< |
$\tau_2$ (fs) & |
| 1039 |
< |
1.136 & |
| 1040 |
< |
1.041 & |
| 1041 |
< |
1.064 & |
| 1042 |
< |
1.109 & |
| 1043 |
< |
1.211 & |
| 1044 |
< |
1.119 & |
| 1045 |
< |
1.039 & |
| 1046 |
< |
1.058 & |
| 1047 |
< |
1.21 & |
| 1048 |
< |
1.15 & |
| 1049 |
< |
1.172 & |
| 1050 |
< |
1.153 & |
| 1051 |
< |
1.125 \\ |
| 1052 |
< |
\end{tabular} |
| 1053 |
< |
\end{table} |
| 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 |
|
|
| 1055 |
– |
|
| 1053 |
|
\section{CONCLUSION} |
| 1054 |
|
In the first paper in this series, we generalized the |
| 1055 |
|
charge-neutralized electrostatic energy originally developed by Wolf |
| 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, |
