trunk/multipole/dielectric.tex

\documentclass[%
 aip,
 jcp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{multirow}
\usepackage{bm}% bold math
\usepackage{natbib}
\usepackage{times}
\usepackage{mathptmx}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}

\begin{document}

\title{Real space electrostatics for multipoles. III. Dielectric Properties}

\author{Madan Lamichhane}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}

\author{Kathie E. Newman}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}

\author{Thomas Parsons}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\author{J. Daniel Gezelter}
 \email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
  Note: This manuscript is a work in progress.   

  We report on the dielectric properties of the shifted potential
  (SP), gradient shifted force (GSF), and Taylor shifted force (TSF)
  real-space methods for multipole interactions that were developed in
  the first two papers in this series.  We find that some subtlety is
  required for computing dielectric properties with the real-space
  methods, particularly when using the common fluctuation formulae.
  Three distinct methods for computing the dielectric constant are
  investigated, including the standard fluctuation formulae,
  potentials of mean force between solvated ions, and direct
  measurement of linear solvent polarization in response to applied
  fields and field gradients.
\end{abstract}

\maketitle

\section{Introduction}

Over the past several years, there has been increasing interest in
pairwise methods for correcting electrostatic interactions in computer
simulations of condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06}
These techniques were developed from the observations and efforts of
Wolf {\it et al.}  and their work towards an $\mathcal{O}(N)$
Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is
able to obtain excellent agreement with Madelung energies in ionic
crystals.\cite{Wolf99} One of the most difficult tests of any new
electrostatic method is the fidelity with which that method can
reproduce the bulk-phase polarizability or equivalently, the
dielectric properties of a fluid.  Of particular interest is the
static dielectric constant, $\epsilon$.  Using the Ewald sum under
tin-foil boundary conditions, $\epsilon$ can be calculated for systems
of non-polarizable substances via
\begin{equation}
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV},
\label{eq:staticDielectric}
\end{equation}
where $\epsilon_0$ is the permittivity of free space and $\langle
M^2\rangle$ is the fluctuation of the system dipole
moment.\cite{Allen:1989fp} The numerator in the fractional term in equation
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box
dipole moment, identical to the quantity calculated in the
finite-system Kirkwood $g$ factor ($G_k$):
\begin{equation}
G_k = \frac{\langle M^2\rangle}{N\mu^2},
\label{eq:KirkwoodFactor}
\end{equation}
where $\mu$ is the dipole moment of a single molecule of the
homogeneous system.\cite{Neumann:1983mz,Neumann:1983yq,Neumann84,Neumann85} The
fluctuation term in both equation (\ref{eq:staticDielectric}) and
(\ref{eq:KirkwoodFactor}) is calculated as follows,
\begin{equation}
\begin{split}
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle 
                        - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\
                   &= \langle M_x^2+M_y^2+M_z^2\rangle
                        - (\langle M_x\rangle^2 + \langle M_x\rangle^2 
                                + \langle M_x\rangle^2).        
\end{split}
\label{eq:fluctBoxDipole}
\end{equation}
This fluctuation term can be accumulated during a simulation; however,
it converges rather slowly, thus requiring multi-nanosecond simulation
times.\cite{Horn04} In the case of tin-foil boundary conditions, the
dielectric/surface term of the Ewald summation is equal to zero.

\section{Dielectric constants for real-space electrostatics}

The new real-space methods require some care when computing dielectric
properties because the interaction between molecular dipoles depends
on the level of approximation being utilized. Consider the cases of
Stockmayer (dipolar) soft spheres that are represented either by two
closely-spaced point charges or by a single point dipole (see
Fig. \ref{fig:stockmayer}).
\begin{figure}
\includegraphics[width=3in]{DielectricFigure}
\caption{With the real-space electrostatic methods, the effective
  dipole tensor, $\mathbf{T}$, governing interactions between
  molecular dipoles is not the same for charge-charge interactions as
  for point dipoles.}
\label{fig:stockmayer}
\end{figure}
In the case where point charges are interacting via an electrostatic
kernel, $v(r)$, the effective {\it molecular} dipole tensor,
$\mathbf{T}$ is obtained from two successive applications of the
gradient operator to the electrostatic kernel,
\begin{equation}
\mathbf{T}_{\alpha \beta}(r) =  \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta}
\left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha}
  r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r}
  v^{\prime}(r) \right)
\end{equation}
where $v(r)$ may be either the bare kernel ($1/r$) or one of the
modified (Wolf or DSF) kernels.  This tensor describes the effective
interaction between molecular dipoles ($\mathbf{D}$) in Gaussian
units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$.

