| 165 |
|
formulation, the total energy for the charge and image were not equal |
| 166 |
|
to the integral of the force expression, and as a result, the total |
| 167 |
|
energy would not be conserved in molecular dynamics (MD) |
| 168 |
< |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
| 168 |
> |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennell and |
| 169 |
|
Gezelter later proposed shifted force variants of the Wolf method with |
| 170 |
|
commensurate force and energy expressions that do not exhibit this |
| 171 |
|
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
| 223 |
|
Even at elevated temperatures, there is local charge balance in an |
| 224 |
|
ionic liquid, where each positive ion has surroundings dominated by |
| 225 |
|
negative ions and vice versa. The reversed-charge images on the |
| 226 |
< |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
| 227 |
< |
the effective multipolar interactions as the charges traverse the |
| 228 |
< |
cutoff boundary. |
| 226 |
> |
cutoff sphere that are integral to the Wolf and damped shifted force |
| 227 |
> |
(DSF) approaches retain the effective multipolar interactions as the |
| 228 |
> |
charges traverse the cutoff boundary. |
| 229 |
|
|
| 230 |
|
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
| 231 |
|
significant orientational averaging that additionally reduces the |
| 232 |
|
effect of long-range multipolar interactions. The image multipoles |
| 233 |
< |
that are introduced in the TSF, GSF, and SP methods mimic this effect |
| 233 |
> |
that are introduced in the Taylor shifted force (TSF), gradient |
| 234 |
> |
shifted force (GSF), and shifted potential (SP) methods mimic this effect |
| 235 |
|
and reduce the effective range of the multipolar interactions as |
| 236 |
|
interacting molecules traverse each other's cutoff boundaries. |
| 237 |
|
|
| 238 |
|
Forces and torques acting on atomic sites are fundamental in driving |
| 239 |
< |
dynamics in molecular simulations, and the damped shifted force (DSF) |
| 240 |
< |
energy kernel provides consistent energies and forces on charged atoms |
| 241 |
< |
within the cutoff sphere. Both the energy and the force go smoothly to |
| 242 |
< |
zero as an atom approaches the cutoff radius. The comparisons of the |
| 243 |
< |
accuracy these quantities between the DSF kernel and SPME was |
| 244 |
< |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
| 245 |
< |
has seen increasing use in molecular systems with relatively uniform |
| 245 |
< |
charge |
| 239 |
> |
dynamics in molecular simulations, and the DSF energy kernel provides |
| 240 |
> |
consistent energies and forces on charged atoms within the cutoff |
| 241 |
> |
sphere. Both the energy and the force go smoothly to zero as an atom |
| 242 |
> |
approaches the cutoff radius. The comparisons of the accuracy these |
| 243 |
> |
quantities between the DSF kernel and SPME was surprisingly |
| 244 |
> |
good.\cite{Fennell:2006lq} As a result, the DSF method has seen |
| 245 |
> |
increasing use in molecular systems with relatively uniform charge |
| 246 |
|
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
| 247 |
|
|
| 248 |
|
\subsection{The damping function} |
| 257 |
|
produce complementary error functions when truncated at a finite |
| 258 |
|
distance. |
| 259 |
|
|
| 260 |
< |
With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method |
| 260 |
> |
With moderate damping coefficients, $\alpha \sim 0.2$ \AA$^{-1}$, the DSF method |
| 261 |
|
produced very good agreement with SPME for interaction energies, |
| 262 |
|
forces and torques for charge-charge |
| 263 |
|
interactions.\cite{Fennell:2006lq} |
| 343 |
|
\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance |
| 344 |
|
dot products: |
| 345 |
|
\begin{equation} |
| 346 |
< |
U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 346 |
> |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 347 |
|
\mathbf{D}_{a} \cdot |
| 348 |
|
\mathbf{D}_{b} \right) v_{21}(r) + |
| 349 |
|
\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right) |
| 372 |
|
which have been projected onto the surface of the cutoff sphere |
| 373 |
|
without changing their relative orientation, |
| 374 |
|
\begin{equation} |
| 375 |
< |
U_{D_{a}D_{b}}(r) = U_{D_{a}D_{b}}(r) - |
| 376 |
< |
U_{D_{a}D_{b}}(r_c) |
| 375 |
> |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r) = U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r) - |
| 376 |
> |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c) |
| 377 |
|
- (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot |
| 378 |
< |
\nabla U_{D_{a}D_{b}}(r_c). |
| 378 |
> |
\nabla U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c). |
| 379 |
|
\end{equation} |
| 380 |
|
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance |
| 381 |
|
(although the signs are reversed for the dipole that has been |
| 482 |
|
in the test cases are given in table~\ref{tab:pars}. |
| 483 |
|
|
| 484 |
|
\begin{table} |
| 485 |
< |
\caption{The parameters used in the systems used to evaluate the new |
| 486 |
< |
real-space methods. The most comprehensive test was a liquid |
| 487 |
< |
composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-} |
| 488 |
< |
ions). This test exercises all orders of the multipolar |
| 489 |
< |
interactions developed in the first paper.\label{tab:pars}} |
| 485 |
> |
\caption{The parameters used in the systems used to evaluate the new |
| 486 |
> |
real-space methods. The most comprehensive test was a liquid |
| 487 |
> |
composed of 2000 soft DQ liquid molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-} |
| 488 |
> |
ions). This test exercises all orders of the multipolar |
| 489 |
> |
interactions developed in the first paper.\label{tab:pars}} |
| 490 |
|
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
| 491 |
|
& \multicolumn{2}{c|}{LJ parameters} & |
| 492 |
|
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
| 499 |
|
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ |
| 500 |
|
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 501 |
|
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ |
| 502 |
< |
SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 502 |
> |
Soft DQ liquid & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 503 |
|
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
| 504 |
|
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
| 505 |
|
\end{tabularx} |
| 511 |
|
charges in addition to the multipolar fluid. The solid-phase |
| 512 |
|
parameters were chosen so that the systems can explore some |
| 513 |
|
orientational freedom for the multipolar sites, while maintaining |
| 514 |
< |
relatively strict translational order. The SSDQ model used here is |
| 515 |
< |
not a particularly accurate water model, but it does test |
| 516 |
< |
dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole |
| 517 |
< |
interactions at roughly the same magnitudes. The last test case, SSDQ |
| 518 |
< |
water with dissolved ions, exercises \textit{all} levels of the |
| 519 |
< |
multipole-multipole interactions we have derived so far and represents |
| 520 |
< |
the most complete test of the new methods. |
| 514 |
> |
relatively strict translational order. The soft DQ liquid model used |
| 515 |
> |
here based loosely on the SSDQO water |
| 516 |
> |
model,\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} but is not itself a |
| 517 |
> |
particularly accurate water model. However, the soft DQ model does |
| 518 |
> |
test dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole |
| 519 |
> |
interactions at roughly the same magnitudes. The last test case, a |
| 520 |
> |
soft DQ liquid with dissolved ions, exercises \textit{all} levels of |
| 521 |
> |
the multipole-multipole interactions we have derived so far and |
| 522 |
> |
represents the most complete test of the new methods. |
| 523 |
|
|
| 524 |
|
In the following section, we present results for the total |
| 525 |
|
electrostatic energy, as well as the electrostatic contributions to |
| 597 |
|
simulations, each system was created with 2,048 randomly-oriented |
| 598 |
|
molecules. These were equilibrated at a temperature of 300K for 1 ns. |
| 599 |
|
Each system was then simulated for 1 ns in the microcanonical (NVE) |
| 600 |
< |
ensemble. We collected 250 different configurations at equal time |
| 601 |
< |
intervals. For the liquid system that included ionic species, we |
| 602 |
< |
converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24 |
| 603 |
< |
\ce{Cl-} ions and re-equilibrated. After equilibration, the system was |
| 604 |
< |
run under the same conditions for 1 ns. A total of 250 configurations |
| 605 |
< |
were collected. In the following comparisons of energies, forces, and |
| 606 |
< |
torques, the Lennard-Jones potentials were turned off and only the |
| 607 |
< |
purely electrostatic quantities were compared with the same values |
| 608 |
< |
obtained via the Ewald sum. |
| 600 |
> |
ensemble with the Dullweber, Leimkuhler, and McLachlan (DLM) |
| 601 |
> |
symplectic splitting integrator using 1 fs |
| 602 |
> |
timesteps.\cite{Dullweber1997} We collected 250 different |
| 603 |
> |
configurations at equal time intervals. For the liquid system that |
| 604 |
> |
included ionic species, we converted 48 randomly-distributed molecules |
| 605 |
> |
into 24 \ce{Na+} and 24 \ce{Cl-} ions and re-equilibrated. After |
| 606 |
> |
equilibration, the system was run under the same conditions for 1 |
| 607 |
> |
ns. A total of 250 configurations were collected. In the following |
| 608 |
> |
comparisons of energies, forces, and torques, the Lennard-Jones |
| 609 |
> |
potentials were turned off and only the purely electrostatic |
| 610 |
> |
quantities were compared with the same values obtained via the Ewald |
| 611 |
> |
sum. |
| 612 |
|
|
| 613 |
|
\subsection{Accuracy of Energy Differences, Forces and Torques} |
| 614 |
|
The pairwise summation techniques (outlined above) were evaluated for |
| 624 |
|
Since none of the real-space methods provide exact energy differences, |
| 625 |
|
we used least square regressions analysis for the six different |
| 626 |
|
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
| 627 |
< |
with the multipolar Ewald reference method. Unitary results for both |
| 628 |
< |
the correlation (slope) and correlation coefficient for these |
| 629 |
< |
regressions indicate perfect agreement between the real-space method |
| 630 |
< |
and the multipolar Ewald sum. |
| 627 |
> |
with the multipolar Ewald reference method. A result of unity for |
| 628 |
> |
both the correlation (slope) and coefficient of determination ($R^2$) |
| 629 |
> |
for these regressions would indicate perfect agreement between the |
| 630 |
> |
real-space method and the multipolar Ewald sum. |
| 631 |
|
|
| 632 |
|
Molecular systems were run long enough to explore independent |
| 633 |
|
configurations and 250 configurations were recorded for comparison. |
| 660 |
|
the forces obtained via the Ewald sum and the real-space methods were |
| 661 |
|
evaluated, |
| 662 |
|
\begin{equation} |
| 663 |
< |
\cos\theta_i = \frac{\vec{f}_i^\mathrm{~Ewald} \cdot |
| 664 |
< |
\vec{f}_i^\mathrm{~GSF}}{\left|\vec{f}_i^\mathrm{~Ewald}\right| \left|\vec{f}_i^\mathrm{~GSF}\right|} |
| 663 |
> |
\cos\theta_i = \frac{\mathbf{f}_i^\mathrm{~Ewald} \cdot |
| 664 |
> |
\mathbf{f}_i^\mathrm{~GSF}}{\left|\mathbf{f}_i^\mathrm{~Ewald}\right| \left|\mathbf{f}_i^\mathrm{~GSF}\right|} |
| 665 |
|
\end{equation} |
| 666 |
|
The total angular displacement of the vectors was calculated as, |
| 667 |
|
\begin{equation} |
| 684 |
|
|
| 685 |
|
\subsection{Energy conservation} |
| 686 |
|
To test conservation the energy for the methods, the mixed molecular |
| 687 |
< |
system of 2000 SSDQ water molecules with 24 \ce{Na+} and 24 \ce{Cl-} |
| 688 |
< |
ions was run for 1 ns in the microcanonical ensemble at an average |
| 689 |
< |
temperature of 300K. Each of the different electrostatic methods |
| 690 |
< |
(Ewald, Hard, SP, GSF, and TSF) was tested for a range of different |
| 691 |
< |
damping values. The molecular system was started with same initial |
| 692 |
< |
positions and velocities for all cutoff methods. The energy drift |
| 693 |
< |
($\delta E_1$) and standard deviation of the energy about the slope |
| 694 |
< |
($\delta E_0$) were evaluated from the total energy of the system as a |
| 695 |
< |
function of time. Although both measures are valuable at |
| 687 |
> |
system of 2000 soft DQ liquid molecules with 24 \ce{Na+} and 24 |
| 688 |
> |
\ce{Cl-} ions was run for 1 ns in the microcanonical ensemble at an |
| 689 |
> |
average temperature of 300K. Each of the different electrostatic |
| 690 |
> |
methods (Ewald, Hard, SP, GSF, and TSF) was tested for a range of |
| 691 |
> |
different damping values. The molecular system was started with same |
| 692 |
> |
initial positions and velocities for all cutoff methods. The energy |
| 693 |
> |
drift ($\delta E_1$) and standard deviation of the energy about the |
| 694 |
> |
slope ($\delta E_0$) were evaluated from the total energy of the |
| 695 |
> |
system as a function of time. Although both measures are valuable at |
| 696 |
|
investigating new methods for molecular dynamics, a useful interaction |
| 697 |
|
model must allow for long simulation times with minimal energy drift. |
| 698 |
|
|
| 701 |
|
|
| 702 |
|
\begin{figure} |
| 703 |
|
\centering |
| 704 |
< |
\includegraphics[width=0.85\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
| 704 |
> |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
| 705 |
|
\caption{Statistical analysis of the quality of configurational |
| 706 |
|
energy differences for the real-space electrostatic methods |
| 707 |
|
compared with the reference Ewald sum. Results with a value equal |
| 708 |
|
to 1 (dashed line) indicate $\Delta E$ values indistinguishable |
| 709 |
|
from those obtained using the multipolar Ewald sum. Different |
| 710 |
|
values of the cutoff radius are indicated with different symbols |
| 711 |
< |
(9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted |
| 711 |
> |
(9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
| 712 |
|
triangles).\label{fig:slopeCorr_energy}} |
| 713 |
|
\end{figure} |
| 714 |
|
|
| 715 |
< |
The combined correlation coefficient and slope for all six systems is |
| 716 |
< |
shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the methods |
| 717 |
< |
reproduce the Ewald configurational energy differences with remarkable |
| 718 |
< |
fidelity. Undamped hard cutoffs introduce a significant amount of |
| 719 |
< |
random scatter in the energy differences which is apparent in the |
| 720 |
< |
reduced value of the correlation coefficient for this method. This |
| 721 |
< |
can be easily understood as configurations which exhibit small |
| 722 |
< |
traversals of a few dipoles or quadrupoles out of the cutoff sphere |
| 723 |
< |
will see large energy jumps when hard cutoffs are used. The |
| 724 |
< |
orientations of the multipoles (particularly in the ordered crystals) |
| 725 |
< |
mean that these energy jumps can go in either direction, producing a |
| 726 |
< |
significant amount of random scatter, but no systematic error. |
| 715 |
> |
The combined coefficient of determination and slope for all six |
| 716 |
> |
systems is shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the |
| 717 |
> |
methods reproduce the Ewald configurational energy differences with |
| 718 |
> |
remarkable fidelity. Undamped hard cutoffs introduce a significant |
| 719 |
> |
amount of random scatter in the energy differences which is apparent |
| 720 |
> |
in the reduced value of $R^2$ for this method. This can be easily |
| 721 |
> |
understood as configurations which exhibit small traversals of a few |
| 722 |
> |
dipoles or quadrupoles out of the cutoff sphere will see large energy |
| 723 |
> |
jumps when hard cutoffs are used. The orientations of the multipoles |
| 724 |
> |
(particularly in the ordered crystals) mean that these energy jumps |
| 725 |
> |
can go in either direction, producing a significant amount of random |
| 726 |
> |
scatter, but no systematic error. |
| 727 |
|
|
| 728 |
|
The TSF method produces energy differences that are highly correlated |
| 729 |
|
with the Ewald results, but it also introduces a significant |
| 734 |
|
effect, particularly for the crystalline systems. |
| 735 |
|
|
| 736 |
|
Both the SP and GSF methods appear to reproduce the Ewald results with |
| 737 |
< |
excellent fidelity, particularly for moderate damping ($\alpha = |
| 738 |
< |
0.1-0.2$\AA$^{-1}$) and with a commonly-used cutoff value ($r_c = |
| 739 |
< |
12$\AA). With the exception of the undamped hard cutoff, and the TSF |
| 737 |
> |
excellent fidelity, particularly for moderate damping ($\alpha \approx |
| 738 |
> |
0.2$~\AA$^{-1}$) and with a commonly-used cutoff value ($r_c = |
| 739 |
> |
12$~\AA). With the exception of the undamped hard cutoff, and the TSF |
| 740 |
|
method with short cutoffs, all of the methods would be appropriate for |
| 741 |
|
use in Monte Carlo simulations. |
| 742 |
|
|
| 766 |
|
energy conservation issues, and this perturbation is evident in the |
| 767 |
|
statistics accumulated for the molecular forces. The GSF |
| 768 |
|
perturbations are minimal, particularly for moderate damping and |
| 769 |
< |
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
| 770 |
< |
reasonable agreement in the correlation coefficient but again the |
| 771 |
< |
systematic error in the forces is concerning if replication of Ewald |
| 767 |
< |
forces is desired. |
| 769 |
> |
commonly-used cutoff values ($r_c = 12$~\AA). The TSF method shows |
| 770 |
> |
reasonable agreement in $R^2$, but again the systematic error in the |
| 771 |
> |
forces is concerning if replication of Ewald forces is desired. |
| 772 |
|
|
| 773 |
|
It is important to note that the forces and torques from the SP and |
| 774 |
|
the Hard cutoffs are not identical. The SP method shifts each |
| 786 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 787 |
|
line) indicate force magnitude values indistinguishable from those |
| 788 |
|
obtained using the multipolar Ewald sum. Different values of the |
| 789 |
< |
cutoff radius are indicated with different symbols (9\AA\ = |
| 790 |
< |
circles, 12\AA\ = squares, and 15\AA\ = inverted |
| 789 |
> |
cutoff radius are indicated with different symbols (9~\AA\ = |
| 790 |
> |
circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
| 791 |
|
triangles).\label{fig:slopeCorr_force}} |
| 792 |
|
\end{figure} |
| 793 |
|
|
| 800 |
|
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 801 |
|
line) indicate force magnitude values indistinguishable from those |
| 802 |
|
obtained using the multipolar Ewald sum. Different values of the |
| 803 |
< |
cutoff radius are indicated with different symbols (9\AA\ = |
| 804 |
< |
circles, 12\AA\ = squares, and 15\AA\ = inverted |
| 803 |
> |
cutoff radius are indicated with different symbols (9~\AA\ = |
| 804 |
> |
circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
| 805 |
|
triangles).\label{fig:slopeCorr_torque}} |
| 806 |
|
\end{figure} |
| 807 |
|
|
| 816 |
|
reproduces the torques in quite good agreement with the Ewald sum. |
| 817 |
|
The other real-space methods can cause some deviations, but excellent |
| 818 |
|
agreement with the Ewald sum torques is recovered at moderate values |
| 819 |
< |
of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1}$) and cutoff |
| 820 |
< |
radius ($r_c \ge 12$\AA). The TSF method exhibits only fair agreement |
| 819 |
> |
of the damping coefficient ($\alpha \approx 0.2$~\AA$^{-1}$) and cutoff |
| 820 |
> |
radius ($r_c \ge 12$~\AA). The TSF method exhibits only fair agreement |
| 821 |
|
in the slope when compared with the Ewald torques even for larger |
| 822 |
|
cutoff radii. It appears that the severity of the perturbations in |
| 823 |
|
the TSF method are most in evidence for the torques. |
| 829 |
|
these quantities. Force and torque vectors for all six systems were |
| 830 |
|
analyzed using Fisher statistics, and the quality of the vector |
| 831 |
|
directionality is shown in terms of circular variance |
| 832 |
< |
($\mathrm{Var}(\theta)$) in figure |
| 833 |
< |
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
| 834 |
< |
from the new real-space methods exhibit nearly-ideal Fisher probability |
| 835 |
< |
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
| 836 |
< |
exhibit the best vectorial agreement with the Ewald sum. The force and |
| 837 |
< |
torque vectors from GSF method also show good agreement with the Ewald |
| 838 |
< |
method, which can also be systematically improved by using moderate |
| 839 |
< |
damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c = |
| 840 |
< |
12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds |
| 841 |
< |
to a distribution with 95\% of force vectors within $6.37^\circ$ of |
| 842 |
< |
the corresponding Ewald forces. The TSF method produces the poorest |
| 843 |
< |
agreement with the Ewald force directions. |
| 832 |
> |
($\mathrm{Var}(\theta)$) in |
| 833 |
> |
Fig. \ref{fig:slopeCorr_circularVariance}. The force and torque |
| 834 |
> |
vectors from the new real-space methods exhibit nearly-ideal Fisher |
| 835 |
> |
probability distributions (Eq.~\ref{eq:pdf}). Both the hard and SP |
| 836 |
> |
cutoff methods exhibit the best vectorial agreement with the Ewald |
| 837 |
> |
sum. The force and torque vectors from GSF method also show good |
| 838 |
> |
agreement with the Ewald method, which can also be systematically |
| 839 |
> |
improved by using moderate damping and a reasonable cutoff radius. For |
| 840 |
> |
$\alpha = 0.2$~\AA$^{-1}$ and $r_c = 12$~\AA, we observe |
| 841 |
> |
$\mathrm{Var}(\theta) = 0.00206$, which corresponds to a distribution |
| 842 |
> |
with 95\% of force vectors within $6.37^\circ$ of the corresponding |
| 843 |
> |
Ewald forces. The TSF method produces the poorest agreement with the |
| 844 |
> |
Ewald force directions. |
| 845 |
|
|
| 846 |
|
Torques are again more perturbed than the forces by the new real-space |
| 847 |
|
methods, but even here the variance is reasonably small. For the same |
| 848 |
< |
method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA), |
| 849 |
< |
the circular variance was 0.01415, corresponds to a distribution which |
| 850 |
< |
has 95\% of torque vectors are within $16.75^\circ$ of the Ewald |
| 851 |
< |
results. Again, the direction of the force and torque vectors can be |
| 852 |
< |
systematically improved by varying $\alpha$ and $r_c$. |
| 848 |
> |
method (GSF) with the same parameters ($\alpha = 0.2$~\AA$^{-1}$, $r_c |
| 849 |
> |
= 12$~\AA), the circular variance was 0.01415, corresponds to a |
| 850 |
> |
distribution which has 95\% of torque vectors are within $16.75^\circ$ |
| 851 |
> |
of the Ewald results. Again, the direction of the force and torque |
| 852 |
> |
vectors can be systematically improved by varying $\alpha$ and $r_c$. |
| 853 |
|
|
| 854 |
|
\begin{figure} |
| 855 |
|
\centering |
| 860 |
|
indicates direction of the force or torque vectors are |
| 861 |
|
indistinguishable from those obtained from the Ewald sum. Here |
| 862 |
|
different symbols represent different values of the cutoff radius |
| 863 |
< |
(9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}} |
| 863 |
> |
(9~\AA\ = circle, 12~\AA\ = square, 15~\AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}} |
| 864 |
|
\end{figure} |
| 865 |
|
|
| 866 |
|
\subsection{Energy conservation\label{sec:conservation}} |
| 867 |
|
|
| 868 |
|
We have tested the conservation of energy one can expect to see with |
| 869 |
< |
the new real-space methods using the SSDQ water model with a small |
| 869 |
> |
the new real-space methods using the soft DQ liquid model with a small |
| 870 |
|
fraction of solvated ions. This is a test system which exercises all |
| 871 |
|
orders of multipole-multipole interactions derived in the first paper |
| 872 |
|
in this series and provides the most comprehensive test of the new |
| 873 |
< |
methods. A liquid-phase system was created with 2000 water molecules |
| 874 |
< |
and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
| 873 |
> |
methods. A liquid-phase system was created with 2000 liquid-phase |
| 874 |
> |
molecules and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
| 875 |
|
temperature of 300K. After equilibration in the canonical (NVT) |
| 876 |
|
ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system |
| 877 |
|
was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald, |
| 878 |
< |
Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA. The |
| 878 |
> |
Hard, SP, GSF, and TSF methods with a cutoff radius of 12~\AA. The |
| 879 |
|
value of the damping coefficient was also varied from the undamped |
| 880 |
< |
case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$ |
| 881 |
< |
\AA$^{-1}$) for all of the real space methods. A sample was also run |
| 882 |
< |
using the multipolar Ewald sum with the same real-space cutoff. |
| 880 |
> |
case ($\alpha = 0$) to a heavily damped case ($\alpha = |
| 881 |
> |
0.3$~\AA$^{-1}$) for all of the real space methods. A sample was also |
| 882 |
> |
run using the multipolar Ewald sum with the same real-space cutoff. |
| 883 |
|
|
| 884 |
|
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
| 885 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
| 901 |
|
\begin{figure} |
| 902 |
|
\centering |
| 903 |
|
\includegraphics[width=\textwidth]{newDrift_12.eps} |
| 904 |
< |
\caption{Energy conservation of the real-space methods for the SSDQ |
| 905 |
< |
water/ion system. $\delta \mathrm{E}_1$ is the linear drift in |
| 906 |
< |
energy over time (in kcal/mol/particle/ns) and $\delta |
| 904 |
> |
\caption{Energy conservation of the real-space methods for the soft |
| 905 |
> |
DQ liauid / ion system. $\delta \mathrm{E}_1$ is the linear drift |
| 906 |
> |
in energy over time (in kcal/mol/particle/ns) and $\delta |
| 907 |
|
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 908 |
|
around this drift (in kcal/mol/particle). Points that appear in |
| 909 |
|
the green region at the bottom exhibit better energy conservation |
| 923 |
|
of the local liquid ordering, one would not expect to see many |
| 924 |
|
differences in $g(r)$. Indeed, the pair distributions are essentially |
| 925 |
|
identical for all of the electrostatic methods studied (for each of |
| 926 |
< |
the different systems under investigation). An example of this |
| 922 |
< |
agreement for the SSDQ water/ion system is shown in |
| 923 |
< |
Fig. \ref{fig:gofr}. |
| 926 |
> |
the different systems under investigation). |
| 927 |
|
|
| 928 |
< |
\begin{figure} |
| 929 |
< |
\centering |
| 927 |
< |
\includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
| 928 |
< |
\caption{The pair distribution functions, $g(r)$, for the SSDQ |
| 929 |
< |
water/ion system obtained using the different real-space methods are |
| 930 |
< |
essentially identical with the result from the Ewald |
| 931 |
< |
treatment.\label{fig:gofr}} |
| 932 |
< |
\end{figure} |
| 928 |
> |
% An example of this agreement for the soft DQ liquid/ion system is |
| 929 |
> |
% shown in Fig. \ref{fig:gofr}. |
| 930 |
|
|
| 931 |
+ |
% \begin{figure} |
| 932 |
+ |
% \centering |
| 933 |
+ |
% \includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
| 934 |
+ |
% \caption{The pair distribution functions, $g(r)$, for the SSDQ |
| 935 |
+ |
% water/ion system obtained using the different real-space methods are |
| 936 |
+ |
% essentially identical with the result from the Ewald |
| 937 |
+ |
% treatment.\label{fig:gofr}} |
| 938 |
+ |
% \end{figure} |
| 939 |
+ |
|
| 940 |
|
There is a minor over-structuring of the first solvation shell when |
| 941 |
|
using TSF or when overdamping with any of the real-space methods. |
| 942 |
|
With moderate damping, GSF and SP produce pair distributions that are |
| 948 |
|
\end{equation} |
| 949 |
|
where $\rho$ is the number density of the site-site pair interactions, |
| 950 |
|
and $a$ is the radial location of the minima following the first peak |
| 951 |
< |
in $g(r)$ ($a = 4.2$ \AA\ for the SSDQ water/ion system). The |
| 951 |
> |
in $g(r)$ ($a = 4.2$~\AA\ for the soft DQ liquid / ion system). The |
| 952 |
|
coordination number is shown as a function of the damping coefficient |
| 953 |
|
for all of the real space methods in Fig. \ref{fig:Props}. |
| 954 |
|
|
| 956 |
|
of the electrostatic energy $\langle U_\mathrm{elect} \rangle / N$ |
| 957 |
|
which is obtained by sampling the liquid-state configurations |
| 958 |
|
experienced by a liquid evolving entirely under the influence of each |
| 959 |
< |
of the methods. In fig \ref{fig:Props} we demonstrate how $\langle |
| 959 |
> |
of the methods. In Fig. \ref{fig:Props} we demonstrate how $\langle |
| 960 |
|
U_\mathrm{elect} \rangle / N$ varies with the damping parameter, |
| 961 |
|
$\alpha$, for each of the methods. |
| 962 |
|
|
| 963 |
|
As in the crystals studied in the first paper, damping is important |
| 964 |
|
for converging the mean electrostatic energy values, particularly for |
| 965 |
|
the two shifted force methods (GSF and TSF). A value of $\alpha |
| 966 |
< |
\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF |
| 966 |
> |
\approx 0.2$~\AA$^{-1}$ is sufficient to converge the SP and GSF |
| 967 |
|
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
| 968 |
< |
dramatic damping ($\alpha \approx 0.28$ \AA$^{-1}$ for $r_c = 9$ \AA). |
| 968 |
> |
dramatic damping ($\alpha \approx 0.28$~\AA$^{-1}$ for $r_c = 9$~\AA). |
| 969 |
|
Overdamping the real-space electrostatic methods occurs with $\alpha > |
| 970 |
< |
0.3$, causing the estimate of the electrostatic energy to drop below |
| 971 |
< |
the Ewald results. |
| 970 |
> |
0.3$~\AA$^{-1}$, causing the estimate of the electrostatic energy to |
| 971 |
> |
drop below the Ewald results. |
| 972 |
|
|
| 973 |
|
These ``optimal'' values of the damping coefficient for structural |
| 974 |
|
features are similar to those observed for DSF electrostatics for |
| 975 |
|
purely point-charge systems, and the range $\alpha= 0.175 \rightarrow |
| 976 |
< |
0.225$ \AA$^{-1}$ for $r_c = 12$\AA\ appears to be an excellent |
| 976 |
> |
0.225$~\AA$^{-1}$ for $r_c = 12$~\AA\ appears to be an excellent |
| 977 |
|
compromise for mixed charge/multipolar systems. |
| 978 |
|
|
| 979 |
|
To test the fidelity of the electrostatic methods at reproducing |
| 987 |
|
which measures long-time behavior and is sensitive to the forces on |
| 988 |
|
the multipoles. The self-diffusion constants (D) were calculated from |
| 989 |
|
linear fits to the long-time portion of the mean square displacement, |
| 990 |
< |
$\langle r^{2}(t) \rangle$.\cite{Allen87} In fig. \ref{fig:Props} we |
| 990 |
> |
$\langle r^{2}(t) \rangle$.\cite{Allen87} In Fig. \ref{fig:Props} we |
| 991 |
|
demonstrate how the diffusion constant depends on the choice of |
| 992 |
|
real-space methods and the damping coefficient. Both the SP and GSF |
| 993 |
|
methods can obtain excellent agreement with Ewald again using moderate |
| 1006 |
|
translational diffusion. Here, $P_l$ is the Legendre polynomial of |
| 1007 |
|
order $l$ and $\hat{\mathbf{A}}_\gamma$ is the unit vector for body |
| 1008 |
|
axis $\gamma$. The reference frame used for our sample dipolar |
| 1009 |
< |
systems has the $z$-axis running along the dipoles, and for the SSDQ |
| 1010 |
< |
water model, the $y$-axis connects the two implied hydrogen atom |
| 1009 |
> |
systems has the $z$-axis running along the dipoles, and for the soft |
| 1010 |
> |
DQ liquid model, the $y$-axis connects the two implied hydrogen-like |
| 1011 |
|
positions. From the orientation autocorrelation functions, we can |
| 1012 |
|
obtain time constants for rotational relaxation either by fitting to a |
| 1013 |
|
multi-exponential model for the orientational relaxation, or by |
| 1014 |
|
integrating the correlation functions. |
| 1015 |
|
|
| 1016 |
< |
In a good water model, the orientational decay times would be |
| 1016 |
> |
In a good model for water, the orientational decay times would be |
| 1017 |
|
comparable to water orientational relaxation times from nuclear |
| 1018 |
|
magnetic resonance (NMR). The relaxation constant obtained from |
| 1019 |
|
$C_2^y(t)$ is normally of experimental interest because it describes |
| 1021 |
|
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
| 1022 |
|
of the dipole-dipole relaxation from a proton NMR signal and can |
| 1023 |
|
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
| 1024 |
< |
In fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$ |
| 1024 |
> |
In Fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$ |
| 1025 |
|
values for the various real-space methods over a range of different |
| 1026 |
|
damping coefficients. The rotational relaxation for the $z$ axis |
| 1027 |
|
primarily probes the torques on the dipoles, while the relaxation for |
| 1028 |
|
the $y$ axis is sensitive primarily to the quadrupolar torques. |
| 1029 |
|
|
| 1030 |
|
\begin{figure} |
| 1031 |
+ |
\includegraphics[width=\textwidth]{properties.eps} |
| 1032 |
|
\caption{Comparison of the structural and dynamic properties for the |
| 1033 |
< |
combined multipolar liquid (SSDQ water + ions) for all of the |
| 1034 |
< |
real-space methods with $r_c = 12$\AA. Electrostatic energies, |
| 1033 |
> |
combined multipolar liquid (soft DQ liquid + ions) for all of the |
| 1034 |
> |
real-space methods with $r_c = 12$~\AA. Electrostatic energies, |
| 1035 |
|
$\langle U_\mathrm{elect} \rangle / N$ (in kcal / mol), |
| 1036 |
|
coordination numbers, $n_C$, diffusion constants (in $10^{-5} |
| 1037 |
|
\mathrm{cm}^2\mathrm{s}^{-1}$), and rotational correlation times |
| 1038 |
|
(in ps) all show excellent agreement with Ewald results for |
| 1039 |
|
damping coefficients in the range $\alpha= 0.175 \rightarrow |
| 1040 |
< |
0.225$ \AA$^{-1}$. \label{fig:Props}} |
| 1034 |
< |
\includegraphics[width=\textwidth]{properties.eps} |
| 1040 |
> |
0.225$~\AA$^{-1}$. \label{fig:Props}} |
| 1041 |
|
\end{figure} |
| 1042 |
|
|
| 1043 |
|
In Fig. \ref{fig:Props} it appears that values for $D$, $\tau_2^y$, |
| 1052 |
|
the SP and GSF methods, as their structural and dynamical properties |
| 1053 |
|
still reproduce the Ewald results even in the completely undamped |
| 1054 |
|
($\alpha = 0$) case. An optimal range for the electrostatic damping |
| 1055 |
< |
parameter appears to be $\alpha= 0.175 \rightarrow 0.225$ \AA$^{-1}$ |
| 1056 |
< |
for $r_c = 12$\AA, which similar to the optimal range found for the |
| 1055 |
> |
parameter appears to be $\alpha= 0.175 \rightarrow 0.225$~\AA$^{-1}$ |
| 1056 |
> |
for $r_c = 12$~\AA, which similar to the optimal range found for the |
| 1057 |
|
damped shifted force potential for point charges.\cite{Fennell:2006lq} |
| 1058 |
|
|
| 1059 |
|
\section{CONCLUSION} |
| 1069 |
|
|
| 1070 |
|
We also developed two natural extensions of the damped shifted-force |
| 1071 |
|
(DSF) model originally proposed by Zahn {\it et al.} and extended by |
| 1072 |
< |
Fennel and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF |
| 1072 |
> |
Fennell and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF |
| 1073 |
|
approaches provide smooth truncation of energies, forces, and torques |
| 1074 |
|
at the real-space cutoff, and both converge to DSF electrostatics for |
| 1075 |
|
point-charge interactions. The TSF model is based on a high-order |
| 1091 |
|
series of the electrostatic kernel at the cutoff radius. The TSF |
| 1092 |
|
method also has the unique property that a large number of derivatives |
| 1093 |
|
can be made to vanish at the cutoff radius. This property has proven |
| 1094 |
< |
useful in past treatments of the corrections to the fluctuation |
| 1095 |
< |
formula for dielectric constants.\cite{Izvekov:2008wo} |
| 1094 |
> |
useful in past treatments of the corrections to the Clausius-Mossotti |
| 1095 |
> |
fluctuation formula for dielectric constants.\cite{Izvekov:2008wo} |
| 1096 |
|
|
| 1097 |
|
Reproduction of both structural and dynamical features in the liquid |
| 1098 |
|
systems is remarkably good for both the SP and GSF models. Pair |
