| 477 |
|
real-space cutoffs. In the first paper of this series, we compared |
| 478 |
|
the dipolar and quadrupolar energy expressions against analytic |
| 479 |
|
expressions for ordered dipolar and quadrupolar |
| 480 |
< |
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} This work uses the |
| 481 |
< |
multipolar Ewald sum as a reference method for comparing energies, |
| 482 |
< |
forces, and torques for molecular models that mimic disordered and |
| 483 |
< |
ordered condensed-phase systems. These test-cases include: |
| 484 |
< |
|
| 480 |
> |
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} In this work, we |
| 481 |
> |
used the multipolar Ewald sum as a reference method for comparing |
| 482 |
> |
energies, forces, and torques for molecular models that mimic |
| 483 |
> |
disordered and ordered condensed-phase systems. These test-cases |
| 484 |
> |
include: |
| 485 |
|
\begin{itemize} |
| 486 |
< |
\item Soft Dipolar fluids ($\sigma = , \epsilon = , |D| = $) |
| 487 |
< |
\item Soft Dipolar solids ($\sigma = , \epsilon = , |D| = $) |
| 488 |
< |
\item Soft Quadrupolar fluids ($\sigma = , \epsilon = , Q_{xx} = ...$) |
| 489 |
< |
\item Soft Quadrupolar solids ($\sigma = , \epsilon = , Q_{xx} = ...$) |
| 490 |
< |
\item A mixed multipole model for water |
| 491 |
< |
\item A mixed multipole models for water with dissolved ions |
| 486 |
> |
\item Soft Dipolar fluids ($\sigma = 3.051$, $\epsilon =0.152$, $|D| = 2.35$) |
| 487 |
> |
\item Soft Dipolar solids ($\sigma = 2.837$, $\epsilon =1.0$, $|D| = 2.35$) |
| 488 |
> |
\item Soft Quadrupolar fluids ($\sigma = 3.051$, $\epsilon =0.152$, $Q_{\alpha\alpha} =\left\{-1,-1,-2.5\right\}$) |
| 489 |
> |
\item Soft Quadrupolar solids ($\sigma = 2.837$, $\epsilon = 1.0$, $Q_{\alpha\alpha} =\left\{-1,-1,-2.5\right\}$) |
| 490 |
> |
\item A mixed multipole model (SSDQ) for water ($\sigma = 3.051$, $\epsilon = 0.152$, $D_z = 2.35$, $Q_{\alpha\alpha} =\left\{-1.35,0,-0.68\right\}$) |
| 491 |
> |
\item A mixed multipole models for water with 48 dissolved ions, 24 |
| 492 |
> |
\ce{Na+}: ($\sigma = 2.579$, $\epsilon =0.118$, $q = 1e$) and 24 |
| 493 |
> |
\ce{Cl-}: ($\sigma = 4.445$, $\epsilon =0.1$l, $q = -1e$) |
| 494 |
|
\end{itemize} |
| 495 |
< |
This last test case exercises all levels of the multipole-multipole |
| 496 |
< |
interactions we have derived so far and represents the most complete |
| 497 |
< |
test of the new methods. |
| 495 |
> |
All Lennard-Jones parameters are in units of \AA\ $(\sigma)$ and kcal |
| 496 |
> |
/ mole $(\epsilon)$. Partial charges are reported in electrons, while |
| 497 |
> |
dipoles are in Debye units, and quadrupoles are in units of Debye-\AA. |
| 498 |
|
|
| 499 |
< |
In the following section, we present results for the total |
| 500 |
< |
electrostatic energy, as well as the electrostatic contributions to |
| 501 |
< |
the force and torque on each molecule. These quantities have been |
| 502 |
< |
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
| 503 |
< |
and have been compared with the values obtaine from the multipolar |
| 504 |
< |
Ewald sum. In Mote Carlo (MC) simulations, the energy differences |
| 505 |
< |
between two configurations is the primary quantity that governs how |
| 506 |
< |
the simulation proceeds. These differences are the most imporant |
| 507 |
< |
indicators of the reliability of a method even if the absolute |
| 508 |
< |
energies are not exact. For each of the multipolar systems listed |
| 509 |
< |
above, we have compared the change in electrostatic potential energy |
| 510 |
< |
($\Delta E$) between 250 statistically-independent configurations. In |
| 511 |
< |
molecular dynamics (MD) simulations, the forces and torques govern the |
| 512 |
< |
behavior of the simulation, so we also compute the electrostatic |
| 513 |
< |
contributions to the forces and torques. |
| 499 |
> |
The last test case exercises all levels of the multipole-multipole |
| 500 |
> |
interactions we have derived so far and represents the most complete |
| 501 |
> |
test of the new methods. In the following section, we present results |
| 502 |
> |
for the total electrostatic energy, as well as the electrostatic |
| 503 |
> |
contributions to the force and torque on each molecule. These |
| 504 |
> |
quantities have been computed using the SP, TSF, and GSF methods, as |
| 505 |
> |
well as a hard cutoff, and have been compared with the values obtaine |
| 506 |
> |
from the multipolar Ewald sum. In Mote Carlo (MC) simulations, the |
| 507 |
> |
energy differences between two configurations is the primary quantity |
| 508 |
> |
that governs how the simulation proceeds. These differences are the |
| 509 |
> |
most imporant indicators of the reliability of a method even if the |
| 510 |
> |
absolute energies are not exact. For each of the multipolar systems |
| 511 |
> |
listed above, we have compared the change in electrostatic potential |
| 512 |
> |
energy ($\Delta E$) between 250 statistically-independent |
| 513 |
> |
configurations. In molecular dynamics (MD) simulations, the forces |
| 514 |
> |
and torques govern the behavior of the simulation, so we also compute |
| 515 |
> |
the electrostatic contributions to the forces and torques. |
| 516 |
|
|
| 517 |
|
\subsection{Model systems} |
| 518 |
|
To sample independent configurations of multipolar crystals, a body |
| 519 |
< |
centered cubic (BCC) crystal which is a minimum energy structure for |
| 519 |
> |
centered cubic (bcc) crystal which is a minimum energy structure for |
| 520 |
|
point dipoles was generated using 3,456 molecules. The multipoles |
| 521 |
|
were translationally locked in their respective crystal sites for |
| 522 |
|
equilibration at a relatively low temperature (50K), so that dipoles |
| 651 |
|
|
| 652 |
|
% \label{fig:barGraph2} |
| 653 |
|
% \end{figure} |
| 654 |
< |
%The correlation coefficient ($R^2$) and slope of the linear regression plots for the energy differences for all six different molecular systems is shown in figure 4a and 4b.The plot shows that the correlation coefficient improves for the SP cutoff method as compared to the undamped hard cutoff method in the case of SSDQC, SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar crystal and liquid, the correlation coefficient is almost unchanged and close to 1. The correlation coefficient is smallest (0.696276 for $r_c$ = 9 $A^o$) for the SSDQC liquid because of the presence of charge-charge and charge-multipole interactions. Since the charge-charge and charge-multipole interaction is long ranged, there is huge deviation of correlation coefficient from 1. Similarly, the quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with compared to interactions in the other multipolar systems, thus the correlation coefficient very close to 1 even for hard cutoff method. The idea of placing image multipole on the surface of the cutoff sphere improves the correlation coefficient and makes it close to 1 for all types of multipolar systems. Similarly the slope is hugely deviated from the correct value for the lower order multipole-multipole interaction and slightly deviated for higher order multipole – multipole interaction. The SP method improves both correlation coefficient ($R^2$) and slope significantly in SSDQC and dipolar systems. The Slope is found to be deviated more in dipolar crystal as compared to liquid which is associated with the large fluctuation in the electrostatic energy in crystal. The GSF also produced better values of correlation coefficient and slope with the proper selection of the damping alpha (Interested reader can consult accompanying supporting material). The TSF method gives good value of correlation coefficient for the dipolar crystal, dipolar liquid, SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the regression slopes are significantly deviated. |
| 654 |
> |
%The correlation coefficient ($R^2$) and slope of the linear |
| 655 |
> |
%regression plots for the energy differences for all six different |
| 656 |
> |
%molecular systems is shown in figure 4a and 4b.The plot shows that |
| 657 |
> |
%the correlation coefficient improves for the SP cutoff method as |
| 658 |
> |
%compared to the undamped hard cutoff method in the case of SSDQC, |
| 659 |
> |
%SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar |
| 660 |
> |
%crystal and liquid, the correlation coefficient is almost unchanged |
| 661 |
> |
%and close to 1. The correlation coefficient is smallest (0.696276 |
| 662 |
> |
%for $r_c$ = 9 $A^\circ$) for the SSDQC liquid because of the presence of |
| 663 |
> |
%charge-charge and charge-multipole interactions. Since the |
| 664 |
> |
%charge-charge and charge-multipole interaction is long ranged, there |
| 665 |
> |
%is huge deviation of correlation coefficient from 1. Similarly, the |
| 666 |
> |
%quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with |
| 667 |
> |
%compared to interactions in the other multipolar systems, thus the |
| 668 |
> |
%correlation coefficient very close to 1 even for hard cutoff |
| 669 |
> |
%method. The idea of placing image multipole on the surface of the |
| 670 |
> |
%cutoff sphere improves the correlation coefficient and makes it close |
| 671 |
> |
%to 1 for all types of multipolar systems. Similarly the slope is |
| 672 |
> |
%hugely deviated from the correct value for the lower order |
| 673 |
> |
%multipole-multipole interaction and slightly deviated for higher |
| 674 |
> |
%order multipole – multipole interaction. The SP method improves both |
| 675 |
> |
%correlation coefficient ($R^2$) and slope significantly in SSDQC and |
| 676 |
> |
%dipolar systems. The Slope is found to be deviated more in dipolar |
| 677 |
> |
%crystal as compared to liquid which is associated with the large |
| 678 |
> |
%fluctuation in the electrostatic energy in crystal. The GSF also |
| 679 |
> |
%produced better values of correlation coefficient and slope with the |
| 680 |
> |
%proper selection of the damping alpha (Interested reader can consult |
| 681 |
> |
%accompanying supporting material). The TSF method gives good value of |
| 682 |
> |
%correlation coefficient for the dipolar crystal, dipolar liquid, |
| 683 |
> |
%SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the |
| 684 |
> |
%regression slopes are significantly deviated. |
| 685 |
> |
|
| 686 |
|
\begin{figure} |
| 687 |
< |
\centering |
| 688 |
< |
\includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 689 |
< |
\caption{The correlation coefficient and regression slope of configurational energy differences for a given method with compared with the reference Ewald method. The value of result equal to 1(dashed line) indicates energy difference is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\ = circle, 12 \AA\ = square 15 \AA\ = inverted triangle)} |
| 690 |
< |
\label{fig:slopeCorr_energy} |
| 691 |
< |
\end{figure} |
| 692 |
< |
The combined correlation coefficient and slope for all six systems is shown in Figure ~\ref{fig:slopeCorr_energy}. The correlation coefficient for the undamped hard cutoff method is does not have good agreement with the Ewald because of the fluctuation of the electrostatic energy in the direct truncation method. This deviation in correlation coefficient is improved by using SP, GSF, and TSF method. But the TSF method worsens the regression slope stating that this method produces statistically more biased result as compared to Ewald. Also the GSF method slightly deviate slope but it can be alleviated by using proper value of damping alpha and cutoff radius. The SP method shows good agreement with Ewald method for all values of damping alpha and radii. |
| 687 |
> |
\centering |
| 688 |
> |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 689 |
> |
\caption{Statistical analysis of the quality of configurational |
| 690 |
> |
energy differences for the real-space electrostatic methods |
| 691 |
> |
compared with the reference Ewald sum. Results with a value equal |
| 692 |
> |
to 1 (dashed line) indicate $\Delta E$ values indistinguishable |
| 693 |
> |
from those obtained using the multipolar Ewald sum. Different |
| 694 |
> |
values of the cutoff radius are indicated with different symbols |
| 695 |
> |
(9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted |
| 696 |
> |
triangles).} |
| 697 |
> |
\label{fig:slopeCorr_energy} |
| 698 |
> |
\end{figure} |
| 699 |
> |
|
| 700 |
> |
The combined correlation coefficient and slope for all six systems is |
| 701 |
> |
shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the methods |
| 702 |
> |
reproduce the Ewald-derived configurational energy differences with |
| 703 |
> |
remarkable fidelity. Undamped hard cutoffs introduce a significant |
| 704 |
> |
amount of random scatter in the energy differences which is apparent |
| 705 |
> |
in the reduced value of the correlation coefficient for this method. |
| 706 |
> |
This can be understood easily as configurations which exhibit only |
| 707 |
> |
small traversals of a few dipoles or quadrupoles out of the cutoff |
| 708 |
> |
sphere will see large energy jumps when hard cutoffs are used. The |
| 709 |
> |
orientations of the multipoles (particularly in the ordered crystals) |
| 710 |
> |
mean that these jumps can go either up or down in energy, producing a |
| 711 |
> |
significant amount of random scatter. |
| 712 |
> |
|
| 713 |
> |
The TSF method produces energy differences that are highly correlated |
| 714 |
> |
with the Ewald results, but it also introduces a significant |
| 715 |
> |
systematic bias in the values of the energies, particularly for |
| 716 |
> |
smaller cutoff values. The TSF method alters the distance dependence |
| 717 |
> |
of different orientational contributions to the energy in a |
| 718 |
> |
non-uniform way, so the size of the cutoff sphere can have a large |
| 719 |
> |
effect on crystalline systems. |
| 720 |
> |
|
| 721 |
> |
Both the SP and GSF methods appear to reproduce the Ewald results with |
| 722 |
> |
excellent fidelity, particularly for moderate damping ($\alpha = |
| 723 |
> |
0.1-0.2$\AA$^{-1}$) and commonly-used cutoff values ($r_c = 12$\AA). |
| 724 |
> |
With the exception of the undamped hard cutoff, and the TSF method |
| 725 |
> |
with short cutoffs, all of the methods would be appropriate for use in |
| 726 |
> |
Monte Carlo simulations. |
| 727 |
> |
|
| 728 |
|
\subsection{Magnitude of the force and torque vectors} |
| 729 |
< |
The comparison of the magnitude of the combined forces and torques for the data accumulated from all system types are shown in Figure ~\ref{fig:slopeCorr_force}. The correlation and slope for the forces agree with the Ewald even for the hard cutoff method. For the system of molecules with higher order multipoles, the interaction is short ranged. Moreover, the force decays more rapidly than the electrostatic energy hence the hard cutoff method also produces good results. Although the pure cutoff gives the good match of the electrostatic force, the discontinuity in the force at the cutoff radius causes problem in the total energy conservation in MD simulations, which will be discussed in detail in subsection D. The correlation coefficient for GSF method also perfectly matches with Ewald but the slope is slightly deviated (due to extra term obtained from the angular differentiation). This deviation in the slope can be alleviated with proper selection of the damping alpha and radii ($\alpha = 0.2$ and $r_c = 12 A^o$ are good choice). The TSF method shows good agreement in the correlation coefficient but the slope is not good as compared to the Ewald. |
| 729 |
> |
|
| 730 |
> |
The comparison of the magnitude of the combined forces and torques for |
| 731 |
> |
the data accumulated from all system types are shown in Figures |
| 732 |
> |
~\ref{fig:slopeCorr_force} and \ref{fig:slopeCorr_torque}. The |
| 733 |
> |
correlation and slope for the forces agree well with the Ewald sum |
| 734 |
> |
even for the hard cutoff method. |
| 735 |
> |
|
| 736 |
> |
For the system of molecules with higher order multipoles, the |
| 737 |
> |
interaction is quite short ranged. Moreover, the force decays more |
| 738 |
> |
rapidly than the electrostatic energy hence the hard cutoff method can |
| 739 |
> |
also produces reasonable agreement. Although the pure cutoff gives |
| 740 |
> |
the good match of the electrostatic force for pairs of molecules |
| 741 |
> |
included within the cutoff sphere, the discontinuity in the force at |
| 742 |
> |
the cutoff radius can potentially cause problems the total energy |
| 743 |
> |
conservation as molecules enter and leave the cutoff sphere. This is |
| 744 |
> |
discussed in detail in section \ref{sec:}. |
| 745 |
> |
|
| 746 |
> |
The two shifted-force methods (GSF and TSF) exhibit a small amount of |
| 747 |
> |
systematic variation and scatter compared with the Ewald forces. The |
| 748 |
> |
shifted-force models intentionally perturb the forces between pairs of |
| 749 |
> |
molecules inside the cutoff sphere in order to correct the energy |
| 750 |
> |
conservation issues, so it is not particularly surprising that this |
| 751 |
> |
perturbation is evident in these same molecular forces. The GSF |
| 752 |
> |
perturbations are minimal, particularly for moderate damping and and |
| 753 |
> |
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
| 754 |
> |
reasonable agreement in the correlation coefficient but again the |
| 755 |
> |
systematic error in the forces is concerning if replication of Ewald |
| 756 |
> |
forces is desired. |
| 757 |
> |
|
| 758 |
|
\begin{figure} |
| 759 |
< |
\centering |
| 760 |
< |
\includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 761 |
< |
\caption{The correlation coefficient and regression slope of the magnitude of the force for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 \AA\ = circle, 12 \AA\ = square 15 \AA\ = inverted triangle). } |
| 762 |
< |
\label{fig:slopeCorr_force} |
| 763 |
< |
\end{figure} |
| 764 |
< |
The torques appears to be very influenced because of extra term generated when the potential energy is modified to get consistent force and torque. The result shows that the torque from the hard cutoff method has good agreement with Ewald. As the potential is modified to make it consistent with the force and torque, the correlation and slope is deviated as shown in Figure~\ref{fig:slopeCorr_torque} for SP, GSF and TSF cutoff methods. But the proper value of the damping alpha and radius can improve the agreement of the GSF with the Ewald method. The TSF method shows worst agreement in the slope as compared to Ewald even for larger cutoff radii. |
| 759 |
> |
\centering |
| 760 |
> |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 761 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 762 |
> |
magnitudes for the real-space electrostatic methods compared with |
| 763 |
> |
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 764 |
> |
line) indicate force magnitude values indistinguishable from those |
| 765 |
> |
obtained using the multipolar Ewald sum. Different values of the |
| 766 |
> |
cutoff radius are indicated with different symbols (9\AA\ = |
| 767 |
> |
circles, 12\AA\ = squares, and 15\AA\ = inverted triangles). } |
| 768 |
> |
\label{fig:slopeCorr_force} |
| 769 |
> |
\end{figure} |
| 770 |
> |
|
| 771 |
> |
|
| 772 |
|
\begin{figure} |
| 773 |
< |
\centering |
| 774 |
< |
\includegraphics[width=0.5 \textwidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 775 |
< |
\caption{The correlation coefficient and regression slope of the magnitude of the torque for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle).} |
| 776 |
< |
\label{fig:slopeCorr_torque} |
| 777 |
< |
\end{figure} |
| 778 |
< |
\subsection{Directionality of the force and torque vectors} |
| 779 |
< |
The accurate evaluation of the direction of the force and torques are also important for the dynamic simulation.In our research, the direction data sets were computed from the purposed method and compared with Ewald using Fisher statistics and results are expressed in terms of circular variance ($Var(\theta$).The force and torque vectors from the purposed method followed Fisher probability distribution function expressed in equation~\ref{eq:pdf}. The circular variance for the force and torque vectors of each molecule in the 250 configurations for all system types is shown in Figure~\ref{fig:slopeCorr_circularVariance}. The direction of the force and torque vectors from hard and SP cutoff methods showed best directional agreement with the Ewald. The force and torque vectors from GSF method also showed good agreement with the Ewald method, which can also be improved by varying damping alpha and cutoff radius.For $\alpha = 0.2$ and $r_c = 12 A^o$, $ Var(\theta) $ for direction of the force was found to be 0.002061 and corresponding value of $\kappa $ was 485.20. Integration of equation ~\ref{eq:pdf} for that corresponding value of $\kappa$ showed that 95\% of force vectors are with in $6.37^o$. The TSF method is the poorest in evaluating accurate direction with compared to Hard, SP, and GSF methods. The circular variance for the direction of the torques is larger as compared to force. For same $\alpha = 0.2, r_c = 12 A^o$ and GSF method, the circular variance was 0.01415, which showed 95\% of torque vectors are within $16.75^o$.The direction of the force and torque vectors can be improved by varying $\alpha$ and $r_c$. |
| 773 |
> |
\centering |
| 774 |
> |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 775 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 776 |
> |
magnitudes for the real-space electrostatic methods compared with |
| 777 |
> |
the reference Ewald sum. Results with a value equal to 1 (dashed |
| 778 |
> |
line) indicate force magnitude values indistinguishable from those |
| 779 |
> |
obtained using the multipolar Ewald sum. Different values of the |
| 780 |
> |
cutoff radius are indicated with different symbols (9\AA\ = |
| 781 |
> |
circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 782 |
> |
\label{fig:slopeCorr_torque} |
| 783 |
> |
\end{figure} |
| 784 |
|
|
| 785 |
+ |
The torques (Fig. \ref{fig:slopeCorr_torque}) appear to be |
| 786 |
+ |
significantly influenced by the choice of real-space method. The |
| 787 |
+ |
torque expressions have the same distance dependence as the energies, |
| 788 |
+ |
which are naturally longer-ranged expressions than the inter-site |
| 789 |
+ |
forces. Torques are also quite sensitive to orientations of |
| 790 |
+ |
neighboring molecules, even those that are near the cutoff distance. |
| 791 |
+ |
|
| 792 |
+ |
The results shows that the torque from the hard cutoff method |
| 793 |
+ |
reproduces the torques in quite good agreement with the Ewald sum. |
| 794 |
+ |
The other real-space methods can cause some significant deviations, |
| 795 |
+ |
but excellent agreement with the Ewald sum torques is recovered at |
| 796 |
+ |
moderate values of the damping coefficient ($\alpha = |
| 797 |
+ |
0.1-0.2$\AA$^{-1}$) and cutoff radius ($r_c \ge 12$\AA). The TSF |
| 798 |
+ |
method exhibits the only fair agreement in the slope as compared to |
| 799 |
+ |
Ewald even for larger cutoff radii. It appears that the severity of |
| 800 |
+ |
the perturbations in the TSF method are most apparent in the torques. |
| 801 |
+ |
|
| 802 |
+ |
\subsection{Directionality of the force and torque vectors} |
| 803 |
+ |
|
| 804 |
+ |
The accurate evaluation of force and torque directions is just as |
| 805 |
+ |
important for molecular dynamics simulations as the magnitudes of |
| 806 |
+ |
these quantities. Force and torque vectors for all six systems were |
| 807 |
+ |
analyzed using Fisher statistics, and the quality of the vector |
| 808 |
+ |
directionality is shown in terms of circular variance |
| 809 |
+ |
($\mathrm{Var}(\theta$) in figure |
| 810 |
+ |
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
| 811 |
+ |
from the new real-space method exhibit nearly-ideal Fisher probability |
| 812 |
+ |
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
| 813 |
+ |
exhibit the best vectorial agreement with the Ewald sum. The force and |
| 814 |
+ |
torque vectors from GSF method also show good agreement with the Ewald |
| 815 |
+ |
method, which can also be systematically improved by using moderate |
| 816 |
+ |
damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c = |
| 817 |
+ |
12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds |
| 818 |
+ |
to a distribution with 95\% of force vectors within $6.37^\circ$ of the |
| 819 |
+ |
corresponding Ewald forces. The TSF method produces the poorest |
| 820 |
+ |
agreement with the Ewald force directions. |
| 821 |
+ |
|
| 822 |
+ |
Torques are again more perturbed by the new real-space methods, than |
| 823 |
+ |
forces, but even here the variance is reasonably small. For the same |
| 824 |
+ |
method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA), |
| 825 |
+ |
the circular variance was 0.01415, corresponds to a distribution which |
| 826 |
+ |
has 95\% of torque vectors are within $16.75^\circ$ of the Ewald |
| 827 |
+ |
results. Again, the direction of the force and torque vectors can be |
| 828 |
+ |
systematically improved by varying $\alpha$ and $r_c$. |
| 829 |
+ |
|
| 830 |
|
\begin{figure} |
| 831 |
< |
\centering |
| 832 |
< |
\includegraphics[width=0.5 \textwidth]{Variance_forceNtorque_modified-crop.pdf} |
| 833 |
< |
\caption{The circular variance of the data sets of the |
| 834 |
< |
direction of the force and torque vectors obtained from a |
| 835 |
< |
given method about reference Ewald method. The result equal |
| 836 |
< |
to 0 (dashed line) indicates direction of the vectors are |
| 837 |
< |
indistinguishable from the Ewald method. Here different |
| 838 |
< |
symbols represent different value of the cutoff radius (9 |
| 839 |
< |
\AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)} |
| 840 |
< |
\label{fig:slopeCorr_circularVariance} |
| 841 |
< |
\end{figure} |
| 688 |
< |
\subsection{Total energy conservation} |
| 689 |
< |
We have tested the conservation of energy in the SSDQC liquid system |
| 690 |
< |
by running system for 1ns in the Hard, SP, GSF and TSF method. The |
| 691 |
< |
Hard cutoff method shows very high energy drifts 433.53 |
| 692 |
< |
KCal/Mol/ns/particle and energy fluctuation 32574.53 Kcal/Mol |
| 693 |
< |
(measured by the SD from the slope) for the undamped case, which makes |
| 694 |
< |
it completely unusable in MD simulations. The SP method also shows |
| 695 |
< |
large value of energy drift 1.289 Kcal/Mol/ns/particle and energy |
| 696 |
< |
fluctuation 0.01953 Kcal/Mol. The energy fluctuation in the SP method |
| 697 |
< |
is due to the non-vanishing nature of the torque and force at the |
| 698 |
< |
cutoff radius. We can improve the energy conservation in some extent |
| 699 |
< |
by the proper selection of the damping alpha but the improvement is |
| 700 |
< |
not good enough, which can be observed in Figure 9a and 9b .The GSF |
| 701 |
< |
and TSF shows very low value of energy drift 0.09016, 0.07371 |
| 702 |
< |
KCal/Mol/ns/particle and fluctuation 3.016E-6, 1.217E-6 Kcal/Mol |
| 703 |
< |
respectively for the undamped case. Since the absolute value of the |
| 704 |
< |
evaluated electrostatic energy, force and torque from TSF method are |
| 705 |
< |
deviated from the Ewald, it does not mimic MD simulations |
| 706 |
< |
appropriately. The electrostatic energy, force and torque from the GSF |
| 707 |
< |
method have very good agreement with the Ewald. In addition, the |
| 708 |
< |
energy drift and energy fluctuation from the GSF method is much better |
| 709 |
< |
than Ewald’s method for reciprocal space vector value ($k_f$) equal to |
| 710 |
< |
7 as shown in Figure~\ref{fig:energyDrift} and |
| 711 |
< |
~\ref{fig:fluctuation}. We can improve the total energy fluctuation |
| 712 |
< |
and drift for the Ewald’s method by increasing size of the reciprocal |
| 713 |
< |
space, which extremely increseses the simulation time. In our current |
| 714 |
< |
simulation, the simulation time for the Hard, SP, and GSF methods are |
| 715 |
< |
about 5.5 times faster than the Ewald method. |
| 831 |
> |
\centering |
| 832 |
> |
\includegraphics[width=0.6 \linewidth]{Variance_forceNtorque_modified-crop.pdf} |
| 833 |
> |
\caption{The circular variance of the direction of the force and |
| 834 |
> |
torque vectors obtained from the real-space methods around the |
| 835 |
> |
reference Ewald vectors. A variance equal to 0 (dashed line) |
| 836 |
> |
indicates direction of the force or torque vectors are |
| 837 |
> |
indistinguishable from those obtained from the Ewald sum. Here |
| 838 |
> |
different symbols represent different values of the cutoff radius |
| 839 |
> |
(9 \AA\ = circle, 12 \AA\ = square, 15 \AA\ = inverted triangle)} |
| 840 |
> |
\label{fig:slopeCorr_circularVariance} |
| 841 |
> |
\end{figure} |
| 842 |
|
|
| 843 |
< |
In Fig.~\ref{fig:energyDrift}, $\delta \mbox{E}_1$ is a measure of the |
| 718 |
< |
linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle |
| 719 |
< |
over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the |
| 720 |
< |
standard deviation of the energy fluctuations in units of $\mbox{kcal |
| 721 |
< |
mol}^{-1}$ per particle. In the bottom plot, it is apparent that the |
| 722 |
< |
energy drift is reduced by a significant amount (2 to 3 orders of |
| 723 |
< |
magnitude improvement at all values of the damping coefficient) by |
| 724 |
< |
chosing either of the shifted-force methods over the hard or SP |
| 725 |
< |
methods. We note that the two shifted-force method can give |
| 726 |
< |
significantly better energy conservation than the multipolar Ewald sum |
| 727 |
< |
with the same choice of real-space cutoffs. |
| 843 |
> |
\subsection{Energy conservation} |
| 844 |
|
|
| 845 |
+ |
We have tested the conservation of energy one can expect to see with |
| 846 |
+ |
the new real-space methods using the SSDQ water model with a small |
| 847 |
+ |
fraction of solvated ions. This is a test system which exercises all |
| 848 |
+ |
orders of multipole-multipole interactions derived in the first paper |
| 849 |
+ |
in this series and provides the most comprehensive test of the new |
| 850 |
+ |
methods. A liquid-phase system was created with 2000 water molecules |
| 851 |
+ |
and 48 dissolved ions at a density of 0.98 g cm${-3}$ and a |
| 852 |
+ |
temperature of 300K. After equilibration, this liquid-phase system |
| 853 |
+ |
was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with |
| 854 |
+ |
a cutoff radius of 9\AA. The value of the damping coefficient was |
| 855 |
+ |
also varied from the undamped case ($\alpha = 0$) to a heavily damped |
| 856 |
+ |
case ($\alpha = 0.3$ \AA$^{-1}$) for the real space methods. A sample |
| 857 |
+ |
was also run using the multipolar Ewald sum. |
| 858 |
+ |
|
| 859 |
+ |
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
| 860 |
+ |
energy over time, $\delta E_1$, and the standard deviation of energy |
| 861 |
+ |
fluctuations around this drift $\delta E_0$. Both of the |
| 862 |
+ |
shifted-force methods (GSF and TSF) provide excellent energy |
| 863 |
+ |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
| 864 |
+ |
while the hard cutoff is essentially unusable for molecular dynamics. |
| 865 |
+ |
SP provides some benefit over the hard cutoff because the energetic |
| 866 |
+ |
jumps that happen as particles leave and enter the cutoff sphere are |
| 867 |
+ |
somewhat reduced. |
| 868 |
+ |
|
| 869 |
+ |
We note that for all tested values of the cutoff radius, the new |
| 870 |
+ |
real-space methods can provide better energy conservation behavior |
| 871 |
+ |
than the multipolar Ewald sum, even when utilizing a relatively large |
| 872 |
+ |
$k$-space cutoff values. |
| 873 |
+ |
|
| 874 |
|
\begin{figure} |
| 875 |
|
\centering |
| 876 |
|
\includegraphics[width=\textwidth]{newDrift.pdf} |
| 877 |
|
\label{fig:energyDrift} |
| 878 |
< |
\caption{Analysis of the energy conservation of the real space |
| 878 |
> |
\caption{Analysis of the energy conservation of the real-space |
| 879 |
|
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
| 880 |
|
energy over time and $\delta \mathrm{E}_0$ is the standard deviation |
| 881 |
|
of energy fluctuations around this drift. All simulations were of a |
| 883 |
|
K starting from the same initial configuration.} |
| 884 |
|
\end{figure} |
| 885 |
|
|
| 886 |
+ |
|
| 887 |
|
\section{CONCLUSION} |
| 888 |
< |
We have generalized the charged neutralized potential energy originally developed by the Wolf et al.\cite{Wolf:1999dn} for the charge-charge interaction to the charge-multipole and multipole-multipole interaction in the SP method for higher order multipoles. Also, we have developed GSF and TSF methods by implementing the modification purposed by Fennel and Gezelter\cite{Fennell:2006lq} for the charge-charge interaction to the higher order multipoles to ensure consistency and smooth truncation of the electrostatic energy, force, and torque for the spherical truncation. The SP methods for multipoles proved its suitability in MC simulations. On the other hand, the results from the GSF method produced good agreement with the Ewald's energy, force, and torque. Also, it shows very good energy conservation in MD simulations. |
| 889 |
< |
The direct truncation of any molecular system without multipole neutralization creates the fluctuation in the electrostatic energy. This fluctuation in the energy is very large for the case of crystal because of long range of multipole ordering (Refer paper I).\cite{PaperI} This is also significant in the case of the liquid because of the local multipole ordering in the molecules. If the net multipole within cutoff radius neutralized within cutoff sphere by placing image multiples on the surface of the sphere, this fluctuation in the energy reduced significantly. Also, the multipole neutralization in the generalized SP method showed very good agreement with the Ewald as compared to direct truncation for the evaluation of the $\triangle E$ between the configurations. |
| 890 |
< |
In MD simulations, the energy conservation is very important. The |
| 891 |
< |
conservation of the total energy can be ensured by i) enforcing the |
| 892 |
< |
smooth truncation of the energy, force and torque in the cutoff radius |
| 893 |
< |
and ii) making the energy, force and torque consistent with each |
| 894 |
< |
other. The GSF and TSF methods ensure the consistency and smooth |
| 895 |
< |
truncation of the energy, force and torque at the cutoff radius, as a |
| 896 |
< |
result show very good total energy conservation. But the TSF method |
| 897 |
< |
does not show good agreement in the absolute value of the |
| 898 |
< |
electrostatic energy, force and torque with the Ewald. The GSF method |
| 899 |
< |
has mimicked Ewald’s force, energy and torque accurately and also |
| 900 |
< |
conserved energy. Therefore, the GSF method is the suitable method for |
| 901 |
< |
evaluating required force field in MD simulations. In addition, the |
| 902 |
< |
energy drift and fluctuation from the GSF method is much better than |
| 903 |
< |
Ewald’s method for finite-sized reciprocal space. |
| 888 |
> |
We have generalized the charged neutralized potential energy |
| 889 |
> |
originally developed by the Wolf et al.\cite{Wolf:1999dn} for the |
| 890 |
> |
charge-charge interaction to the charge-multipole and |
| 891 |
> |
multipole-multipole interaction in the SP method for higher order |
| 892 |
> |
multipoles. Also, we have developed GSF and TSF methods by |
| 893 |
> |
implementing the modification purposed by Fennel and |
| 894 |
> |
Gezelter\cite{Fennell:2006lq} for the charge-charge interaction to the |
| 895 |
> |
higher order multipoles to ensure consistency and smooth truncation of |
| 896 |
> |
the electrostatic energy, force, and torque for the spherical |
| 897 |
> |
truncation. The SP methods for multipoles proved its suitability in MC |
| 898 |
> |
simulations. On the other hand, the results from the GSF method |
| 899 |
> |
produced good agreement with the Ewald's energy, force, and |
| 900 |
> |
torque. Also, it shows very good energy conservation in MD |
| 901 |
> |
simulations. The direct truncation of any molecular system without |
| 902 |
> |
multipole neutralization creates the fluctuation in the electrostatic |
| 903 |
> |
energy. This fluctuation in the energy is very large for the case of |
| 904 |
> |
crystal because of long range of multipole ordering (Refer paper |
| 905 |
> |
I).\cite{PaperI} This is also significant in the case of the liquid |
| 906 |
> |
because of the local multipole ordering in the molecules. If the net |
| 907 |
> |
multipole within cutoff radius neutralized within cutoff sphere by |
| 908 |
> |
placing image multiples on the surface of the sphere, this fluctuation |
| 909 |
> |
in the energy reduced significantly. Also, the multipole |
| 910 |
> |
neutralization in the generalized SP method showed very good agreement |
| 911 |
> |
with the Ewald as compared to direct truncation for the evaluation of |
| 912 |
> |
the $\triangle E$ between the configurations. In MD simulations, the |
| 913 |
> |
energy conservation is very important. The conservation of the total |
| 914 |
> |
energy can be ensured by i) enforcing the smooth truncation of the |
| 915 |
> |
energy, force and torque in the cutoff radius and ii) making the |
| 916 |
> |
energy, force and torque consistent with each other. The GSF and TSF |
| 917 |
> |
methods ensure the consistency and smooth truncation of the energy, |
| 918 |
> |
force and torque at the cutoff radius, as a result show very good |
| 919 |
> |
total energy conservation. But the TSF method does not show good |
| 920 |
> |
agreement in the absolute value of the electrostatic energy, force and |
| 921 |
> |
torque with the Ewald. The GSF method has mimicked Ewald’s force, |
| 922 |
> |
energy and torque accurately and also conserved energy. Therefore, the |
| 923 |
> |
GSF method is the suitable method for evaluating required force field |
| 924 |
> |
in MD simulations. In addition, the energy drift and fluctuation from |
| 925 |
> |
the GSF method is much better than Ewald’s method for finite-sized |
| 926 |
> |
reciprocal space. |
| 927 |
|
|
| 928 |
|
Note that the TSF, GSF, and SP models are $\mathcal{O}(N)$ methods |
| 929 |
|
that can be made extremely efficient using spline interpolations of |