When utilizing the new real-space methods for point dipoles, the
tensor is explicitly constructed,
\begin{equation}
\mathbf{T}_{\alpha \beta}(r)  =  \delta_{\alpha \beta} v_{21}(r) +
\frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r) 
\end{equation}
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of
the approximation. Although the Taylor-shifted (TSF) and
gradient-shifted (GSF) models produce to the same $v(r)$ function for
point charges, they have distinct forms for the dipole-dipole
interactions.

Neumann and Steinhauser showed~\cite{Neumann:1983mz,Neumann:1983yq}
that the relative dielectric permittivity $\epsilon$ is given by the
general fluctuation formula,
\begin{equation}
\frac{\epsilon - 1}{\epsilon +2} = \frac{4 \pi}{3} \frac{\left<M^2\right>}{3 V k_B
  T} \left( 1 - \frac{3}{4 \pi} \frac{\epsilon -1}{\epsilon + 2}
  \tilde{T}(0) \right)
\end{equation}
where $\left<M^2\right>$ is the mean fluctuation of the square of the
net dipole moment of the simulation cell, $V$ is the volume of the
cell, and $\tilde{T}(0) = \tilde{T̃}_{xx}(0) = \tilde{T̃}_{yy}(0) =
\tilde{T̃}_{zz}(0)$ is the $\mathbf{k} = 0$ limit of Fourier transform
of the diagonal term of the effective molecular dipolar tensor,
\begin{equation}
\tilde{T}(0) = \int_V \mathbf{T}(\mathbf{r})  d\mathbf{r}
\end{equation}
where the integral is carried out over the relevant geometry for the
interaction.  For the real-space methods, the integration can be done
easily up to the imposed cutoff radius, while for the Ewald sum, the
results also depend on details of the $k$-space
calculation.\cite{Neumann:1983yq}

For molecules composed of point charges interacting via the bare
$(1/r)$ kernel, the simple cutoff (SC) and minimum image (MI)
approaches give $\tilde{T}(0) = 0$ which means that the
Clausius-Mosotti equation,
\begin{equation}
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon
  -1}{\epsilon+2} 
\end{equation}
is the relevant fluctuation formula for the relative permittivity.

Zahn {\it et al.}\cite{Zahn02} showed that for the damped shifted
force charge-charge kernel, $\tilde{T}(0) = 4 \pi / 3$. This was later
generalized by Izvekov {\it et al.} for all point-charge kernels which
have forces that go to zero at a cutoff radius and which maintain a
pole of first order at $r=0$.\cite{Izvekov:2008wo} When $\tilde{T}(0)
= 4 \pi/ 3$, the expression for the dielectric constant reduces to the
widely-used {\it conducting boundary} formula,
\begin{equation}
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon -
  1}{3}
\end{equation}
It is convenient to define a quantity $Q = \frac{3}{4 \pi}
\tilde{T}(0)$, and to combine all of the fluctuation formulae
together:
\begin{equation}
\frac{4 \pi}{3} \frac{\left< M^2 \right>}{3 V k_B T} =
\left\{\begin{array}{ll}
\frac{\epsilon-1}{\epsilon+2} \left[1- \frac{\epsilon-1}{\epsilon+2}
  Q\right]^{-1} & \mathrm{General~Case} \\
\frac{\epsilon-1}{\epsilon+2} & Q \rightarrow 0 \mathrm{~limit} \\
\frac{\epsilon-1}{3} & Q \rightarrow 1 \mathrm{~limit}
\end{array}\right.
\end{equation}

The Clausius-Mossotti $(Q\rightarrow 0)$ approach is subject to noise
and error magnification,\cite{Allen:1989fp} and the conducting
boundary approach $(Q \rightarrow 1)$ has become widely used as a
result.  If the electrostatic method being used has $Q < 1$, the
relative dielectric permittivity $\epsilon$ can be estimated once the
conducting boundary fluctuation result has been found,
\begin{equation}
\epsilon = \frac{(Q+2) (\epsilon_\text{CB}-1)+3} {(Q-1) (\epsilon_\text{CB}-1)+3}
\end{equation}
Note that this expression becomes quite numerically sensitive when the
value of $Q$ deviates significantly from 1.

We have derived a set of expressions for the value of $Q$ for the new
real space methods that are shown in Table \ref{tab:Q}.

\begin{table}
  \caption{Expressions for the dielectric correction factor ($Q$) for the
    real-space electrostatic methods in terms of the damping parameter
    ($\alpha$) and the cutoff radius ($r_c$).  The Ewald-Kornfeld result 
    derived in Refs. \onlinecite{Adams:1980rt,Adams:1981fr,Neumann:1983yq} is shown for comparison using the Ewald 
    convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). }
\label{tab:Q}
\begin{tabular}{l|c|c|c|c|}
       & \multicolumn{2}{c|}{damped} & \multicolumn{2}{c|}{undamped} \\ \cline{2-3} \cline{4-5}
