| 35 |
|
%\linenumbers\relax % Commence numbering lines |
| 36 |
|
\usepackage{amsmath} |
| 37 |
|
\usepackage{times} |
| 38 |
< |
\usepackage{mathptm} |
| 38 |
> |
\usepackage{mathptmx} |
| 39 |
|
\usepackage{tabularx} |
| 40 |
|
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 41 |
|
\usepackage{url} |
| 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. |
| 73 |
> |
against the multipolar Ewald sum as a reference method. The tests |
| 74 |
> |
were carried out in a variety of condensed-phase environments which |
| 75 |
> |
were designed to test all levels of the multipole-multipole |
| 76 |
> |
interactions. Comparisons of the energy differences between |
| 77 |
> |
configurations, molecular forces, and torques were used to analyze |
| 78 |
> |
how well the real-space models perform relative to the more |
| 79 |
> |
computationally expensive Ewald treatment. We have also investigated the |
| 80 |
> |
energy conservation properties of the new methods in molecular |
| 81 |
> |
dynamics simulations using all of these methods. The SP method shows |
| 82 |
> |
excellent agreement with configurational energy differences, forces, |
| 83 |
> |
and torques, and would be suitable for use in Monte Carlo |
| 84 |
> |
calculations. Of the two new shifted-force methods, the GSF |
| 85 |
> |
approach shows the best agreement with Ewald-derived energies, |
| 86 |
> |
forces, and torques and exhibits energy conservation properties that |
| 87 |
> |
make it an excellent choice for efficiently computing electrostatic |
| 88 |
> |
interactions in molecular dynamics simulations. |
| 89 |
|
\end{abstract} |
| 90 |
|
|
| 91 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 92 |
|
% Classification Scheme. |
| 93 |
< |
\keywords{Electrostatics, Multipoles, Real-space} |
| 93 |
> |
%\keywords{Electrostatics, Multipoles, Real-space} |
| 94 |
|
|
| 95 |
|
\maketitle |
| 96 |
|
|
| 120 |
|
method may require modification to compute interactions for |
| 121 |
|
interfacial molecular systems such as membranes and liquid-vapor |
| 122 |
|
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
| 123 |
< |
To simulate interfacial systems, Parry’s extension of the 3D Ewald sum |
| 123 |
> |
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
| 124 |
|
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
| 125 |
< |
periodicity in the Ewald’s method can also be problematic for |
| 126 |
< |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
| 127 |
< |
methods that were developed to handle two-dimensional (2D) |
| 128 |
< |
electrostatic interactions in interfacial systems have not had similar |
| 129 |
< |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 130 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
| 125 |
> |
periodicity in the Ewald method can also be problematic for molecular |
| 126 |
> |
interfaces.\cite{Fennell:2006lq} Modified Ewald methods that were |
| 127 |
> |
developed to handle two-dimensional (2D) electrostatic interactions in |
| 128 |
> |
interfacial systems have not seen similar particle-mesh |
| 129 |
> |
treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
| 130 |
> |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
| 131 |
> |
with system size. |
| 132 |
|
|
| 133 |
|
\subsection{Real-space methods} |
| 134 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
| 135 |
|
method for calculating electrostatic interactions between point |
| 136 |
|
charges. They argued that the effective Coulomb interaction in |
| 137 |
< |
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
| 138 |
< |
an ordered lattice (e.g. when computing the Madelung constant of an |
| 139 |
< |
ionic solid), the material can be considered as a set of ions |
| 137 |
> |
condensed phase systems is actually short ranged.\cite{Wolf92,Wolf95} |
| 138 |
> |
For an ordered lattice (e.g., when computing the Madelung constant of |
| 139 |
> |
an ionic solid), the material can be considered as a set of ions |
| 140 |
|
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
| 141 |
|
an effective distance dependence for the electrostatic interactions of |
| 142 |
< |
$r^{-5}$ (see figure \ref{fig:NaCl}. For this reason, careful |
| 143 |
< |
applications of Wolf's method are able to obtain accurate estimates of |
| 144 |
< |
Madelung constants using relatively short cutoff radii. Recently, |
| 145 |
< |
Fukuda used neutralization of the higher order moments for the |
| 146 |
< |
calculation of the electrostatic interaction of the point charges |
| 147 |
< |
system.\cite{Fukuda:2013sf} |
| 142 |
> |
$r^{-5}$ (see figure \ref{fig:schematic}). For this reason, careful |
| 143 |
> |
application of Wolf's method can obtain accurate estimates of Madelung |
| 144 |
> |
constants using relatively short cutoff radii. Recently, Fukuda used |
| 145 |
> |
neutralization of the higher order moments for calculation of the |
| 146 |
> |
electrostatic interactions in point charge |
| 147 |
> |
systems.\cite{Fukuda:2013sf} |
| 148 |
|
|
| 149 |
< |
\begin{figure}[h!] |
| 149 |
> |
\begin{figure} |
| 150 |
|
\centering |
| 151 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 152 |
< |
\caption{Top: NaCl crystal showing how spherical truncation can |
| 153 |
< |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
| 154 |
< |
molecules interact with the central ion. Bottom: A dipolar |
| 155 |
< |
crystal exhibiting similar behavior and illustrating how the |
| 156 |
< |
effective dipole-octupole interactions can be disrupted by |
| 157 |
< |
spherical truncation.} |
| 158 |
< |
\label{fig:NaCl} |
| 151 |
> |
\includegraphics[width=\linewidth]{schematic.pdf} |
| 152 |
> |
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
| 153 |
> |
charges (in the smaller grey circle), so interactions are |
| 154 |
> |
effectively charge-multipole at longer distances. With hard |
| 155 |
> |
cutoffs, motion of individual charges in and out of the cutoff |
| 156 |
> |
sphere can break the effective multipolar ordering. Bottom: |
| 157 |
> |
dipolar crystals and fluids have a similar effective |
| 158 |
> |
\textit{quadrupolar} ordering (in the smaller grey circles), and |
| 159 |
> |
orientational averaging helps to reduce the effective range of the |
| 160 |
> |
interactions in the fluid. Placement of reversed image multipoles |
| 161 |
> |
on the surface of the cutoff sphere recovers the effective |
| 162 |
> |
higher-order multipole behavior.} |
| 163 |
> |
\label{fig:schematic} |
| 164 |
|
\end{figure} |
| 165 |
|
|
| 166 |
|
The direct truncation of interactions at a cutoff radius creates |
| 167 |
< |
truncation defects. Wolf \textit{et al.} further argued that |
| 167 |
> |
truncation defects. Wolf \textit{et al.} argued that |
| 168 |
|
truncation errors are due to net charge remaining inside the cutoff |
| 169 |
|
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
| 170 |
|
placing an image charge on the surface of the cutoff sphere for every |
| 184 |
|
|
| 185 |
|
Considering the interaction of one central ion in an ionic crystal |
| 186 |
|
with a portion of the crystal at some distance, the effective Columbic |
| 187 |
< |
potential is found to be decreasing as $r^{-5}$. If one views the |
| 188 |
< |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
| 187 |
> |
potential is found to decrease as $r^{-5}$. If one views the \ce{NaCl} |
| 188 |
> |
crystal as a simple cubic (SC) structure with an octupolar |
| 189 |
|
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
| 190 |
|
rapidly to the Madelung energy than the dipolar |
| 191 |
|
approximation.\cite{Wolf92} To find the correct Madelung constant, |
| 192 |
|
Lacman suggested that the NaCl structure could be constructed in a way |
| 193 |
|
that the finite crystal terminates with complete \ce{(NaCl)4} |
| 194 |
< |
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded |
| 195 |
< |
by opposite charges. Similarly for each pair of charges, there is an |
| 196 |
< |
opposite pair of charge adjacent to it. The central ion sees what is |
| 197 |
< |
effectively a set of octupoles at large distances. These facts suggest |
| 192 |
< |
that the Madelung constants are relatively short ranged for perfect |
| 193 |
< |
ionic crystals.\cite{Wolf:1999dn} |
| 194 |
> |
molecules.\cite{Lacman65} The central ion sees what is effectively a |
| 195 |
> |
set of octupoles at large distances. These facts suggest that the |
| 196 |
> |
Madelung constants are relatively short ranged for perfect ionic |
| 197 |
> |
crystals.\cite{Wolf:1999dn} |
| 198 |
|
|
| 199 |
|
One can make a similar argument for crystals of point multipoles. The |
| 200 |
|
Luttinger and Tisza treatment of energy constants for dipolar lattices |
| 212 |
|
unstable. |
| 213 |
|
|
| 214 |
|
In ionic crystals, real-space truncation can break the effective |
| 215 |
< |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
| 216 |
< |
swings in the electrostatic energy as the cutoff radius is increased |
| 217 |
< |
(or as individual ions move back and forth across the boundary). This |
| 218 |
< |
is why the image charges were necessary for the Wolf sum to exhibit |
| 219 |
< |
rapid convergence. Similarly, the real-space truncation of point |
| 220 |
< |
multipole interactions breaks higher order multipole arrangements, and |
| 221 |
< |
image multipoles are required for real-space treatments of |
| 218 |
< |
electrostatic energies. |
| 215 |
> |
multipolar arrangements (see Fig. \ref{fig:schematic}), causing |
| 216 |
> |
significant swings in the electrostatic energy as individual ions move |
| 217 |
> |
back and forth across the boundary. This is why the image charges are |
| 218 |
> |
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
| 219 |
> |
the real-space truncation of point multipole interactions breaks |
| 220 |
> |
higher order multipole arrangements, and image multipoles are required |
| 221 |
> |
for real-space treatments of electrostatic energies. |
| 222 |
|
|
| 223 |
+ |
The shorter effective range of electrostatic interactions is not |
| 224 |
+ |
limited to perfect crystals, but can also apply in disordered fluids. |
| 225 |
+ |
Even at elevated temperatures, there is, on average, local charge |
| 226 |
+ |
balance in an ionic liquid, where each positive ion has surroundings |
| 227 |
+ |
dominated by negaitve ions and vice versa. The reversed-charge images |
| 228 |
+ |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
| 229 |
+ |
retain the effective multipolar interactions as the charges traverse |
| 230 |
+ |
the cutoff boundary. |
| 231 |
+ |
|
| 232 |
+ |
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
| 233 |
+ |
significant orientational averaging that additionally reduces the |
| 234 |
+ |
effect of long-range multipolar interactions. The image multipoles |
| 235 |
+ |
that are introduced in the TSF, GSF, and SP methods mimic this effect |
| 236 |
+ |
and reduce the effective range of the multipolar interactions as |
| 237 |
+ |
interacting molecules traverse each other's cutoff boundaries. |
| 238 |
+ |
|
| 239 |
|
% Because of this reason, although the nature of electrostatic |
| 240 |
|
% interaction short ranged, the hard cutoff sphere creates very large |
| 241 |
|
% fluctuation in the electrostatic energy for the perfect crystal. In |
| 292 |
|
rough approximation. Atomic sites can also be assigned point |
| 293 |
|
multipoles and polarizabilities to increase the accuracy of the |
| 294 |
|
molecular model. Recently, water has been modeled with point |
| 295 |
< |
multipoles up to octupolar |
| 296 |
< |
order.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 295 |
> |
multipoles up to octupolar order using the soft sticky |
| 296 |
> |
dipole-quadrupole-octupole (SSDQO) |
| 297 |
> |
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
| 298 |
|
multipoles up to quadrupolar order have also been coupled with point |
| 299 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
| 300 |
< |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk}. But |
| 300 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But |
| 301 |
|
using point multipole with the real space truncation without |
| 302 |
|
accounting for multipolar neutrality will create energy conservation |
| 303 |
|
issues in molecular dynamics (MD) simulations. |
| 337 |
|
\begin{equation} |
| 338 |
|
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
| 339 |
|
\end{equation} |
| 340 |
< |
where the multipole operator for site $\bf a$, |
| 341 |
< |
\begin{equation} |
| 342 |
< |
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} |
| 343 |
< |
+ Q_{{\bf a}\alpha\beta} |
| 324 |
< |
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
| 325 |
< |
\end{equation} |
| 326 |
< |
is expressed in terms of the point charge, $C_{\bf a}$, dipole, |
| 327 |
< |
$D_{{\bf a}\alpha}$, and quadrupole, $Q_{{\bf a}\alpha\beta}$, for |
| 328 |
< |
object $\bf a$. Note that in this work, we use the primitive |
| 329 |
< |
quadrupole tensor, $Q_{a\alpha,\beta}=\frac{1}{2}\sum_{k\;in\;a}q_k |
| 330 |
< |
r_{k\alpha}r_{k\beta}$ to represent point quadrupoles on a site. |
| 340 |
> |
where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is |
| 341 |
> |
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
| 342 |
> |
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
| 343 |
> |
$\bf a$, etc. |
| 344 |
|
|
| 345 |
< |
Interactions between multipoles can be expressed as higher derivatives |
| 346 |
< |
of the bare Coulomb potential, so one way of ensuring that the forces |
| 347 |
< |
and torques vanish at the cutoff distance is to include a larger |
| 348 |
< |
number of terms in the truncated Taylor expansion, e.g., |
| 349 |
< |
% |
| 350 |
< |
\begin{equation} |
| 351 |
< |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
| 352 |
< |
\end{equation} |
| 353 |
< |
% |
| 354 |
< |
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 355 |
< |
Thus, for $f(r)=1/r$, we find |
| 356 |
< |
% |
| 357 |
< |
\begin{equation} |
| 358 |
< |
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} . |
| 359 |
< |
\end{equation} |
| 360 |
< |
This function is an approximate electrostatic potential that has |
| 361 |
< |
vanishing second derivatives at the cutoff radius, making it suitable |
| 362 |
< |
for shifting the forces and torques of charge-dipole interactions. |
| 345 |
> |
% Interactions between multipoles can be expressed as higher derivatives |
| 346 |
> |
% of the bare Coulomb potential, so one way of ensuring that the forces |
| 347 |
> |
% and torques vanish at the cutoff distance is to include a larger |
| 348 |
> |
% number of terms in the truncated Taylor expansion, e.g., |
| 349 |
> |
% % |
| 350 |
> |
% \begin{equation} |
| 351 |
> |
% f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert _{r_c} . |
| 352 |
> |
% \end{equation} |
| 353 |
> |
% % |
| 354 |
> |
% The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
| 355 |
> |
% Thus, for $f(r)=1/r$, we find |
| 356 |
> |
% % |
| 357 |
> |
% \begin{equation} |
| 358 |
> |
% 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} . |
| 359 |
> |
% \end{equation} |
| 360 |
> |
% This function is an approximate electrostatic potential that has |
| 361 |
> |
% vanishing second derivatives at the cutoff radius, making it suitable |
| 362 |
> |
% for shifting the forces and torques of charge-dipole interactions. |
| 363 |
|
|
| 364 |
< |
In general, the TSF potential for any multipole-multipole interaction |
| 365 |
< |
can be written |
| 364 |
> |
The TSF potential for any multipole-multipole interaction can be |
| 365 |
> |
written |
| 366 |
|
\begin{equation} |
| 367 |
|
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
| 368 |
|
\label{generic} |
| 369 |
|
\end{equation} |
| 370 |
< |
with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for |
| 371 |
< |
charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and |
| 372 |
< |
$n=4$ for quadrupole-quadrupole. To ensure smooth convergence of the |
| 373 |
< |
energy, force, and torques, the required number of terms from Taylor |
| 374 |
< |
series expansion in $f_n(r)$ must be performed for different |
| 375 |
< |
multipole-multipole interactions. |
| 370 |
> |
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
| 371 |
> |
of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for |
| 372 |
> |
charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole |
| 373 |
> |
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for |
| 374 |
> |
quadrupole-quadrupole. To ensure smooth convergence of the energy, |
| 375 |
> |
force, and torques, a Taylor expansion with $n$ terms must be |
| 376 |
> |
performed at cutoff radius ($r_c$) to obtain $f_n(r)$. |
| 377 |
|
|
| 378 |
< |
To carry out the same procedure for a damped electrostatic kernel, we |
| 379 |
< |
replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
| 380 |
< |
Many of the derivatives of the damped kernel are well known from |
| 381 |
< |
Smith's early work on multipoles for the Ewald |
| 382 |
< |
summation.\cite{Smith82,Smith98} |
| 378 |
> |
% To carry out the same procedure for a damped electrostatic kernel, we |
| 379 |
> |
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
| 380 |
> |
% Many of the derivatives of the damped kernel are well known from |
| 381 |
> |
% Smith's early work on multipoles for the Ewald |
| 382 |
> |
% summation.\cite{Smith82,Smith98} |
| 383 |
|
|
| 384 |
< |
Note that increasing the value of $n$ will add additional terms to the |
| 385 |
< |
electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 386 |
< |
$(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 387 |
< |
functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 388 |
< |
f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 389 |
< |
for computing multipole energies, forces, and torques, and smooth |
| 390 |
< |
cutoffs of these quantities can be guaranteed as long as the number of |
| 391 |
< |
terms in the Taylor series exceeds the derivative order required. |
| 384 |
> |
% Note that increasing the value of $n$ will add additional terms to the |
| 385 |
> |
% electrostatic potential, e.g., $f_2(r)$ includes orders up to |
| 386 |
> |
% $(r-r_c)^3/r_c^4$, and so on. Successive derivatives of the $f_n(r)$ |
| 387 |
> |
% functions are denoted $g_2(r) = f^\prime_2(r)$, $h_2(r) = |
| 388 |
> |
% f^{\prime\prime}_2(r)$, etc. These higher derivatives are required |
| 389 |
> |
% for computing multipole energies, forces, and torques, and smooth |
| 390 |
> |
% cutoffs of these quantities can be guaranteed as long as the number of |
| 391 |
> |
% terms in the Taylor series exceeds the derivative order required. |
| 392 |
|
|
| 393 |
|
For multipole-multipole interactions, following this procedure results |
| 394 |
< |
in separate radial functions for each distinct orientational |
| 395 |
< |
contribution to the potential, and ensures that the forces and torques |
| 396 |
< |
from {\it each} of these contributions will vanish at the cutoff |
| 397 |
< |
radius. For example, the direct dipole dot product ($\mathbf{D}_{i} |
| 398 |
< |
\cdot \mathbf{D}_{j}$) is treated differently than the dipole-distance |
| 394 |
> |
in separate radial functions for each of the distinct orientational |
| 395 |
> |
contributions to the potential, and ensures that the forces and |
| 396 |
> |
torques from each of these contributions will vanish at the cutoff |
| 397 |
> |
radius. For example, the direct dipole dot product |
| 398 |
> |
($\mathbf{D}_{\bf a} |
| 399 |
> |
\cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance |
| 400 |
|
dot products: |
| 401 |
|
\begin{equation} |
| 402 |
< |
U_{D_{i}D_{j}}(r)= -\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{i} \cdot |
| 403 |
< |
\mathbf{D}_{j} \right) \frac{g_2(r)}{r} |
| 404 |
< |
-\frac{1}{4\pi \epsilon_0} |
| 405 |
< |
\left( \mathbf{D}_{i} \cdot \hat{r} \right) |
| 406 |
< |
\left( \mathbf{D}_{j} \cdot \hat{r} \right) \left(h_2(r) - |
| 392 |
< |
\frac{g_2(r)}{r} \right) |
| 402 |
> |
U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
| 403 |
> |
\mathbf{D}_{\bf a} \cdot |
| 404 |
> |
\mathbf{D}_{\bf b} \right) v_{21}(r) + |
| 405 |
> |
\left( \mathbf{D}_{\bf a} \cdot \hat{r} \right) |
| 406 |
> |
\left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right] |
| 407 |
|
\end{equation} |
| 408 |
|
|
| 409 |
< |
The electrostatic forces and torques acting on the central multipole |
| 410 |
< |
site due to another site within cutoff sphere are derived from |
| 409 |
> |
For the Taylor shifted (TSF) method with the undamped kernel, |
| 410 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} + |
| 411 |
> |
\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4} |
| 412 |
> |
- \frac{6}{r r_c^2}$. In these functions, one can easily see the |
| 413 |
> |
connection to unmodified electrostatics as well as the smooth |
| 414 |
> |
transition to zero in both these functions as $r\rightarrow r_c$. The |
| 415 |
> |
electrostatic forces and torques acting on the central multipole due |
| 416 |
> |
to another site within the cutoff sphere are derived from |
| 417 |
|
Eq.~\ref{generic}, accounting for the appropriate number of |
| 418 |
|
derivatives. Complete energy, force, and torque expressions are |
| 419 |
|
presented in the first paper in this series (Reference |
| 421 |
|
|
| 422 |
|
\subsection{Gradient-shifted force (GSF)} |
| 423 |
|
|
| 424 |
< |
A second (and significantly simpler) method involves shifting the |
| 425 |
< |
gradient of the raw coulomb potential for each particular multipole |
| 424 |
> |
A second (and conceptually simpler) method involves shifting the |
| 425 |
> |
gradient of the raw Coulomb potential for each particular multipole |
| 426 |
|
order. For example, the raw dipole-dipole potential energy may be |
| 427 |
|
shifted smoothly by finding the gradient for two interacting dipoles |
| 428 |
|
which have been projected onto the surface of the cutoff sphere |
| 429 |
|
without changing their relative orientation, |
| 430 |
< |
\begin{displaymath} |
| 431 |
< |
U_{D_{i}D_{j}}(r_{ij}) = U_{D_{i}D_{j}}(r_{ij}) - U_{D_{i}D_{j}}(r_c) |
| 432 |
< |
- (r_{ij}-r_c) \hat{r}_{ij} \cdot |
| 433 |
< |
\vec{\nabla} V_{D_{i}D_{j}}(r) \Big \lvert _{r_c} |
| 434 |
< |
\end{displaymath} |
| 435 |
< |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{i}$ |
| 436 |
< |
and $\mathbf{D}_{j}$, are retained at the cutoff distance (although |
| 437 |
< |
the signs are reversed for the dipole that has been projected onto the |
| 438 |
< |
cutoff sphere). In many ways, this simpler approach is closer in |
| 439 |
< |
spirit to the original shifted force method, in that it projects a |
| 440 |
< |
neutralizing multipole (and the resulting forces from this multipole) |
| 441 |
< |
onto a cutoff sphere. The resulting functional forms for the |
| 442 |
< |
potentials, forces, and torques turn out to be quite similar in form |
| 443 |
< |
to the Taylor-shifted approach, although the radial contributions are |
| 444 |
< |
significantly less perturbed by the Gradient-shifted approach than |
| 445 |
< |
they are in the Taylor-shifted method. |
| 430 |
> |
\begin{equation} |
| 431 |
> |
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
| 432 |
> |
U_{D_{\bf a} D_{\bf b}}(r_c) |
| 433 |
> |
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
| 434 |
> |
\nabla U_{D_{\bf a}D_{\bf b}}(r_c). |
| 435 |
> |
\end{equation} |
| 436 |
> |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf |
| 437 |
> |
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
| 438 |
> |
(although the signs are reversed for the dipole that has been |
| 439 |
> |
projected onto the cutoff sphere). In many ways, this simpler |
| 440 |
> |
approach is closer in spirit to the original shifted force method, in |
| 441 |
> |
that it projects a neutralizing multipole (and the resulting forces |
| 442 |
> |
from this multipole) onto a cutoff sphere. The resulting functional |
| 443 |
> |
forms for the potentials, forces, and torques turn out to be quite |
| 444 |
> |
similar in form to the Taylor-shifted approach, although the radial |
| 445 |
> |
contributions are significantly less perturbed by the gradient-shifted |
| 446 |
> |
approach than they are in the Taylor-shifted method. |
| 447 |
|
|
| 448 |
+ |
For the gradient shifted (GSF) method with the undamped kernel, |
| 449 |
+ |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
| 450 |
+ |
$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. |
| 451 |
+ |
Again, these functions go smoothly to zero as $r\rightarrow r_c$, and |
| 452 |
+ |
because the Taylor expansion retains only one term, they are |
| 453 |
+ |
significantly less perturbed than the TSF functions. |
| 454 |
+ |
|
| 455 |
|
In general, the gradient shifted potential between a central multipole |
| 456 |
|
and any multipolar site inside the cutoff radius is given by, |
| 457 |
|
\begin{equation} |
| 458 |
< |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 459 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
| 460 |
< |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
| 458 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 459 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} |
| 460 |
> |
\cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 461 |
|
\label{generic2} |
| 462 |
|
\end{equation} |
| 463 |
|
where the sum describes a separate force-shifting that is applied to |
| 464 |
< |
each orientational contribution to the energy. |
| 464 |
> |
each orientational contribution to the energy. In this expression, |
| 465 |
> |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
| 466 |
> |
($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ |
| 467 |
> |
represent the orientations the multipoles. |
| 468 |
|
|
| 469 |
|
The third term converges more rapidly than the first two terms as a |
| 470 |
|
function of radius, hence the contribution of the third term is very |
| 471 |
|
small for large cutoff radii. The force and torque derived from |
| 472 |
< |
equation \ref{generic2} are consistent with the energy expression and |
| 472 |
> |
Eq. \ref{generic2} are consistent with the energy expression and |
| 473 |
|
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
| 474 |
|
can be considered generalizations of the original DSF method for |
| 475 |
|
higher order multipole interactions. GSF and TSF are also identical up |
| 477 |
|
the energy, force and torque for higher order multipole-multipole |
| 478 |
|
interactions. Complete energy, force, and torque expressions for the |
| 479 |
|
GSF potential are presented in the first paper in this series |
| 480 |
< |
(Reference~\onlinecite{PaperI}) |
| 480 |
> |
(Reference~\onlinecite{PaperI}). |
| 481 |
|
|
| 482 |
|
|
| 483 |
|
\subsection{Shifted potential (SP) } |
| 490 |
|
interactions with the central multipole and the image. This |
| 491 |
|
effectively shifts the total potential to zero at the cutoff radius, |
| 492 |
|
\begin{equation} |
| 493 |
< |
U_{SP}(\vec r)=\sum U(\vec r) - U(\vec r_c) |
| 493 |
> |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 494 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 495 |
|
\label{eq:SP} |
| 496 |
|
\end{equation} |
| 497 |
|
where the sum describes separate potential shifting that is done for |
| 503 |
|
multipoles that reorient after leaving the cutoff sphere can re-enter |
| 504 |
|
the cutoff sphere without perturbing the total energy. |
| 505 |
|
|
| 506 |
< |
The potential energy between a central multipole and other multipolar |
| 507 |
< |
sites then goes smoothly to zero as $r \rightarrow r_c$. However, the |
| 508 |
< |
force and torque obtained from the shifted potential (SP) are |
| 509 |
< |
discontinuous at $r_c$. Therefore, MD simulations will still |
| 510 |
< |
experience energy drift while operating under the SP potential, but it |
| 511 |
< |
may be suitable for Monte Carlo approaches where the configurational |
| 512 |
< |
energy differences are the primary quantity of interest. |
| 506 |
> |
For the shifted potential (SP) method with the undamped kernel, |
| 507 |
> |
$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) = |
| 508 |
> |
\frac{3}{r^3} - \frac{3}{r_c^3}$. The potential energy between a |
| 509 |
> |
central multipole and other multipolar sites goes smoothly to zero as |
| 510 |
> |
$r \rightarrow r_c$. However, the force and torque obtained from the |
| 511 |
> |
shifted potential (SP) are discontinuous at $r_c$. MD simulations |
| 512 |
> |
will still experience energy drift while operating under the SP |
| 513 |
> |
potential, but it may be suitable for Monte Carlo approaches where the |
| 514 |
> |
configurational energy differences are the primary quantity of |
| 515 |
> |
interest. |
| 516 |
|
|
| 517 |
< |
\subsection{The Self term} |
| 517 |
> |
\subsection{The Self Term} |
| 518 |
|
In the TSF, GSF, and SP methods, a self-interaction is retained for |
| 519 |
|
the central multipole interacting with its own image on the surface of |
| 520 |
|
the cutoff sphere. This self interaction is nearly identical with the |
| 536 |
|
used the multipolar Ewald sum as a reference method for comparing |
| 537 |
|
energies, forces, and torques for molecular models that mimic |
| 538 |
|
disordered and ordered condensed-phase systems. The parameters used |
| 539 |
< |
in the test-cases are given in table~\ref{tab:pars}. |
| 539 |
> |
in the test cases are given in table~\ref{tab:pars}. |
| 540 |
|
|
| 541 |
|
\begin{table} |
| 542 |
|
\label{tab:pars} |
| 554 |
|
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
| 555 |
|
\AA\textsuperscript{2})} \\ \hline |
| 556 |
|
Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ |
| 557 |
< |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & 10,000 & 17.6 &17.6 & 0 \\ |
| 557 |
> |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ |
| 558 |
|
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 559 |
< |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & 10,000 & 17.6&17.6&0 \\ |
| 559 |
> |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ |
| 560 |
|
SSDQ water & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
| 561 |
|
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
| 562 |
|
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
| 581 |
|
electrostatic energy, as well as the electrostatic contributions to |
| 582 |
|
the force and torque on each molecule. These quantities have been |
| 583 |
|
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
| 584 |
< |
and have been compared with the values obtaine from the multipolar |
| 585 |
< |
Ewald sum. In Mote Carlo (MC) simulations, the energy differences |
| 584 |
> |
and have been compared with the values obtained from the multipolar |
| 585 |
> |
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
| 586 |
|
between two configurations is the primary quantity that governs how |
| 587 |
|
the simulation proceeds. These differences are the most imporant |
| 588 |
|
indicators of the reliability of a method even if the absolute |
| 629 |
|
recomputed at each time step. |
| 630 |
|
|
| 631 |
|
\subsection{Model systems} |
| 632 |
< |
To sample independent configurations of multipolar crystals, a body |
| 633 |
< |
centered cubic (bcc) crystal which is a minimum energy structure for |
| 634 |
< |
point dipoles was generated using 3,456 molecules. The multipoles |
| 635 |
< |
were translationally locked in their respective crystal sites for |
| 636 |
< |
equilibration at a relatively low temperature (50K), so that dipoles |
| 637 |
< |
or quadrupoles could freely explore all accessible orientations. The |
| 638 |
< |
translational constraints were removed, and the crystals were |
| 639 |
< |
simulated for 10 ps in the microcanonical (NVE) ensemble with an |
| 640 |
< |
average temperature of 50 K. Configurations were sampled at equal |
| 641 |
< |
time intervals for the comparison of the configurational energy |
| 642 |
< |
differences. The crystals were not simulated close to the melting |
| 643 |
< |
points in order to avoid translational deformation away of the ideal |
| 644 |
< |
lattice geometry. |
| 645 |
< |
|
| 646 |
< |
For dipolar, quadrupolar, and mixed-multipole liquid simulations, each |
| 647 |
< |
system was created with 2048 molecules oriented randomly. These were |
| 632 |
> |
To sample independent configurations of the multipolar crystals, body |
| 633 |
> |
centered cubic (bcc) crystals, which exhibit the minimum energy |
| 634 |
> |
structures for point dipoles, were generated using 3,456 molecules. |
| 635 |
> |
The multipoles were translationally locked in their respective crystal |
| 636 |
> |
sites for equilibration at a relatively low temperature (50K) so that |
| 637 |
> |
dipoles or quadrupoles could freely explore all accessible |
| 638 |
> |
orientations. The translational constraints were then removed, the |
| 639 |
> |
systems were re-equilibrated, and the crystals were simulated for an |
| 640 |
> |
additional 10 ps in the microcanonical (NVE) ensemble with an average |
| 641 |
> |
temperature of 50 K. The balance between moments of inertia and |
| 642 |
> |
particle mass were chosen to allow orientational sampling without |
| 643 |
> |
significant translational motion. Configurations were sampled at |
| 644 |
> |
equal time intervals in order to compare configurational energy |
| 645 |
> |
differences. The crystals were simulated far from the melting point |
| 646 |
> |
in order to avoid translational deformation away of the ideal lattice |
| 647 |
> |
geometry. |
| 648 |
|
|
| 649 |
< |
system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K |
| 650 |
< |
temperature after equilibration. We collected 250 different |
| 651 |
< |
configurations in equal interval of time. For the ions mixed liquid |
| 652 |
< |
system, we converted 48 different molecules into 24 \ce{Na+} and 24 |
| 653 |
< |
\ce{Cl-} ions and equilibrated. After equilibration, the system was run |
| 654 |
< |
at the same environment for 1ns and 250 configurations were |
| 655 |
< |
collected. While comparing energies, forces, and torques with Ewald |
| 656 |
< |
method, Lennard-Jones potentials were turned off and purely |
| 657 |
< |
electrostatic interaction had been compared. |
| 649 |
> |
For dipolar, quadrupolar, and mixed-multipole \textit{liquid} |
| 650 |
> |
simulations, each system was created with 2,048 randomly-oriented |
| 651 |
> |
molecules. These were equilibrated at a temperature of 300K for 1 ns. |
| 652 |
> |
Each system was then simulated for 1 ns in the microcanonical (NVE) |
| 653 |
> |
ensemble. We collected 250 different configurations at equal time |
| 654 |
> |
intervals. For the liquid system that included ionic species, we |
| 655 |
> |
converted 48 randomly-distributed molecules into 24 \ce{Na+} and 24 |
| 656 |
> |
\ce{Cl-} ions and re-equilibrated. After equilibration, the system was |
| 657 |
> |
run under the same conditions for 1 ns. A total of 250 configurations |
| 658 |
> |
were collected. In the following comparisons of energies, forces, and |
| 659 |
> |
torques, the Lennard-Jones potentials were turned off and only the |
| 660 |
> |
purely electrostatic quantities were compared with the same values |
| 661 |
> |
obtained via the Ewald sum. |
| 662 |
|
|
| 663 |
|
\subsection{Accuracy of Energy Differences, Forces and Torques} |
| 664 |
|
The pairwise summation techniques (outlined above) were evaluated for |
| 672 |
|
should be identical for all methods. |
| 673 |
|
|
| 674 |
|
Since none of the real-space methods provide exact energy differences, |
| 675 |
< |
we used least square regressions analysiss for the six different |
| 675 |
> |
we used least square regressions analysis for the six different |
| 676 |
|
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
| 677 |
|
with the multipolar Ewald reference method. Unitary results for both |
| 678 |
|
the correlation (slope) and correlation coefficient for these |
| 683 |
|
configurations and 250 configurations were recorded for comparison. |
| 684 |
|
Each system provided 31,125 energy differences for a total of 186,750 |
| 685 |
|
data points. Similarly, the magnitudes of the forces and torques have |
| 686 |
< |
also been compared by using least squares regression analyses. In the |
| 686 |
> |
also been compared using least squares regression analysis. In the |
| 687 |
|
forces and torques comparison, the magnitudes of the forces acting in |
| 688 |
|
each molecule for each configuration were evaluated. For example, our |
| 689 |
|
dipolar liquid simulation contains 2048 molecules and there are 250 |
| 868 |
|
molecules inside each other's cutoff spheres in order to correct the |
| 869 |
|
energy conservation issues, and this perturbation is evident in the |
| 870 |
|
statistics accumulated for the molecular forces. The GSF |
| 871 |
< |
perturbations are minimal, particularly for moderate damping and and |
| 871 |
> |
perturbations are minimal, particularly for moderate damping and |
| 872 |
|
commonly-used cutoff values ($r_c = 12$\AA). The TSF method shows |
| 873 |
|
reasonable agreement in the correlation coefficient but again the |
| 874 |
|
systematic error in the forces is concerning if replication of Ewald |
| 925 |
|
these quantities. Force and torque vectors for all six systems were |
| 926 |
|
analyzed using Fisher statistics, and the quality of the vector |
| 927 |
|
directionality is shown in terms of circular variance |
| 928 |
< |
($\mathrm{Var}(\theta$) in figure |
| 928 |
> |
($\mathrm{Var}(\theta)$) in figure |
| 929 |
|
\ref{fig:slopeCorr_circularVariance}. The force and torque vectors |
| 930 |
|
from the new real-space methods exhibit nearly-ideal Fisher probability |
| 931 |
|
distributions (Eq.~\ref{eq:pdf}). Both the hard and SP cutoff methods |
| 980 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
| 981 |
|
fluctuations around this drift $\delta E_0$. Both of the |
| 982 |
|
shifted-force methods (GSF and TSF) provide excellent energy |
| 983 |
< |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
| 983 |
> |
conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), |
| 984 |
|
while the hard cutoff is essentially unusable for molecular dynamics. |
| 985 |
|
SP provides some benefit over the hard cutoff because the energetic |
| 986 |
|
jumps that happen as particles leave and enter the cutoff sphere are |
| 995 |
|
|
| 996 |
|
\begin{figure} |
| 997 |
|
\centering |
| 998 |
< |
\includegraphics[width=\textwidth]{newDrift.pdf} |
| 998 |
> |
\includegraphics[width=\textwidth]{newDrift_12.pdf} |
| 999 |
|
\label{fig:energyDrift} |
| 1000 |
|
\caption{Analysis of the energy conservation of the real-space |
| 1001 |
|
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |
| 1002 |
< |
energy over time and $\delta \mathrm{E}_0$ is the standard deviation |
| 1003 |
< |
of energy fluctuations around this drift. All simulations were of a |
| 1004 |
< |
2000-molecule simulation of SSDQ water with 48 ionic charges at 300 |
| 1005 |
< |
K starting from the same initial configuration. All runs utilized |
| 1006 |
< |
the same real-space cutoff, $r_c = 12$\AA.} |
| 1002 |
> |
energy over time (in kcal / mol / particle / ns) and $\delta |
| 1003 |
> |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
| 1004 |
> |
around this drift (in kcal / mol / particle). All simulations were |
| 1005 |
> |
of a 2000-molecule simulation of SSDQ water with 48 ionic charges at |
| 1006 |
> |
300 K starting from the same initial configuration. All runs |
| 1007 |
> |
utilized the same real-space cutoff, $r_c = 12$\AA.} |
| 1008 |
|
\end{figure} |
| 1009 |
|
|
| 1010 |
|
|
| 1084 |
|
handling of energies, forces, and torques as multipoles cross the |
| 1085 |
|
real-space cutoff boundary. |
| 1086 |
|
|
| 1087 |
+ |
\begin{acknowledgments} |
| 1088 |
+ |
JDG acknowledges helpful discussions with Christopher |
| 1089 |
+ |
Fennell. Support for this project was provided by the National |
| 1090 |
+ |
Science Foundation under grant CHE-1362211. Computational time was |
| 1091 |
+ |
provided by the Center for Research Computing (CRC) at the |
| 1092 |
+ |
University of Notre Dame. |
| 1093 |
+ |
\end{acknowledgments} |
| 1094 |
+ |
|
| 1095 |
|
%\bibliographystyle{aip} |
| 1096 |
|
\newpage |
| 1097 |
|
\bibliography{references} |
