1 |
\documentclass[aip,jcp,amsmath,amssymb,preprint]{revtex4-1} |
2 |
|
3 |
\usepackage{graphicx}% Include figure files |
4 |
\usepackage{dcolumn}% Align table columns on decimal point |
5 |
\usepackage{multirow} |
6 |
\usepackage{bm} |
7 |
\usepackage{natbib} |
8 |
\usepackage{times} |
9 |
\usepackage{mathptmx} |
10 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
11 |
\usepackage{url} |
12 |
\usepackage{braket} |
13 |
\usepackage{tabularx} |
14 |
\usepackage{rotating} |
15 |
\usepackage{multirow} |
16 |
\usepackage{booktabs} |
17 |
|
18 |
\newcolumntype{Y}{>{\centering\arraybackslash}X} |
19 |
|
20 |
\begin{document} |
21 |
|
22 |
\title{Real space electrostatics for multipoles. III. Dielectric Properties} |
23 |
|
24 |
\author{Madan Lamichhane} |
25 |
\affiliation{Department of Physics, University |
26 |
of Notre Dame, Notre Dame, IN 46556} |
27 |
\author{Thomas Parsons} |
28 |
\affiliation{Department of Chemistry and Biochemistry, University |
29 |
of Notre Dame, Notre Dame, IN 46556} |
30 |
\author{Kathie E. Newman} |
31 |
\affiliation{Department of Physics, University |
32 |
of Notre Dame, Notre Dame, IN 46556} |
33 |
\author{J. Daniel Gezelter} |
34 |
\email{gezelter@nd.edu.} |
35 |
\affiliation{Department of Chemistry and Biochemistry, University |
36 |
of Notre Dame, Notre Dame, IN 46556} |
37 |
|
38 |
\date{\today}% It is always \today, today, |
39 |
% but any date may be explicitly specified |
40 |
|
41 |
\begin{abstract} |
42 |
In the first two papers in this series, we developed new shifted |
43 |
potential (SP), gradient shifted force (GSF), and Taylor shifted |
44 |
force (TSF) real-space methods for multipole interactions in |
45 |
condensed phase simulations. Here, we discuss the dielectric |
46 |
properties of fluids that emerge from simulations using these |
47 |
methods. Most electrostatic methods (including the Ewald sum) |
48 |
require correction to the conducting boundary fluctuation formula |
49 |
for the static dielectric constants, and we discuss the derivation |
50 |
of these corrections for the new real space methods. For quadrupolar |
51 |
fluids, the analogous material property is the quadrupolar |
52 |
susceptibility. As in the dipolar case, the fluctuation formula for |
53 |
the quadrupolar susceptibility has corrections that depend on the |
54 |
electrostatic method being utillized. One of the most important |
55 |
effects measured by both the static dielectric and quadrupolar |
56 |
susceptibility is the ability to screen charges embedded in the |
57 |
fluid. We use potentials of mean force between solvated ions to |
58 |
discuss how geometric factors can lead to distance-dependent |
59 |
screening in both quadrupolar and dipolar fluids. |
60 |
\end{abstract} |
61 |
|
62 |
\maketitle |
63 |
|
64 |
\section{Introduction} |
65 |
|
66 |
Over the past several years, there has been increasing interest in |
67 |
pairwise or ``real space'' methods for computing electrostatic |
68 |
interactions in condensed phase |
69 |
simulations.\cite{Wolf99,Zahn02,Kast03,Beckd05,Ma05,Wu:2005nr,Fennell06,Fukuda:2013uq,Wang:2016kx} |
70 |
These techniques were initially developed by Wolf {\it et al.} in |
71 |
their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99} |
72 |
Wolf's method of using cutoff neutralization and electrostatic damping |
73 |
is able to obtain excellent agreement with Madelung energies in ionic |
74 |
crystals.\cite{Wolf99} |
75 |
|
76 |
Zahn \textit{et al.}\cite{Zahn02} and Fennell and Gezelter extended |
77 |
this method using shifted force approximations at the cutoff distance |
78 |
in order to conserve total energy in molecular dynamics |
79 |
simulations.\cite{Fennell06} Other recent advances in real-space |
80 |
methods for systems of point charges have included explicit |
81 |
elimination of the net multipole moments inside the cutoff sphere |
82 |
around each charge site.\cite{Fukuda:2013uq,Wang:2016kx} |
83 |
|
84 |
In the previous two papers in this series, we developed three |
85 |
generalized real space methods: shifted potential (SP), gradient |
86 |
shifted force (GSF), and Taylor shifted force (TSF).\cite{PaperI, |
87 |
PaperII} These methods evaluate electrostatic interactions for |
88 |
charges and higher order multipoles using a finite-radius cutoff |
89 |
sphere. The neutralization and damping of local moments within the |
90 |
cutoff sphere is a multipolar generalization of Wolf's sum. In the |
91 |
GSF and TSF methods, additional terms are added to the potential |
92 |
energy so that forces and torques also vanish smoothly at the cutoff |
93 |
radius. This ensures that the total energy is conserved in molecular |
94 |
dynamics simulations. |
95 |
|
96 |
One of the most stringent tests of any new electrostatic method is the |
97 |
fidelity with which that method can reproduce the bulk-phase |
98 |
polarizability or equivalently, the dielectric properties of a |
99 |
fluid. Before the advent of computer simulations, Kirkwood and Onsager |
100 |
developed fluctuation formulae for the dielectric properties of |
101 |
dipolar fluids.\cite{Kirkwood39,Onsagar36} Along with projections of |
102 |
the frequency-dependent dielectric to zero frequency, these |
103 |
fluctuation formulae are now widely used to predict the static |
104 |
dielectric constant of simulated materials. |
105 |
|
106 |
If we consider a system of dipolar or quadrupolar molecules under the |
107 |
influence of an external field or field gradient, the net polarization |
108 |
of the system will largely be proportional to the applied |
109 |
perturbation.\cite{Chitanvis96, Stern-Feller03, SalvchovI14, |
110 |
SalvchovII14} In simulations, the net polarization of the system is |
111 |
also determined by the interactions \textit{between} the |
112 |
molecules. Therefore the macroscopic polarizablity obtained from a |
113 |
simulation depends on the details of the electrostatic interaction |
114 |
methods that were employed in the simulation. To determine the |
115 |
relevant physical properties of the multipolar fluid from the system |
116 |
fluctuations, the interactions between molecules must be incorporated |
117 |
into the formalism for the bulk properties. |
118 |
|
119 |
In most simulations, bulk materials are treated using periodic |
120 |
replicas of small regions, and this level of approximation requires |
121 |
corrections to the fluctuation formulae that were derived for the bulk |
122 |
fluids. In 1983 Neumann proposed a general formula for evaluating |
123 |
dielectric properties of dipolar fluids using both Ewald and |
124 |
real-space cutoff methods.\cite{NeumannI83} Steinhauser and Neumann |
125 |
used this formula to evaluate the corrected dielectric constant for |
126 |
the Stockmayer fluid using two different methods: Ewald-Kornfeld (EK) |
127 |
and reaction field (RF) methods.\cite{NeumannII83} |
128 |
|
129 |
Zahn \textit{et al.}\cite{Zahn02} utilized this approach and evaluated |
130 |
the correction factor for using damped shifted charge-charge |
131 |
kernel. This was later generalized by Izvekov \textit{et |
132 |
al.},\cite{Izvekov08} who noted that the expression for the |
133 |
dielectric constant reduces to widely-used conducting boundary formula |
134 |
for real-space methods that have first derivatives that vanish at the |
135 |
cutoff radius. |
136 |
|
137 |
One of the primary topics of this paper is the derivation of |
138 |
correction factors for the three new real space methods. We find that |
139 |
the correction formulae for dipolar molecules depends not only on the |
140 |
methodology being used, but also on whether the molecular dipoles are |
141 |
treated using point charges or point dipoles. We derive correction |
142 |
factors for both cases. |
143 |
|
144 |
In quadrupolar fluids, the relationship between quadrupolar |
145 |
susceptibility and the dielectric constant is not as straightforward |
146 |
as in the dipolar case. The effective dielectric constant depends on |
147 |
the geometry of the external (or internal) field |
148 |
perturbation.\cite{Ernst92} Significant efforts have been made to |
149 |
increase our understanding the dielectric properties of these |
150 |
fluids,\cite{JeonI03,JeonII03,Chitanvis96} although a general |
151 |
correction formula has not yet been developed. |
152 |
|
153 |
In this paper we derive general formulae for calculating the |
154 |
quadrupolar susceptibility of quadrupolar fluids. We also evaluate the |
155 |
correction factor for SP, GSF, and TSF methods for quadrupolar fluids |
156 |
interacting via point charges, point dipoles or directly through |
157 |
quadrupole-quadrupole interactions. |
158 |
|
159 |
We also calculate the screening behavior for two ions immersed in |
160 |
multipolar fluids to estimate the distance dependence of charge |
161 |
screening in both dipolar and quadrupolar fluids. We use three |
162 |
distinct methods to compare our analytical results with computer |
163 |
simulations: |
164 |
\begin{enumerate} |
165 |
\item responses of the fluid to external perturbations, |
166 |
\item fluctuations of system multipole moments, and |
167 |
\item potentials of mean force between solvated ions, |
168 |
\end{enumerate} |
169 |
|
170 |
\begin{figure} |
171 |
\includegraphics[width=\linewidth]{Schematic} |
172 |
\caption{Dielectric properties of a fluid measure the response to |
173 |
external electric fields and gradients (left panel), or internal |
174 |
fields and gradients generated by the molecules themselves (central |
175 |
panel), or fields produced by embedded ions (right panel). The |
176 |
dielectric constant measures all three responses in dipolar fluids. |
177 |
In quadrupolar liquids, the relevant bulk property is the |
178 |
quadrupolar susceptibility, and the geometry of the perturbation |
179 |
determines the effective dielectric screening.} |
180 |
\label{fig:schematic} |
181 |
\end{figure} |
182 |
|
183 |
Under the influence of weak external fields, the bulk polarization of |
184 |
the system is primarily a linear response to the perturbation, where |
185 |
proportionality constant depends on the electrostatic interactions |
186 |
between the multipoles. The fluctuation formulae connect bulk |
187 |
properties of the fluid to equilibrium fluctuations in the system |
188 |
multipolar moments during a simulation. These fluctuations also depend |
189 |
on the form of the electrostatic interactions between molecules. |
190 |
Therefore, the connections between the actual bulk properties and both |
191 |
the computed fluctuation and external field responses must be modified |
192 |
accordingly. |
193 |
|
194 |
The potential of mean force (PMF) allows calculation of an effective |
195 |
dielectric constant or screening factor from the potential energy |
196 |
between ions before and after dielectric material is introduced. |
197 |
Computing the PMF between embedded point charges is an additional |
198 |
check on the bulk properties computed via the other two methods. |
199 |
|
200 |
\section{The Real-Space Methods} |
201 |
|
202 |
In the first paper in this series, we derived interaction energies, as |
203 |
well as expressions for the forces and torques for point multipoles |
204 |
interacting via three new real-space methods.\cite{PaperI} The Taylor |
205 |
shifted-force (TSF) method modifies the electrostatic kernel, |
206 |
$f(r) = 1/r$, so that all forces and torques go smoothly to zero at |
207 |
the cutoff radius, |
208 |
\begin{equation} |
209 |
U^{\text{TSF}}_{ab} = M_{a} M_{b} f_n(r). |
210 |
\label{TSF} |
211 |
\end{equation} |
212 |
Here the multipole operator for site $a$, $M_{a}$, is expressed in |
213 |
terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and |
214 |
quadrupole, $\mathsf{Q}_{a}$, for object $a$, etc. Because each of |
215 |
the multipole operators includes a set of implied gradient operators, |
216 |
an approximate electrostatic kernel, $f_n(r)$ is Taylor-expanded |
217 |
around the cutoff radius, so that $n + 1$ derivatives vanish as |
218 |
$r \rightarrow r_c$. This ensures smooth convergence of the energy, |
219 |
forces, and torques as molecules leave and reenter each others cutoff |
220 |
spheres. The order of the Taylor expansion is determined by the |
221 |
multipolar order of the interaction. That is, smooth |
222 |
quadrupole-quadrupole forces require the fifth derivative to vanish at |
223 |
the cutoff radius, so the appropriate function Taylor expansion will |
224 |
be of fifth order. |
225 |
|
226 |
Following this procedure results in separate radial functions for each |
227 |
of the distinct orientational contributions to the potential. For |
228 |
example, in dipole-dipole interactions, the direct dipole dot product |
229 |
($\mathbf{D}_{a} \cdot \mathbf{D}_{b}$) is treated differently than |
230 |
the dipole-distance dot products: |
231 |
\begin{equation} |
232 |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
233 |
\mathbf{D}_{a} \cdot |
234 |
\mathbf{D}_{b} \right) v_{21}(r) + |
235 |
\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right) |
236 |
\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right] |
237 |
\end{equation} |
238 |
The two radial functions, $v_{21}(r)$ and $v_{22}(r)$ are both $1/r^3$ |
239 |
in standard electrostatics, but they become distinct in the TSF |
240 |
method, the relative orientations for interacting molecules so that |
241 |
forces and torques vanish smoothly as the molecules drift beyond the |
242 |
cutoff radius. |
243 |
|
244 |
A second and somewhat simpler approach involves shifting the gradient |
245 |
of the Coulomb potential for each particular multipole order, |
246 |
\begin{equation} |
247 |
U^{\text{GSF}}_{ab} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
248 |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c) |
249 |
\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, |
250 |
\mathsf{B}) \right] |
251 |
\label{generic2} |
252 |
\end{equation} |
253 |
where the sum describes a separate force-shifting that is applied to |
254 |
each orientational contribution to the energy, i.e. $v_{21}$ and |
255 |
$v_{22}$ are shifted separately. In this expression, |
256 |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
257 |
($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$ represent |
258 |
the orientations the multipoles. Because this procedure is equivalent |
259 |
to using the gradient of an image multipole placed at the cutoff |
260 |
sphere for shifting the force, this method is called the gradient |
261 |
shifted-force (GSF) approach. |
262 |
|
263 |
Both the TSF and GSF approaches can be thought of as multipolar |
264 |
extensions of the original damped shifted-force (DSF) approach that |
265 |
was developed for point charges. There is also a multipolar extension |
266 |
of the Wolf sum that is obtained by projecting an image multipole onto |
267 |
the surface of the cutoff sphere, and including the interactions with |
268 |
the central multipole and the image. This effectively shifts only the |
269 |
total potential to zero at the cutoff radius, |
270 |
\begin{equation} |
271 |
U^{\text{SP}}_{ab} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
272 |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
273 |
\label{eq:SP} |
274 |
\end{equation} |
275 |
where the sum again describes separate potential shifting that is done |
276 |
for each orientational contribution to the energy. The potential |
277 |
energy between a central multipole and other multipolar sites goes |
278 |
smoothly to zero as $r \rightarrow r_c$, but the forces and torques |
279 |
obtained from this shifted potential (SP) approach are discontinuous |
280 |
at $r_c$. |
281 |
|
282 |
All three of the new real space methods share a common structure: the |
283 |
various orientational contributions to multipolar interaction energies |
284 |
require separate treatment of their radial functions, and these are |
285 |
tabulated for both the raw Coulombic kernel ($1/r$) as well as the |
286 |
damped kernel ($\mathrm{erfc}(\alpha r)/r$), in the first paper of this |
287 |
series.\cite{PaperI} The second paper in this series evaluated the |
288 |
fidelity with which the three new methods reproduced Ewald-based |
289 |
results for a number of model systems.\cite{PaperII} One of the major |
290 |
findings was that moderately-damped GSF simulations produced nearly |
291 |
identical behavior with Ewald-based simulations, but the real-space |
292 |
methods scale linearly with system size. |
293 |
|
294 |
\section{Dipolar Fluids and the Dielectric Constant} |
295 |
|
296 |
Dielectric properties of a fluid arise mainly from responses of the |
297 |
fluid to either applied fields or transient fields internal to the |
298 |
fluid. In response to an applied field, the molecules have electronic |
299 |
polarizabilities, changes to internal bond lengths, and reorientations |
300 |
towards the direction of the applied field. There is an added |
301 |
complication that in the presence of external field, the perturbation |
302 |
experienced by any single molecule is not only due to external field |
303 |
but also to the fields produced by the all other molecules in the |
304 |
system. |
305 |
|
306 |
\subsection{Response to External Perturbations} |
307 |
|
308 |
In the presence of uniform electric field $\mathbf{E}$, an individual |
309 |
molecule with a permanent dipole moment $p_o$ will realign along the |
310 |
direction of the field with an average polarization given by |
311 |
\begin{equation} |
312 |
\braket{\mathbf{p}} = \epsilon_0 \alpha_p \mathbf{E}, |
313 |
\end{equation} |
314 |
where $\alpha_p = {p_o}^2 / 3 \epsilon_0 k_B T$ is the contribution to |
315 |
molecular polarizability due solely to reorientation dynamics. |
316 |
Because the applied field must overcome thermal motion, the |
317 |
orientational polarization depends inversely on the temperature. |
318 |
|
319 |
Likewise, a condensed phase system of permanent dipoles will also |
320 |
polarize along the direction of an applied field. The polarization |
321 |
density of the system is |
322 |
\begin{equation} |
323 |
\textbf{P} = \epsilon_o \alpha_{D} \mathbf{E}. |
324 |
\label{pertDipole} |
325 |
\end{equation} |
326 |
The constant $\alpha_D$ is the macroscopic polarizability, which is an |
327 |
emergent property of the dipolar fluid. Note that the polarizability, |
328 |
$\alpha_D$ is distinct from the dipole susceptibility, $\chi_D$, |
329 |
which is the quantity most directly related to the static dielectric |
330 |
constant, $\epsilon = 1 + \chi_D$. |
331 |
|
332 |
\subsection{Fluctuation Formula} |
333 |
% Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can |
334 |
% be calculated for systems of orientationally-polarized molecules, |
335 |
% \begin{equation} |
336 |
% \epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
337 |
% \label{eq:staticDielectric} |
338 |
% \end{equation} |
339 |
% where $\epsilon_0$ is the permittivity of free space and |
340 |
% $\langle M^2\rangle$ is the fluctuation of the system dipole |
341 |
% moment.\cite{Allen89} The numerator in the fractional term in |
342 |
% equation (\ref{eq:staticDielectric}) is identical to the quantity |
343 |
% calculated in the finite-system Kirkwood $g$ factor ($G_k$): |
344 |
% \begin{equation} |
345 |
% G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
346 |
% \label{eq:KirkwoodFactor} |
347 |
% \end{equation} |
348 |
% where $\mu$ is the dipole moment of a single molecule of the |
349 |
% homogeneous system.\cite{NeumannI83,NeumannII83,Neumann84,Neumann85} The |
350 |
% fluctuation term in both equation (\ref{eq:staticDielectric}) and |
351 |
% (\ref{eq:KirkwoodFactor}) is calculated as follows, |
352 |
% \begin{equation} |
353 |
% \begin{split} |
354 |
% \langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
355 |
% - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
356 |
% &= \langle M_x^2+M_y^2+M_z^2\rangle |
357 |
% - (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
358 |
% + \langle M_x\rangle^2). |
359 |
% \end{split} |
360 |
% \label{eq:fluctBoxDipole} |
361 |
% \end{equation} |
362 |
% This fluctuation term can be accumulated during a simulation; however, |
363 |
% it converges rather slowly, thus requiring multi-nanosecond simulation |
364 |
% times.\cite{Horn04} In the case of tin-foil boundary conditions, the |
365 |
% dielectric/surface term of the Ewald summation is equal to zero. |
366 |
|
367 |
For a system of dipolar molecules at thermal equilibrium, we can |
368 |
define both a system dipole moment, $\mathbf{M} = \sum_i \mathbf{p}_i$ |
369 |
as well as a dipole polarization density, |
370 |
$\mathbf{P} = \braket{\mathbf{M}}/V$. |
371 |
|
372 |
In the presence of applied field $\mathbf{E}$, the polarization |
373 |
density can be expressed in terms of fluctuations in the net dipole |
374 |
moment, |
375 |
\begin{equation} |
376 |
\mathbf{P} = \epsilon_o \frac{\braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E} |
377 |
\label{flucDipole} |
378 |
\end{equation} |
379 |
This has structural similarity with the Boltzmann average for the |
380 |
polarization of a single molecule. Here |
381 |
$ \braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2$ measures fluctuations |
382 |
in the net dipole moment, |
383 |
\begin{equation} |
384 |
\langle \mathbf{M}^2 \rangle - \langle \mathbf{M} \rangle^2 = |
385 |
\langle M_x^2+M_y^2+M_z^2 \rangle - \left( \langle M_x \rangle^2 + |
386 |
\langle M_y \rangle^2 + |
387 |
\langle M_z \rangle^2 \right). |
388 |
\label{eq:flucDip} |
389 |
\end{equation} |
390 |
For the limiting case $\textbf{E} \rightarrow 0 $, the ensemble |
391 |
average of both the net dipole moment $\braket{\mathbf{M}}$ and |
392 |
dipolar polarization $\bf{P}$ tends to vanish but |
393 |
$\braket{\mathbf{M}^2}$ does not. The macroscopic dipole |
394 |
polarizability can therefore be written as, |
395 |
\begin{equation} |
396 |
\alpha_D = \frac{\braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2}{3 \epsilon_o V k_B T} |
397 |
\end{equation} |
398 |
This relationship between a macroscopic property of dipolar material |
399 |
and microscopic fluctuations is true even when the applied field |
400 |
$ \textbf{E} \rightarrow 0 $. |
401 |
|
402 |
\subsection{Correction Factors} |
403 |
\label{sec:corrFactor} |
404 |
In the presence of a uniform external field $ \mathbf{E}^\circ$, the |
405 |
total electric field at $\mathbf{r}$ depends on the polarization |
406 |
density at all other points in the system,\cite{NeumannI83} |
407 |
\begin{equation} |
408 |
\mathbf{E}(\mathbf{r}) = \mathbf{E}^\circ(\mathbf{r}) + |
409 |
\frac{1}{4\pi\epsilon_o} \int d\mathbf{r}^\prime |
410 |
\mathbf{T}(\mathbf{r}-\mathbf{r}^\prime)\cdot |
411 |
{\mathbf{P}(\mathbf{r}^\prime)}. |
412 |
\label{eq:localField} |
413 |
\end{equation} |
414 |
$\mathbf{T}$ is the dipole interaction tensor connecting dipoles at |
415 |
$\mathbf{r}^\prime$ with the point of interest ($\mathbf{r}$). |
416 |
|
417 |
In simulations of dipolar fluids, the molecular dipoles may be |
418 |
represented either by closely-spaced point charges or by |
419 |
point dipoles (see Fig. \ref{fig:tensor}). |
420 |
\begin{figure} |
421 |
\includegraphics[width=\linewidth]{Tensors} |
422 |
\caption{In the real-space electrostatic methods, the molecular dipole |
423 |
tensor, $\mathbf{T}_{\alpha\beta}(r)$, is not the same for |
424 |
charge-charge interactions as for point dipoles (left panel). The |
425 |
same holds true for the molecular quadrupole tensor (right panel), |
426 |
$\mathbf{T}_{\alpha\beta\gamma\delta}(r)$, which can have distinct |
427 |
forms if the molecule is represented by charges, dipoles, or point |
428 |
quadrupoles.} |
429 |
\label{fig:tensor} |
430 |
\end{figure} |
431 |
In the case where point charges are interacting via an electrostatic |
432 |
kernel, $v(r)$, the effective {\it molecular} dipole tensor, |
433 |
$\mathbf{T}$ is obtained from two successive applications of the |
434 |
gradient operator to the electrostatic kernel, |
435 |
\begin{align} |
436 |
\mathbf{T}_{\alpha \beta}(r) =& \nabla_\alpha \nabla_\beta |
437 |
\left(v(r)\right) \\ |
438 |
=& \delta_{\alpha \beta} |
439 |
\left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha} |
440 |
r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r} |
441 |
v^{\prime}(r) \right) |
442 |
\label{dipole-chargeTensor} |
443 |
\end{align} |
444 |
where $v(r)$ may be either the bare kernel ($1/r$) or one of the |
445 |
modified (Wolf or DSF) kernels. This tensor describes the effective |
446 |
interaction between molecular dipoles ($\mathbf{D}$) in Gaussian units |
447 |
as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$. |
448 |
|
449 |
When utilizing any of the three new real-space methods for point |
450 |
\textit{dipoles}, the tensor is explicitly constructed, |
451 |
\begin{equation} |
452 |
\mathbf{T}_{\alpha \beta}(r) = \delta_{\alpha \beta} v_{21}(r) + |
453 |
\frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r) |
454 |
\label{dipole-diopleTensor} |
455 |
\end{equation} |
456 |
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of |
457 |
the approximation.\cite{PaperI,PaperII} Although the Taylor-shifted |
458 |
(TSF) and gradient-shifted (GSF) models produce to the same $v(r)$ |
459 |
function for point charges, they have distinct forms for the |
460 |
dipole-dipole interaction. |
461 |
|
462 |
Using a constitutive relation, the polarization density |
463 |
$\mathbf{P}(\mathbf{r})$ is given by, |
464 |
\begin{equation} |
465 |
\mathbf{P}(\mathbf{r}) = \epsilon_o \chi_D |
466 |
\left(\mathbf{E}^\circ(\mathbf{r}) + \frac{1}{4\pi\epsilon_o} \int |
467 |
d\mathbf{r}^\prime \mathbf{T}(\mathbf{r}-\mathbf{r}^\prime ) \cdot \mathbf{P}(\mathbf{r}^\prime)\right). |
468 |
\label{constDipole} |
469 |
\end{equation} |
470 |
Note that $\chi_D$ explicitly depends on the details of the dipole |
471 |
interaction tensor. Neumann \textit{et al.} |
472 |
\cite{NeumannI83,NeumannII83,Neumann84,Neumann85} derived an elegant |
473 |
way to modify the fluctuation formula to correct for approximate |
474 |
interaction tensors. This correction was derived using a Fourier |
475 |
representation of the interaction tensor, |
476 |
$\tilde{\mathbf{T}}(\mathbf{k})$, and involves the quantity, |
477 |
\begin{equation} |
478 |
A = \frac{3}{4\pi}\tilde{\mathbf{T}}(0) = \frac{3}{4\pi} \int_V |
479 |
d\mathbf{r} \mathbf{T}(\mathbf{r}) |
480 |
\end{equation} |
481 |
which is the $k \rightarrow 0$ limit of $\tilde{\mathbf{T}}$. Using |
482 |
this quantity (originally called $Q$ in |
483 |
refs. \onlinecite{NeumannI83,NeumannII83,Neumann84,Neumann85}), the |
484 |
dielectric constant can be computed |
485 |
\begin{equation} |
486 |
\epsilon = \frac{3+(A + 2)(\epsilon_{CB}-1)}{3 + (A -1)(\epsilon_{CB}-1)} |
487 |
\label{correctionFormula} |
488 |
\end{equation} |
489 |
where $\epsilon_{CB}$ is the widely-used conducting boundary |
490 |
expression for the dielectric constant, |
491 |
\begin{equation} |
492 |
\epsilon_{CB} = 1 + \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 |
493 |
\epsilon_o V k_B T} |
494 |
\label{conductingBoundary} |
495 |
\end{equation} |
496 |
Eqs. (\ref{correctionFormula}) and (\ref{conductingBoundary}) allows |
497 |
estimation of the static dielectric constant from fluctuations easily |
498 |
computed directly from simulations. |
499 |
|
500 |
% Here $\chi^*_D$ is a dipolar susceptibility can be expressed in terms of dielectric constant as $ \chi^*_D = \epsilon - 1$ which different than macroscopic dipolar polarizability $\alpha_D$ in the sections \ref{subsec:perturbation} and \ref{subsec:fluctuation}. We can split integral into two parts: singular part i.e $|\textbf{r}-\textbf{r}'|\rightarrow 0 $ and non-singular part i.e $|\textbf{r}-\textbf{r}'| > 0 $ . The singular part of the integral can be written as,\cite{NeumannI83, Jackson98} |
501 |
% \begin{equation} |
502 |
% \frac{1}{4\pi\epsilon_o} \int_{|\textbf{r}-\textbf{r}'| \rightarrow 0} d^3r'\; \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')} = - \frac{\textbf{P}(\textbf{r})}{3\epsilon_o} |
503 |
% \label{singular} |
504 |
% \end{equation} |
505 |
% Substituting equation (\ref{singular}) in the equation (\ref{constDipole}) we get, |
506 |
% \begin{equation} |
507 |
% \textbf{P}(\textbf{r}) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(\textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int_{|\textbf{r}-\textbf{r}'| > 0} d^3r'\; \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}\right). |
508 |
% \end{equation} |
509 |
% For both polarization and electric field homogeneous, this can be easily solved using Fourier transformation, |
510 |
% \begin{equation} |
511 |
% \textbf{P}(\mathbf{k}) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(1- \frac{3}{4\pi}\;\frac{\chi^*_D}{\chi^*_D + 3}\; \textbf{T}(\mathbf{k})\right)^{-1}\textbf{E}^o(\mathbf{k}). |
512 |
% \end{equation} |
513 |
% For homogeneous applied field Fourier component is non-zero only if $\kappa = 0$. Therefore, |
514 |
% \begin{equation} |
515 |
% \textbf{P}(0) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(1- \frac{\chi^*_D}{\chi^*_D + 3}\; A_{dipole})\right)^{-1}\textbf{E}^o({0}). |
516 |
% \label{fourierDipole} |
517 |
% \end{equation} |
518 |
% where $A_{dipole}=\frac{3}{4\pi}T(0) = \frac{3}{4\pi} \int_V d^3r\;T(r)$. Now equation (\ref{fourierDipole}) can be compared with equation (\ref{flucDipole}). Therefore, |
519 |
% \begin{equation} |
520 |
% \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T} = \frac{3\;\chi^*_D}{\chi^*_D + 3} \left(1- \frac{\chi^*_D}{\chi^*_D + 3}\; A_{dipole})\right)^{-1} |
521 |
% \end{equation} |
522 |
% Substituting $\chi^*_D = \epsilon-1$ and $ \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T} = \epsilon_{CB}-1 = \alpha_D$ in above equation we get, |
523 |
% \begin{equation} |
524 |
% \epsilon = \frac{3+(A_{dipole} + 2)(\epsilon_{CB}-1)}{3+(A_{dipole} -1)(\epsilon_{CB}-1)} = \frac{3+(A_{dipole} + 2)\alpha_D}{3+(A_{dipole} -1)\alpha_D} |
525 |
% \label{correctionFormula} |
526 |
% \end{equation} |
527 |
|
528 |
We have utilized the Neumann \textit{et al.} approach for the three |
529 |
new real-space methods, and obtain method-dependent correction |
530 |
factors. The expression for the correction factor also depends on |
531 |
whether the simulation involves point charges or point dipoles to |
532 |
represent the molecular dipoles. These corrections factors are |
533 |
listed in Table \ref{tab:A}. |
534 |
\begin{table} |
535 |
\caption{Expressions for the dipolar correction factor ($A$) for the |
536 |
real-space electrostatic methods in terms of the damping parameter |
537 |
($\alpha$) and the cutoff radius ($r_c$). The Ewald-Kornfeld result |
538 |
derived in Refs. \onlinecite{Adams80,Adams81,NeumannI83} is shown |
539 |
for comparison using the Ewald convergence parameter ($\kappa$) |
540 |
and the real-space cutoff value ($r_c$). } |
541 |
\label{tab:A} |
542 |
\begin{tabular}{l|c|c} |
543 |
\toprule |
544 |
Method & point charges & point dipoles \\ |
545 |
\colrule |
546 |
Spherical Cutoff (SC) & $\mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\ |
547 |
Shifted Potental (SP) & $ \mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(r_c \alpha) -\frac{2 \alpha r_c}{\sqrt{\pi}}\left(1+\frac{2\alpha^2 {r_c}^2}{3} \right)e^{-\alpha^2{r_c}^2} $\\ |
548 |
Gradient-shifted (GSF) & 1 & $\mathrm{erf}(\alpha r_c)-\frac{2 \alpha r_c}{\sqrt{\pi}} \left(1 + \frac{2 \alpha^2 r_c^2}{3} + \frac{\alpha^4 r_c^4}{3}\right)e^{-\alpha^2 r_c^2} $ \\ |
549 |
Taylor-shifted (TSF) & 1 & 1 \\ \colrule |
550 |
Ewald-Kornfeld (EK) & \multicolumn{2}{c}{$\mathrm{erf}(r_c \kappa) - |
551 |
\frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 |
552 |
r_c^2}$} \\ |
553 |
\botrule |
554 |
\end{tabular} |
555 |
\end{table} |
556 |
Note that for point charges, the GSF and TSF methods produce estimates |
557 |
of the dielectric that need no correction, and the TSF method likewise |
558 |
needs no correction for point dipoles. |
559 |
|
560 |
\section{Quadrupolar Fluids and the Quadrupolar Susceptibility} |
561 |
\subsection{Response to External Perturbations} |
562 |
|
563 |
A molecule with a permanent quadrupole, $\mathsf{q}$, will align in |
564 |
the presence of an electric field gradient $\nabla\mathbf{E}$. The |
565 |
anisotropic polarization of the quadrupole is given |
566 |
by,\cite{AduGyamfi78,AduGyamfi81} |
567 |
\begin{equation} |
568 |
\braket{\mathsf{q}} - \frac{\mathbf{I}}{3} |
569 |
\mathrm{Tr}(\mathsf{q}) = \epsilon_o \alpha_q \nabla\mathbf{E}, |
570 |
\end{equation} |
571 |
where $\alpha_q = q_o^2 / 15 \epsilon_o k_B T $ is a molecular quadrupole |
572 |
polarizablity and $q_o$ is an effective quadrupole moment for the molecule, |
573 |
\begin{equation} |
574 |
q_o^2 = 3 \mathsf{q}:\mathsf{q} - \mathrm{Tr}(\mathsf{q})^2. |
575 |
\end{equation} |
576 |
|
577 |
In the presence of an external field gradient, a system of quadrupolar |
578 |
molecules also organizes with an anisotropic polarization, |
579 |
\begin{equation} |
580 |
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q}) = \epsilon_o |
581 |
\alpha_Q \nabla\mathbf{E} |
582 |
\end{equation} |
583 |
where $\mathsf{Q}$ is the traced quadrupole density of the system and |
584 |
$\alpha_Q$ is a macroscopic quadrupole polarizability which has |
585 |
dimensions of $\mathrm{length}^{-2}$. Equivalently, the traceless form |
586 |
may be used, |
587 |
\begin{equation} |
588 |
\mathsf{\Theta} = 3 \epsilon_o \alpha_Q \nabla\mathbf{E}, |
589 |
\label{pertQuad} |
590 |
\end{equation} |
591 |
where |
592 |
$\mathsf{\Theta} = 3\mathsf{Q} - \mathbf{I} \mathrm{Tr}(\mathsf{Q})$ |
593 |
is the traceless tensor that also describes the system quadrupole |
594 |
density. It is this tensor that will be utilized to derive correction |
595 |
factors below. |
596 |
|
597 |
\subsection{Fluctuation Formula} |
598 |
As in the dipolar case, we may define a system quadrupole moment, |
599 |
$\mathsf{M}_Q = \sum_i \mathsf{q}_i$ and the traced quadrupolar |
600 |
density, $\mathsf{Q} = \mathsf{M}_Q / V$. A fluctuation formula can |
601 |
be written for a system comprising quadrupolar |
602 |
molecules,\cite{LoganI81,LoganII82,LoganIII82} |
603 |
\begin{equation} |
604 |
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q}) = \epsilon_o |
605 |
\frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o |
606 |
V k_B T} \cdot \nabla\mathbf{E}. |
607 |
\label{flucQuad} |
608 |
\end{equation} |
609 |
Some care is needed in the definitions of the averaged quantities. These |
610 |
refer to the effective quadrupole moment of the system, and they are |
611 |
computed as follows, |
612 |
\begin{align} |
613 |
\braket{\mathsf{M}_Q^2} &= \braket{3 \mathsf{M}_Q:\mathsf{M}_Q - |
614 |
\mathrm{Tr}(\mathsf{M}_Q)^2 }\\ |
615 |
\braket{\mathsf{M}_Q}^2 &= 3 \braket{\mathsf{M}_Q}:\braket{\mathsf{M}_Q} - |
616 |
\mathrm{Tr}(\braket{\mathsf{M}_Q})^2 |
617 |
\label{eq:flucQuad} |
618 |
\end{align} |
619 |
The macroscopic quadrupolarizability is given by, |
620 |
\begin{equation} |
621 |
\alpha_Q = \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o V k_B T} |
622 |
\label{propConstQuad} |
623 |
\end{equation} |
624 |
This relationship between a macroscopic property of a quadrupolar |
625 |
fluid and microscopic fluctuations should be true even when the |
626 |
applied field gradient $\nabla \mathbf{E} \rightarrow 0$. |
627 |
|
628 |
\subsection{Correction Factors} |
629 |
In this section we generalize the treatment of Neumann \textit{et al.} |
630 |
for quadrupolar fluids. Interactions involving multiple quadrupoles |
631 |
are rank 4 tensors, and we therefore describe quantities in this |
632 |
section using Einstein notation. |
633 |
|
634 |
In the presence of a uniform external field gradient, |
635 |
$\partial_\alpha {E}^\circ_\beta $, the total field gradient at |
636 |
$\mathbf{r}$ depends on the quadrupole polarization density at all |
637 |
other points in the system, |
638 |
\begin{equation} |
639 |
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha |
640 |
{E}^\circ_\beta(\mathbf{r}) + \frac{1}{4\pi \epsilon_o}\int |
641 |
T_{\alpha\beta\gamma\delta}({\mathbf{r}-\mathbf{r}^\prime}) |
642 |
Q_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime |
643 |
\label{gradMaxwell} |
644 |
\end{equation} |
645 |
where and $T_{\alpha\beta\gamma\delta}$ is the quadrupole interaction |
646 |
tensor connecting quadrupoles at $\mathbf{r}^\prime$ with the point of |
647 |
interest ($\mathbf{r}$). |
648 |
|
649 |
% \begin{figure} |
650 |
% \includegraphics[width=3in]{QuadrupoleFigure} |
651 |
% \caption{With the real-space electrostatic methods, the molecular |
652 |
% quadrupole tensor, $\mathbf{T}_{\alpha\beta\gamma\delta}(r)$, |
653 |
% governing interactions between molecules is not the same for |
654 |
% quadrupoles represented via sets of charges, point dipoles, or by |
655 |
% single point quadrupoles.} |
656 |
% \label{fig:quadrupolarFluid} |
657 |
% \end{figure} |
658 |
|
659 |
In simulations of quadrupolar fluids, the molecular quadrupoles may be |
660 |
represented by closely-spaced point charges, by multiple point |
661 |
dipoles, or by a single point quadrupole (see |
662 |
Fig. \ref{fig:tensor}). In the case where point charges are |
663 |
interacting via an electrostatic kernel, $v(r)$, the effective |
664 |
molecular quadrupole tensor can obtained from four successive |
665 |
applications of the gradient operator to the electrostatic kernel, |
666 |
\begin{eqnarray} |
667 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}) &=& \nabla_\alpha \nabla_\beta |
668 |
\nabla_\gamma \nabla_\delta v(r) \\ |
669 |
&=& \left(\delta_{\alpha\beta}\delta_{\gamma\delta} + |
670 |
\delta_{\alpha\gamma}\delta_{\beta\delta}+ |
671 |
\delta_{\alpha\delta}\delta_{\beta\gamma}\right)\left(-\frac{v^\prime(r)}{r^3} |
672 |
+ \frac{v^{\prime \prime}(r)}{r^2}\right) \nonumber \\ |
673 |
& &+ \left(\delta_{\alpha\beta} r_\gamma r_\delta + 5 \mathrm{~permutations} |
674 |
\right) \left(\frac{3v^\prime(r)}{r^5}-\frac{3v^{\prime \prime}(r)}{r^4} + |
675 |
\frac{v^{\prime \prime \prime}(r)}{r^3}\right) \nonumber \\ |
676 |
& &+ r_\alpha r_\beta r_\gamma r_\delta |
677 |
\left(-\frac{15v^\prime(r)}{r^7}+\frac{15v^{\prime \prime}(r)}{r^6}-\frac{6v^{\prime |
678 |
\prime \prime}(r)}{r^5} + \frac{v^{\prime \prime \prime \prime}(r)}{r^4}\right), |
679 |
\label{quadCharge} |
680 |
\end{eqnarray} |
681 |
where $v(r)$ can either be the electrostatic kernel ($1/r$) or one of |
682 |
the modified (Wolf or DSF) kernels. |
683 |
|
684 |
Similarly, when representing quadrupolar molecules with multiple point |
685 |
dipoles, the molecular quadrupole interaction tensor can be obtained |
686 |
using two successive applications of the gradient operator to the |
687 |
dipole interaction tensor, |
688 |
\begin{eqnarray} |
689 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}) &=& \nabla_\alpha \nabla_\beta |
690 |
v_{\gamma\delta}(\mathbf{r}) \\ |
691 |
& = & \delta_{\alpha\beta}\delta_{\gamma\delta} \frac{v^\prime_{21}(r)}{r} + |
692 |
\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+ |
693 |
\delta_{\alpha\delta}\delta_{\beta\gamma}\right)\frac{v_{22}(r)}{r^2} |
694 |
\nonumber\\ |
695 |
& &+ \delta_{\gamma\delta} r_\alpha r_\beta |
696 |
\left(\frac{v^{\prime \prime}_{21}(r)}{r^2}-\frac{v^\prime_{21}(r)}{r^3} \right)\nonumber \\ |
697 |
& &+\left(\delta_{\alpha\beta} r_\gamma r_\delta + |
698 |
\delta_{\alpha\gamma} r_\beta r_\delta +\delta_{\alpha\delta} |
699 |
r_\gamma r_\beta + \delta_{\beta\gamma} r_\alpha r_\delta |
700 |
+\delta_{\beta\delta} r_\alpha r_\gamma \right) |
701 |
\left(\frac{v^\prime_{22}(r)}{r^3}-\frac{2v_{22}(r)}{r^4}\right) |
702 |
\nonumber \\ |
703 |
& &+ r_\alpha r_\beta r_\gamma r_\delta |
704 |
\left(\frac{v^{\prime |
705 |
\prime}_{22}(r)}{r^4}-\frac{5v^\prime_{22}(r)}{r^5}+\frac{8v_{22}(r)}{r^6}\right), |
706 |
\label{quadDip} |
707 |
\end{eqnarray} |
708 |
where $v_{\gamma\delta}(\mathbf{r})$ is a dipole-dipole interaction |
709 |
tensor that depends on the level of the |
710 |
approximation.\cite{PaperI,PaperII} Similarly $v_{21}(r)$ and |
711 |
$v_{22}(r)$ are the radial function for different real space cutoff |
712 |
methods defined in the first paper in this series.\cite{PaperI} |
713 |
|
714 |
For quadrupolar liquids modeled using point quadrupoles, the |
715 |
interaction tensor can be constructed as, |
716 |
\begin{eqnarray} |
717 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}) & = & |
718 |
\left(\delta_{\alpha\beta}\delta_{\gamma\delta} |
719 |
+ |
720 |
\delta_{\alpha\gamma}\delta_{\beta\delta}+ |
721 |
\delta_{\alpha\delta}\delta_{\beta\gamma}\right)v_{41}(r) |
722 |
+ \delta_{\gamma\delta} |
723 |
r_\alpha r_\beta |
724 |
\frac{v_{42}(r)}{r^2} |
725 |
\nonumber \\ |
726 |
& & + r_\alpha r_\beta r_\gamma r_\delta \left(\frac{v_{43}(r)}{r^4}\right), |
727 |
\label{quadRadial} |
728 |
\end{eqnarray} |
729 |
where again $v_{41}(r)$, $v_{42}(r)$, and $v_{43}(r)$ are radial |
730 |
functions defined in Paper I of the series. \cite{PaperI} Note that |
731 |
these radial functions have different functional forms depending on |
732 |
the level of approximation being employed. |
733 |
|
734 |
The integral in equation (\ref{gradMaxwell}) can be divided into two |
735 |
parts, $|\mathbf{r}-\mathbf{r}^\prime|\rightarrow 0 $ and |
736 |
$|\mathbf{r}-\mathbf{r}^\prime|> 0$. Since the self-contribution to |
737 |
the field gradient vanishes at the singularity (see the supporting |
738 |
information), equation (\ref{gradMaxwell}) can be written as, |
739 |
\begin{equation} |
740 |
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha {E}^\circ_\beta(\mathbf{r}) + |
741 |
\frac{1}{4\pi \epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 } |
742 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime) |
743 |
{Q}_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime. |
744 |
\end{equation} |
745 |
If $\mathbf{r} = \mathbf{r}^\prime$ is excluded from the integration, |
746 |
the total gradient can be most easily expressed in terms of |
747 |
traceless quadrupole density as below,\cite{LoganI81} |
748 |
\begin{equation} |
749 |
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha |
750 |
{E}^\circ_\beta(\mathbf{r}) + \frac{1}{12\pi |
751 |
\epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 } |
752 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime) \Theta_{\gamma\delta}(\mathbf{r}') d\mathbf{r}^\prime, |
753 |
\end{equation} |
754 |
where |
755 |
$\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$ |
756 |
is the traceless quadrupole density. In analogy to |
757 |
Eq. (\ref{pertQuad}) above, the quadrupolarization density may now be |
758 |
related to the actual quadrupolar susceptibility, $\chi_Q$, |
759 |
\begin{equation} |
760 |
\frac{1}{3}{\Theta}_{\alpha\beta}(\mathbf{r}) = \epsilon_o {\chi}_Q |
761 |
\left(\partial_\alpha {E}^\circ_\beta(\mathbf{r}) + \frac{1}{12\pi |
762 |
\epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 } |
763 |
T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime) |
764 |
\Theta_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime \right). |
765 |
\end{equation} |
766 |
For toroidal boundaries and with a uniform imposed |
767 |
$\partial_\alpha E^\circ_\beta$, the quadrupole density |
768 |
${\Theta}_{\alpha\beta}$ will be uniform over the entire space. After |
769 |
performing a Fourier transform (see the Appendix in |
770 |
ref. \onlinecite{NeumannI83}) we obtain, |
771 |
\begin{equation} |
772 |
\frac{1}{3}\tilde{\Theta}_{\alpha\beta}(\mathbf{k})= |
773 |
\epsilon_o {\chi}_Q \left[{\partial_\alpha |
774 |
\tilde{E}^\circ_\beta}(\mathbf{k})+ \frac{1}{12\pi |
775 |
\epsilon_o} \tilde{T}_{\alpha\beta\gamma\delta}(\mathbf{k}) |
776 |
\tilde{\Theta}_{\gamma\delta}(\mathbf{k})\right]. |
777 |
\end{equation} |
778 |
|
779 |
The only surviving components of the quadrupole polarization will be |
780 |
in the direction of the applied field gradient, so we may simplify |
781 |
this expression considerably. That is, the only relevant polarization |
782 |
term is $\tilde{\Theta}_{\alpha\beta}(\mathbf{k})$ and the only relevant |
783 |
tensor contributions will have doubly-repeated indices, e.g. |
784 |
$\tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k})$. Therefore we may |
785 |
consider each component of the interaction tensor separately in a |
786 |
scalar expression, |
787 |
\begin{eqnarray} |
788 |
\frac{1}{3}\tilde{\Theta}_{\alpha\beta}(\mathbf{k}) &=& \epsilon_o |
789 |
{\chi}_Q |
790 |
\left[{\partial_\alpha |
791 |
\tilde{E}^\circ_\beta}(\mathbf{k})+ |
792 |
\frac{1}{12\pi |
793 |
\epsilon_o}\tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k}) |
794 |
{\tilde{\Theta}}_{\alpha\beta}(\mathbf{k})\right] |
795 |
\nonumber \\ |
796 |
&=& \epsilon_o {\chi}_Q \left(1-\frac{1}{4\pi} {\chi}_Q \tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k})\right)^{-1} |
797 |
{\partial_\alpha \tilde{E}^\circ_\beta}(\mathbf{k}) |
798 |
\label{fourierQuad} |
799 |
\end{eqnarray} |
800 |
If the applied field gradient is homogeneous over the |
801 |
entire volume, ${\partial_ \alpha \tilde{E}^\circ_\beta}(\mathbf{k}) = 0 $ except at |
802 |
$ \mathbf{k} = 0$. As in the dipolar case, the only relevant |
803 |
contribution of the interaction tensor will also be when $\mathbf{k} = 0$. |
804 |
Therefore equation (\ref{fourierQuad}) can be written as, |
805 |
\begin{equation} |
806 |
\frac{1}{3}{\tilde{\Theta}}_{\alpha\beta}({0}) = \epsilon_o {\chi}_Q |
807 |
\left(1-\frac{1}{4\pi} {\chi}_Q \tilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1} \partial_\alpha \tilde{E}^\circ_\beta({0}) |
808 |
\label{fourierQuad2} |
809 |
\end{equation} |
810 |
where $\tilde{T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as, |
811 |
\begin{equation} |
812 |
\tilde{T}_{\alpha\beta\alpha\beta}({0}) = \int_V {T}_{\alpha\beta\alpha\beta}(\mathbf{r})d\mathbf{r} |
813 |
\label{realTensorQaud} |
814 |
\end{equation} |
815 |
We may now define a quantity to represent the contribution from the |
816 |
$\mathbf{k} = 0$ portion of the interaction tensor, |
817 |
\begin{equation} |
818 |
B = \frac{1}{4 \pi} \int_V T_{\alpha \beta \alpha \beta}(\mathbf{r}) |
819 |
d\mathbf{r} |
820 |
\end{equation} |
821 |
which has been integrated over the interaction volume $V$ and has |
822 |
units of $\mathrm{length}^{-2}$. |
823 |
|
824 |
In terms of traced quadrupole moment, equation (\ref{fourierQuad2}) |
825 |
can be written as, |
826 |
\begin{equation} |
827 |
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q}) |
828 |
= \epsilon_o {\chi}_Q\left(1- {\chi}_Q B \right)^{-1} \nabla \mathbf{E}^\circ |
829 |
\label{tracedConstQuad} |
830 |
\end{equation} |
831 |
Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we obtain, |
832 |
\begin{equation} |
833 |
\frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o |
834 |
V k_B T} = \chi_Q \left(1 - \chi_Q B \right)^{-1}, |
835 |
\end{equation} |
836 |
or equivalently, |
837 |
\begin{equation} |
838 |
\chi_Q = \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o |
839 |
V k_B T} \left(1 + B \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o |
840 |
V k_B T} \right)^{-1} |
841 |
\label{eq:finalForm} |
842 |
\end{equation} |
843 |
Eq. (\ref{eq:finalForm}) now expresses a bulk property (the |
844 |
quadrupolar susceptibility, $\chi_Q$) in terms of a straightforward |
845 |
fluctuation in the system quadrupole moment and a quadrupolar |
846 |
correction factor ($B$). The correction factors depend on the cutoff |
847 |
method being employed in the simulation, and these are listed in Table |
848 |
\ref{tab:B}. |
849 |
|
850 |
In terms of the macroscopic quadrupole polarizability, $\alpha_Q$, |
851 |
which may be thought of as the ``conducting boundary'' version of the |
852 |
susceptibility, |
853 |
\begin{equation} |
854 |
\chi_Q = \frac{\alpha_Q}{1 + B \alpha_Q}. |
855 |
\label{eq:quadrupolarSusceptiblity} |
856 |
\end{equation} |
857 |
If an electrostatic method produces $B \rightarrow 0$, the computed |
858 |
quadrupole polarizability and quadrupole susceptibility converge to |
859 |
the same value. |
860 |
|
861 |
\begin{sidewaystable} |
862 |
\caption{Expressions for the quadrupolar correction factor ($B$) for |
863 |
the real-space electrostatic methods in terms of the damping |
864 |
parameter ($\alpha$) and the cutoff radius ($r_c$). The units |
865 |
of the correction factor are $ \mathrm{length}^{-2}$ for quadrupolar |
866 |
fluids.} |
867 |
\label{tab:B} |
868 |
\begin{tabular}{l|c|c|c} |
869 |
\toprule |
870 |
Method & charges & dipoles & quadrupoles \\\colrule |
871 |
Spherical Cutoff (SC) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $ -\frac{8 {\alpha}^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\ |
872 |
Shifted Potental (SP) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $- \frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$& $ -\frac{16 \alpha^7 {r_c}^5}{9\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\ |
873 |
Gradient-shifted (GSF) & $- \frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & 0 & $-\frac{4{\alpha}^7{r_c}^5 }{9\sqrt{\pi}}e^{-\alpha^2 r_c^2}(-1+2\alpha ^2 r_c^2)$\\ |
874 |
Taylor-shifted (TSF) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} |
875 |
e^{-\alpha^2 r_c^2}$ & |
876 |
$\frac{4\;\mathrm{erfc(\alpha |
877 |
r_c)}}{{r_c}^2} + |
878 |
\frac{8 |
879 |
\alpha}{3\sqrt{\pi}r_c}e^{-\alpha^2 |
880 |
{r_c}^2}\left(3+ 2 |
881 |
\alpha^2 {r_c}^2 + |
882 |
\alpha^4 {r_c}^4\right) |
883 |
$ & |
884 |
$\frac{10\;\mathrm{erfc}(\alpha r_c )}{{r_c}^2} + \frac{4{\alpha}}{9\sqrt{\pi}{r_c}}e^{-\alpha^2 r_c^2}\left(45 + 30\alpha ^2 {r_c}^2 + 12\alpha^4 {r_c}^4 + 3\alpha^6 {r_c}^6 + 2 \alpha^8 {r_c}^8\right)$ \\ |
885 |
\colrule |
886 |
Ewald-Kornfeld (EK) & \multicolumn{3}{c}{ Nothing Yet} \\ |
887 |
\botrule |
888 |
\end{tabular} |
889 |
\end{sidewaystable} |
890 |
|
891 |
\section{Screening of Charges by Multipolar Fluids} |
892 |
\label{sec:PMF} |
893 |
In a dipolar fluid, the static dielectric constant is also a measure |
894 |
of the ability of the fluid to screen charges from one another. A set |
895 |
of point charges creates an inhomogeneous field in the fluid, and the |
896 |
fluid responds to this field as if it was created externally or via |
897 |
local polarization fluctuations. For this reason, the dielectric |
898 |
constant can be used to estimate an effective potential between two |
899 |
point charges ($C_i$ and $C_j$) embedded in the fluid, |
900 |
\begin{equation} |
901 |
U_\mathrm{effective} = \frac{C_i C_j}{4 \pi \epsilon_0 \epsilon |
902 |
r_{ij}}. |
903 |
\label{eq:effectivePot} |
904 |
\end{equation} |
905 |
|
906 |
The same set of point charges can also create an inhomogeneous field |
907 |
\textit{gradient}, and this will cause a response in a quadrupolar |
908 |
fluid that will also cause an effective screening. As discussed above, |
909 |
however, the relevant phyiscal property in quadrupolar fluids is the |
910 |
susceptibility, $\chi_Q$. The screening dielectric associated with |
911 |
the quadrupolar susceptibility is defined as,\cite{Ernst92} |
912 |
\begin{equation} |
913 |
\epsilon = 1 + \chi_Q G = 1 + G \frac{\alpha_Q}{1 + \alpha_Q B} |
914 |
\label{eq:dielectricFromQuadrupoles} |
915 |
\end{equation} |
916 |
where $G$ is a geometrical factor that depends on the geometry of the |
917 |
field perturbation, |
918 |
\begin{equation} |
919 |
G = \frac{\int_V \left| \nabla \mathbf{E}^\circ \right|^2 d\mathbf{r}} |
920 |
{\int_V \left|\mathbf{E}^\circ\right|^2 d\mathbf{r}} |
921 |
\end{equation} |
922 |
integrated over the interaction volume. Note that this geometrical |
923 |
factor is also required to compute effective dielectric constants even |
924 |
when the field gradient is homogeneous over the entire sample. |
925 |
|
926 |
To measure effective screening in a multipolar fluid, we compute an |
927 |
effective interaction potential, the potential of mean force (PMF), |
928 |
between positively and negatively charged ions when they screened by |
929 |
the intervening fluid. The PMF is obtained from a sequence of |
930 |
simulations in which two ions are constrained to a fixed distance, and |
931 |
the average constraint force to hold them at a fixed distance $r$ is |
932 |
collected during a long simulation,\cite{Wilfred07} |
933 |
\begin{equation} |
934 |
w(r) = \int_{r_o}^{r}\left\langle \frac{\partial f}{\partial r^\prime} |
935 |
\right\rangle_{r^\prime} dr^\prime + 2k_BT \ln\left(\frac{r}{r_o}\right) + w(r_o), |
936 |
\label{eq:pmf} |
937 |
\end{equation} |
938 |
where $\braket{\partial f/\partial r^\prime}_{r^\prime}$ is the mean |
939 |
constraint force required to hold the ions at distance $r^\prime$, |
940 |
$2k_BT \log(r/r_o)$ is the Fixman factor,\cite{Fixman:1974fk} and |
941 |
$r_o$ is a reference position (usually taken as a large separation |
942 |
betwen the ions). If the dielectric constant is a good measure of the |
943 |
screening at all inter-ion separations, we would expect $w(r)$ to have |
944 |
the form in Eq. (\ref{eq:effectivePot}). Because real fluids are not |
945 |
continuum dielectrics, the effective dielectric constant is a function |
946 |
of the interionic separation, |
947 |
\begin{equation} |
948 |
\epsilon(r) = \frac{u_\mathrm{raw}(r) }{w(r)} |
949 |
\end{equation} |
950 |
where $u_\mathrm{raw}(r)$ is the direct charge-charge interaction |
951 |
potential that is in use during the simulation. $\epsilon(r)$ may |
952 |
vary considerably from the bulk estimates at short distances, although |
953 |
it should converge to the bulk value as the separation between the |
954 |
ions increases. |
955 |
|
956 |
\section{Simulation Methodology} |
957 |
To test the formalism developed in the preceding sections, we have |
958 |
carried out computer simulations using three different techniques: i) |
959 |
simulations in the presence of external fields, ii) equilibrium |
960 |
calculations of box moment fluctuations, and iii) potentials of mean |
961 |
force (PMF) between embedded ions. In all cases, the fluids were |
962 |
composed of point multipoles protected by a Lennard-Jones potential. |
963 |
The parameters used in the test systems are given in table |
964 |
\ref{Tab:C}. |
965 |
|
966 |
\begin{table} |
967 |
\caption{\label{Tab:C}The parameters used in simulations to evaluate |
968 |
the dielectric response of the new real-space methods.} |
969 |
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
970 |
& \multicolumn{2}{c|}{LJ parameters} & |
971 |
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
972 |
Test system & $\sigma$& $\epsilon$ & $C$ & $D$ & |
973 |
$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass & $I_{xx}$ & $I_{yy}$ & |
974 |
$I_{zz}$ \\ \cline{6-8}\cline{10-12} |
975 |
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
976 |
\AA\textsuperscript{2})} \\ \hline |
977 |
Dipolar fluid & 3.41 & 0.2381 & - & 1.4026 &-&-&-& 39.948 & 11.613 & 11.613 & 0.0 \\ |
978 |
Quadrupolar fluid & 2.985 & 0.265 & - & - & 0.0 & 0.0 &-2.139 & 18.0153 & 43.0565 & 43.0565 & 0.0 \\ |
979 |
\ce{q+} & 1.0 & 0.1 & +1 & - & - & - & - & 22.98 & - & - & - \\ |
980 |
\ce{q-} & 1.0 & 0.1 & -1 & - & - & - & - & 22.98 & - & - & - \\ \hline |
981 |
\end{tabularx} |
982 |
\end{table} |
983 |
|
984 |
The first of the test systems consists entirely of fluids of point |
985 |
dipolar or quadrupolar molecules in the presence of constant field or |
986 |
field gradients. Since there are no isolated charges within the |
987 |
system, the divergence of the field will be zero, \textit{i.e.} |
988 |
$\vec{\nabla} \cdot \mathbf{E} = 0$. This condition can be satisfied |
989 |
by using the relatively simple applied potential as described in the |
990 |
supporting information. |
991 |
|
992 |
When a constant electric field or field gradient is applied to the |
993 |
system, the molecules align along the direction of the applied field, |
994 |
and polarize to a degree determined both by the strength of the field |
995 |
and the fluid's polarizibility. We have calculated ensemble averages |
996 |
of the box dipole and quadrupole moments as a function of the strength |
997 |
of the applied fields. If the fields are sufficiently weak, the |
998 |
response is linear in the field strength, and one can easily compute |
999 |
the polarizability directly from the simulations. |
1000 |
|
1001 |
The second set of test systems consists of equilibrium simulations of |
1002 |
fluids of point dipolar or quadrupolar molecules simulated in the |
1003 |
absence of any external perturbation. The fluctuation of the ensemble |
1004 |
averages of the box multipolar moment was calculated for each of the |
1005 |
multipolar fluids. The box multipolar moments were computed as simple |
1006 |
sums over the instantaneous molecular moments, and fluctuations in |
1007 |
these quantities were obtained from Eqs. (\ref{eq:flucDip}) and |
1008 |
(\ref{eq:flucQuad}). The macroscopic polarizabilities of the system at |
1009 |
a were derived using Eqs.(\ref{flucDipole}) and (\ref{flucQuad}). |
1010 |
|
1011 |
The final system consists of dipolar or quadrupolar fluids with two |
1012 |
oppositely charged ions embedded within the fluid. These ions are |
1013 |
constrained to be at fixed distance throughout a simulation, although |
1014 |
they are allowed to move freely throughout the fluid while satisfying |
1015 |
that constraint. Separate simulations were run at a range of |
1016 |
constraint distances. A dielectric screening factor was computed using |
1017 |
the ratio between the potential between the two ions in the absence of |
1018 |
the fluid medium and the PMF obtained from the simulations. |
1019 |
|
1020 |
We carried out these simulations for all three of the new real-space |
1021 |
electrostatic methods (SP, GSF, and TSF) that were developed in the |
1022 |
first paper (Ref. \onlinecite{PaperI}) in the series. The radius of |
1023 |
the cutoff sphere was taken to be 12~\AA. Each of the real space |
1024 |
methods also depends on an adjustable damping parameter $\alpha$ (in |
1025 |
units of $\mathrm{length}^{-1}$). We have selected ten different |
1026 |
values of damping parameter: 0.0, 0.05, 0.1, 0.15, 0.175, 0.2, 0.225, |
1027 |
0.25, 0.3, and 0.35~\AA$^{-1}$ in our simulations of the dipolar |
1028 |
liquids, while four values were chosen for the quadrupolar fluids: |
1029 |
0.0, 0.1, 0.2, and 0.3~\AA$^{-1}$. |
1030 |
|
1031 |
For each of the methods and systems listed above, a reference |
1032 |
simulation was carried out using a multipolar implementation of the |
1033 |
Ewald sum.\cite{Smith82,Smith98} A default tolerance of |
1034 |
$1 \times 10^{-8}$~kcal/mol was used in all Ewald calculations, |
1035 |
resulting in Ewald coefficient 0.3119~\AA$^{-1}$ for a cutoff radius |
1036 |
of 12~\AA. All of the electrostatics and constraint methods were |
1037 |
implemented in our group's open source molecular simulation program, |
1038 |
OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
1039 |
this work. |
1040 |
|
1041 |
Dipolar systems contained 2048 Lennard-Jones-protected point dipolar |
1042 |
(Stockmayer) molecules with reduced density $\rho^* = 0.822$, |
1043 |
temperature $T^* = 1.15$, moment of inertia $I^* = 0.025$, and dipole |
1044 |
moment $\mu^* = \sqrt{3.0}$. These systems were equilibrated for |
1045 |
0.5~ns in the canonical (NVT) ensemble. Data collection was carried |
1046 |
out over a 1~ns simulation in the microcanonical (NVE) ensemble. Box |
1047 |
dipole moments were sampled every fs. For simulations with external |
1048 |
perturbations, field strengths ranging from $0 - 10 \times |
1049 |
10^{-4}$~V/\AA with increments of $ 10^{-4}$~V/\AA were carried out |
1050 |
for each system. |
1051 |
|
1052 |
Quadrupolar systems contained 4000 linear point quadrupoles with a |
1053 |
density $2.338 \mathrm{~g/cm}^3$ at a temperature of 500~K. These |
1054 |
systems were equilibrated for 200~ps in a canonical (NVT) ensemble. |
1055 |
Data collection was carried out over a 500~ps simulation in the |
1056 |
microcanonical (NVE) ensemble. Components of box quadrupole moments |
1057 |
were sampled every 100 fs. For quadrupolar simulations with external |
1058 |
field gradients, field strengths ranging from |
1059 |
$0 - 9 \times 10^{-2}$~V/\AA$^2$ with increments of |
1060 |
$10^{-2}$~V/\AA$^2$ were carried out for each system. |
1061 |
|
1062 |
To carry out the PMF simulations, two of the multipolar molecules in |
1063 |
the test system were converted into \ce{q+} and \ce{q-} ions and |
1064 |
constrained to remain at a fixed distance for the duration of the |
1065 |
simulation. The constrained distance was then varied from 5--12~\AA. |
1066 |
In the PMF calculations, all simulations were equilibrated for 500 ps |
1067 |
in the NVT ensemble and run for 5 ns in the microcanonical (NVE) |
1068 |
ensemble. Constraint forces were sampled every 20~fs. |
1069 |
|
1070 |
\section{Results} |
1071 |
\subsection{Dipolar fluid} |
1072 |
\begin{figure} |
1073 |
\includegraphics[width=5in]{dielectricFinal_Dipole.pdf} |
1074 |
\caption{The polarizability ($\alpha_D$), correction factor ($A$), and |
1075 |
dielectric constant ($\epsilon$) for the test dipolar fluid. The |
1076 |
left panels were computed using external fields, and those on the |
1077 |
right are the result of equilibrium fluctuations. In the GSF and SP |
1078 |
methods, the corrections are large in with small values of $\alpha$, |
1079 |
and a optimal damping coefficient is evident around 0.25 \AA$^{-1}$. |
1080 |
The dashed lines in the upper panel indicate back-calculation of the |
1081 |
polarizability using the Ewald estimate (Refs. \onlinecite{Adams81} |
1082 |
and \onlinecite{NeumannI83}) for the dielectric constant.} |
1083 |
\label{fig:dielectricDipole} |
1084 |
\end{figure} |
1085 |
The macroscopic polarizability ($\alpha_D$) for the dipolar fluid is |
1086 |
shown in the upper panels in Fig. \ref{fig:dielectricDipole}. The |
1087 |
polarizability obtained from the both perturbation and fluctuation |
1088 |
approaches are in excellent agreement with each other. The data also |
1089 |
show a stong dependence on the damping parameter for both the Shifted |
1090 |
Potential (SP) and Gradient Shifted force (GSF) methods, while Taylor |
1091 |
shifted force (TSF) is largely independent of the damping |
1092 |
parameter. |
1093 |
|
1094 |
The calculated correction factors discussed in section |
1095 |
\ref{sec:corrFactor} are shown in the middle panels. Because the TSF |
1096 |
method has $A = 1$ for all values of the damping parameter, the |
1097 |
computed polarizabilities need no correction for the dielectric |
1098 |
calculation. The value of $A$ varies with the damping parameter in |
1099 |
both the SP and GSF methods, and inclusion of the correction yields |
1100 |
dielectric estimates (shown in the lower panel) that are generally too |
1101 |
large until the damping reaches $\sim$~0.25~\AA$^{-1}$. Above this |
1102 |
value, the dielectric constants are generally in good agreement with |
1103 |
previous simulation results.\cite{NeumannI83} |
1104 |
|
1105 |
Figure \ref{fig:dielectricDipole} also contains back-calculations of |
1106 |
the polarizability using the reference (Ewald) simulation |
1107 |
results.\cite{NeumannI83} These are indicated with dashed lines in the |
1108 |
upper panels. It is clear that the expected polarizability for the SP |
1109 |
and GSF methods are quite close to results obtained from the |
1110 |
simulations. This indicates that the correction formula for the |
1111 |
dipolar fluid (Eq. \ref{correctionFormula}) is quite sensitive when |
1112 |
the value of $\mathrm{A}$ deviates significantly from unity. |
1113 |
|
1114 |
These results also suggest an optimal value for the damping parameter |
1115 |
of ($\alpha \sim 0.2-0.3$~\AA$^{-1}$ when evaluating dielectric |
1116 |
constants of point dipolar fluids using either the perturbation and |
1117 |
fluctuation approaches for the new real-space methods. |
1118 |
|
1119 |
\begin{figure} |
1120 |
\includegraphics[width=4in]{ScreeningFactor_Dipole.pdf} |
1121 |
\caption{The distance-dependent screening factor, $\epsilon(r)$, |
1122 |
between two ions immersed in the dipolar fluid. The new methods are |
1123 |
shown in separate panels, and different values of the damping |
1124 |
parameter ($\alpha$) are indicated with different symbols. All of |
1125 |
the methods appear to be converging to the bulk dielectric constant |
1126 |
($\sim 65$) at large ion separations. CHECK PLOT} |
1127 |
\label{fig:ScreeningFactor_Dipole} |
1128 |
\end{figure} |
1129 |
We have also evaluated the distance-dependent screening factor, |
1130 |
$\epsilon(r)$, between two oppositely charged ions when they are |
1131 |
placed in the dipolar fluid. These results were computed using |
1132 |
Eq. \ref{eq:pmf} and are shown in Fig. |
1133 |
\ref{fig:ScreeningFactor_Dipole}. |
1134 |
|
1135 |
The screening factor is similar to the dielectric constant, but |
1136 |
measures a local property of the ions in the fluid and depends on both |
1137 |
ion-dipole and dipole-dipole interactions. These interactions depend |
1138 |
on the distance between ions as well as the electrostatic interaction |
1139 |
methods utilized in the simulations. The screening should converge to |
1140 |
the dielectric constant when the field due to ions is small. This |
1141 |
occurs when the ions are separated (or when the damping parameter is |
1142 |
large). In Fig. \ref{fig:ScreeningFactor_Dipole} we observe that for |
1143 |
the higher value of damping alpha \textit{i.e.} |
1144 |
$\alpha = 0.2$~\AA$^{-1}$ and $0.3$~\AA$^{-1}$ and large separation |
1145 |
between ions, the screening factor does indeed approach the correct |
1146 |
dielectric constant. |
1147 |
|
1148 |
REVISIT THESE THREE PARAGRAPHS: |
1149 |
|
1150 |
It is also notable that the TSF method again displays smaller |
1151 |
perturbations away from the correct dielectric screening behavior. We |
1152 |
also observe that for TSF method yields high dielectric screening even |
1153 |
for lower values of $\alpha$. |
1154 |
|
1155 |
At short distances, the presence of the ions creates a strong local |
1156 |
field that acts to align nearby dipoles nearly perfectly in opposition |
1157 |
to the field from the ions. This has the effect of increasing the |
1158 |
effective screening when the ions are brought close to one another. |
1159 |
|
1160 |
This is due to the absence of corrections for dipole-dipole |
1161 |
interactions as discussed above. However, the TSF remains quite |
1162 |
perturbative for ion-dipole interactions, an effect that is most |
1163 |
apparent at small ion separations. |
1164 |
|
1165 |
\subsection{Quadrupolar fluid} |
1166 |
\begin{figure} |
1167 |
\includegraphics[width=4in]{polarizabilityFinal_Quad.pdf} |
1168 |
\caption{The quadrupole polarizability ($\alpha_Q$), correction factor |
1169 |
($B$), and susceptibility ($\chi_Q$) for the test quadrupolar |
1170 |
fluid. The left panels were computed using external field gradients, |
1171 |
and those on the right are the result of equilibrium fluctuations. |
1172 |
The GSF and SP methods allow nearly unmodified use of the |
1173 |
``conducting boundary'' or polarizability results in place of the |
1174 |
bulk susceptibility.} |
1175 |
\label{fig:dielectricQuad} |
1176 |
\end{figure} |
1177 |
The polarizability ($\alpha_Q$), correction factor ($B$), and |
1178 |
susceptibility ($\chi_Q$) for the quadrupolar fluid is plotted against |
1179 |
damping parameter Fig. \ref{fig:dielectricQuad}. In quadrupolar |
1180 |
fluids, both the polarizability and susceptibility have units of |
1181 |
$\mathrm{length}^2$ Although the susceptibility has dimensionality, it |
1182 |
is the relevant measure of macroscopic quadrupolar |
1183 |
properties.\cite{JeonI03, JeonII03}. The left panel in |
1184 |
Fig. \ref{fig:dielectricQuad} shows results obtained from the applied |
1185 |
field gradient simulations whereas the results from the equilibrium |
1186 |
fluctuation formula are plotted in the right panels. |
1187 |
|
1188 |
The susceptibility for the quadrupolar fluid is obtained from |
1189 |
quadrupolarizability and a correction factor using |
1190 |
Eq. (\ref{eq:quadrupolarSusceptiblity}). The susceptibilities are |
1191 |
shown in the bottom panels of Fig. \ref{fig:dielectricQuad}. All three |
1192 |
methods: (SP, GSF, and TSF) produce similar susceptibilities over the |
1193 |
range of damping parameters. This shows that susceptibility derived |
1194 |
using the quadrupolarizability and the correction factors are |
1195 |
essentially independent of the electrostatic method utilized in the |
1196 |
simulation. |
1197 |
|
1198 |
\begin{figure} |
1199 |
\includegraphics[width=5in]{ScreeningFactor_Quad.pdf} |
1200 |
\caption{\label{fig:screenQuad} The distance-dependent screening |
1201 |
factor, $\epsilon(r)$, between two ions immersed in the quadrupolar |
1202 |
fluid. Results from the perturbation and fluctuation methods are |
1203 |
shown in right and central panels. Here the susceptibility is |
1204 |
calculated from the simulation and the geometrical factor is |
1205 |
evaluated analytically, using the field and field-gradient produced |
1206 |
by ions. The right hand panel shows the screening factor obtained |
1207 |
from the PMF calculations.} |
1208 |
\end{figure} |
1209 |
A more difficult test of the quadrupolar susceptibility is made by |
1210 |
comparing with direct calculation of the electrostatic screening using |
1211 |
the potential of mean force (PMF). Since the effective dielectric |
1212 |
constant for a quadrupolar fluid depends on the geometry of the field |
1213 |
and field gradient, this is not a physical property of the quadrupolar |
1214 |
fluid. |
1215 |
|
1216 |
The geometrical factor for embedded ions changes with the ion |
1217 |
separation distance. It is therefore reasonable to treat the |
1218 |
dielectric constant as a distance-dependent screening factor. Since |
1219 |
the quadrupolar molecules couple with the gradient of the field, the |
1220 |
distribution of the quadrupoles will be inhomogeneously distributed |
1221 |
around the point charges. Hence the distribution of quadrupolar |
1222 |
molecules should be taken into account when computing the geometrical |
1223 |
factors in the presence of this perturbation, |
1224 |
\begin{eqnarray} |
1225 |
G &=& \frac{\int_V g(\mathbf{r}) \left|\nabla \mathbf{E}^\circ |
1226 |
\right|^2 d\mathbf{r}}{\int_V\left|\mathbf{E}^\circ\right|^2 |
1227 |
d\mathbf{r}} \nonumber \\ \nonumber \\ |
1228 |
&=& \frac{ 2 \pi \int_{-1}^{1} \int_{0}^{R} r^2 g(r, |
1229 |
\cos\theta) \left|\nabla \mathbf{E}^\circ \right|^2 dr d\cos\theta }{\int_V\left|\mathbf{E}^\circ\right|^2 d\mathbf{r}} |
1230 |
\label{eq:geometricalFactor} |
1231 |
\end{eqnarray} |
1232 |
where $g(r,\cos\theta)$ is a distribution function for the quadrupoles |
1233 |
with respect to an origin at midpoint of a line joining the two probe |
1234 |
charges. |
1235 |
|
1236 |
The effective screening factor is plotted against ion separation |
1237 |
distance in Fig. \ref{fig:screenQuad}. The screening evaluated from |
1238 |
the perturbation and fluctuation methods are shown in right and |
1239 |
central panels. Here the susceptibility is calculated from the |
1240 |
simulation and the geometrical factor is evaluated analytically, using |
1241 |
the field and field-gradient produced by ions. The right hand panel |
1242 |
shows the screening factor obtained from the PMF calculations. |
1243 |
|
1244 |
We note that the screening factor obtained from both the perturbation |
1245 |
and fluctuation formula show good agreement with the screening factor |
1246 |
calculated using PMF method. As there are no large differences in |
1247 |
quadrupole-quadruople interactions for various real-space |
1248 |
methods,\cite{PaperI,PaperII} we generally good agreement for the |
1249 |
screening factors using any of the real space methods. |
1250 |
|
1251 |
\section{Conclusions} |
1252 |
We have used both perturbation and fluctuation approaches to evaluate |
1253 |
dielectric properties for dipolar and quadrupolar fluids. The static |
1254 |
dielectric constant is the relevant bulk property for dipolar fluids, |
1255 |
while the quadrupolar susceptibility plays a similar role for |
1256 |
quadrupoles. Corrections to both the static dielectric constant and |
1257 |
the quadrupolar susceptibility were derived for three new real space |
1258 |
electrostatic methods, and these corrections were tested against a |
1259 |
third measure of dielectric screening, the potential of mean force |
1260 |
between two ions immersed in the fluids. |
1261 |
|
1262 |
For the dipolar fluids, we find that the polarizability evaluated |
1263 |
using the perturbation and fluctuation methods show excellent |
1264 |
agreement, indicating that equilibrium calculations of the dipole |
1265 |
fluctuations are good measures of bulk polarizability. One of the |
1266 |
findings of the second paper in this series is that the moderately |
1267 |
damped GSF and SP methods were most suitable for molecular dynamics |
1268 |
and Monte Carlo simulations, respectively.\cite{PaperII} |
1269 |
|
1270 |
The dielectric constant evaluated using the computed polarizability |
1271 |
and correction factors agrees well with the previous Ewald-based |
1272 |
simulation results \cite{Adams81,NeumannI83} for moderate damping |
1273 |
parameters in the range 0.25--0.3~\AA~$^{-1}$. |
1274 |
|
1275 |
Although the TSF method alters many dynamic and structural features in |
1276 |
multipolar liquids,\cite{PaperII} it is surprisingly good at computing |
1277 |
bulk dielectric properties at nearly all ranges of the damping |
1278 |
parameter. In fact, the correction factor, $A = 1$, for the TSF |
1279 |
method so the conducting boundary formula is essentially correct when |
1280 |
using this method for point dipolar fluids. |
1281 |
|
1282 |
The dielectric correction formula (equation ~\ref{correctionFormula}) |
1283 |
is extremely sensitive when the correction factor ($A$) deviates from |
1284 |
1, and this behavior is also seen in the effective screening of ions |
1285 |
embedded in the fluid. |
1286 |
|
1287 |
As in the dipolar case, the quadpole polarizability evaluated from |
1288 |
both perturbation and fluctuation simulations show excellent |
1289 |
agreement, again confirming that equilibrium fluctuation calculations |
1290 |
are sufficient to reproduce bulk dielectric properties in these |
1291 |
fluids. The quadrupolar susceptibility calculated via our derived |
1292 |
correction factors tends to produce the same result for all three real |
1293 |
space methods. Similarly, the screening factor calculated using the |
1294 |
susceptibility and a weighted geometric factor provides good agreement |
1295 |
with results obtained directly via potentials of mean force. For |
1296 |
quadrupolar fluids, the distance dependence of the electrostatic |
1297 |
interaction is significantly reduced and the correction factors are |
1298 |
all small. These points suggest that how an electrostatic method |
1299 |
treats the cutoff radius become less consequential for higher order |
1300 |
multipoles. |
1301 |
|
1302 |
For this reason, we renew our recommendation that the |
1303 |
moderately-damped GSF method is an excellent choice for molecular |
1304 |
dynamics simulation where point-multipole interactions are being |
1305 |
utilized to compute bulk properties of fluids. |
1306 |
|
1307 |
\newpage |
1308 |
|
1309 |
\bibliography{dielectric_new} |
1310 |
|
1311 |
\end{document} |