Method & $Q_\mathrm{charges}$ & $Q_\mathrm{dipoles}$ & $Q_\mathrm{charges}$ & $Q_\mathrm{dipoles}$ \\
\hline
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\\
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\\
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\\
Taylor-shifted  (TSF) & 1 & 1 & 1 & 1\\ 
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
\end{tabular}
\end{table}

\newpage

\appendix
\section{Point-multipolar interactions with a spatially-varying electric field}

We can treat objects $a$, $b$, and $c$ containing embedded collections
of charges. When we define the primitive moments, we sum over that
collections of charges using a local coordinate system within each
object.  The point charge, dipole, and quadrupole for object $a$ are
given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively.
These are the primitive multipoles which can be expressed as a
distribution of charges,
\begin{align}
C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\
D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\
Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in }  a} q_k
r_{k\alpha}  r_{k\beta} . \label{eq:quadrupole}
\end{align}
Note that the definition of the primitive quadrupole here differs from
the standard traceless form, and contains an additional Taylor-series
based factor of $1/2$.  In Paper 1, we derived the forces and torques
each object exerts on the others.

Here we must also consider an external electric field that varies in
space: $\mathbf E(\mathbf r)$.  Each of the local charges $q_k$ in
object $a$ will then experience a slightly different field.  This
electric field can be expanded in a Taylor series around the local
origin of each object.  A different Taylor series expansion is carried
out for each object.

For a particular charge $q_k$, the electric field at that site's
position is given by:
\begin{equation}
E_\gamma + \nabla_\delta E_\gamma r_{k \delta} 
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
r_{k \varepsilon} + ... 
\end{equation}
Note that the electric field is always evaluated at the origin of the
objects, and treating each object using point multipoles simplifies
this greatly.

To find the force exerted on object $a$ by the electric field, one
takes the electric field expression, and multiplies it by $q_k$, and
then sum over all charges in $a$:

\begin{align}
F_\gamma &=  \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta} 
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
r_{k \varepsilon} + ...  \rbrace  \\
 &= C_a E_\gamma + D_{a  \delta} \nabla_\delta E_\gamma 
+ Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma +
... 
\end{align}

Similarly, the torque exerted by the field on $a$ can be expressed as
\begin{align}
\tau_\alpha &=  \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\
 & =  \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k
 r_{k\beta} E_\gamma(\mathbf r_k) \\
 & = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma 
+ 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta
E_\gamma + ...
\end{align}

The last term is essentially identical with form derived by Torres del
Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation
utilized a traceless form of the quadrupole that is different than the
primitive definition in use here.  We note that the Levi-Civita symbol
can be eliminated by utilizing the matrix cross product in an
identical form as in Ref. \onlinecite{Smith98}:
\begin{equation}
\left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta
\left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta}
  -\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta} 
\right]
\label{eq:matrixCross}
\end{equation}
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of
the matrix indices.  In table \ref{tab:UFT} we give compact
expressions for how the multipole sites interact with an external
field that has exhibits spatial variations.

\begin{table}
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
  $(\mathbf{\tau})$ expressions for a multipolar site embedded in an
  electric field with spatial variations, $\mathbf{E}(\mathbf{r})$.
\label{tab:UFT}}
\begin{tabular}{r|ccc}
  & Charge & Dipole & Quadrupole \\ \hline
$U$ &  $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\
$\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ &  $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\
$\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$
\end{tabular}
\end{table}


\newpage

\bibliography{multipole}

\end{document}
Revision:	4218
Committed:	Thu Sep 11 18:08:53 2014 UTC (10 years ago) by gezelter
Content type:	application/x-tex
File size:	15816 byte(s)
Log Message:	Changes
#	User	Rev	Content
1	gezelter	4195	\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	gezelter	4218	\title{Real space electrostatics for multipoles. III. Dielectric Properties}
24	gezelter	4195
25	gezelter	4218	\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	gezelter	4195	\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	gezelter	4218	\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	gezelter	4195	\maketitle
62
63	gezelter	4218	\section{Introduction}
64	gezelter	4195
65	gezelter	4218	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	gezelter	4195
113	gezelter	4218	\section{Dielectric constants for real-space electrostatics}
114
115	gezelter	4195	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	gezelter	4196	interaction between molecular dipoles ($\mathbf{D}$) in Gaussian
142	gezelter	4195	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	gezelter	4218
246	gezelter	4195	\newpage
247
248	gezelter	4218	\appendix
249			\section{Point-multipolar interactions with a spatially-varying electric field}
250	gezelter	4195
251	gezelter	4218	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	gezelter	4195	\bibliography{multipole}
345
346			\end{document}