| 36 |
|
\usepackage{amsmath} |
| 37 |
|
\usepackage{times} |
| 38 |
|
\usepackage{mathptm} |
| 39 |
+ |
\usepackage{tabularx} |
| 40 |
|
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 41 |
|
\usepackage{url} |
| 42 |
|
\usepackage[english]{babel} |
| 43 |
|
|
| 44 |
+ |
\newcolumntype{Y}{>{\centering\arraybackslash}X} |
| 45 |
|
|
| 46 |
|
\begin{document} |
| 47 |
|
|
| 48 |
< |
\preprint{AIP/123-QED} |
| 48 |
> |
%\preprint{AIP/123-QED} |
| 49 |
|
|
| 50 |
< |
\title[Efficient electrostatics for condensed-phase multipoles]{Real space alternatives to the Ewald |
| 51 |
< |
Sum. II. Comparison of Simulation Methodologies} % Force line breaks with \\ |
| 50 |
> |
\title{Real space alternatives to the Ewald |
| 51 |
> |
Sum. II. Comparison of Methods} % Force line breaks with \\ |
| 52 |
|
|
| 53 |
|
\author{Madan Lamichhane} |
| 54 |
|
\affiliation{Department of Physics, University |
| 67 |
|
% but any date may be explicitly specified |
| 68 |
|
|
| 69 |
|
\begin{abstract} |
| 70 |
< |
We have tested our recently developed shifted potential, gradient-shifted force, and Taylor-shifted force methods for the higher-order multipoles against Ewald’s method in different types of liquid and crystalline system. In this paper, we have also investigated the conservation of total energy in the molecular dynamic simulation using all of these methods. The shifted potential method shows better agreement with the Ewald in the energy differences between different configurations as compared to the direct truncation. Both the gradient shifted force and Taylor-shifted force methods reproduce very good energy conservation. But the absolute energy, force and torque evaluated from the gradient shifted force method shows better result as compared to taylor-shifted force method. Hence the gradient-shifted force method suitably mimics the electrostatic interaction in the molecular dynamic simulation. |
| 70 |
> |
We have tested the real-space shifted potential (SP), |
| 71 |
> |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
| 72 |
> |
for multipoles that were developed in the first paper in this series |
| 73 |
> |
against a reference method. The tests were carried out in a variety |
| 74 |
> |
of condensed-phase environments which were designed to test all |
| 75 |
> |
levels of the multipole-multipole interactions. Comparisons of the |
| 76 |
> |
energy differences between configurations, molecular forces, and |
| 77 |
> |
torques were used to analyze how well the real-space models perform |
| 78 |
> |
relative to the more computationally expensive Ewald sum. We have |
| 79 |
> |
also investigated the energy conservation properties of the new |
| 80 |
> |
methods in molecular dynamics simulations using all of these |
| 81 |
> |
methods. The SP method shows excellent agreement with |
| 82 |
> |
configurational energy differences, forces, and torques, and would |
| 83 |
> |
be suitable for use in Monte Carlo calculations. Of the two new |
| 84 |
> |
shifted-force methods, the GSF approach shows the best agreement |
| 85 |
> |
with Ewald-derived energies, forces, and torques and exhibits energy |
| 86 |
> |
conservation properties that make it an excellent choice for |
| 87 |
> |
efficiently computing electrostatic interactions in molecular |
| 88 |
> |
dynamics simulations. |
| 89 |
|
\end{abstract} |
| 90 |
|
|
| 91 |
< |
\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 91 |
> |
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 92 |
|
% Classification Scheme. |
| 93 |
|
\keywords{Electrostatics, Multipoles, Real-space} |
| 94 |
|
|
| 335 |
|
number of terms in the truncated Taylor expansion, e.g., |
| 336 |
|
% |
| 337 |
|
\begin{equation} |
| 338 |
< |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-R_c)^m}{m!} f^{(m)} \Big \lvert _{R_c} . |
| 338 |
> |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
| 339 |
|
\end{equation} |
| 340 |
|
% |
| 341 |
|
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 342 |
|
Thus, for $f(r)=1/r$, we find |
| 343 |
|
% |
| 344 |
|
\begin{equation} |
| 345 |
< |
f_1(r)=\frac{1}{r}- \frac{1}{R_c} + (r - R_c) \frac{1}{R_c^2} - \frac{(r-R_c)^2}{R_c^3} . |
| 345 |
> |
f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . |
| 346 |
|
\end{equation} |
| 347 |
|
This function is an approximate electrostatic potential that has |
| 348 |
|
vanishing second derivatives at the cutoff radius, making it suitable |
| 369 |
|
|
| 370 |
|
Note that increasing the value of $n$ will add additional terms to the |
| 371 |
|
electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 372 |
< |
$(r-R_c)^3/R_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 372 |
> |
$(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 373 |
|
functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 374 |
|
f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 375 |
|
for computing multipole energies, forces, and torques, and smooth |
| 397 |
|
Eq.~\ref{generic}, accounting for the appropriate number of |
| 398 |
|
derivatives. Complete energy, force, and torque expressions are |
| 399 |
|
presented in the first paper in this series (Reference |
| 400 |
< |
\citep{PaperI}). |
| 400 |
> |
\onlinecite{PaperI}). |
| 401 |
|
|
| 402 |
|
\subsection{Gradient-shifted force (GSF)} |
| 403 |
|
|
| 408 |
|
which have been projected onto the surface of the cutoff sphere |
| 409 |
|
without changing their relative orientation, |
| 410 |
|
\begin{displaymath} |
| 411 |
< |
U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(R_c) |
| 412 |
< |
- (r_{ij}-R_c) \hat{r}_{ij} \cdot |
| 413 |
< |
\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{R_c} |
| 411 |
> |
U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(r_c) |
| 412 |
> |
- (r_{ij}-r_c) \hat{r}_{ij} \cdot |
| 413 |
> |
\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{r_c} |
| 414 |
|
\end{displaymath} |
| 415 |
|
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$ |
| 416 |
|
and $\mathbf{D}_{j}$, are retained at the cutoff distance (although |
| 439 |
|
function of radius, hence the contribution of the third term is very |
| 440 |
|
small for large cutoff radii. The force and torque derived from |
| 441 |
|
equation \ref{generic2} are consistent with the energy expression and |
| 442 |
< |
approach zero as $r \rightarrow R_c$. Both the GSF and TSF methods |
| 442 |
> |
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
| 443 |
|
can be considered generalizations of the original DSF method for |
| 444 |
|
higher order multipole interactions. GSF and TSF are also identical up |
| 445 |
|
to the charge-dipole interaction but generate different expressions in |
| 446 |
|
the energy, force and torque for higher order multipole-multipole |
| 447 |
|
interactions. Complete energy, force, and torque expressions for the |
| 448 |
|
GSF potential are presented in the first paper in this series |
| 449 |
< |
(Reference \citep{PaperI}) |
| 449 |
> |
(Reference~\onlinecite{PaperI}) |
| 450 |
|
|
| 451 |
|
|
| 452 |
|
\subsection{Shifted potential (SP) } |
| 466 |
|
each orientational contribution to the energy (e.g. the direct dipole |
| 467 |
|
product contribution is shifted {\it separately} from the |
| 468 |
|
dipole-distance terms in dipole-dipole interactions). Note that this |
| 469 |
< |
is not a simple shifting of the total potential at $R_c$. Each radial |
| 469 |
> |
is not a simple shifting of the total potential at $r_c$. Each radial |
| 470 |
|
contribution is shifted separately. One consequence of this is that |
| 471 |
|
multipoles that reorient after leaving the cutoff sphere can re-enter |
| 472 |
|
the cutoff sphere without perturbing the total energy. |
| 473 |
|
|
| 474 |
|
The potential energy between a central multipole and other multipolar |
| 475 |
< |
sites then goes smoothly to zero as $r \rightarrow R_c$. However, the |
| 475 |
> |
sites then goes smoothly to zero as $r \rightarrow r_c$. However, the |
| 476 |
|
force and torque obtained from the shifted potential (SP) are |
| 477 |
< |
discontinuous at $R_c$. Therefore, MD simulations will still |
| 477 |
> |
discontinuous at $r_c$. Therefore, MD simulations will still |
| 478 |
|
experience energy drift while operating under the SP potential, but it |
| 479 |
|
may be suitable for Monte Carlo approaches where the configurational |
| 480 |
|
energy differences are the primary quantity of interest. |
| 485 |
|
the cutoff sphere. This self interaction is nearly identical with the |
| 486 |
|
self-terms that arise in the Ewald sum for multipoles. Complete |
| 487 |
|
expressions for the self terms are presented in the first paper in |
| 488 |
< |
this series (Reference \citep{PaperI}) |
| 488 |
> |
this series (Reference \onlinecite{PaperI}). |
| 489 |
|
|
| 490 |
|
|
| 491 |
|
\section{\label{sec:methodology}Methodology} |
| 500 |
|
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} In this work, we |
| 501 |
|
used the multipolar Ewald sum as a reference method for comparing |
| 502 |
|
energies, forces, and torques for molecular models that mimic |
| 503 |
< |
disordered and ordered condensed-phase systems. These test-cases |
| 504 |
< |
include: |
| 485 |
< |
\begin{itemize} |
| 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 |
< |
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. |
| 503 |
> |
disordered and ordered condensed-phase systems. The parameters used |
| 504 |
> |
in the test-cases are given in table~\ref{tab:pars}. |
| 505 |
|
|
| 506 |
< |
The last test case exercises all levels of the multipole-multipole |
| 507 |
< |
interactions we have derived so far and represents the most complete |
| 508 |
< |
test of the new methods. In the following section, we present results |
| 509 |
< |
for the total electrostatic energy, as well as the electrostatic |
| 510 |
< |
contributions to the force and torque on each molecule. These |
| 511 |
< |
quantities have been computed using the SP, TSF, and GSF methods, as |
| 512 |
< |
well as a hard cutoff, and have been compared with the values obtaine |
| 513 |
< |
from the multipolar Ewald sum. In Mote Carlo (MC) simulations, the |
| 514 |
< |
energy differences between two configurations is the primary quantity |
| 515 |
< |
that governs how the simulation proceeds. These differences are the |
| 516 |
< |
most imporant indicators of the reliability of a method even if the |
| 517 |
< |
absolute energies are not exact. For each of the multipolar systems |
| 518 |
< |
listed above, we have compared the change in electrostatic potential |
| 519 |
< |
energy ($\Delta E$) between 250 statistically-independent |
| 520 |
< |
configurations. In molecular dynamics (MD) simulations, the forces |
| 521 |
< |
and torques govern the behavior of the simulation, so we also compute |
| 522 |
< |
the electrostatic contributions to the forces and torques. |
| 506 |
> |
\begin{table} |
| 507 |
> |
\label{tab:pars} |
| 508 |
> |
\caption{The parameters used in the systems used to evaluate the new |
| 509 |
> |
real-space methods. The most comprehensive test was a liquid |
| 510 |
> |
composed of 2000 SSDQ molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-} |
| 511 |
> |
ions). This test excercises all orders of the multipolar |
| 512 |
> |
interactions developed in the first paper.} |
| 513 |
> |
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
| 514 |
> |
& \multicolumn{2}{c|}{LJ parameters} & |
| 515 |
> |
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
| 516 |
> |
Test system & $\sigma$& $\epsilon$ & $C$ & $D$ & |
| 517 |
> |
$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass & $I_{xx}$ & $I_{yy}$ & |
| 518 |
> |
$I_{zz}$ \\ \cline{6-8}\cline{10-12} |
| 519 |
> |
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
| 520 |
> |
\AA\textsuperscript{2})} \\ \hline |
| 521 |
> |
Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ |
| 522 |
> |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & 10,000 & 17.6 &17.6 & 0 \\ |
| 523 |
> |
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 524 |
> |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & 10,000 & 17.6&17.6&0 \\ |
| 525 |
> |
SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 526 |
> |
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
| 527 |
> |
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
| 528 |
> |
\end{tabularx} |
| 529 |
> |
\end{table} |
| 530 |
> |
The systems consist of pure multipolar solids (both dipole and |
| 531 |
> |
quadrupole), pure multipolar liquids (both dipole and quadrupole), a |
| 532 |
> |
fluid composed of sites containing both dipoles and quadrupoles |
| 533 |
> |
simultaneously, and a final test case that includes ions with point |
| 534 |
> |
charges in addition to the multipolar fluid. The solid-phase |
| 535 |
> |
parameters were chosen so that the systems can explore some |
| 536 |
> |
orientational freedom for the multipolar sites, while maintaining |
| 537 |
> |
relatively strict translational order. The SSDQ model used here is |
| 538 |
> |
not a particularly accurate water model, but it does test |
| 539 |
> |
dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole |
| 540 |
> |
interactions at roughly the same magnitudes. The last test case, SSDQ |
| 541 |
> |
water with dissolved ions, exercises \textit{all} levels of the |
| 542 |
> |
multipole-multipole interactions we have derived so far and represents |
| 543 |
> |
the most complete test of the new methods. |
| 544 |
|
|
| 545 |
+ |
In the following section, we present results for the total |
| 546 |
+ |
electrostatic energy, as well as the electrostatic contributions to |
| 547 |
+ |
the force and torque on each molecule. These quantities have been |
| 548 |
+ |
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
| 549 |
+ |
and have been compared with the values obtaine from the multipolar |
| 550 |
+ |
Ewald sum. In Mote Carlo (MC) simulations, the energy differences |
| 551 |
+ |
between two configurations is the primary quantity that governs how |
| 552 |
+ |
the simulation proceeds. These differences are the most imporant |
| 553 |
+ |
indicators of the reliability of a method even if the absolute |
| 554 |
+ |
energies are not exact. For each of the multipolar systems listed |
| 555 |
+ |
above, we have compared the change in electrostatic potential energy |
| 556 |
+ |
($\Delta E$) between 250 statistically-independent configurations. In |
| 557 |
+ |
molecular dynamics (MD) simulations, the forces and torques govern the |
| 558 |
+ |
behavior of the simulation, so we also compute the electrostatic |
| 559 |
+ |
contributions to the forces and torques. |
| 560 |
+ |
|
| 561 |
+ |
\subsection{Implementation} |
| 562 |
+ |
The real-space methods developed in the first paper in this series |
| 563 |
+ |
have been implemented in our group's open source molecular simulation |
| 564 |
+ |
program, OpenMD,\cite{openmd} which was used for all calculations in |
| 565 |
+ |
this work. The complementary error function can be a relatively slow |
| 566 |
+ |
function on some processors, so all of the radial functions are |
| 567 |
+ |
precomputed on a fine grid and are spline-interpolated to provide |
| 568 |
+ |
values when required. |
| 569 |
+ |
|
| 570 |
+ |
Using the same simulation code, we compare to a multipolar Ewald sum |
| 571 |
+ |
with a reciprocal space cutoff, $k_\mathrm{max} = 7$. Our version of |
| 572 |
+ |
the Ewald sum is a re-implementation of the algorithm originally |
| 573 |
+ |
proposed by Smith that does not use the particle mesh or smoothing |
| 574 |
+ |
approximations.\cite{Smith82,Smith98} In all cases, the quantities |
| 575 |
+ |
being compared are the electrostatic contributions to energies, force, |
| 576 |
+ |
and torques. All other contributions to these quantities (i.e. from |
| 577 |
+ |
Lennard-Jones interactions) are removed prior to the comparisons. |
| 578 |
+ |
|
| 579 |
+ |
The convergence parameter ($\alpha$) also plays a role in the balance |
| 580 |
+ |
of the real-space and reciprocal-space portions of the Ewald |
| 581 |
+ |
calculation. Typical molecular mechanics packages set this to a value |
| 582 |
+ |
that depends on the cutoff radius and a tolerance (typically less than |
| 583 |
+ |
$1 \times 10^{-4}$ kcal/mol). Smaller tolerances are typically |
| 584 |
+ |
associated with increasing accuracy at the expense of computational |
| 585 |
+ |
time spent on the reciprocal-space portion of the |
| 586 |
+ |
summation.\cite{Perram88,Essmann:1995pb} A default tolerance of $1 \times |
| 587 |
+ |
10^{-8}$ kcal/mol was used in all Ewald calculations, resulting in |
| 588 |
+ |
Ewald coefficient 0.3119 \AA$^{-1}$ for a cutoff radius of 12 \AA. |
| 589 |
+ |
|
| 590 |
+ |
The real-space models have self-interactions that provide |
| 591 |
+ |
contributions to the energies only. Although the self interaction is |
| 592 |
+ |
a rapid calculation, we note that in systems with fluctuating charges |
| 593 |
+ |
or point polarizabilities, the self-term is not static and must be |
| 594 |
+ |
recomputed at each time step. |
| 595 |
+ |
|
| 596 |
|
\subsection{Model systems} |
| 597 |
|
To sample independent configurations of multipolar crystals, a body |
| 598 |
|
centered cubic (bcc) crystal which is a minimum energy structure for |
| 778 |
|
|
| 779 |
|
The combined correlation coefficient and slope for all six systems is |
| 780 |
|
shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the methods |
| 781 |
< |
reproduce the Ewald-derived configurational energy differences with |
| 782 |
< |
remarkable fidelity. Undamped hard cutoffs introduce a significant |
| 783 |
< |
amount of random scatter in the energy differences which is apparent |
| 784 |
< |
in the reduced value of the correlation coefficient for this method. |
| 785 |
< |
This can be understood easily as configurations which exhibit only |
| 786 |
< |
small traversals of a few dipoles or quadrupoles out of the cutoff |
| 787 |
< |
sphere will see large energy jumps when hard cutoffs are used. The |
| 781 |
> |
reproduce the Ewald configurational energy differences with remarkable |
| 782 |
> |
fidelity. Undamped hard cutoffs introduce a significant amount of |
| 783 |
> |
random scatter in the energy differences which is apparent in the |
| 784 |
> |
reduced value of the correlation coefficient for this method. This |
| 785 |
> |
can be easily understood as configurations which exhibit small |
| 786 |
> |
traversals of a few dipoles or quadrupoles out of the cutoff sphere |
| 787 |
> |
will see large energy jumps when hard cutoffs are used. The |
| 788 |
|
orientations of the multipoles (particularly in the ordered crystals) |
| 789 |
< |
mean that these jumps can go either up or down in energy, producing a |
| 790 |
< |
significant amount of random scatter. |
| 789 |
> |
mean that these energy jumps can go in either direction, producing a |
| 790 |
> |
significant amount of random scatter, but no systematic error. |
| 791 |
|
|
| 792 |
|
The TSF method produces energy differences that are highly correlated |
| 793 |
|
with the Ewald results, but it also introduces a significant |
| 795 |
|
smaller cutoff values. The TSF method alters the distance dependence |
| 796 |
|
of different orientational contributions to the energy in a |
| 797 |
|
non-uniform way, so the size of the cutoff sphere can have a large |
| 798 |
< |
effect on crystalline systems. |
| 798 |
> |
effect, particularly for the crystalline systems. |
| 799 |
|
|
| 800 |
|
Both the SP and GSF methods appear to reproduce the Ewald results with |
| 801 |
|
excellent fidelity, particularly for moderate damping ($\alpha = |
| 802 |
< |
0.1-0.2$\AA$^{-1}$) and commonly-used cutoff values ($r_c = 12$\AA). |
| 803 |
< |
With the exception of the undamped hard cutoff, and the TSF method |
| 804 |
< |
with short cutoffs, all of the methods would be appropriate for use in |
| 805 |
< |
Monte Carlo simulations. |
| 802 |
> |
0.1-0.2$\AA$^{-1}$) and with a commonly-used cutoff value ($r_c = |
| 803 |
> |
12$\AA). With the exception of the undamped hard cutoff, and the TSF |
| 804 |
> |
method with short cutoffs, all of the methods would be appropriate for |
| 805 |
> |
use in Monte Carlo simulations. |
| 806 |
|
|
| 807 |
|
\subsection{Magnitude of the force and torque vectors} |
| 808 |
|
|
| 809 |
< |
The comparison of the magnitude of the combined forces and torques for |
| 810 |
< |
the data accumulated from all system types are shown in Figures |
| 809 |
> |
The comparisons of the magnitudes of the forces and torques for the |
| 810 |
> |
data accumulated from all six systems are shown in Figures |
| 811 |
|
~\ref{fig:slopeCorr_force} and \ref{fig:slopeCorr_torque}. The |
| 812 |
|
correlation and slope for the forces agree well with the Ewald sum |
| 813 |
< |
even for the hard cutoff method. |
| 813 |
> |
even for the hard cutoffs. |
| 814 |
|
|
| 815 |
< |
For the system of molecules with higher order multipoles, the |
| 816 |
< |
interaction is quite short ranged. Moreover, the force decays more |
| 817 |
< |
rapidly than the electrostatic energy hence the hard cutoff method can |
| 818 |
< |
also produces reasonable agreement. Although the pure cutoff gives |
| 819 |
< |
the good match of the electrostatic force for pairs of molecules |
| 820 |
< |
included within the cutoff sphere, the discontinuity in the force at |
| 821 |
< |
the cutoff radius can potentially cause problems the total energy |
| 822 |
< |
conservation as molecules enter and leave the cutoff sphere. This is |
| 823 |
< |
discussed in detail in section \ref{sec:}. |
| 815 |
> |
For systems of molecules with only multipolar interactions, the pair |
| 816 |
> |
energy contributions are quite short ranged. Moreover, the force |
| 817 |
> |
decays more rapidly than the electrostatic energy, hence the hard |
| 818 |
> |
cutoff method can also produce reasonable agreement for this quantity. |
| 819 |
> |
Although the pure cutoff gives reasonably good electrostatic forces |
| 820 |
> |
for pairs of molecules included within each other's cutoff spheres, |
| 821 |
> |
the discontinuity in the force at the cutoff radius can potentially |
| 822 |
> |
cause energy conservation problems as molecules enter and leave the |
| 823 |
> |
cutoff spheres. This is discussed in detail in section |
| 824 |
> |
\ref{sec:conservation}. |
| 825 |
|
|
| 826 |
|
The two shifted-force methods (GSF and TSF) exhibit a small amount of |
| 827 |
|
systematic variation and scatter compared with the Ewald forces. The |
| 828 |
|
shifted-force models intentionally perturb the forces between pairs of |
| 829 |
< |
molecules inside the cutoff sphere in order to correct the energy |
| 830 |
< |
conservation issues, so it is not particularly surprising that this |
| 831 |
< |
perturbation is evident in these same molecular forces. The GSF |
| 829 |
> |
molecules inside each other's cutoff spheres in order to correct the |
| 830 |
> |
energy conservation issues, and this perturbation is evident in the |
| 831 |
> |
statistics accumulated for the molecular forces. The GSF |
| 832 |
|
perturbations are minimal, particularly for moderate damping and and |
| 833 |
|
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
| 834 |
|
reasonable agreement in the correlation coefficient but again the |
| 871 |
|
|
| 872 |
|
The results shows that the torque from the hard cutoff method |
| 873 |
|
reproduces the torques in quite good agreement with the Ewald sum. |
| 874 |
< |
The other real-space methods can cause some significant deviations, |
| 875 |
< |
but excellent agreement with the Ewald sum torques is recovered at |
| 876 |
< |
moderate values of the damping coefficient ($\alpha = |
| 877 |
< |
0.1-0.2$\AA$^{-1}$) and cutoff radius ($r_c \ge 12$\AA). The TSF |
| 878 |
< |
method exhibits the only fair agreement in the slope as compared to |
| 879 |
< |
Ewald even for larger cutoff radii. It appears that the severity of |
| 880 |
< |
the perturbations in the TSF method are most apparent in the torques. |
| 874 |
> |
The other real-space methods can cause some deviations, but excellent |
| 875 |
> |
agreement with the Ewald sum torques is recovered at moderate values |
| 876 |
> |
of the damping coefficient ($\alpha = 0.1-0.2$\AA$^{-1}$) and cutoff |
| 877 |
> |
radius ($r_c \ge 12$\AA). The TSF method exhibits only fair agreement |
| 878 |
> |
in the slope when compared with the Ewald torques even for larger |
| 879 |
> |
cutoff radii. It appears that the severity of the perturbations in |
| 880 |
> |
the TSF method are most in evidence for the torques. |
| 881 |
|
|
| 882 |
|
\subsection{Directionality of the force and torque vectors} |
| 883 |
|
|
| 888 |
|
directionality is shown in terms of circular variance |
| 889 |
|
($\mathrm{Var}(\theta$) in figure |
| 890 |
|
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
| 891 |
< |
from the new real-space method exhibit nearly-ideal Fisher probability |
| 891 |
> |
from the new real-space methods exhibit nearly-ideal Fisher probability |
| 892 |
|
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
| 893 |
|
exhibit the best vectorial agreement with the Ewald sum. The force and |
| 894 |
|
torque vectors from GSF method also show good agreement with the Ewald |
| 895 |
|
method, which can also be systematically improved by using moderate |
| 896 |
< |
damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c = |
| 896 |
> |
damping and a reasonable cutoff radius. For $\alpha = 0.2$ and $r_c = |
| 897 |
|
12$\AA, we observe $\mathrm{Var}(\theta) = 0.00206$, which corresponds |
| 898 |
< |
to a distribution with 95\% of force vectors within $6.37^\circ$ of the |
| 899 |
< |
corresponding Ewald forces. The TSF method produces the poorest |
| 898 |
> |
to a distribution with 95\% of force vectors within $6.37^\circ$ of |
| 899 |
> |
the corresponding Ewald forces. The TSF method produces the poorest |
| 900 |
|
agreement with the Ewald force directions. |
| 901 |
|
|
| 902 |
< |
Torques are again more perturbed by the new real-space methods, than |
| 903 |
< |
forces, but even here the variance is reasonably small. For the same |
| 902 |
> |
Torques are again more perturbed than the forces by the new real-space |
| 903 |
> |
methods, but even here the variance is reasonably small. For the same |
| 904 |
|
method (GSF) with the same parameters ($\alpha = 0.2, r_c = 12$\AA), |
| 905 |
|
the circular variance was 0.01415, corresponds to a distribution which |
| 906 |
|
has 95\% of torque vectors are within $16.75^\circ$ of the Ewald |
| 920 |
|
\label{fig:slopeCorr_circularVariance} |
| 921 |
|
\end{figure} |
| 922 |
|
|
| 923 |
< |
\subsection{Energy conservation} |
| 923 |
> |
\subsection{Energy conservation\label{sec:conservation}} |
| 924 |
|
|
| 925 |
|
We have tested the conservation of energy one can expect to see with |
| 926 |
|
the new real-space methods using the SSDQ water model with a small |
| 928 |
|
orders of multipole-multipole interactions derived in the first paper |
| 929 |
|
in this series and provides the most comprehensive test of the new |
| 930 |
|
methods. A liquid-phase system was created with 2000 water molecules |
| 931 |
< |
and 48 dissolved ions at a density of 0.98 g cm${-3}$ and a |
| 931 |
> |
and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
| 932 |
|
temperature of 300K. After equilibration, this liquid-phase system |
| 933 |
|
was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with |
| 934 |
< |
a cutoff radius of 9\AA. The value of the damping coefficient was |
| 934 |
> |
a cutoff radius of 12\AA. The value of the damping coefficient was |
| 935 |
|
also varied from the undamped case ($\alpha = 0$) to a heavily damped |
| 936 |
< |
case ($\alpha = 0.3$ \AA$^{-1}$) for the real space methods. A sample |
| 937 |
< |
was also run using the multipolar Ewald sum. |
| 936 |
> |
case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods. A |
| 937 |
> |
sample was also run using the multipolar Ewald sum with the same |
| 938 |
> |
real-space cutoff. |
| 939 |
|
|
| 940 |
|
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
| 941 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
| 945 |
|
while the hard cutoff is essentially unusable for molecular dynamics. |
| 946 |
|
SP provides some benefit over the hard cutoff because the energetic |
| 947 |
|
jumps that happen as particles leave and enter the cutoff sphere are |
| 948 |
< |
somewhat reduced. |
| 948 |
> |
somewhat reduced, but like the Wolf method for charges, the SP method |
| 949 |
> |
would not be as useful for molecular dynamics as either of the |
| 950 |
> |
shifted-force methods. |
| 951 |
|
|
| 952 |
|
We note that for all tested values of the cutoff radius, the new |
| 953 |
|
real-space methods can provide better energy conservation behavior |
| 963 |
|
energy over time and $\delta \mathrm{E}_0$ is the standard deviation |
| 964 |
|
of energy fluctuations around this drift. All simulations were of a |
| 965 |
|
2000-molecule simulation of SSDQ water with 48 ionic charges at 300 |
| 966 |
< |
K starting from the same initial configuration.} |
| 966 |
> |
K starting from the same initial configuration. All runs utilized |
| 967 |
> |
the same real-space cutoff, $r_c = 12$\AA.} |
| 968 |
|
\end{figure} |
| 969 |
|
|
| 970 |
|
|
| 971 |
|
\section{CONCLUSION} |
| 972 |
< |
We have generalized the charged neutralized potential energy |
| 973 |
< |
originally developed by the Wolf et al.\cite{Wolf:1999dn} for the |
| 974 |
< |
charge-charge interaction to the charge-multipole and |
| 975 |
< |
multipole-multipole interaction in the SP method for higher order |
| 976 |
< |
multipoles. Also, we have developed GSF and TSF methods by |
| 977 |
< |
implementing the modification purposed by Fennel and |
| 978 |
< |
Gezelter\cite{Fennell:2006lq} for the charge-charge interaction to the |
| 979 |
< |
higher order multipoles to ensure consistency and smooth truncation of |
| 980 |
< |
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. |
| 972 |
> |
In the first paper in this series, we generalized the |
| 973 |
> |
charge-neutralized electrostatic energy originally developed by Wolf |
| 974 |
> |
\textit{et al.}\cite{Wolf:1999dn} to multipole-multipole interactions |
| 975 |
> |
up to quadrupolar order. The SP method is essentially a |
| 976 |
> |
multipole-capable version of the Wolf model. The SP method for |
| 977 |
> |
multipoles provides excellent agreement with Ewald-derived energies, |
| 978 |
> |
forces and torques, and is suitable for Monte Carlo simulations, |
| 979 |
> |
although the forces and torques retain discontinuities at the cutoff |
| 980 |
> |
distance that prevents its use in molecular dynamics. |
| 981 |
|
|
| 982 |
< |
Note that the TSF, GSF, and SP models are $\mathcal{O}(N)$ methods |
| 983 |
< |
that can be made extremely efficient using spline interpolations of |
| 984 |
< |
the radial functions. They require no Fourier transforms or $k$-space |
| 985 |
< |
sums, and guarantee the smooth handling of energies, forces, and |
| 986 |
< |
torques as multipoles cross the real-space cutoff boundary. |
| 982 |
> |
We also developed two natural extensions of the damped shifted-force |
| 983 |
> |
(DSF) model originally proposed by Fennel and |
| 984 |
> |
Gezelter.\cite{Fennell:2006lq} The GSF and TSF approaches provide |
| 985 |
> |
smooth truncation of energies, forces, and torques at the real-space |
| 986 |
> |
cutoff, and both converge to DSF electrostatics for point-charge |
| 987 |
> |
interactions. The TSF model is based on a high-order truncated Taylor |
| 988 |
> |
expansion which can be relatively perturbative inside the cutoff |
| 989 |
> |
sphere. The GSF model takes the gradient from an images of the |
| 990 |
> |
interacting multipole that has been projected onto the cutoff sphere |
| 991 |
> |
to derive shifted force and torque expressions, and is a significantly |
| 992 |
> |
more gentle approach. |
| 993 |
|
|
| 994 |
+ |
Of the two newly-developed shifted force models, the GSF method |
| 995 |
+ |
produced quantitative agreement with Ewald energy, force, and torques. |
| 996 |
+ |
It also performs well in conserving energy in MD simulations. The |
| 997 |
+ |
Taylor-shifted (TSF) model provides smooth dynamics, but these take |
| 998 |
+ |
place on a potential energy surface that is significantly perturbed |
| 999 |
+ |
from Ewald-based electrostatics. |
| 1000 |
+ |
|
| 1001 |
+ |
% The direct truncation of any electrostatic potential energy without |
| 1002 |
+ |
% multipole neutralization creates large fluctuations in molecular |
| 1003 |
+ |
% simulations. This fluctuation in the energy is very large for the case |
| 1004 |
+ |
% of crystal because of long range of multipole ordering (Refer paper |
| 1005 |
+ |
% I).\cite{PaperI} This is also significant in the case of the liquid |
| 1006 |
+ |
% because of the local multipole ordering in the molecules. If the net |
| 1007 |
+ |
% multipole within cutoff radius neutralized within cutoff sphere by |
| 1008 |
+ |
% placing image multiples on the surface of the sphere, this fluctuation |
| 1009 |
+ |
% in the energy reduced significantly. Also, the multipole |
| 1010 |
+ |
% neutralization in the generalized SP method showed very good agreement |
| 1011 |
+ |
% with the Ewald as compared to direct truncation for the evaluation of |
| 1012 |
+ |
% the $\triangle E$ between the configurations. In MD simulations, the |
| 1013 |
+ |
% energy conservation is very important. The conservation of the total |
| 1014 |
+ |
% energy can be ensured by i) enforcing the smooth truncation of the |
| 1015 |
+ |
% energy, force and torque in the cutoff radius and ii) making the |
| 1016 |
+ |
% energy, force and torque consistent with each other. The GSF and TSF |
| 1017 |
+ |
% methods ensure the consistency and smooth truncation of the energy, |
| 1018 |
+ |
% force and torque at the cutoff radius, as a result show very good |
| 1019 |
+ |
% total energy conservation. But the TSF method does not show good |
| 1020 |
+ |
% agreement in the absolute value of the electrostatic energy, force and |
| 1021 |
+ |
% torque with the Ewald. The GSF method has mimicked Ewald’s force, |
| 1022 |
+ |
% energy and torque accurately and also conserved energy. |
| 1023 |
+ |
|
| 1024 |
+ |
The only cases we have found where the new GSF and SP real-space |
| 1025 |
+ |
methods can be problematic are those which retain a bulk dipole moment |
| 1026 |
+ |
at large distances (e.g. the $Z_1$ dipolar lattice). In ferroelectric |
| 1027 |
+ |
materials, uniform weighting of the orientational contributions can be |
| 1028 |
+ |
important for converging the total energy. In these cases, the |
| 1029 |
+ |
damping function which causes the non-uniform weighting can be |
| 1030 |
+ |
replaced by the bare electrostatic kernel, and the energies return to |
| 1031 |
+ |
the expected converged values. |
| 1032 |
+ |
|
| 1033 |
+ |
Based on the results of this work, the GSF method is a suitable and |
| 1034 |
+ |
efficient replacement for the Ewald sum for evaluating electrostatic |
| 1035 |
+ |
interactions in MD simulations. Both methods retain excellent |
| 1036 |
+ |
fidelity to the Ewald energies, forces and torques. Additionally, the |
| 1037 |
+ |
energy drift and fluctuations from the GSF electrostatics are better |
| 1038 |
+ |
than a multipolar Ewald sum for finite-sized reciprocal spaces. |
| 1039 |
+ |
Because they use real-space cutoffs with moderate cutoff radii, the |
| 1040 |
+ |
GSF and SP models approach $\mathcal{O}(N)$ scaling as the system size |
| 1041 |
+ |
increases. Additionally, they can be made extremely efficient using |
| 1042 |
+ |
spline interpolations of the radial functions. They require no |
| 1043 |
+ |
Fourier transforms or $k$-space sums, and guarantee the smooth |
| 1044 |
+ |
handling of energies, forces, and torques as multipoles cross the |
| 1045 |
+ |
real-space cutoff boundary. |
| 1046 |
+ |
|
| 1047 |
|
%\bibliographystyle{aip} |
| 1048 |
|
\newpage |
| 1049 |
|
\bibliography{references} |
