83 |
|
energy conservation properties that make it an excellent choice for |
84 |
|
efficient computation of electrostatic interactions in molecular |
85 |
|
dynamics simulations. Both SP and GSF are able to reproduce |
86 |
< |
structural and dyanamical properties in the liquid models with |
86 |
> |
structural and dynamical properties in the liquid models with |
87 |
|
excellent fidelity. |
88 |
|
\end{abstract} |
89 |
|
|
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 |
172 |
|
were also proposed by Chen \textit{et |
173 |
|
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
174 |
< |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
174 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfully |
175 |
|
used additional neutralization of higher order moments for systems of |
176 |
|
point charges.\cite{Fukuda:2013sf} |
177 |
|
|
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 aproaches 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 excercises 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. |
646 |
|
simulations. Because the real space methods reweight the different |
647 |
|
orientational contributions to the energies, it is also important to |
648 |
|
understand how the methods impact the \textit{directionality} of the |
649 |
< |
force and torque vectors. Fisher developed a probablity density |
649 |
> |
force and torque vectors. Fisher developed a probability density |
650 |
|
function to analyse directional data sets, |
651 |
|
\begin{equation} |
652 |
|
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta e^{\kappa \cos\theta} |
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 |
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. |
876 |
> |
ensemble using a Nos\'e-Hoover thermostat, six |
877 |
> |
statistically-independent replicas of this liquid-phase system were |
878 |
> |
run in the microcanonical (NVE) ensemble under the Ewald, Hard, SP, |
879 |
> |
GSF, and TSF methods with a cutoff radius of 12~\AA. The value of the |
880 |
> |
damping coefficient was also varied from the undamped case ($\alpha = |
881 |
> |
0$) to a heavily damped case ($\alpha = 0.3$~\AA$^{-1}$) for all of |
882 |
> |
the real space methods. A sample was also run using the multipolar |
883 |
> |
Ewald sum with the same real-space cutoff. |
884 |
|
|
885 |
|
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
886 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
901 |
|
|
902 |
|
\begin{figure} |
903 |
|
\centering |
904 |
< |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
905 |
< |
\caption{Energy conservation of the real-space methods for the SSDQ |
906 |
< |
water/ion system. $\delta \mathrm{E}_1$ is the linear drift in |
907 |
< |
energy over time (in kcal/mol/particle/ns) and $\delta |
904 |
> |
\includegraphics[width=\textwidth]{finalDrift.eps} |
905 |
> |
\caption{Energy conservation of the real-space methods for the soft |
906 |
> |
DQ liquid / ion system. $\delta \mathrm{E}_1$ is the linear drift |
907 |
> |
in energy over time (in kcal/mol/particle/ns) and $\delta |
908 |
|
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
909 |
|
around this drift (in kcal/mol/particle). Points that appear in |
910 |
|
the green region at the bottom exhibit better energy conservation |
924 |
|
of the local liquid ordering, one would not expect to see many |
925 |
|
differences in $g(r)$. Indeed, the pair distributions are essentially |
926 |
|
identical for all of the electrostatic methods studied (for each of |
927 |
< |
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}. |
927 |
> |
the different systems under investigation). |
928 |
|
|
929 |
< |
\begin{figure} |
930 |
< |
\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} |
929 |
> |
% An example of this agreement for the soft DQ liquid/ion system is |
930 |
> |
% shown in Fig. \ref{fig:gofr}. |
931 |
|
|
932 |
< |
There is a minor overstructuring of the first solvation shell when |
932 |
> |
% \begin{figure} |
933 |
> |
% \centering |
934 |
> |
% \includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
935 |
> |
% \caption{The pair distribution functions, $g(r)$, for the SSDQ |
936 |
> |
% water/ion system obtained using the different real-space methods are |
937 |
> |
% essentially identical with the result from the Ewald |
938 |
> |
% treatment.\label{fig:gofr}} |
939 |
> |
% \end{figure} |
940 |
> |
|
941 |
> |
There is a minor over-structuring of the first solvation shell when |
942 |
|
using TSF or when overdamping with any of the real-space methods. |
943 |
|
With moderate damping, GSF and SP produce pair distributions that are |
944 |
|
identical (within numerical noise) to their Ewald counterparts. The |
945 |
< |
degree of overstructuring can be measured most easily using the |
945 |
> |
degree of over-structuring can be measured most easily using the |
946 |
|
coordination number, |
947 |
|
\begin{equation} |
948 |
|
n_C = 4\pi\rho \int_{0}^{a}r^2\text{g}(r)dr, |
949 |
|
\end{equation} |
950 |
|
where $\rho$ is the number density of the site-site pair interactions, |
951 |
< |
$a$ and is the radial location of the minima following the first peak |
952 |
< |
in $g(r)$ ($a = 4.2$ \AA for the SSDQ water/ion system). The |
951 |
> |
and $a$ is the radial location of the minima following the first peak |
952 |
> |
in $g(r)$ ($a = 4.2$~\AA\ for the soft DQ liquid / ion system). The |
953 |
|
coordination number is shown as a function of the damping coefficient |
954 |
< |
for all of the real space methods in Fig. \ref{fig:Props}. |
954 |
> |
for all of the real space methods in Fig. \ref{fig:Props}. |
955 |
|
|
956 |
|
A more demanding test of modified electrostatics is the average value |
957 |
|
of the electrostatic energy $\langle U_\mathrm{elect} \rangle / N$ |
958 |
|
which is obtained by sampling the liquid-state configurations |
959 |
|
experienced by a liquid evolving entirely under the influence of each |
960 |
< |
of the methods. In fig \ref{fig:Props} we demonstrate how $\langle |
960 |
> |
of the methods. In Fig. \ref{fig:Props} we demonstrate how $\langle |
961 |
|
U_\mathrm{elect} \rangle / N$ varies with the damping parameter, |
962 |
< |
$\alpha$, for each of the methods. |
962 |
> |
$\alpha$, for each of the methods. |
963 |
|
|
964 |
|
As in the crystals studied in the first paper, damping is important |
965 |
|
for converging the mean electrostatic energy values, particularly for |
966 |
|
the two shifted force methods (GSF and TSF). A value of $\alpha |
967 |
< |
\approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF |
967 |
> |
\approx 0.2$~\AA$^{-1}$ is sufficient to converge the SP and GSF |
968 |
|
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
969 |
< |
dramatic damping ($\alpha \approx 0.28$ \AA$^{-1}$ for $r_c = 9$ \AA). |
969 |
> |
dramatic damping ($\alpha \approx 0.28$~\AA$^{-1}$ for $r_c = 9$~\AA). |
970 |
|
Overdamping the real-space electrostatic methods occurs with $\alpha > |
971 |
< |
0.3$, causing the estimate of the electrostatic energy to drop below |
972 |
< |
the Ewald results. |
971 |
> |
0.3$~\AA$^{-1}$, causing the estimate of the electrostatic energy to |
972 |
> |
drop below the Ewald results. |
973 |
|
|
974 |
< |
These ``optimal'' values of the damping coefficient are slightly |
975 |
< |
larger than what were observed for DSF electrostatics for purely |
976 |
< |
point-charge systems, although the range $\alpha= 0.175 \rightarrow |
977 |
< |
0.225$ \AA$^{-1}$ for $r_c = 12$\AA\ appears to be an excellent |
974 |
> |
These ``optimal'' values of the damping coefficient for structural |
975 |
> |
features are similar to those observed for DSF electrostatics for |
976 |
> |
purely point-charge systems, and the range $\alpha= 0.175 \rightarrow |
977 |
> |
0.225$~\AA$^{-1}$ for $r_c = 12$~\AA\ appears to be an excellent |
978 |
|
compromise for mixed charge/multipolar systems. |
979 |
|
|
980 |
|
To test the fidelity of the electrostatic methods at reproducing |
986 |
|
\label{eq:diff} |
987 |
|
\end{equation} |
988 |
|
which measures long-time behavior and is sensitive to the forces on |
989 |
< |
the multipoles. The self-diffusion constants (D) were calculated from |
989 |
> |
the multipoles. The self-diffusion constants (D) were calculated from |
990 |
|
linear fits to the long-time portion of the mean square displacement, |
991 |
< |
$\langle r^{2}(t) \rangle$.\cite{Allen87} In fig. \ref{fig:Props} we |
991 |
> |
$\langle r^{2}(t) \rangle$.\cite{Allen87} In Fig. \ref{fig:Props} we |
992 |
|
demonstrate how the diffusion constant depends on the choice of |
993 |
|
real-space methods and the damping coefficient. Both the SP and GSF |
994 |
|
methods can obtain excellent agreement with Ewald again using moderate |
996 |
|
|
997 |
|
In addition to translational diffusion, orientational relaxation times |
998 |
|
were calculated for comparisons with the Ewald simulations and with |
999 |
< |
experiments. These values were determined from the same 1~ns |
1000 |
< |
microcanonical trajectories used for translational diffusion by |
994 |
< |
calculating the orientational time correlation function, |
999 |
> |
experiments. These values were determined by calculating the |
1000 |
> |
orientational time correlation function, |
1001 |
|
\begin{equation} |
1002 |
|
C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t) |
1003 |
|
\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle, |
1004 |
|
\label{eq:OrientCorr} |
1005 |
|
\end{equation} |
1006 |
< |
where $P_l$ is the Legendre polynomial of order $l$ and |
1007 |
< |
$\hat{\mathbf{A}}_\gamma$ is the unit vector for body axis $\gamma$. |
1008 |
< |
The reference frame used for our sample dipolar systems has the |
1009 |
< |
$z$-axis running along the dipoles, and for the SSDQ water model, the |
1010 |
< |
$y$-axis connects the two implied hydrogen atom positions. From the |
1011 |
< |
orientation autocorrelation functions, we can obtain time constants |
1012 |
< |
for rotational relaxation either by fitting an exponential function or |
1013 |
< |
by integrating the entire correlation function. In a good water |
1014 |
< |
model, these decay times would be comparable to water orientational |
1015 |
< |
relaxation times from nuclear magnetic resonance (NMR). The relaxation |
1010 |
< |
constant obtained from $C_2^y(t)$ is normally of experimental interest |
1011 |
< |
because it describes the relaxation of the principle axis connecting |
1012 |
< |
the hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the |
1013 |
< |
intermolecular portion of the dipole-dipole relaxation from a proton |
1014 |
< |
NMR signal and should provide an estimate of the NMR relaxation time |
1015 |
< |
constant.\cite{Impey82} |
1006 |
> |
from the same 350 ps microcanonical trajectories that were used for |
1007 |
> |
translational diffusion. Here, $P_l$ is the Legendre polynomial of |
1008 |
> |
order $l$ and $\hat{\mathbf{A}}_\gamma$ is the unit vector for body |
1009 |
> |
axis $\gamma$. The reference frame used for our sample dipolar |
1010 |
> |
systems has the $z$-axis running along the dipoles, and for the soft |
1011 |
> |
DQ liquid model, the $y$-axis connects the two implied hydrogen-like |
1012 |
> |
positions. From the orientation autocorrelation functions, we can |
1013 |
> |
obtain time constants for rotational relaxation either by fitting to a |
1014 |
> |
multi-exponential model for the orientational relaxation, or by |
1015 |
> |
integrating the correlation functions. |
1016 |
|
|
1017 |
< |
Results for the diffusion constants and orientational relaxation times |
1018 |
< |
are shown in figure \ref{fig:Props}. From this data, it is apparent |
1019 |
< |
that the values for both $D$ and $\tau_2$ using the Ewald sum are |
1020 |
< |
reproduced with reasonable fidelity by the GSF method. |
1017 |
> |
In a good model for water, the orientational decay times would be |
1018 |
> |
comparable to water orientational relaxation times from nuclear |
1019 |
> |
magnetic resonance (NMR). The relaxation constant obtained from |
1020 |
> |
$C_2^y(t)$ is normally of experimental interest because it describes |
1021 |
> |
the relaxation of the principle axis connecting the hydrogen |
1022 |
> |
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
1023 |
> |
of the dipole-dipole relaxation from a proton NMR signal and can |
1024 |
> |
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
1025 |
> |
In Fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$ |
1026 |
> |
values for the various real-space methods over a range of different |
1027 |
> |
damping coefficients. The rotational relaxation for the $z$ axis |
1028 |
> |
primarily probes the torques on the dipoles, while the relaxation for |
1029 |
> |
the $y$ axis is sensitive primarily to the quadrupolar torques. |
1030 |
|
|
1031 |
|
\begin{figure} |
1032 |
+ |
\includegraphics[width=\textwidth]{properties.eps} |
1033 |
|
\caption{Comparison of the structural and dynamic properties for the |
1034 |
< |
combined multipolar liquid (SSDQ water + ions) for all of the |
1035 |
< |
real-space methods with $r_c = 12$\AA. Electrostatic energies, |
1034 |
> |
combined multipolar liquid (soft DQ liquid + ions) for all of the |
1035 |
> |
real-space methods with $r_c = 12$~\AA. Electrostatic energies, |
1036 |
|
$\langle U_\mathrm{elect} \rangle / N$ (in kcal / mol), |
1037 |
< |
coordination numbers, $n_C$, diffusion constants (in cm$^2$ |
1038 |
< |
s$^{-1}$), and rotational correlation times (in fs) all show |
1039 |
< |
excellent agreement with Ewald results for damping coefficients in |
1040 |
< |
the range $\alpha= 0.175 \rightarrow 0.225$ |
1041 |
< |
\AA$^{-1}$. \label{fig:Props}} |
1032 |
< |
\includegraphics[width=\textwidth]{properties.eps} |
1037 |
> |
coordination numbers, $n_C$, diffusion constants (in $10^{-5} |
1038 |
> |
\mathrm{cm}^2\mathrm{s}^{-1}$), and rotational correlation times |
1039 |
> |
(in ps) all show excellent agreement with Ewald results for |
1040 |
> |
damping coefficients in the range $\alpha= 0.175 \rightarrow |
1041 |
> |
0.225$~\AA$^{-1}$. \label{fig:Props}} |
1042 |
|
\end{figure} |
1043 |
|
|
1044 |
+ |
In Fig. \ref{fig:Props} it appears that values for $D$, $\tau_2^y$, |
1045 |
+ |
and $\tau_2^z$ using the Ewald sum are reproduced with excellent |
1046 |
+ |
fidelity by the GSF and SP methods. All of the real space methods can |
1047 |
+ |
be \textit{overdamped}, which reduces the effective range of multipole |
1048 |
+ |
interactions, causing structural and dynamical changes from the |
1049 |
+ |
correct behavior. Because overdamping weakens orientational |
1050 |
+ |
preferences between adjacent molecules, it manifests as too-rapid |
1051 |
+ |
orientational decay coupled with faster diffusion and |
1052 |
+ |
over-coordination of the liquid. Underdamping is less problematic for |
1053 |
+ |
the SP and GSF methods, as their structural and dynamical properties |
1054 |
+ |
still reproduce the Ewald results even in the completely undamped |
1055 |
+ |
($\alpha = 0$) case. An optimal range for the electrostatic damping |
1056 |
+ |
parameter appears to be $\alpha= 0.175 \rightarrow 0.225$~\AA$^{-1}$ |
1057 |
+ |
for $r_c = 12$~\AA, which similar to the optimal range found for the |
1058 |
+ |
damped shifted force potential for point charges.\cite{Fennell:2006lq} |
1059 |
|
|
1060 |
|
\section{CONCLUSION} |
1061 |
|
In the first paper in this series, we generalized the |
1070 |
|
|
1071 |
|
We also developed two natural extensions of the damped shifted-force |
1072 |
|
(DSF) model originally proposed by Zahn {\it et al.} and extended by |
1073 |
< |
Fennel and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF |
1073 |
> |
Fennell and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF |
1074 |
|
approaches provide smooth truncation of energies, forces, and torques |
1075 |
|
at the real-space cutoff, and both converge to DSF electrostatics for |
1076 |
|
point-charge interactions. The TSF model is based on a high-order |
1080 |
|
sphere to derive shifted force and torque expressions, and is a |
1081 |
|
significantly more gentle approach. |
1082 |
|
|
1083 |
< |
The GSF method produced quantitative agreement with Ewald energy, |
1084 |
< |
force, and torques. It also performs well in conserving energy in MD |
1083 |
> |
The GSF method produces quantitative agreement with Ewald energies, |
1084 |
> |
forces, and torques. It also performs well in conserving energy in MD |
1085 |
|
simulations. The Taylor-shifted (TSF) model provides smooth dynamics, |
1086 |
|
but these take place on a potential energy surface that is |
1087 |
|
significantly perturbed from Ewald-based electrostatics. Because it |
1092 |
|
series of the electrostatic kernel at the cutoff radius. The TSF |
1093 |
|
method also has the unique property that a large number of derivatives |
1094 |
|
can be made to vanish at the cutoff radius. This property has proven |
1095 |
< |
useful in past treatments of the corrections to the fluctuation |
1096 |
< |
formula for dielectric constants.\cite{Izvekov:2008wo} |
1095 |
> |
useful in past treatments of the corrections to the Clausius-Mossotti |
1096 |
> |
fluctuation formula for dielectric constants.\cite{Izvekov:2008wo} |
1097 |
|
|
1098 |
|
Reproduction of both structural and dynamical features in the liquid |
1099 |
|
systems is remarkably good for both the SP and GSF models. Pair |
1119 |
|
Based on the results of this work, we can conclude that the GSF method |
1120 |
|
is a suitable and efficient replacement for the Ewald sum for |
1121 |
|
evaluating electrostatic interactions in modern MD simulations, and |
1122 |
< |
the SP meethod would be an excellent choice for Monte Carlo |
1122 |
> |
the SP method would be an excellent choice for Monte Carlo |
1123 |
|
simulations where smooth forces and energy conservation are not |
1124 |
|
important. Both the SP and GSF methods retain excellent fidelity to |
1125 |
|
the Ewald energies, forces and torques. Additionally, the energy |
1126 |
|
drift and fluctuations from the GSF electrostatics are significantly |
1127 |
< |
better than a multipolar Ewald sum for finite-sized reciprocal spaces. |
1127 |
> |
better than a multipolar Ewald sum for finite-sized reciprocal spaces, |
1128 |
> |
and physical properties are reproduced accurately. |
1129 |
|
|
1130 |
|
As in all purely pairwise cutoff methods, the SP, GSF and TSF methods |
1131 |
|
are expected to scale approximately {\it linearly} with system size, |