| 498 |
|
|
| 499 |
|
\begin{figure} |
| 500 |
|
\centering |
| 501 |
< |
\includegraphics[width = \linewidth]{./figures/dualLinear.pdf} |
| 501 |
> |
\includegraphics[width = 3.5in]{./figures/dualLinear.pdf} |
| 502 |
|
\caption{Example least squares regressions of the configuration energy |
| 503 |
|
differences for SPC/E water systems. The upper plot shows a data set |
| 504 |
|
with a poor correlation coefficient ($R^2$), while the lower plot |
| 2352 |
|
|
| 2353 |
|
\begin{figure} |
| 2354 |
|
\centering |
| 2355 |
< |
\includegraphics[width=5.5in]{./figures/t5peDynamics.pdf} |
| 2355 |
> |
\includegraphics[width=3.5in]{./figures/t5peDynamics.pdf} |
| 2356 |
|
\caption{Diffusion constants ({\it upper}) and reorientational time |
| 2357 |
|
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf} |
| 2358 |
|
technique compared with experiment. Data at temperatures less that |
| 2466 |
|
standard charge-dipole potential (Eq. (\ref{eq:chargeDipole})). Note |
| 2467 |
|
that this damping term is dependent upon distance and not upon |
| 2468 |
|
orientation, and that it is acting on what was originally an |
| 2469 |
< |
$r^{-3}_{ij}$ function. By writing the damped form in this manner, we |
| 2469 |
> |
$r^{-3}$ function. By writing the damped form in this manner, we |
| 2470 |
|
can collect the damping into one function and apply it to the original |
| 2471 |
|
potential when damping is desired. This works well for potentials that |
| 2472 |
|
have only one $r^{-n}$ term (where $n$ is an odd positive integer); |
| 2520 |
|
generating function, |
| 2521 |
|
\begin{equation} |
| 2522 |
|
c_n(r_{ij}) = \frac{2^n(\alpha r_{ij})^{2n-1}e^{-\alpha^2r^2_{ij}}} |
| 2523 |
< |
{\sqrt{\pi}(2n-1)!!} + c_{n-1}(r_{ij}), |
| 2523 |
> |
{(2n-1)!!\sqrt{\pi}} + c_{n-1}(r_{ij}), |
| 2524 |
|
\label{eq:dampingGeneratingFunc} |
| 2525 |
|
\end{equation} |
| 2526 |
|
where, |
| 2586 |
|
In order to find these optimal values, we mapped out the static |
| 2587 |
|
dielectric constant as a function of both the damping parameter and |
| 2588 |
|
cutoff radius for several different water models. To calculate the |
| 2589 |
< |
static dielectric constant, we performed 5ns $NPT$ calculations at 9, |
| 2590 |
< |
10, 11, and 12\AA cutoff radii, each with damping parameter values |
| 2591 |
< |
ranging from 0 to 0.35\AA$^{-1}$ using the TIP5P-E, TIP4P-Ew, SPC/E, |
| 2589 |
> |
static dielectric constant, we performed 5ns $NPT$ calculations on |
| 2590 |
> |
systems of 512 water molecules, using the TIP5P-E, TIP4P-Ew, SPC/E, |
| 2591 |
|
and SSD/RF water models. TIP4P-Ew is a reparametrized version of the |
| 2592 |
|
four-point transferable intermolecular potential (TIP4P) for water |
| 2593 |
|
targeted for use with the Ewald summation.\cite{Horn04} SSD/RF is the |
| 2594 |
|
reaction field modified variant of the soft sticky dipole (SSD) model |
| 2595 |
< |
for water, and this model is discussed in more detail in the next |
| 2596 |
< |
chapter. One thing to note about it, electrostatic interactions are |
| 2597 |
< |
handled via dipole-dipole interactions rather than charge-charge |
| 2598 |
< |
interactions like the other three models. Damping of the dipole-dipole |
| 2599 |
< |
interaction was handled as described in section |
| 2600 |
< |
\ref{sec:dampingMultipoles}. |
| 2595 |
> |
for water\cite{Fennell04} This model is discussed in more detail in |
| 2596 |
> |
the next chapter. One thing to note about it, electrostatic |
| 2597 |
> |
interactions are handled via dipole-dipole interactions rather than |
| 2598 |
> |
charge-charge interactions like the other three models. Damping of the |
| 2599 |
> |
dipole-dipole interaction was handled as described in section |
| 2600 |
> |
\ref{sec:dampingMultipoles}. Each of these systems were studied with |
| 2601 |
> |
cutoff radii of 9, 10, 11, and 12\AA\ and with damping parameter values |
| 2602 |
> |
ranging from 0 to 0.35\AA$^{-1}$. |
| 2603 |
|
\begin{figure} |
| 2604 |
|
\centering |
| 2605 |
< |
\includegraphics[width=3.5in]{./figures/dielectricMap.pdf} |
| 2606 |
< |
\caption{The static dielectric constant for the A: TIP5P-E, B: TIP4P-Ew, |
| 2607 |
< |
C: SPC/E, and D: SSD/RF water models as a function of cutoff radius |
| 2608 |
< |
($R_\textrm{c}$) and damping coefficient ($\alpha$).} |
| 2605 |
> |
\includegraphics[width=\linewidth]{./figures/dielectricMap.pdf} |
| 2606 |
> |
\caption{The static dielectric constant for the TIP5P-E (A), TIP4P-Ew |
| 2607 |
> |
(B), SPC/E (C), and SSD/RF (D) water models as a function of cutoff |
| 2608 |
> |
radius ($R_\textrm{c}$) and damping coefficient ($\alpha$).} |
| 2609 |
|
\label{fig:dielectricMap} |
| 2610 |
|
\end{figure} |
| 2611 |
|
|
| 2613 |
|
\ref{fig:dielectricMap} in the form of shaded contour plots. An |
| 2614 |
|
interesting aspect of all four contour plots is that the dielectric |
| 2615 |
|
constant is effectively linear with respect to $\alpha$ and |
| 2616 |
< |
$R_\textrm{c}$ in the low to moderate damping regions. Another point |
| 2617 |
< |
to note is that choosing $\alpha$ and $R_\textrm{c}$ identical to |
| 2618 |
< |
those used in studies with the Ewald summation results in the same |
| 2619 |
< |
calculated dielectric constant. As an example, in the paper outlining |
| 2620 |
< |
the development of TIP5P-E, the real-space cutoff and Ewald |
| 2621 |
< |
coefficient were tethered to the system size, and for a 512 molecule |
| 2622 |
< |
system are approximately 12\AA and 0.25\AA$^{-1}$ |
| 2623 |
< |
respectively.\cite{Rick04} These parameters resulted in a dielectric |
| 2624 |
< |
constant of 92$\pm$14, while with {\sc sf} these parameters give a |
| 2625 |
< |
dielectric constant of 90.8$\pm$0.9. Another example comes from the |
| 2626 |
< |
TIP4P-Ew paper where $\alpha$ and $R_\textrm{c}$ were chosen to be |
| 2627 |
< |
9.5\AA and 0.35\AA$^{-1}$, and these parameters resulted in a |
| 2628 |
< |
$\epsilon_0$ equal to 63$\pm$1.\cite{Horn04} We did not perform |
| 2629 |
< |
calculations with these exact parameters, but interpolating between |
| 2630 |
< |
surrounding values gives a $\epsilon_0$ of 61$\pm$1. Seeing a |
| 2631 |
< |
dependence of the dielectric constant on $\alpha$ and $R_\textrm{c}$ |
| 2632 |
< |
with the {\sc sf} technique, it might be interesting to investigate |
| 2633 |
< |
the dielectric dependence when using the Ewald summation. |
| 2634 |
< |
|
| 2634 |
< |
|
| 2616 |
> |
$R_\textrm{c}$ in the low to moderate damping regions, and the slope |
| 2617 |
> |
is 0.025\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is |
| 2618 |
> |
that choosing $\alpha$ and $R_\textrm{c}$ identical to those used in |
| 2619 |
> |
studies with the Ewald summation results in the same calculated |
| 2620 |
> |
dielectric constant. As an example, in the paper outlining the |
| 2621 |
> |
development of TIP5P-E, the real-space cutoff and Ewald coefficient |
| 2622 |
> |
were tethered to the system size, and for a 512 molecule system are |
| 2623 |
> |
approximately 12\AA\ and 0.25\AA$^{-1}$ respectively.\cite{Rick04} |
| 2624 |
> |
These parameters resulted in a dielectric constant of 92$\pm$14, while |
| 2625 |
> |
with {\sc sf} these parameters give a dielectric constant of |
| 2626 |
> |
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where |
| 2627 |
> |
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5\AA\ and |
| 2628 |
> |
0.35\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal |
| 2629 |
> |
to 63$\pm$1.\cite{Horn04} We did not perform calculations with these |
| 2630 |
> |
exact parameters, but interpolating between surrounding values gives a |
| 2631 |
> |
$\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric |
| 2632 |
> |
constant on $\alpha$ and $R_\textrm{c}$ with the {\sc sf} technique, |
| 2633 |
> |
it might be interesting to investigate the dielectric dependence of |
| 2634 |
> |
the real-space Ewald parameters. |
| 2635 |
|
|
| 2636 |
+ |
Although it is tempting to choose damping parameters equivalent to |
| 2637 |
+ |
these Ewald examples, the results discussed in sections |
| 2638 |
+ |
\ref{sec:EnergyResults} through \ref{sec:IndividualResults} indicate |
| 2639 |
+ |
that values this high are destructive to both the energetics and |
| 2640 |
+ |
dynamics. Ideally, $\alpha$ should not exceed 0.3\AA$^{-1}$ for any of |
| 2641 |
+ |
the cutoff values in this range. If the optimal damping parameter is |
| 2642 |
+ |
chosen to be midway between 0.275 and 0.3\AA$^{-1}$ (0.2875\AA$^{-1}$) |
| 2643 |
+ |
at the 9\AA\ cutoff, then the linear trend with $R_\textrm{c}$ will |
| 2644 |
+ |
always keep $\alpha$ below 0.3\AA$^{-1}$. This linear progression |
| 2645 |
+ |
would give values of 0.2875, 0.2625, 0.2375, and 0.2125\AA$^{-1}$ for |
| 2646 |
+ |
cutoff radii of 9, 10, 11, and 12\AA. Setting this to be the default |
| 2647 |
+ |
behavior for the damped {\sc sf} technique will result in consistent |
| 2648 |
+ |
dielectric behavior for these and other condensed molecular systems, |
| 2649 |
+ |
regardless of the chosen cutoff radius. The static dielectric |
| 2650 |
+ |
constants for TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be |
| 2651 |
+ |
approximately fixed at 74, 52, 58, and 89 respectively. These values |
| 2652 |
+ |
are generally lower than the values reported in the literature; |
| 2653 |
+ |
however, the relative dielectric behavior scales as expected when |
| 2654 |
+ |
comparing the models to one another. |
| 2655 |
|
|
| 2656 |
|
\section{Conclusions}\label{sec:PairwiseConclusions} |
| 2657 |
|
|
| 2665 |
|
energetic and dynamic characteristics exhibited by simulations |
| 2666 |
|
employing lattice summation techniques. The cumulative energy |
| 2667 |
|
difference results showed the undamped {\sc sf} and moderately damped |
| 2668 |
< |
{\sc sp} methods produced results nearly identical to {\sc spme}. |
| 2669 |
< |
Similarly for the dynamic features, the undamped or moderately damped |
| 2670 |
< |
{\sc sf} and moderately damped {\sc sp} methods produce force and |
| 2671 |
< |
torque vector magnitude and directions very similar to the expected |
| 2672 |
< |
values. These results translate into long-time dynamic behavior |
| 2673 |
< |
equivalent to that produced in simulations using {\sc spme}. |
| 2668 |
> |
{\sc sp} methods produced results nearly identical to the Ewald |
| 2669 |
> |
summation. Similarly for the dynamic features, the undamped or |
| 2670 |
> |
moderately damped {\sc sf} and moderately damped {\sc sp} methods |
| 2671 |
> |
produce force and torque vector magnitude and directions very similar |
| 2672 |
> |
to the expected values. These results translate into long-time |
| 2673 |
> |
dynamic behavior equivalent to that produced in simulations using the |
| 2674 |
> |
Ewald summation. A detailed study of water simulations showed that |
| 2675 |
> |
liquid properties calculated when using {\sc sf} will also be |
| 2676 |
> |
equivalent to those obtained using the Ewald summation. |
| 2677 |
|
|
| 2678 |
|
As in all purely-pairwise cutoff methods, these methods are expected |
| 2679 |
|
to scale approximately {\it linearly} with system size, and they are |
| 2701 |
|
data of interest is either structural or short-time dynamical |
| 2702 |
|
quantities. For long-time dynamics and collective motions, the safest |
| 2703 |
|
pairwise method we have evaluated is the {\sc sf} method with an |
| 2704 |
< |
electrostatic damping between 0.2 and 0.25\AA$^{-1}$. |
| 2704 |
> |
electrostatic damping between 0.2 and 0.25\AA$^{-1}$. It is also |
| 2705 |
> |
important to note that the static dielectric constant in water |
| 2706 |
> |
simulations is highly dependent on both $\alpha$ and |
| 2707 |
> |
$R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc |
| 2708 |
> |
sf} method should use an $\alpha$ of 0.2175\AA$^{-1}$ for an |
| 2709 |
> |
$R_\textrm{c}$ of 12\AA, and $\alpha$ should decrease by |
| 2710 |
> |
0.025\AA$^{-1}$ for every 1\AA\ increase in cutoff radius. |
| 2711 |
|
|
| 2712 |
|
We are not suggesting that there is any flaw with the Ewald sum; in |
| 2713 |
|
fact, it is the standard by which these simple pairwise sums have been |
