| 47 |
|
|
| 48 |
|
%\preprint{AIP/123-QED} |
| 49 |
|
|
| 50 |
< |
\title{Real space alternatives to the Ewald Sum. II. Comparison of Methods} |
| 50 |
> |
\title{Real space electrostatics for multipoles. II. Comparisons with |
| 51 |
> |
the Ewald Sum} |
| 52 |
|
|
| 53 |
|
\author{Madan Lamichhane} |
| 54 |
|
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
| 995 |
|
of a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 996 |
|
300 K starting from the same initial configuration. All runs |
| 997 |
|
utilized the same real-space cutoff, $r_c = 12$\AA.} |
| 998 |
+ |
\end{figure} |
| 999 |
+ |
|
| 1000 |
+ |
\subsection{Reproduction of Structural Features\label{sec:structure}} |
| 1001 |
+ |
One of the best tests of modified interaction potentials is the |
| 1002 |
+ |
fidelity with which they can reproduce structural features in a |
| 1003 |
+ |
liquid. One commonly-utilized measure of structural ordering is the |
| 1004 |
+ |
pair distribution function, $g(r)$, which measures local density |
| 1005 |
+ |
deviations in relation to the bulk density. In the electrostatic |
| 1006 |
+ |
approaches studied here, the short-range repulsion from the |
| 1007 |
+ |
Lennard-Jones potential is identical for the various electrostatic |
| 1008 |
+ |
methods, and since short range repulsion determines much of the local |
| 1009 |
+ |
liquid ordering, one would not expect to see any differences in |
| 1010 |
+ |
$g(r)$. Indeed, the pair distributions are essentially identical for |
| 1011 |
+ |
all of the electrostatic methods studied (for each of the different |
| 1012 |
+ |
systems under investigation). Interested readers may consult the |
| 1013 |
+ |
supplementary information for plots of these pair distribution |
| 1014 |
+ |
functions. |
| 1015 |
+ |
|
| 1016 |
+ |
A direct measure of the structural features that is a more |
| 1017 |
+ |
enlightening test of the modified electrostatic methods is the average |
| 1018 |
+ |
value of the electrostatic energy $\langle U_\mathrm{elect} \rangle$ |
| 1019 |
+ |
which is obtained by sampling the liquid-state configurations |
| 1020 |
+ |
experienced by a liquid evolving entirely under the influence of the |
| 1021 |
+ |
methods being investigated. In figure \ref{fig:Uelect} we show how |
| 1022 |
+ |
$\langle U_\mathrm{elect} \rangle$ for varies with the damping parameter, |
| 1023 |
+ |
$\alpha$, for each of the methods. |
| 1024 |
+ |
|
| 1025 |
+ |
\begin{figure} |
| 1026 |
+ |
\centering |
| 1027 |
+ |
\includegraphics[width=\textwidth]{averagePotentialEnergy_r9_12.eps} |
| 1028 |
+ |
\label{fig:Uelect} |
| 1029 |
+ |
\caption{The average electrostatic potential energy, |
| 1030 |
+ |
$\langle U_\mathrm{elect} \rangle$ for the SSDQ water with ions as a function |
| 1031 |
+ |
of the damping parameter, $\alpha$, for each of the real-space |
| 1032 |
+ |
electrostatic methods. Top panel: simulations run with a real-space |
| 1033 |
+ |
cutoff, $r_c = 9$\AA. Bottom panel: the same quantity, but with a |
| 1034 |
+ |
larger cutoff, $r_c = 12$\AA.} |
| 1035 |
|
\end{figure} |
| 1036 |
+ |
|
| 1037 |
+ |
It is clear that moderate damping is important for converging the mean |
| 1038 |
+ |
potential energy values, particularly for the two shifted force |
| 1039 |
+ |
methods (GSF and TSF). A value of $\alpha \approx 0.18$ \AA$^{-1}$ is |
| 1040 |
+ |
sufficient to converge the SP and GSF energies with a cutoff of 12 |
| 1041 |
+ |
\AA, while shorter cutoffs require more dramatic damping ($\alpha |
| 1042 |
+ |
\approx 0.36$ \AA$^{-1}$ for $r_c = 9$ \AA). It is also clear from |
| 1043 |
+ |
fig. \ref{fig:Uelect} that it is possible to overdamp the real-space |
| 1044 |
+ |
electrostatic methods, causing the estimate of the energy to drop |
| 1045 |
+ |
below the Ewald results. |
| 1046 |
+ |
|
| 1047 |
+ |
These ``optimal'' values of the damping coefficient are slightly |
| 1048 |
+ |
larger than what were observed for DSF electrostatics for purely |
| 1049 |
+ |
point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for |
| 1050 |
+ |
$r_c = 12$\AA appears to be an excellent compromise for mixed charge |
| 1051 |
+ |
multipole systems. |
| 1052 |
|
|
| 1053 |
+ |
\subsection{Reproduction of Dynamic Properties\label{sec:structure}} |
| 1054 |
+ |
To test the fidelity of the electrostatic methods at reproducing |
| 1055 |
+ |
dynamics in a multipolar liquid, it is also useful to look at |
| 1056 |
+ |
transport properties, particularly the diffusion constant, |
| 1057 |
+ |
\begin{equation} |
| 1058 |
+ |
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
| 1059 |
+ |
\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle |
| 1060 |
+ |
\label{eq:diff} |
| 1061 |
+ |
\end{equation} |
| 1062 |
+ |
which measures long-time behavior and is sensitive to the forces on |
| 1063 |
+ |
the multipoles. For the soft dipolar fluid, and the SSDQ liquid |
| 1064 |
+ |
systems, the self-diffusion constants (D) were calculated from linear |
| 1065 |
+ |
fits to the long-time portion of the mean square displacement |
| 1066 |
+ |
($\langle r^{2}(t) \rangle$).\cite{Allen87} |
| 1067 |
|
|
| 1068 |
+ |
In addition to translational diffusion, orientational relaxation times |
| 1069 |
+ |
were calculated for comparisons with the Ewald simulations and with |
| 1070 |
+ |
experiments. These values were determined from the same 1~ns $NVE$ |
| 1071 |
+ |
trajectories used for translational diffusion by calculating the |
| 1072 |
+ |
orientational time correlation function, |
| 1073 |
+ |
\begin{equation} |
| 1074 |
+ |
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t) |
| 1075 |
+ |
\cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle, |
| 1076 |
+ |
\label{eq:OrientCorr} |
| 1077 |
+ |
\end{equation} |
| 1078 |
+ |
where $P_l$ is the Legendre polynomial of order $l$ and |
| 1079 |
+ |
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along |
| 1080 |
+ |
axis $\gamma$. The body-fixed reference frame used for our |
| 1081 |
+ |
orientational correlation functions has the $z$-axis running along the |
| 1082 |
+ |
dipoles, and for the SSDQ water model, the $y$-axis connects the two |
| 1083 |
+ |
implied hydrogen atoms. |
| 1084 |
+ |
|
| 1085 |
+ |
From the orientation autocorrelation functions, we can obtain time |
| 1086 |
+ |
constants for rotational relaxation either by fitting an exponential |
| 1087 |
+ |
function or by integrating the entire correlation function. These |
| 1088 |
+ |
decay times are directly comparable to water orientational relaxation |
| 1089 |
+ |
times from nuclear magnetic resonance (NMR). The relaxation constant |
| 1090 |
+ |
obtained from $C_2^y(t)$ is normally of experimental interest because |
| 1091 |
+ |
it describes the relaxation of the principle axis connecting the |
| 1092 |
+ |
hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular |
| 1093 |
+ |
portion of the dipole-dipole relaxation from a proton NMR signal and |
| 1094 |
+ |
should provide an estimate of the NMR relaxation time |
| 1095 |
+ |
constant.\cite{Impey82} |
| 1096 |
+ |
|
| 1097 |
+ |
Results for the diffusion constants and orientational relaxation times |
| 1098 |
+ |
are shown in figure \ref{fig:dynamics}. From this data, it is apparent |
| 1099 |
+ |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
| 1100 |
+ |
reproduced with high fidelity by the GSF method. |
| 1101 |
+ |
|
| 1102 |
+ |
The $\tau_2$ results in \ref{fig:dynamics} show a much greater |
| 1103 |
+ |
difference between the real-space and the Ewald results. |
| 1104 |
+ |
|
| 1105 |
+ |
|
| 1106 |
|
\section{CONCLUSION} |
| 1107 |
|
In the first paper in this series, we generalized the |
| 1108 |
|
charge-neutralized electrostatic energy originally developed by Wolf |
