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 |