| 1 |
\documentclass[% |
| 2 |
aip, |
| 3 |
jcp, |
| 4 |
amsmath,amssymb, |
| 5 |
preprint,% |
| 6 |
% reprint,% |
| 7 |
%author-year,% |
| 8 |
%author-numerical,% |
| 9 |
]{revtex4-1} |
| 10 |
|
| 11 |
\usepackage{graphicx}% Include figure files |
| 12 |
\usepackage{dcolumn}% Align table columns on decimal point |
| 13 |
\usepackage{multirow} |
| 14 |
\usepackage{bm}% bold math |
| 15 |
\usepackage{natbib} |
| 16 |
\usepackage{times} |
| 17 |
\usepackage{mathptmx} |
| 18 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 19 |
\usepackage{url} |
| 20 |
|
| 21 |
\begin{document} |
| 22 |
|
| 23 |
\title{Real space electrostatics for multipoles. III. Dielectric Properties} |
| 24 |
|
| 25 |
\author{Madan Lamichhane} |
| 26 |
\affiliation{Department of Physics, University |
| 27 |
of Notre Dame, Notre Dame, IN 46556} |
| 28 |
|
| 29 |
\author{Kathie E. Newman} |
| 30 |
\affiliation{Department of Physics, University |
| 31 |
of Notre Dame, Notre Dame, IN 46556} |
| 32 |
|
| 33 |
\author{Thomas Parsons} |
| 34 |
\affiliation{Department of Chemistry and Biochemistry, University |
| 35 |
of Notre Dame, Notre Dame, IN 46556} |
| 36 |
|
| 37 |
\author{J. Daniel Gezelter} |
| 38 |
\email{gezelter@nd.edu.} |
| 39 |
\affiliation{Department of Chemistry and Biochemistry, University |
| 40 |
of Notre Dame, Notre Dame, IN 46556} |
| 41 |
|
| 42 |
\date{\today}% It is always \today, today, |
| 43 |
% but any date may be explicitly specified |
| 44 |
|
| 45 |
\begin{abstract} |
| 46 |
Note: This manuscript is a work in progress. |
| 47 |
|
| 48 |
We report on the dielectric properties of the shifted potential |
| 49 |
(SP), gradient shifted force (GSF), and Taylor shifted force (TSF) |
| 50 |
real-space methods for multipole interactions that were developed in |
| 51 |
the first two papers in this series. We find that some subtlety is |
| 52 |
required for computing dielectric properties with the real-space |
| 53 |
methods, particularly when using the common fluctuation formulae. |
| 54 |
Three distinct methods for computing the dielectric constant are |
| 55 |
investigated, including the standard fluctuation formulae, |
| 56 |
potentials of mean force between solvated ions, and direct |
| 57 |
measurement of linear solvent polarization in response to applied |
| 58 |
fields and field gradients. |
| 59 |
\end{abstract} |
| 60 |
|
| 61 |
\maketitle |
| 62 |
|
| 63 |
\section{Introduction} |
| 64 |
|
| 65 |
Over the past several years, there has been increasing interest in |
| 66 |
pairwise methods for correcting electrostatic interactions in computer |
| 67 |
simulations of condensed molecular |
| 68 |
systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06} |
| 69 |
These techniques were developed from the observations and efforts of |
| 70 |
Wolf {\it et al.} and their work towards an $\mathcal{O}(N)$ |
| 71 |
Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is |
| 72 |
able to obtain excellent agreement with Madelung energies in ionic |
| 73 |
crystals.\cite{Wolf99} One of the most difficult tests of any new |
| 74 |
electrostatic method is the fidelity with which that method can |
| 75 |
reproduce the bulk-phase polarizability or equivalently, the |
| 76 |
dielectric properties of a fluid. Of particular interest is the |
| 77 |
static dielectric constant, $\epsilon$. Using the Ewald sum under |
| 78 |
tin-foil boundary conditions, $\epsilon$ can be calculated for systems |
| 79 |
of non-polarizable substances via |
| 80 |
\begin{equation} |
| 81 |
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
| 82 |
\label{eq:staticDielectric} |
| 83 |
\end{equation} |
| 84 |
where $\epsilon_0$ is the permittivity of free space and $\langle |
| 85 |
M^2\rangle$ is the fluctuation of the system dipole |
| 86 |
moment.\cite{Allen:1989fp} The numerator in the fractional term in equation |
| 87 |
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box |
| 88 |
dipole moment, identical to the quantity calculated in the |
| 89 |
finite-system Kirkwood $g$ factor ($G_k$): |
| 90 |
\begin{equation} |
| 91 |
G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
| 92 |
\label{eq:KirkwoodFactor} |
| 93 |
\end{equation} |
| 94 |
where $\mu$ is the dipole moment of a single molecule of the |
| 95 |
homogeneous system.\cite{Neumann:1983mz,Neumann:1983yq,Neumann84,Neumann85} The |
| 96 |
fluctuation term in both equation (\ref{eq:staticDielectric}) and |
| 97 |
(\ref{eq:KirkwoodFactor}) is calculated as follows, |
| 98 |
\begin{equation} |
| 99 |
\begin{split} |
| 100 |
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
| 101 |
- \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
| 102 |
&= \langle M_x^2+M_y^2+M_z^2\rangle |
| 103 |
- (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
| 104 |
+ \langle M_x\rangle^2). |
| 105 |
\end{split} |
| 106 |
\label{eq:fluctBoxDipole} |
| 107 |
\end{equation} |
| 108 |
This fluctuation term can be accumulated during a simulation; however, |
| 109 |
it converges rather slowly, thus requiring multi-nanosecond simulation |
| 110 |
times.\cite{Horn04} In the case of tin-foil boundary conditions, the |
| 111 |
dielectric/surface term of the Ewald summation is equal to zero. |
| 112 |
|
| 113 |
\section{Dielectric constants for real-space electrostatics} |
| 114 |
|
| 115 |
The new real-space methods require some care when computing dielectric |
| 116 |
properties because the interaction between molecular dipoles depends |
| 117 |
on the level of approximation being utilized. Consider the cases of |
| 118 |
Stockmayer (dipolar) soft spheres that are represented either by two |
| 119 |
closely-spaced point charges or by a single point dipole (see |
| 120 |
Fig. \ref{fig:stockmayer}). |
| 121 |
\begin{figure} |
| 122 |
\includegraphics[width=3in]{DielectricFigure} |
| 123 |
\caption{With the real-space electrostatic methods, the effective |
| 124 |
dipole tensor, $\mathbf{T}$, governing interactions between |
| 125 |
molecular dipoles is not the same for charge-charge interactions as |
| 126 |
for point dipoles.} |
| 127 |
\label{fig:stockmayer} |
| 128 |
\end{figure} |
| 129 |
In the case where point charges are interacting via an electrostatic |
| 130 |
kernel, $v(r)$, the effective {\it molecular} dipole tensor, |
| 131 |
$\mathbf{T}$ is obtained from two successive applications of the |
| 132 |
gradient operator to the electrostatic kernel, |
| 133 |
\begin{equation} |
| 134 |
\mathbf{T}_{\alpha \beta}(r) = \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta} |
| 135 |
\left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha} |
| 136 |
r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r} |
| 137 |
v^{\prime}(r) \right) |
| 138 |
\end{equation} |
| 139 |
where $v(r)$ may be either the bare kernel ($1/r$) or one of the |
| 140 |
modified (Wolf or DSF) kernels. This tensor describes the effective |
| 141 |
interaction between molecular dipoles ($\mathbf{D}$) in Gaussian |
| 142 |
units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$. |
| 143 |
|
| 144 |
When utilizing the new real-space methods for point dipoles, the |
| 145 |
tensor is explicitly constructed, |
| 146 |
\begin{equation} |
| 147 |
\mathbf{T}_{\alpha \beta}(r) = \delta_{\alpha \beta} v_{21}(r) + |
| 148 |
\frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r) |
| 149 |
\end{equation} |
| 150 |
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of |
| 151 |
the approximation. Although the Taylor-shifted (TSF) and |
| 152 |
gradient-shifted (GSF) models produce to the same $v(r)$ function for |
| 153 |
point charges, they have distinct forms for the dipole-dipole |
| 154 |
interactions. |
| 155 |
|
| 156 |
Neumann and Steinhauser showed~\cite{Neumann:1983mz,Neumann:1983yq} |
| 157 |
that the relative dielectric permittivity $\epsilon$ is given by the |
| 158 |
general fluctuation formula, |
| 159 |
\begin{equation} |
| 160 |
\frac{\epsilon - 1}{\epsilon +2} = \frac{4 \pi}{3} \frac{\left<M^2\right>}{3 V k_B |
| 161 |
T} \left( 1 - \frac{3}{4 \pi} \frac{\epsilon -1}{\epsilon + 2} |
| 162 |
\tilde{T}(0) \right) |
| 163 |
\end{equation} |
| 164 |
where $\left<M^2\right>$ is the mean fluctuation of the square of the |
| 165 |
net dipole moment of the simulation cell, $V$ is the volume of the |
| 166 |
cell, and $\tilde{T}(0) = \tilde{T̃}_{xx}(0) = \tilde{T̃}_{yy}(0) = |
| 167 |
\tilde{T̃}_{zz}(0)$ is the $\mathbf{k} = 0$ limit of Fourier transform |
| 168 |
of the diagonal term of the effective molecular dipolar tensor, |
| 169 |
\begin{equation} |
| 170 |
\tilde{T}(0) = \int_V \mathbf{T}(\mathbf{r}) d\mathbf{r} |
| 171 |
\end{equation} |
| 172 |
where the integral is carried out over the relevant geometry for the |
| 173 |
interaction. For the real-space methods, the integration can be done |
| 174 |
easily up to the imposed cutoff radius, while for the Ewald sum, the |
| 175 |
results also depend on details of the $k$-space |
| 176 |
calculation.\cite{Neumann:1983yq} |
| 177 |
|
| 178 |
For molecules composed of point charges interacting via the bare |
| 179 |
$(1/r)$ kernel, the simple cutoff (SC) and minimum image (MI) |
| 180 |
approaches give $\tilde{T}(0) = 0$ which means that the |
| 181 |
Clausius-Mosotti equation, |
| 182 |
\begin{equation} |
| 183 |
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon |
| 184 |
-1}{\epsilon+2} |
| 185 |
\end{equation} |
| 186 |
is the relevant fluctuation formula for the relative permittivity. |
| 187 |
|
| 188 |
Zahn {\it et al.}\cite{Zahn02} showed that for the damped shifted |
| 189 |
force charge-charge kernel, $\tilde{T}(0) = 4 \pi / 3$. This was later |
| 190 |
generalized by Izvekov {\it et al.} for all point-charge kernels which |
| 191 |
have forces that go to zero at a cutoff radius and which maintain a |
| 192 |
pole of first order at $r=0$.\cite{Izvekov:2008wo} When $\tilde{T}(0) |
| 193 |
= 4 \pi/ 3$, the expression for the dielectric constant reduces to the |
| 194 |
widely-used {\it conducting boundary} formula, |
| 195 |
\begin{equation} |
| 196 |
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon - |
| 197 |
1}{3} |
| 198 |
\end{equation} |
| 199 |
It is convenient to define a quantity $Q = \frac{3}{4 \pi} |
| 200 |
\tilde{T}(0)$, and to combine all of the fluctuation formulae |
| 201 |
together: |
| 202 |
\begin{equation} |
| 203 |
\frac{4 \pi}{3} \frac{\left< M^2 \right>}{3 V k_B T} = |
| 204 |
\left\{\begin{array}{ll} |
| 205 |
\frac{\epsilon-1}{\epsilon+2} \left[1- \frac{\epsilon-1}{\epsilon+2} |
| 206 |
Q\right]^{-1} & \mathrm{General~Case} \\ |
| 207 |
\frac{\epsilon-1}{\epsilon+2} & Q \rightarrow 0 \mathrm{~limit} \\ |
| 208 |
\frac{\epsilon-1}{3} & Q \rightarrow 1 \mathrm{~limit} |
| 209 |
\end{array}\right. |
| 210 |
\end{equation} |
| 211 |
|
| 212 |
The Clausius-Mossotti $(Q\rightarrow 0)$ approach is subject to noise |
| 213 |
and error magnification,\cite{Allen:1989fp} and the conducting |
| 214 |
boundary approach $(Q \rightarrow 1)$ has become widely used as a |
| 215 |
result. If the electrostatic method being used has $Q < 1$, the |
| 216 |
relative dielectric permittivity $\epsilon$ can be estimated once the |
| 217 |
conducting boundary fluctuation result has been found, |
| 218 |
\begin{equation} |
| 219 |
\epsilon = \frac{(Q+2) (\epsilon_\text{CB}-1)+3} {(Q-1) (\epsilon_\text{CB}-1)+3} |
| 220 |
\end{equation} |
| 221 |
Note that this expression becomes quite numerically sensitive when the |
| 222 |
value of $Q$ deviates significantly from 1. |
| 223 |
|
| 224 |
We have derived a set of expressions for the value of $Q$ for the new |
| 225 |
real space methods that are shown in Table \ref{tab:Q}. |
| 226 |
|
| 227 |
\begin{table} |
| 228 |
\caption{Expressions for the dielectric correction factor ($Q$) for the |
| 229 |
real-space electrostatic methods in terms of the damping parameter |
| 230 |
($\alpha$) and the cutoff radius ($r_c$). The Ewald-Kornfeld result |
| 231 |
derived in Refs. \onlinecite{Adams:1980rt,Adams:1981fr,Neumann:1983yq} is shown for comparison using the Ewald |
| 232 |
convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). } |
| 233 |
\label{tab:Q} |
| 234 |
\begin{tabular}{l|c|c|c|c|} |
| 235 |
& \multicolumn{2}{c|}{damped} & \multicolumn{2}{c|}{undamped} \\ \cline{2-3} \cline{4-5} |
| 236 |
Method & $Q_\mathrm{charges}$ & $Q_\mathrm{dipoles}$ & $Q_\mathrm{charges}$ & $Q_\mathrm{dipoles}$ \\ |
| 237 |
\hline |
| 238 |
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}$ &0 &0\\ |
| 239 |
Shifted Potental (SP) & $ \mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(\alpha r_c)-\frac{2 \alpha r_c }{\sqrt{\pi }} \left(1 + \frac{2 \alpha^2 r_c^2}{3} \right) e^{-\alpha ^2 r_c^2} $ &0 &0\\ |
| 240 |
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} $ &1 &1\\ |
| 241 |
Taylor-shifted (TSF) & 1 & 1 & 1 & 1\\ |
| 242 |
Ewald-Kornfeld (EK) & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ & - & - \\\hline |
| 243 |
\end{tabular} |
| 244 |
\end{table} |
| 245 |
|
| 246 |
\newpage |
| 247 |
|
| 248 |
\appendix |
| 249 |
\section{Point-multipolar interactions with a spatially-varying electric field} |
| 250 |
|
| 251 |
We can treat objects $a$, $b$, and $c$ containing embedded collections |
| 252 |
of charges. When we define the primitive moments, we sum over that |
| 253 |
collections of charges using a local coordinate system within each |
| 254 |
object. The point charge, dipole, and quadrupole for object $a$ are |
| 255 |
given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively. |
| 256 |
These are the primitive multipoles which can be expressed as a |
| 257 |
distribution of charges, |
| 258 |
\begin{align} |
| 259 |
C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\ |
| 260 |
D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\ |
| 261 |
Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in } a} q_k |
| 262 |
r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} |
| 263 |
\end{align} |
| 264 |
Note that the definition of the primitive quadrupole here differs from |
| 265 |
the standard traceless form, and contains an additional Taylor-series |
| 266 |
based factor of $1/2$. In Paper 1, we derived the forces and torques |
| 267 |
each object exerts on the others. |
| 268 |
|
| 269 |
Here we must also consider an external electric field that varies in |
| 270 |
space: $\mathbf E(\mathbf r)$. Each of the local charges $q_k$ in |
| 271 |
object $a$ will then experience a slightly different field. This |
| 272 |
electric field can be expanded in a Taylor series around the local |
| 273 |
origin of each object. A different Taylor series expansion is carried |
| 274 |
out for each object. |
| 275 |
|
| 276 |
For a particular charge $q_k$, the electric field at that site's |
| 277 |
position is given by: |
| 278 |
\begin{equation} |
| 279 |
E_\gamma + \nabla_\delta E_\gamma r_{k \delta} |
| 280 |
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta} |
| 281 |
r_{k \varepsilon} + ... |
| 282 |
\end{equation} |
| 283 |
Note that the electric field is always evaluated at the origin of the |
| 284 |
objects, and treating each object using point multipoles simplifies |
| 285 |
this greatly. |
| 286 |
|
| 287 |
To find the force exerted on object $a$ by the electric field, one |
| 288 |
takes the electric field expression, and multiplies it by $q_k$, and |
| 289 |
then sum over all charges in $a$: |
| 290 |
|
| 291 |
\begin{align} |
| 292 |
F_\gamma &= \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta} |
| 293 |
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta} |
| 294 |
r_{k \varepsilon} + ... \rbrace \\ |
| 295 |
&= C_a E_\gamma + D_{a \delta} \nabla_\delta E_\gamma |
| 296 |
+ Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma + |
| 297 |
... |
| 298 |
\end{align} |
| 299 |
|
| 300 |
Similarly, the torque exerted by the field on $a$ can be expressed as |
| 301 |
\begin{align} |
| 302 |
\tau_\alpha &= \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\ |
| 303 |
& = \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k |
| 304 |
r_{k\beta} E_\gamma(\mathbf r_k) \\ |
| 305 |
& = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma |
| 306 |
+ 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta |
| 307 |
E_\gamma + ... |
| 308 |
\end{align} |
| 309 |
|
| 310 |
The last term is essentially identical with form derived by Torres del |
| 311 |
Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation |
| 312 |
utilized a traceless form of the quadrupole that is different than the |
| 313 |
primitive definition in use here. We note that the Levi-Civita symbol |
| 314 |
can be eliminated by utilizing the matrix cross product in an |
| 315 |
identical form as in Ref. \onlinecite{Smith98}: |
| 316 |
\begin{equation} |
| 317 |
\left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta |
| 318 |
\left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta} |
| 319 |
-\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta} |
| 320 |
\right] |
| 321 |
\label{eq:matrixCross} |
| 322 |
\end{equation} |
| 323 |
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of |
| 324 |
the matrix indices. In table \ref{tab:UFT} we give compact |
| 325 |
expressions for how the multipole sites interact with an external |
| 326 |
field that has exhibits spatial variations. |
| 327 |
|
| 328 |
\begin{table} |
| 329 |
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque |
| 330 |
$(\mathbf{\tau})$ expressions for a multipolar site embedded in an |
| 331 |
electric field with spatial variations, $\mathbf{E}(\mathbf{r})$. |
| 332 |
\label{tab:UFT}} |
| 333 |
\begin{tabular}{r|ccc} |
| 334 |
& Charge & Dipole & Quadrupole \\ \hline |
| 335 |
$U$ & $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\ |
| 336 |
$\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ & $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\ |
| 337 |
$\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$ |
| 338 |
\end{tabular} |
| 339 |
\end{table} |
| 340 |
|
| 341 |
|
| 342 |
\newpage |
| 343 |
|
| 344 |
\bibliography{multipole} |
| 345 |
|
| 346 |
\end{document} |
