| 1 |
< |
\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES} |
| 1 |
> |
\chapter[ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES]{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION TECHNIQUES} |
| 2 |
|
|
| 3 |
|
In molecular simulations, proper accumulation of electrostatic |
| 4 |
|
interactions is essential and one of the most |
| 1067 |
|
long-range electrostatic |
| 1068 |
|
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
| 1069 |
|
Without this correction, the pressure term on the central particle |
| 1070 |
< |
from the surroundings is missing. Because they expand to compensate |
| 1071 |
< |
for this added pressure term when this correction is included, systems |
| 1072 |
< |
composed of these particles tend to under-predict the density of water |
| 1073 |
< |
under standard conditions. When using any form of long-range |
| 1074 |
< |
electrostatic correction, it has become common practice to develop or |
| 1075 |
< |
utilize a reparametrized water model that corrects for this |
| 1070 |
> |
from the surroundings is missing. When this correction is included, |
| 1071 |
> |
systems of these particles expand to compensate for this added |
| 1072 |
> |
pressure term and under-predict the density of water under standard |
| 1073 |
> |
conditions. When using any form of long-range electrostatic |
| 1074 |
> |
correction, it has become common practice to develop or utilize a |
| 1075 |
> |
reparametrized water model that corrects for this |
| 1076 |
|
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 1077 |
< |
this practice and was optimized specifically for use with the Ewald |
| 1077 |
> |
this practice and was optimized for use with the Ewald |
| 1078 |
|
summation.\cite{Rick04} In his publication, Rick preserved the |
| 1079 |
|
geometry and point charge magnitudes in TIP5P and focused on altering |
| 1080 |
< |
the Lennard-Jones parameters to correct the density at |
| 1081 |
< |
298K.\cite{Rick04} With the density corrected, he compared common |
| 1082 |
< |
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1083 |
< |
9~\AA\ cutoff. |
| 1080 |
> |
the Lennard-Jones parameters to correct the density at 298~K. With the |
| 1081 |
> |
density corrected, he compared common water properties for TIP5P-E |
| 1082 |
> |
using the Ewald sum with TIP5P using a 9~\AA\ cutoff. |
| 1083 |
|
|
| 1084 |
|
In the following sections, we compared these same water properties |
| 1085 |
|
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1112 |
|
brackets of equation |
| 1113 |
|
\ref{eq:MolecularPressure}) is directly dependent on the interatomic |
| 1114 |
|
forces. Since the {\sc sp} method does not modify the forces (see |
| 1115 |
< |
section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} |
| 1115 |
> |
section \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} |
| 1116 |
|
will be identical to that obtained without an electrostatic |
| 1117 |
|
correction. The {\sc sf} method does alter the virial component and, |
| 1118 |
|
by way of the modified pressures, should provide densities more in |
| 1126 |
|
temperatures. The average densities were calculated from the later |
| 1127 |
|
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 1128 |
|
method for accumulating statistics, these sequences were spliced into |
| 1129 |
< |
200 segments to calculate the average density and standard deviation |
| 1130 |
< |
at each temperature.\cite{Mahoney00} |
| 1129 |
> |
200 segments, each providing an average density. These 200 density |
| 1130 |
> |
values were used to calculate the average and standard deviation of |
| 1131 |
> |
the density at each temperature.\cite{Mahoney00} |
| 1132 |
|
|
| 1133 |
|
\begin{figure} |
| 1134 |
|
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 1148 |
|
TIP5P-E using differing electrostatic corrections overlaid on the |
| 1149 |
|
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 1150 |
|
technique are close to, though typically lower than, those calculated |
| 1151 |
< |
while using the Ewald summation. These slightly reduced densities |
| 1152 |
< |
indicate that the pressure component from the image charges at |
| 1153 |
< |
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 1154 |
< |
portion of the Ewald summation. Bringing the image charges closer to |
| 1155 |
< |
the central particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than |
| 1156 |
< |
the preferred 12~\AA\ R$_\textrm{c}$) increases the strength of their |
| 1157 |
< |
interactions, resulting in a further reduction of the densities. |
| 1151 |
> |
using the Ewald summation. These slightly reduced densities indicate |
| 1152 |
> |
that the pressure component from the image charges at R$_\textrm{c}$ |
| 1153 |
> |
is larger than that exerted by the reciprocal-space portion of the |
| 1154 |
> |
Ewald summation. Bringing the image charges closer to the central |
| 1155 |
> |
particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the |
| 1156 |
> |
preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image |
| 1157 |
> |
charge interactions on the central particle and results in a further |
| 1158 |
> |
reduction of the densities. |
| 1159 |
|
|
| 1160 |
|
Because the strength of the image charge interactions has a noticeable |
| 1161 |
|
effect on the density, we would expect the use of electrostatic |
| 1211 |
|
non-polarizable models.\cite{Sorenson00} This excellent agreement with |
| 1212 |
|
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To |
| 1213 |
|
check whether the choice of using the Ewald summation or the {\sc sf} |
| 1214 |
< |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K |
| 1215 |
< |
and 1atm were determined for the systems compared in the previous |
| 1216 |
< |
section. |
| 1214 |
> |
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at |
| 1215 |
> |
298~K and 1~atm were determined for the systems compared in the |
| 1216 |
> |
previous section. |
| 1217 |
|
|
| 1218 |
|
\begin{figure} |
| 1219 |
|
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 1226 |
|
\end{figure} |
| 1227 |
|
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc |
| 1228 |
|
sf} technique with a various parameters are overlaid on the |
| 1229 |
< |
$g_\textrm{OO}(r)$ while using the Ewald summation in figure |
| 1230 |
< |
\ref{fig:t5peGofRs}. The differences in density do not appear to have |
| 1231 |
< |
any effect on the liquid structure as the $g_\textrm{OO}(r)$s are |
| 1232 |
< |
indistinguishable. These results indicate that the $g_\textrm{OO}(r)$ |
| 1233 |
< |
is insensitive to the choice of electrostatic correction. |
| 1229 |
> |
$g_\textrm{OO}(r)$ while using the Ewald summation in |
| 1230 |
> |
figure~\ref{fig:t5peGofRs}. The differences in density do not appear |
| 1231 |
> |
to have any effect on the liquid structure as the $g_\textrm{OO}(r)$s |
| 1232 |
> |
are indistinguishable. These results indicate that the |
| 1233 |
> |
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic |
| 1234 |
> |
correction. |
| 1235 |
|
|
| 1236 |
|
\subsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 1237 |
|
|
| 1246 |
|
good set for comparisons involving the {\sc sf} technique. |
| 1247 |
|
|
| 1248 |
|
The $\Delta H_\textrm{vap}$ is the enthalpy change required to |
| 1249 |
< |
transform one mol of substance from the liquid phase to the gas |
| 1249 |
> |
transform one mole of substance from the liquid phase to the gas |
| 1250 |
|
phase.\cite{Berry00} In molecular simulations, this quantity can be |
| 1251 |
|
determined via |
| 1252 |
|
\begin{equation} |
| 1361 |
|
|
| 1362 |
|
As observed for the density in section \ref{sec:t5peDensity}, the |
| 1363 |
|
property trends with temperature seen when using the Ewald summation |
| 1364 |
< |
are reproduced with the {\sc sf} technique. One noticable difference |
| 1364 |
> |
are reproduced with the {\sc sf} technique. One noticeable difference |
| 1365 |
|
between the properties calculated using the two methods are the lower |
| 1366 |
|
$\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be |
| 1367 |
|
expected due to the direct weakening of the electrostatic interaction |
| 1372 |
|
steeper than the experimental trend, indirectly resulting in the |
| 1373 |
|
inflated $C_p$ values at all temperatures. |
| 1374 |
|
|
| 1375 |
< |
Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ |
| 1376 |
< |
values all overlap within error. As indicated for the $\Delta |
| 1377 |
< |
H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph, |
| 1378 |
< |
the deviations between experiment and simulation in this region are |
| 1379 |
< |
not the fault of the electrostatic summation methods but are due to |
| 1380 |
< |
the TIP5P class model itself. Like most rigid, non-polarizable, |
| 1381 |
< |
point-charge water models, the density decreases with temperature at a |
| 1382 |
< |
much faster rate than experiment (see figure |
| 1383 |
< |
\ref{fig:t5peDensities}). The reduced density leads to the inflated |
| 1375 |
> |
Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values |
| 1376 |
> |
all overlap within error. As indicated for the $\Delta H_\textrm{vap}$ |
| 1377 |
> |
and $C_p$ results discussed in the previous paragraph, the deviations |
| 1378 |
> |
between experiment and simulation in this region are not the fault of |
| 1379 |
> |
the electrostatic summation methods but are due to the geometry and |
| 1380 |
> |
parameters of the TIP5P class of water models. Like most rigid, |
| 1381 |
> |
non-polarizable, point-charge water models, the density decreases with |
| 1382 |
> |
temperature at a much faster rate than experiment (see figure |
| 1383 |
> |
\ref{fig:t5peDensities}). This reduced density leads to the inflated |
| 1384 |
|
compressibility and expansivity values at higher temperatures seen |
| 1385 |
|
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability |
| 1386 |
< |
and many-body effects are required in order for simulation to overcome |
| 1387 |
< |
these differences with experiment.\cite{Laasonen93,Donchev06} |
| 1386 |
> |
and many-body effects are required in order for water models to |
| 1387 |
> |
overcome differences between simulation-based and experimentally |
| 1388 |
> |
determined densities at these higher |
| 1389 |
> |
temperatures.\cite{Laasonen93,Donchev06} |
| 1390 |
|
|
| 1391 |
|
At temperatures below the freezing point for experimental water, the |
| 1392 |
|
differences between {\sc sf} and the Ewald summation results are more |
| 1395 |
|
particularly in the case of {\sc sf} without damping. This points to |
| 1396 |
|
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 1397 |
|
temperatures below 250~K in the {\sc sf} simulations, indicating that |
| 1398 |
< |
disorder in the reciprical-space term of the Ewald summation might act |
| 1398 |
> |
disorder in the reciprocal-space term of the Ewald summation might act |
| 1399 |
|
to loosen up the local structure more than the image-charges in {\sc |
| 1400 |
< |
sf}. Because the systems are locked in different regions of |
| 1401 |
< |
phase-space, comparisons between properties at these temperatures are |
| 1402 |
< |
not exactly fair. This observation is explored in more detail in |
| 1403 |
< |
section \ref{sec:t5peDynamics}. |
| 1400 |
> |
sf}. The damped {\sc sf} actually makes a better comparison with |
| 1401 |
> |
experiment in this region, particularly for the $\alpha_p$ values. The |
| 1402 |
> |
local interactions in the undamped {\sc sf} technique appear to be too |
| 1403 |
> |
strong since the property change is much more dramatic than the damped |
| 1404 |
> |
forms, while the Ewald summation appears to weight the |
| 1405 |
> |
reciprocal-space interactions at the expense the local interactions, |
| 1406 |
> |
disagreeing with the experimental results. This observation is |
| 1407 |
> |
explored in more detail in section \ref{sec:t5peDynamics}. |
| 1408 |
|
|
| 1409 |
|
The final thermodynamic property displayed in figure |
| 1410 |
|
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
| 1414 |
|
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the |
| 1415 |
|
converged $\epsilon$ values accumulated for the {\sc sf} |
| 1416 |
|
simulations. Lack of a damping function results in dielectric |
| 1417 |
< |
constants significantly smaller than that obtained using the Ewald |
| 1417 |
> |
constants significantly smaller than those obtained using the Ewald |
| 1418 |
|
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the |
| 1419 |
|
agreement considerably. It should be noted that the choice of the |
| 1420 |
|
``Ewald coefficient'' value also has a significant effect on the |
| 1426 |
|
sf}; however, the choice of cutoff radius also plays an important |
| 1427 |
|
role. In section \ref{sec:dampingDielectric}, this connection is |
| 1428 |
|
further explored as optimal damping coefficients for different choices |
| 1429 |
< |
of $R_\textrm{c}$ are determined for {\sc sf} for capturing the |
| 1430 |
< |
dielectric behavior. |
| 1429 |
> |
of $R_\textrm{c}$ are determined for {\sc sf} in order to best capture |
| 1430 |
> |
the dielectric behavior. |
| 1431 |
|
|
| 1432 |
|
\subsection{Dynamic Properties}\label{sec:t5peDynamics} |
| 1433 |
|
|
| 1434 |
|
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 1435 |
< |
method, 200~ps $NVE$ simulations were performed for each temperature at |
| 1436 |
< |
the average density reported by the $NPT$ simulations. The |
| 1437 |
< |
self-diffusion constants ($D$) were calculated with the Einstein |
| 1438 |
< |
relation using the mean square displacement (MSD), |
| 1435 |
> |
method, 200~ps $NVE$ simulations were performed for each temperature |
| 1436 |
> |
at the average density reported by the $NPT$ simulations. The |
| 1437 |
> |
self-diffusion constants ($D$) were calculated using the mean square |
| 1438 |
> |
displacement (MSD) form of the Einstein relation, |
| 1439 |
|
\begin{equation} |
| 1440 |
|
D = \lim_{t\rightarrow\infty} |
| 1441 |
|
\frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t}, |
| 1451 |
|
\item linear diffusive regime, and |
| 1452 |
|
\item a region with poor statistics. |
| 1453 |
|
\end{enumerate} |
| 1454 |
< |
The slope from the linear region (region 2) is used to calculate $D$. |
| 1454 |
> |
The slope from the linear regime (region 2) is used to calculate $D$. |
| 1455 |
|
\begin{figure} |
| 1456 |
|
\centering |
| 1457 |
|
\includegraphics[width=3.5in]{./figures/ExampleMSD.pdf} |
| 1470 |
|
labeled frame axes.} |
| 1471 |
|
\label{fig:waterFrame} |
| 1472 |
|
\end{figure} |
| 1473 |
< |
In addition to translational diffusion, reorientational time constants |
| 1473 |
> |
In addition to translational diffusion, orientational relaxation times |
| 1474 |
|
were calculated for comparisons with the Ewald simulations and with |
| 1475 |
< |
experiments. These values were determined from 25~ps $NVE$ trajectories |
| 1476 |
< |
through calculation of the orientational time correlation function, |
| 1475 |
> |
experiments. These values were determined from 25~ps $NVE$ |
| 1476 |
> |
trajectories through calculation of the orientational time correlation |
| 1477 |
> |
function, |
| 1478 |
|
\begin{equation} |
| 1479 |
|
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t) |
| 1480 |
|
\cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle, |
| 1541 |
|
relaxes faster than experiment with the Ewald sum while tracking |
| 1542 |
|
experiment fairly well when using the {\sc sf} technique, independent |
| 1543 |
|
of the choice of damping constant. Their are several possible reasons |
| 1544 |
< |
for this deviation between techniques. The Ewald results were taken |
| 1545 |
< |
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick |
| 1546 |
< |
calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of |
| 1547 |
< |
0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, |
| 1548 |
< |
placing the result more in line with that obtained using the Ewald |
| 1549 |
< |
sum. These results support this explanation; however, recomputing the |
| 1550 |
< |
results to meet a poorer statistical standard is |
| 1544 |
> |
for this deviation between techniques. The Ewald results were |
| 1545 |
> |
calculated using shorter (10ps) trajectories than the {\sc sf} results |
| 1546 |
> |
(25ps). A quick calculation from a 10~ps trajectory with {\sc sf} with |
| 1547 |
> |
an $\alpha$ of 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in |
| 1548 |
> |
$\tau_2^y$, placing the result more in line with that obtained using |
| 1549 |
> |
the Ewald sum. This example supports this explanation; however, |
| 1550 |
> |
recomputing the results to meet a poorer statistical standard is |
| 1551 |
|
counter-productive. Assuming the Ewald results are not the product of |
| 1552 |
|
poor statistics, differences in techniques to integrate the |
| 1553 |
|
orientational motion could also play a role. {\sc shake} is the most |
| 1554 |
|
commonly used technique for approximating rigid-body orientational |
| 1555 |
< |
motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and |
| 1555 |
> |
motion,\cite{Ryckaert77} whereas in {\sc oopse}, we maintain and |
| 1556 |
|
integrate the entire rotation matrix using the {\sc dlm} |
| 1557 |
|
method.\cite{Meineke05} Since {\sc shake} is an iterative constraint |
| 1558 |
|
technique, if the convergence tolerances are raised for increased |
| 1580 |
|
multipoles. In a mixed system of monopoles and multipoles, the |
| 1581 |
|
undamped {\sc sf} potential needs only to shift the force terms of the |
| 1582 |
|
monopole (and use the monopole potential of equation (\ref{eq:SFPot})) |
| 1583 |
< |
and smoothly cutoff the multipole interactions with a switching |
| 1583 |
> |
and smoothly truncate the multipole interactions with a switching |
| 1584 |
|
function. The switching function is required in order to conserve |
| 1585 |
< |
energy, because a discontinuity will exist at $R_\textrm{c}$ in the |
| 1586 |
< |
absence of shifting terms. |
| 1585 |
> |
energy, because a discontinuity will exist in both the potential and |
| 1586 |
> |
forces at $R_\textrm{c}$ in the absence of shifting terms. |
| 1587 |
|
|
| 1588 |
|
If we consider damping the {\sc sf} potential (Eq. (\ref{eq:DSFPot})), |
| 1589 |
|
then we need to incorporate the complimentary error function term into |
| 1674 |
|
\end{equation} |
| 1675 |
|
Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional |
| 1676 |
|
term. Continuing with higher rank tensors, we can obtain the damping |
| 1677 |
< |
functions for higher multipoles as well as the forces. Each subsequent |
| 1677 |
> |
functions for higher multipole potentials and forces. Each subsequent |
| 1678 |
|
damping function includes one additional term, and we can simplify the |
| 1679 |
|
procedure for obtaining these terms by writing out the following |
| 1680 |
|
generating function, |
| 1708 |
|
c_1(r_{ij}), |
| 1709 |
|
\label{eq:dampDipoleDipole} |
| 1710 |
|
\end{equation} |
| 1711 |
< |
$c_2(r_{ij})$ and $c_1(r_{ij})$ respectively dampen these two |
| 1712 |
< |
parts. The forces for the damped dipole-dipole interaction, |
| 1711 |
> |
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts |
| 1712 |
> |
respectively. The forces for the damped dipole-dipole interaction, |
| 1713 |
|
\begin{equation} |
| 1714 |
|
\begin{split} |
| 1715 |
|
F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
| 1739 |
|
constant is calculated from the long-time fluctuations of the system's |
| 1740 |
|
accumulated dipole moment (Eq. (\ref{eq:staticDielectric})), so it is |
| 1741 |
|
going to be quite sensitive to the choice of damping parameter. We |
| 1742 |
< |
would like to choose an optimal damping constant for any particular |
| 1743 |
< |
cutoff radius choice that would properly capture the dielectric |
| 1744 |
< |
behavior of the liquid. |
| 1742 |
> |
would like to choose optimal damping constants such that any arbitrary |
| 1743 |
> |
choice of cutoff radius will properly capture the dielectric behavior |
| 1744 |
> |
of the liquid. |
| 1745 |
|
|
| 1746 |
|
In order to find these optimal values, we mapped out the static |
| 1747 |
|
dielectric constant as a function of both the damping parameter and |
| 1752 |
|
four-point transferable intermolecular potential (TIP4P) for water |
| 1753 |
|
targeted for use with the Ewald summation.\cite{Horn04} SSD/RF is the |
| 1754 |
|
reaction field modified variant of the soft sticky dipole (SSD) model |
| 1755 |
< |
for water\cite{Fennell04} This model is discussed in more detail in |
| 1755 |
> |
for water.\cite{Fennell04} This model is discussed in more detail in |
| 1756 |
|
the next chapter. One thing to note about it, electrostatic |
| 1757 |
|
interactions are handled via dipole-dipole interactions rather than |
| 1758 |
|
charge-charge interactions like the other three models. Damping of the |
| 1785 |
|
with {\sc sf} these parameters give a dielectric constant of |
| 1786 |
|
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where |
| 1787 |
|
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and |
| 1788 |
< |
0.35~\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal |
| 1789 |
< |
to 63$\pm$1.\cite{Horn04} We did not perform calculations with these |
| 1790 |
< |
exact parameters, but interpolating between surrounding values gives a |
| 1791 |
< |
$\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric |
| 1792 |
< |
constant on $\alpha$ and $R_\textrm{c}$ with the {\sc sf} technique, |
| 1793 |
< |
it might be interesting to investigate the dielectric dependence of |
| 1794 |
< |
the real-space Ewald parameters. |
| 1788 |
> |
0.35~\AA$^{-1}$, and these parameters resulted in a dielectric |
| 1789 |
> |
constant equal to 63$\pm$1.\cite{Horn04} We did not perform |
| 1790 |
> |
calculations with these exact parameters, but interpolating between |
| 1791 |
> |
surrounding values gives a dielectric constant of 61$\pm$1. Since the |
| 1792 |
> |
dielectric constant is dependent on $\alpha$ and $R_\textrm{c}$ with |
| 1793 |
> |
the {\sc sf} technique, it might be interesting to investigate the |
| 1794 |
> |
dielectric dependence of the real-space Ewald parameters. |
| 1795 |
|
|
| 1796 |
|
Although it is tempting to choose damping parameters equivalent to |
| 1797 |
|
these Ewald examples, the results discussed in sections |
| 1825 |
|
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1826 |
|
energetic and dynamic characteristics exhibited by simulations |
| 1827 |
|
employing lattice summation techniques. The cumulative energy |
| 1828 |
< |
difference results showed the undamped {\sc sf} and moderately damped |
| 1829 |
< |
{\sc sp} methods produced results nearly identical to the Ewald |
| 1828 |
> |
difference results showed that the undamped {\sc sf} and moderately |
| 1829 |
> |
damped {\sc sp} methods produce results nearly identical to the Ewald |
| 1830 |
|
summation. Similarly for the dynamic features, the undamped or |
| 1831 |
|
moderately damped {\sc sf} and moderately damped {\sc sp} methods |
| 1832 |
|
produce force and torque vector magnitude and directions very similar |
| 1839 |
|
As in all purely-pairwise cutoff methods, these methods are expected |
| 1840 |
|
to scale approximately {\it linearly} with system size, and they are |
| 1841 |
|
easily parallelizable. This should result in substantial reductions |
| 1842 |
< |
in the computational cost of performing large simulations. |
| 1842 |
> |
in the computational cost associated with large-scale simulations. |
| 1843 |
|
|
| 1844 |
|
Aside from the computational cost benefit, these techniques have |
| 1845 |
|
applicability in situations where the use of the Ewald sum can prove |
| 1858 |
|
systems containing point charges, most structural features will be |
| 1859 |
|
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1860 |
|
method with an electrostatic damping of 0.2~\AA$^{-1}$. These methods |
| 1861 |
< |
would also be appropriate for molecular dynamics simulations where the |
| 1862 |
< |
data of interest is either structural or short-time dynamical |
| 1861 |
> |
would also be appropriate in molecular dynamics simulations where the |
| 1862 |
> |
data of interest are either structural or short-time dynamical |
| 1863 |
|
quantities. For long-time dynamics and collective motions, the safest |
| 1864 |
|
pairwise method we have evaluated is the {\sc sf} method with an |
| 1865 |
|
electrostatic damping between 0.2 and 0.25~\AA$^{-1}$. It is also |
| 1868 |
|
$R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc |
| 1869 |
|
sf} method should use an $\alpha$ of 0.2175~\AA$^{-1}$ for an |
| 1870 |
|
$R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease by |
| 1871 |
< |
0.025~\AA$^{-1}$ for every 1~\AA\ increase in cutoff radius. |
| 1871 |
> |
0.025~\AA$^{-1}$ for every 1~\AA\ increase in the cutoff radius. |
| 1872 |
|
|
| 1873 |
|
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1874 |
< |
fact, it is the standard by which these simple pairwise sums have been |
| 1875 |
< |
judged. However, these results do suggest that in the typical |
| 1874 |
> |
fact, it is the standard by which these simple pairwise methods have |
| 1875 |
> |
been judged. However, these results do suggest that in the typical |
| 1876 |
|
simulations performed today, the Ewald summation may no longer be |
| 1877 |
|
required to obtain the level of accuracy most researchers have come to |
| 1878 |
|
expect. |
