trunk/multipole/dielectric_new.tex

\documentclass[aip,jcp,amsmath,amssymb,preprint]{revtex4-1}

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

\newcolumntype{Y}{>{\centering\arraybackslash}X}

\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{Thomas Parsons}
\affiliation{Department of Chemistry and Biochemistry, 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{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}
  In the first two papers in this series, we developed new shifted
  potential (SP), gradient shifted force (GSF), and Taylor shifted
  force (TSF) real-space methods for multipole interactions in
  condensed phase simulations.  Here, we discuss the dielectric
  properties of fluids that emerge from simulations using these
  methods.  Most electrostatic methods (including the Ewald sum)
  require correction to the conducting boundary fluctuation formula
  for the static dielectric constants, and we discuss the derivation
  of these corrections for the new real space methods. For quadrupolar
  fluids, the analogous material property is the quadrupolar
  susceptibility. As in the dipolar case, the fluctuation formula for
  the quadrupolar susceptibility has corrections that depend on the
  electrostatic method being utillized. One of the most important
  effects measured by both the static dielectric and quadrupolar
  susceptibility is the ability to screen charges embedded in the
  fluid.  We use potentials of mean force between solvated ions to
  discuss how geometric factors can lead to distance-dependent
  screening in both quadrupolar and dipolar fluids.
\end{abstract}

\maketitle

\section{Introduction}

Over the past several years, there has been increasing interest in
pairwise or ``real space'' methods for computing electrostatic
interactions in condensed phase
simulations.\cite{Wolf99,Zahn02,Kast03,Beckd05,Ma05,Wu:2005nr,Fennell06,Fukuda:2013uq,Wang:2016kx}
These techniques were initially developed by Wolf {\it et al.}  in
their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99}
Wolf's method of using cutoff neutralization and electrostatic damping
is able to obtain excellent agreement with Madelung energies in ionic
crystals.\cite{Wolf99}

Zahn \textit{et al.}\cite{Zahn02} and Fennell and Gezelter extended
this method using shifted force approximations at the cutoff distance
in order to conserve total energy in molecular dynamics
simulations.\cite{Fennell06} Other recent advances in real-space
methods for systems of point charges have included explicit
elimination of the net multipole moments inside the cutoff sphere
around each charge site.\cite{Fukuda:2013uq,Wang:2016kx}

In the previous two papers in this series, we developed three
generalized real space methods: shifted potential (SP), gradient
shifted force (GSF), and Taylor shifted force (TSF).\cite{PaperI,
  PaperII} These methods evaluate electrostatic interactions for
charges and higher order multipoles using a finite-radius cutoff
sphere.  The neutralization and damping of local moments within the
cutoff sphere is a multipolar generalization of Wolf's sum.  In the
GSF and TSF methods, additional terms are added to the potential
energy so that forces and torques also vanish smoothly at the cutoff
radius.  This ensures that the total energy is conserved in molecular
dynamics simulations.

One of the most stringent 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. Before the advent of computer simulations, Kirkwood and Onsager
developed fluctuation formulae for the dielectric properties of
dipolar fluids.\cite{Kirkwood39,Onsagar36} Along with projections of
the frequency-dependent dielectric to zero frequency, these
fluctuation formulae are now widely used to predict the static
dielectric constant of simulated materials.

If we consider a system of dipolar or quadrupolar molecules under the
influence of an external field or field gradient, the net polarization
of the system will largely be proportional to the applied
perturbation.\cite{Chitanvis96, Stern-Feller03, SalvchovI14,
  SalvchovII14} In simulations, the net polarization of the system is
also determined by the interactions \textit{between} the
molecules. Therefore the macroscopic polarizablity obtained from a
simulation depends on the details of the electrostatic interaction
methods that were employed in the simulation. To determine the
relevant physical properties of the multipolar fluid from the system
fluctuations, the interactions between molecules must be incorporated
into the formalism for the bulk properties.

In most simulations, bulk materials are treated using periodic
replicas of small regions, and this level of approximation requires
corrections to the fluctuation formulae that were derived for the bulk
fluids. In 1983 Neumann proposed a general formula for evaluating
dielectric properties of dipolar fluids using both Ewald and
real-space cutoff methods.\cite{NeumannI83} Steinhauser and Neumann
used this formula to evaluate the corrected dielectric constant for
the Stockmayer fluid using two different methods: Ewald-Kornfield (EK)
and reaction field (RF) methods.\cite{NeumannII83}

Zahn \textit{et al.}\cite{Zahn02} utilized this approach and evaluated
the correction factor for using damped shifted charge-charge
kernel. This was later generalized by Izvekov \textit{et
  al.},\cite{Izvekov08} who noted that the expression for the
dielectric constant reduces to widely-used conducting boundary formula
for real-space methods that have first derivatives that vanish at the
cutoff radius.

One of the primary topics of this paper is the derivation of
correction factors for the three new real space methods.  We find that
the correction formulae for dipolar molecules depends not only on the
methodology being used, but also on whether the molecular dipoles are
treated using point charges or point dipoles.  We derive correction
factors for both cases.

In quadrupolar fluids, the relationship between quadrupolar
susceptibility and the dielectric constant is not as straightforward
as in the dipolar case.  The effective dielectric constant depends on
the geometry of the external (or internal) field
perturbation.\cite{Ernst92} Significant efforts have been made to
increase our understanding the dielectric properties of these
fluids,\cite{JeonI03,JeonII03,Chitanvis96} although a general
correction formula has not yet been developed.

In this paper we derive general formulae for calculating the
quadrupolar susceptibility of quadrupolar fluids. We also evaluate the
correction factor for SP, GSF, and TSF methods for quadrupolar fluids
interacting via point charges, point dipoles or directly through
quadrupole-quadrupole interactions.

We also calculate the screening behavior for two ions immersed in
multipolar fluids to estimate the distance dependence of charge
screening in both dipolar and quadrupolar fluids.  We use three
distinct methods to compare our analytical results with computer
simulations:
\begin{enumerate}
\item responses of the fluid to external perturbations,
\item fluctuations of system multipole moments, and
\item potentials of mean force between solvated ions,
\end{enumerate}

Under the influence of weak external fields, the bulk polarization of
the system is primarily a linear response to the perturbation, where
proportionality constant depends on the electrostatic interactions
between the multipoles. The fluctuation formulae connect bulk
properties of the fluid to equilibrium fluctuations in the system
multipolar moments during a simulation. These fluctuations also depend
on the form of the electrostatic interactions between molecules.
Therefore, the connections between the actual bulk properties and both
the computed fluctuation and external field responses must be modified
accordingly.

The potential of mean force (PMF) allows calculation of an effective
dielectric constant or screening factor from the potential energy
between ions before and after dielectric material is introduced.
Computing the PMF between embedded point charges is an additional
check on the bulk properties computed via the other two methods.

\section{The Real-Space Methods}

\section{Dipolar Fluids and the Dielectric Constant}

Dielectric properties of a fluid arise mainly from responses of the
fluid to either applied fields or transient fields internal to the
fluid. In response to an applied field, the molecules have electronic
polarizabilities, changes to internal bond lengths, and reorientations
towards the direction of the applied field. There is an added
complication that in the presence of external field, the perturbation
experienced by any single molecule is not only due to external field
but also to the fields produced by the all other molecules in the
system.

\subsection{Response to External Perturbations}

In the presence of uniform electric field $\mathbf{E}$, an individual
molecule with a permanent dipole moment $p_o$ will realign along the
direction of the field with an average polarization given by
\begin{equation}
\braket{\mathbf{p}} = \epsilon_0 \alpha_p \mathbf{E},
\end{equation}
where $\alpha_p = {p_o}^2 / 3 \epsilon_0 k_B T$ is the contribution to
molecular polarizability due solely to reorientation dynamics.
Because the applied field must overcome thermal motion, the
orientational polarization depends inversely on the temperature.

Likewise, a condensed phase system of permanent dipoles will also
polarize along the direction of an applied field. The polarization
density of the system is
\begin{equation}
\textbf{P} = \epsilon_o \alpha_{D} \mathbf{E}. 
\label{pertDipole}
\end{equation} 
The constant $\alpha_D$ is the macroscopic polarizability, which is an
emergent property of the dipolar fluid.  Note that the polarizability,
$\alpha_D$ is distinct from the dipole susceptibility, $\chi_D$,
which is the quantity most directly related to the static dielectric
constant, $\epsilon = 1 + \chi_D$.

\subsection{Fluctuation Formula}
% Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can
% be calculated for systems of orientationally-polarized molecules,
% \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{Allen89} The numerator in the fractional term in
% equation (\ref{eq:staticDielectric}) is 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{NeumannI83,NeumannII83,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.

For a system of dipolar molecules at thermal equilibrium, we can
define both a system dipole moment, $\mathbf{M} = \sum_i \mathbf{p}_i$
as well as a dipole polarization density,
$\mathbf{P} = \braket{\mathbf{M}}/V$.

In the presence of applied field $\mathbf{E}$, the polarization
density can be expressed in terms of fluctuations in the net dipole
moment,
\begin{equation}
\mathbf{P} = \epsilon_o \frac{\braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E}
\label{flucDipole}
\end{equation}
This has structural similarity with the Boltzmann average for the
polarization of a single molecule. Here
$ \braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2$ measures fluctuations
in the net dipole moment,
\begin{equation}
 \langle \mathbf{M}^2 \rangle - \langle \mathbf{M} \rangle^2 =
 \langle M_x^2+M_y^2+M_z^2 \rangle - \left( \langle M_x \rangle^2 +
   \langle M_y \rangle^2 +
   \langle M_z \rangle^2 \right).
\label{eq:flucDip}
\end{equation}
For the limiting case $\textbf{E} \rightarrow 0 $, the ensemble
average of both the net dipole moment $\braket{\mathbf{M}}$ and
dipolar polarization $\bf{P}$ tends to vanish but
$\braket{\mathbf{M}^2}$ does not.  The macroscopic dipole
polarizability can therefore be written as,
\begin{equation}
  \alpha_D = \frac{\braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2}{3 \epsilon_o V k_B T}
\end{equation}
This relationship between a macroscopic property of dipolar material
and microscopic fluctuations is true even when the applied field
$ \textbf{E} \rightarrow 0 $.

\subsection{Correction Factors}

In the presence of a uniform external field $ \mathbf{E}^\circ$, the
total electric field at $\mathbf{r}$ depends on the polarization
density at all other points in the system,\cite{NeumannI83}
\begin{equation}
\mathbf{E}(\mathbf{r}) = \mathbf{E}^\circ(\mathbf{r}) +
\frac{1}{4\pi\epsilon_o} \int d\mathbf{r}^\prime
\mathbf{T}(\mathbf{r}-\mathbf{r}^\prime)\cdot
{\mathbf{P}(\mathbf{r}^\prime)}.
\label{eq:localField}
\end{equation}
$\mathbf{T}$ is the dipole interaction tensor connecting dipoles at
$\mathbf{r}^\prime$ with the point of interest ($\mathbf{r}$).

In simulations of dipolar fluids, the molecular dipoles may be
represented either by closely-spaced point charges or by 
point dipoles (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{align}
\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)
\label{dipole-chargeTensor}
\end{align}
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 any of the three new real-space methods for point
\textit{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) 
\label{dipole-diopleTensor}
\end{equation}
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of
the approximation.\cite{PaperI,PaperII} 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 interaction.
 
Using a constitutive relation, the polarization density
$\mathbf{P}(\mathbf{r})$ is given by,
\begin{equation}
\mathbf{P}(\mathbf{r}) = \epsilon_o \chi_D
\left(\mathbf{E}^\circ(\mathbf{r}) + \frac{1}{4\pi\epsilon_o} \int
  d\mathbf{r}^\prime \mathbf{T}(\mathbf{r}-\mathbf{r}^\prime ) \cdot \mathbf{P}(\mathbf{r}^\prime)\right).
\label{constDipole}
\end{equation}
Note that $\chi_D$ explicitly depends on the details of the dipole
interaction tensor.  Neumann \textit{et al.}
\cite{NeumannI83,NeumannII83,Neumann84,Neumann85} derived an elegant
way to modify the fluctuation formula to correct for approximate
interaction tensors. This correction was derived using a Fourier
representation of the interaction tensor,
$\tilde{\mathbf{T}}(\mathbf{k})$, and involves the quantity,
\begin{equation}
 A = \frac{3}{4\pi}\tilde{\mathbf{T}}(0) = \frac{3}{4\pi} \int_V
 d\mathbf{r} \mathbf{T}(\mathbf{r})
\end{equation}
which is the $k \rightarrow 0$ limit of $\tilde{\mathbf{T}}$.  Using
this quantity (originally called $Q$ in
refs. \onlinecite{NeumannI83,NeumannII83,Neumann84,Neumann85}), the
dielectric constant can be computed
\begin{equation}
\epsilon = \frac{3+(A + 2)(\epsilon_{CB}-1)}{3 + (A -1)(\epsilon_{CB}-1)}
\label{correctionFormula}
\end{equation}
where $\epsilon_{CB}$ is the widely-used conducting boundary
expression for the dielectric constant,
\begin{equation}
\epsilon_{CB} = 1 + \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3
  \epsilon_o V k_B T}
\label{correctionFormula}
\end{equation}
Equation (\ref{correctionFormula}) allows estimation of the static
dielectric constant from fluctuations easily computed directly from
simulations.


% 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}
% \begin{equation}
% \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}
% \label{singular}
% \end{equation}
% Substituting equation (\ref{singular}) in the equation (\ref{constDipole}) we get,
% \begin{equation}
% \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).
% \end{equation}
% For both polarization and electric field homogeneous, this can be easily solved using Fourier transformation, 
% \begin{equation}
% \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}).
% \end{equation}
% For homogeneous applied field Fourier component is non-zero only if $\kappa = 0$. Therefore,
% \begin{equation}
% \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}).
% \label{fourierDipole}
% \end{equation}
% 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,
% \begin{equation}
% \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}
% \end{equation}
% 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,
% \begin{equation}
% \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}
% \label{correctionFormula}
% \end{equation}


We have utilized the Neumann \textit{et al.} approach for thre three
new real-space methods, and obtain method-dependent correction
factors.  The expression for the correction factor also depends on
whether the simulation involves point charges or point dipoles to
represent the molecular dipoles.  These corrections factors are 
listed in Table \ref{tab:A}.
\begin{table}
  \caption{Expressions for the dipolar correction factor ($A$) 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{Adams80,Adams81,NeumannI83} is shown
    for comparison using the Ewald convergence parameter ($\kappa$)
    and the real-space cutoff value ($r_c$). } 
\label{tab:A}
{%
\begin{tabular}{l|c|c|c}
       
Method & point charges & point 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}$ \\
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} $\\
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} $ \\
Taylor-shifted  (TSF) & 1 & 1 \\ 
Ewald-Kornfeld (EK) & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa
                      r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ & \\
\end{tabular}%
}
\end{table}
Note that for point charges, the GSF and TSF methods produce estimates
of the dielectric that need no correction, and the TSF method likewise
needs no correction for point dipoles. 

\section{Quadrupolar Fluids and the Quadrupolar Susceptibility}
\subsection{Response to External Perturbations}

A molecule with a permanent quadrupole, $\mathsf{q}$, will align in
the presence of an electric field gradient $\nabla\mathbf{E}$.  The
anisotropic polarization of the quadrupole is given
by,\cite{AduGyamfi78,AduGyamfi81}
\begin{equation}
\braket{\mathsf{q}} - \frac{\mathbf{I}}{3}
\mathrm{Tr}(\mathsf{q}) = \epsilon_o \alpha_q \nabla\mathbf{E},
\end{equation}
where $\alpha_q = q_o^2 / 15 \epsilon_o k_B T $ is a molecular quadrupole
polarizablity and $q_o$ is an effective quadrupole moment for the molecule,
\begin{equation}
 q_o^2 = 3 \mathsf{q}:\mathsf{q} - \mathrm{Tr}(\mathsf{q})^2.
\end{equation}

In the presence of an external field gradient, a system of quadrupolar
molecules also organizes with an anisotropic polarization,
\begin{equation}
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q}) =  \epsilon_o
\alpha_Q  \nabla\mathbf{E}
\end{equation}
where $\mathsf{Q}$ is the traced quadrupole density of the system and
$\alpha_Q$ is a macroscopic quadrupole polarizability which has
dimensions of $\mathrm{length}^{-2}$. Equivalently, the traceless form
may be used,
\begin{equation}
\mathsf{\Theta} = 3 \epsilon_o  \alpha_Q \nabla\mathbf{E},
\label{pertQuad}
\end{equation} 
where
$\mathsf{\Theta} = 3\mathsf{Q} - \mathbf{I} \mathrm{Tr}(\mathsf{Q})$
is the traceless tensor that also describes the system quadrupole
density.  It is this tensor that will be utilized to derive correction
factors below.

\subsection{Fluctuation Formula}
As in the dipolar case, we may define a system quadrupole moment,
$\mathsf{M}_Q = \sum_i \mathsf{q}_i$ and the traced quadrupolar
density, $\mathsf{Q} = \mathsf{M}_Q / V$.  A fluctuation formula can
be written for a system comprising quadrupolar
molecules,\cite{LoganI81,LoganII82,LoganIII82}
\begin{equation}
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q}) = \epsilon_o
\frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o
  V k_B T} \cdot \nabla\mathbf{E}.
\label{flucQuad}
\end{equation}
Some care is needed in the definitions of the averaged quantities.  These
refer to the effective quadrupole moment of the system, and they are
computed as follows,
\begin{align}
\braket{\mathsf{M}_Q^2} &= \braket{3 \mathsf{M}_Q:\mathsf{M}_Q -
  \mathrm{Tr}(\mathsf{M}_Q)^2 }\\
\braket{\mathsf{M}_Q}^2 &= 3 \braket{\mathsf{M}_Q}:\braket{\mathsf{M}_Q} -
\mathrm{Tr}(\braket{\mathsf{M}_Q})^2
\label{eq:flucQuad}
\end{align}
The macroscopic quadrupolarizability is given by,
\begin{equation}
  \alpha_Q = \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o V k_B T}
\label{propConstQuad}
\end{equation}
This relationship between a macroscopic property of a quadrupolar
fluid and microscopic fluctuations should be true even when the
applied field gradient $\nabla \mathbf{E} \rightarrow 0$.

\subsection{Correction Factors}
In this section we generalize the treatment of Neumann \textit{et al.}
for quadrupolar fluids. Interactions involving multiple quadrupoles
are rank 4 tensors, and we therefore describe quantities in this
section using Einstein notation.

In the presence of a uniform external field gradient,
$\partial_\alpha {E}^\circ_\beta $, the total field gradient at
$\mathbf{r}$ depends on the quadrupole polarization density at all
other points in the system,
\begin{equation}
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha
{E}^\circ_\beta(\mathbf{r}) + \frac{1}{4\pi \epsilon_o}\int
T_{\alpha\beta\gamma\delta}({\mathbf{r}-\mathbf{r}^\prime})
Q_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime 
\label{gradMaxwell}
\end{equation}
where and $T_{\alpha\beta\gamma\delta}$ is the quadrupole interaction
tensor connecting quadrupoles at $\mathbf{r}^\prime$ with the point of
interest ($\mathbf{r}$).

\begin{figure}
\includegraphics[width=3in]{QuadrupoleFigure}
\caption{With the real-space electrostatic methods, the molecular
  quadrupole tensor, $\mathbf{T}_{\alpha\beta\gamma\delta}(r)$,
  governing interactions between molecules is not the same for
  quadrupoles represented via sets of charges, point dipoles, or by
  single point quadrupoles.}
\label{fig:quadrupolarFluid}
\end{figure}

In simulations of quadrupolar fluids, the molecular quadrupoles may be
represented by closely-spaced point charges, by multiple point
dipoles, or by a single point quadrupole (see
Fig. \ref{fig:quadrupolarFluid}).  In the case where point charges are
interacting via an electrostatic kernel, $v(r)$, the effective
molecular quadrupole tensor can obtained from four successive
applications of the gradient operator to the electrostatic kernel,
\begin{eqnarray}
T_{\alpha\beta\gamma\delta}(\mathbf{r}) &=& \nabla_\alpha \nabla_\beta
                                   \nabla_\gamma \nabla_\delta v(r) \\
 &=& \left(\delta_{\alpha\beta}\delta_{\gamma\delta} +
     \delta_{\alpha\gamma}\delta_{\beta\delta}+
     \delta_{\alpha\delta}\delta_{\beta\gamma}\right)\left(-\frac{v^\prime(r)}{r^3}
     + \frac{v^{\prime \prime}(r)}{r^2}\right) \nonumber \\
 & &+ \left(\delta_{\alpha\beta} r_\gamma r_\delta + 5  \mathrm{~permutations}
     \right) \left(\frac{3v^\prime(r)}{r^5}-\frac{3v^{\prime \prime}(r)}{r^4} +
     \frac{v^{\prime \prime \prime}(r)}{r^3}\right) \nonumber \\
 & &+ r_\alpha r_\beta r_\gamma r_\delta
     \left(-\frac{15v^\prime(r)}{r^7}+\frac{15v^{\prime \prime}(r)}{r^6}-\frac{6v^{\prime
     \prime \prime}(r)}{r^5} + \frac{v^{\prime \prime \prime \prime}(r)}{r^4}\right),
\label{quadCharge}
\end{eqnarray}
where $v(r)$ can either be the electrostatic kernel ($1/r$) or one of
the modified (Wolf or DSF) kernels.  

Similarly, when representing quadrupolar molecules with multiple point
dipoles, the molecular quadrupole interaction tensor can be obtained
using two successive applications of the gradient operator to the
dipole interaction tensor,
\begin{eqnarray}
T_{\alpha\beta\gamma\delta}(\mathbf{r}) &=& \nabla_\alpha \nabla_\beta
                                            v_{\gamma\delta}(\mathbf{r}) \\ 
& = & \delta_{\alpha\beta}\delta_{\gamma\delta} \frac{v^\prime_{21}(r)}{r} +
      \left(\delta_{\alpha\gamma}\delta_{\beta\delta}+
      \delta_{\alpha\delta}\delta_{\beta\gamma}\right)\frac{v_{22}(r)}{r^2}
      \nonumber\\ 
 & &+ \delta_{\gamma\delta} r_\alpha r_\beta
     \left(\frac{v^{\prime \prime}_{21}(r)}{r^2}-\frac{v^\prime_{21}(r)}{r^3} \right)\nonumber \\
 & &+\left(\delta_{\alpha\beta} r_\gamma r_\delta +
     \delta_{\alpha\gamma} r_\beta r_\delta  +\delta_{\alpha\delta}
     r_\gamma r_\beta + \delta_{\beta\gamma} r_\alpha r_\delta
     +\delta_{\beta\delta} r_\alpha r_\gamma  \right)
     \left(\frac{v^\prime_{22}(r)}{r^3}-\frac{2v_{22}(r)}{r^4}\right)
     \nonumber \\
 & &+ r_\alpha r_\beta r_\gamma r_\delta
     \left(\frac{v^{\prime
     \prime}_{22}(r)}{r^4}-\frac{5v^\prime_{22}(r)}{r^5}+\frac{8v_{22}(r)}{r^6}\right),
\label{quadDip}
\end{eqnarray}
where $v_{\gamma\delta}(\mathbf{r})$ is a dipole-dipole interaction
tensor that depends on the level of the
approximation.\cite{PaperI,PaperII} Similarly $v_{21}(r)$ and
$v_{22}(r)$ are the radial function for different real space cutoff
methods defined in the first paper in this series.\cite{PaperI}

For quadrupolar liquids modeled using point quadrupoles, the
interaction tensor can be constructed as,
\begin{eqnarray}
T_{\alpha\beta\gamma\delta}(\mathbf{r}) & = &
                                              \left(\delta_{\alpha\beta}\delta_{\gamma\delta}
                                              +
                                              \delta_{\alpha\gamma}\delta_{\beta\delta}+
                                              \delta_{\alpha\delta}\delta_{\beta\gamma}\right)v_{41}(r)
                                              + \delta_{\gamma\delta}
                                              r_\alpha r_\beta
                                              \frac{v_{42}(r)}{r^2}
                                              \nonumber \\  
& & + r_\alpha r_\beta r_\gamma r_\delta  \left(\frac{v_{43}(r)}{r^4}\right), 
\label{quadRadial}
\end{eqnarray}
where again $v_{41}(r)$, $v_{42}(r)$, and $v_{43}(r)$ are radial
functions defined in Paper I of the series. \cite{PaperI} Note that
these radial functions have different functional forms depending on
the level of approximation being employed.

The integral in equation (\ref{gradMaxwell}) can be divided into two
parts, $|\mathbf{r}-\mathbf{r}^\prime|\rightarrow 0 $ and
$|\mathbf{r}-\mathbf{r}^\prime|> 0$. Since the self-contribution to
the field gradient vanishes at the singularity (see Appendix
\ref{singularQuad}), equation (\ref{gradMaxwell}) can be written as,
\begin{equation}
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha {E}^\circ_\beta(\mathbf{r}) +
  \frac{1}{4\pi \epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 }
  T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime)
  {Q}_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime.
\end{equation}
If $\mathbf{r} = \mathbf{r}^\prime$ is excluded from the integration,
the total gradient can be most easily expressed in terms of
traceless quadrupole density as below,\cite{LoganI81}
\begin{equation}
\partial_\alpha E_\beta(\mathbf{r}) = \partial_\alpha
{E}^\circ_\beta(\mathbf{r}) + \frac{1}{12\pi
  \epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 }
T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime) \Theta_{\gamma\delta}(\mathbf{r}') d\mathbf{r}^\prime,
\end{equation}
where
$\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$
is the traceless quadrupole density. In analogy to
Eq. (\ref{pertQuad}) above, the quadrupolarization density may now be
related to the actual quadrupolar susceptibility, $\chi_Q$,
\begin{equation}
\frac{1}{3}{\Theta}_{\alpha\beta}(\mathbf{r}) = \epsilon_o {\chi}_Q
\left(\partial_\alpha {E}^\circ_\beta(\mathbf{r}) + \frac{1}{12\pi
    \epsilon_o}\int\limits_{|\mathbf{r}-\mathbf{r}^\prime|> 0 }
  T_{\alpha\beta\gamma\delta}(\mathbf{r}-\mathbf{r}^\prime)
  \Theta_{\gamma\delta}(\mathbf{r}^\prime) d\mathbf{r}^\prime \right).
\end{equation}
For toroidal boundaries and with a uniform imposed
$\partial_\alpha E^\circ_\beta$, the quadrupole density
${\Theta}_{\alpha\beta}$ will be uniform over the entire space. After
performing a Fourier transform (see the Appendix in
ref. \onlinecite{NeumannI83}) we obtain,
\begin{equation}
\frac{1}{3}\tilde{\Theta}_{\alpha\beta}(\mathbf{k})=
\epsilon_o {\chi}_Q \left[{\partial_\alpha
    \tilde{E}^\circ_\beta}(\mathbf{k})+ \frac{1}{12\pi
    \epsilon_o} \tilde{T}_{\alpha\beta\gamma\delta}(\mathbf{k})
 \tilde{\Theta}_{\gamma\delta}(\mathbf{k})\right].
\end{equation} 

The only surviving components of the quadrupole polarization will be
in the direction of the applied field gradient, so we may simplify
this expression considerably.  That is, the only relevant polarization
term is $\tilde{\Theta}_{\alpha\beta}(\mathbf{k})$ and the only relevant
tensor contributions will have doubly-repeated indices, e.g.
$\tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k})$. Therefore we may
consider each component of the interaction tensor separately in a
scalar expression,
\begin{eqnarray}
\frac{1}{3}\tilde{\Theta}_{\alpha\beta}(\mathbf{k}) &=& \epsilon_o
                                                      {\chi}_Q
                                                      \left[{\partial_\alpha
                                                      \tilde{E}^\circ_\beta}(\mathbf{k})+
                                                      \frac{1}{12\pi
                                                      \epsilon_o}\tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k})
                                                      {\tilde{\Theta}}_{\alpha\beta}(\mathbf{k})\right]
                                                      \nonumber \\
&=& \epsilon_o {\chi}_Q \left(1-\frac{1}{4\pi} {\chi}_Q \tilde{T}_{\alpha\beta\alpha\beta}(\mathbf{k})\right)^{-1}
{\partial_\alpha \tilde{E}^\circ_\beta}(\mathbf{k}) 
\label{fourierQuad}
\end{eqnarray} 
If the applied field gradient is homogeneous over the
entire volume, ${\partial_ \alpha \tilde{E}^\circ_\beta}(\mathbf{k}) = 0 $ except at
$ \mathbf{k} = 0$.  As in the dipolar case, the only relevant
contribution of the interaction tensor will also be when $\mathbf{k} = 0$.
Therefore equation (\ref{fourierQuad}) can be written as,
\begin{equation}
\frac{1}{3}{\tilde{\Theta}}_{\alpha\beta}({0}) = \epsilon_o {\chi}_Q
\left(1-\frac{1}{4\pi} {\chi}_Q \tilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1} \partial_\alpha \tilde{E}^\circ_\beta({0})
\label{fourierQuad2}
\end{equation}
where $\tilde{T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as,
\begin{equation}
\tilde{T}_{\alpha\beta\alpha\beta}({0}) = \int_V {T}_{\alpha\beta\alpha\beta}(\mathbf{r})d\mathbf{r}
\label{realTensorQaud}
\end{equation}
We may now define a quantity to represent the contribution from the
$\mathbf{k} = 0$ portion of the interaction tensor,
\begin{equation}
B = \frac{1}{4 \pi} \int_V T_{\alpha \beta \alpha \beta}(\mathbf{r})
d\mathbf{r}
\end{equation}
which has been integrated over the interaction volume $V$ and has
units of $\mathrm{length}^{-2}$.

In terms of traced quadrupole moment, equation (\ref{fourierQuad2})
can be written as,
\begin{equation}
\mathsf{Q} - \frac{\mathbf{I}}{3} \mathrm{Tr}(\mathsf{Q})
= \epsilon_o {\chi}_Q\left(1-  {\chi}_Q B \right)^{-1} \nabla \mathbf{E}^\circ
\label{tracedConstQuad}
\end{equation}
Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we obtain,
\begin{equation}
\frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o
  V k_B T} = \chi_Q \left(1 - \chi_Q B \right)^{-1},
\end{equation}
or equivalently,
\begin{equation}
\chi_Q = \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o
  V k_B T} \left(1 + B \frac{\braket{\mathsf{M}_Q^2}-{\braket{\mathsf{M}_Q}}^2}{15 \epsilon_o
  V k_B T} \right)^{-1}
\label{eq:finalForm}
\end{equation}
Eq. (\ref{eq:finalForm}) now expresses a bulk property (the
quadrupolar susceptibility, $\chi_Q$) in terms of a straightforward
fluctuation in the system quadrupole moment and a quadrupolar
correction factor ($B$).  The correction factors depend on the cutoff
method being employed in the simulation, and these are listed in Table
\ref{tab:B}.  

In terms of the macroscopic quadrupole polarizability, $\alpha_Q$,
which may be thought of as the ``conducting boundary'' version of the
susceptibility,
\begin{equation}
\chi_Q = \frac{\alpha_Q}{1 + B \alpha_Q}.
\end{equation}
If an electrostatic method produces $B \rightarrow 0$, the computed
quadrupole polarizability and quadrupole susceptibility converge to
the same value.

\begin{sidewaystable}
  \caption{Expressions for the quadrupolar correction factor ($B$) for
    the real-space electrostatic methods in terms of the damping
    parameter ($\alpha$) and the cutoff radius ($r_c$). The units
    of the correction factor are $ \mathrm{length}^{-2}$ for quadrupolar
    fluids.} 
\label{tab:B}
\begin{tabular}{l|c|c|c|c}
       
Method & charges & dipoles & quadrupoles \\\hline
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}$ \\
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}$  \\
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)$\\
Taylor-shifted  (TSF) &  $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}}
                        e^{-\alpha^2 r_c^2}$ &
                                               $\frac{4\;\mathrm{erfc(\alpha r_c)}}{{r_c}^2} + \frac{8 \alpha}{3\sqrt{\pi}r_c}e^{-\alpha^2 {r_c}^2}\left(3+ 2 \alpha^2 {r_c}^2 + \alpha^4 {r_c}^4\right)  $ & $\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)$ \\
Ewald-Kornfeld (EK) & & &
\end{tabular}
\end{sidewaystable}

\section{Screening of Charges by Multipolar Fluids}
\label{sec:PMF}
In a dipolar fluid, the static dielectric constant is also a measure
of the ability of the fluid to screen charges from one another.  A set
of point charges creates an inhomogeneous field in the fluid, and the
fluid responds to this field as if it was created externally or via
local polarization fluctuations. For this reason, the dielectric
constant can be used to estimate an effective potential between two
point charges ($C_i$ and $C_j$) embedded in the fluid,
\begin{equation}
U_\mathrm{effective} = \frac{C_i C_j}{4 \pi \epsilon_0 \epsilon
  r_{ij}}.
\label{eq:effectivePot}
\end{equation}

The same set of point charges can also create an inhomogeneous field
\textit{gradient}, and this will cause a response in a quadrupolar
fluid that will also cause an effective screening. As discussed above,
however, the relevant phyiscal property in quadrupolar fluids is the
susceptibility, $\chi_Q$.  The screening dielectric associated with
the quadrupolar susceptibility is defined as,\cite{Ernst92}
\begin{equation}
\epsilon = 1 + \chi_Q G = 1 + G \frac{\alpha_Q}{1 + \alpha_Q  B}
\label{eq:dielectricFromQuadrupoles}
\end{equation}
where $G$ is a geometrical factor that depends on the geometry of the
field perturbation,
\begin{equation}
G = \frac{\int_V \left| \nabla \mathbf{E}^\circ \right|^2 d\mathbf{r}}
{\int_V \left|\mathbf{E}^\circ\right|^2 d\mathbf{r}}
\end{equation}
integrated over the interaction volume. Note that this geometrical
factor is also required to compute effective dielectric constants even
when the field gradient is homogeneous over the entire sample.

To measure effective screening in a multipolar fluid, we compute an
effective interaction potential, the potential of mean force (PMF),
between positively and negatively charged ions when they screened by
the intervening fluid.  The PMF is obtained from a sequence of
simulations in which two ions are constrained to a fixed distance, and
the average constraint force to hold them at a fixed distance $r$ is
collected during a long simulation,\cite{Wilfred07}
\begin{equation}
w(r) = \int_{r_o}^{r}\left\langle \frac{\partial f}{\partial r^\prime}
\right\rangle_{r^\prime} dr^\prime + 2k_BT \ln\left(\frac{r}{r_o}\right) + w(r_o),
\label{eq:pmf}
\end{equation}
where $\braket{\partial f/\partial r^\prime}_{r^\prime}$ is the mean
constraint force required to hold the ions at distance $r^\prime$,
$2k_BT \log(r/r_o)$ is the Fixman factor,\cite{Fixman:1974fk} and
$r_o$ is a reference position (usually taken as a large separation
betwen the ions). If the dielectric constant is a good measure of the
screening at all inter-ion separations, we would expect $w(r)$ to have
the form in Eq. (\ref{eq:effectivePot}).  Because real fluids are not
continuum dielectrics, the effective dielectric constant is a function
of the interionic separation,
\begin{equation}
\epsilon(r) = \frac{u_\mathrm{raw}(r) }{w(r)} 
\end{equation}
where $u_\mathrm{raw}(r)$ is the direct charge-charge interaction
potential that is in use during the simulation.  $\epsilon(r)$ may
vary considerably from the bulk estimates at short distances, although
it should converge to the bulk value as the separation between the
ions increases.

\section{Simulation Methodology}
To test the formalism developed in the preceding sections, we have
carried out computer simulations using three different techniques: i)
simulations in the presence of external fields, ii) equilibrium
calculations of box moment fluctuations, and iii) potentials of mean
force (PMF) between embedded ions. In all cases, the fluids were
composed of point multipoles protected by a Lennard-Jones potential.
The parameters used in the test systems are given in table
\ref{Tab:C}.

\begin{table}
  \caption{\label{Tab:C}The parameters used in simulations to evaluate
    the dielectric response of the new real-space methods.} 
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline
             & \multicolumn{2}{c|}{LJ parameters} &
             \multicolumn{5}{c|}{Electrostatic moments} & & & & \\
 Test system & $\sigma$& $\epsilon$ & $C$ & $D$  &
 $Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass  & $I_{xx}$ & $I_{yy}$ &
 $I_{zz}$ \\ \cline{6-8}\cline{10-12}
 & (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu
 \AA\textsuperscript{2})} \\ \hline
    Dipolar fluid & 3.41 & 0.2381 & - & 1.4026 &-&-&-& 39.948 & 11.613 & 11.613 & 0.0 \\
Quadrupolar fluid & 2.985 & 0.265 & - & - & 0.0 & 0.0 &-2.139 & 18.0153 & 43.0565 & 43.0565 & 0.0  \\
              \ce{q+} & 1.0 & 0.1 & +1 & - & - & - & - & 22.98 & - & - & - \\
              \ce{q-} & 1.0 & 0.1 & -1 & - & - & - & - & 22.98 & - & - & - \\ \hline
\end{tabularx}
\end{table}

The first of the test systems consists entirely of fluids of point
dipolar or quadrupolar molecules in the presence of constant field or
field gradients.  Since there are no isolated charges within the
system, the divergence of the field will be zero, \textit{i.e.}
$\vec{\nabla} \cdot \mathbf{E} = 0$. This condition can be satisfied
by using the relatively simple applied potential as described in the
supporting information.

When a constant electric field or field gradient is applied to the
system, the molecules align along the direction of the applied field,
and polarize to a degree determined both by the strength of the field
and the fluid's polarizibility.  We have calculated ensemble averages
of the box dipole and quadrupole moments as a function of the strength
of the applied fields.  If the fields are sufficiently weak, the
response is linear in the field strength, and one can easily compute
the polarizability directly from the simulations.

The second set of test systems consists of equilibrium simulations of
fluids of point dipolar or quadrupolar molecules simulated in the
absence of any external perturbation. The fluctuation of the ensemble
averages of the box multipolar moment was calculated for each of the
multipolar fluids. The box multipolar moments were computed as simple
sums over the instantaneous molecular moments, and fluctuations in
these quantities were obtained from Eqs. (\ref{eq:flucDip}) and
(\ref{eq:flucQuad}). The macroscopic polarizabilities of the system at
a were derived using Eqs.(\ref{flucDipole}) and (\ref{flucQuad}).

The final system consists of dipolar or quadrupolar fluids with two
oppositely charged ions embedded within the fluid. These ions are
constrained to be at fixed distance throughout a simulation, although
they are allowed to move freely throughout the fluid while satisfying
that constraint. Separate simulations were run at a range of
constraint distances. A dielectric screening factor was computed using
the ratio between the potential between the two ions in the absence of
the fluid medium and the PMF obtained from the simulations.

We carried out these simulations for all three of the new real-space
electrostatic methods (SP, GSF, and TSF) that were developed in the
first paper (Ref. \onlinecite{PaperI}) in the series. The radius of
the cutoff sphere was taken to be 12~\AA. Each of the real space
methods also depends on an adjustable damping parameter $\alpha$ (in
units of $\mathrm{length}^{-1}$).  We have selected ten different
values of damping parameter: 0.0, 0.05, 0.1, 0.15, 0.175, 0.2, 0.225,
0.25, 0.3, and 0.35~\AA$^{-1}$ in our simulations of the dipolar
liquids, while four values were chosen for the quadrupolar fluids:
0.0, 0.1, 0.2, and 0.3~\AA$^{-1}$.

For each of the methods and systems listed above, a reference
simulation was carried out using a multipolar implementation of the
Ewald sum.\cite{Smith82,Smith98} A default tolerance of
$1 \times 10^{-8}$~kcal/mol was used in all Ewald calculations,
resulting in Ewald coefficient 0.3119~\AA$^{-1}$ for a cutoff radius
of 12~\AA.  All of the electrostatics and constraint methods were
implemented in our group's open source molecular simulation program,
OpenMD,\cite{Meineke05,openmd} which was used for all calculations in
this work.

Dipolar systems contained 2048 Lennard-Jones-protected point dipolar
(Stockmayer) molecules with reduced density $\rho^* = 0.822$,
temperature $T^* = 1.15$, moment of inertia $I^* = 0.025$, and dipole
moment $\mu^* = \sqrt{3.0}$.  These systems were equilibrated for
0.5~ns in the canonical (NVT) ensemble.  Data collection was carried
out over a 1~ns simulation in the microcanonical (NVE) ensemble.  Box
dipole moments were sampled every fs.  For simulations with external
perturbations, field strengths ranging from
$0 - 10 \times 10^{-4} \mathrm{~V/\r{A}}$ with increments of
$ 10^{-4} \mathrm{~V/\r{A}}$ were carried out for each system.

Quadrupolar systems contained 4000 linear point quadrupoles with a
density $2.338 \mathrm{~g/cm}^3$ at a temperature of 500~K.  These
systems were equilibrated for 200~ps in a canonical (NVT) ensemble.
Data collection was carried out over a 500~ps simulation in the
microcanonical (NVE) ensemble.  Components of box quadrupole moments
were sampled every 100 fs.  For quadrupolar simulations with external
field gradients, field strengths ranging from
$0 - 9 \times 10^{-2} \mathrm{~V/\r{A}}^2$ with increments of
$10^{-2} \mathrm{~V/\r{A}^2}$ were carried out for each system.

To carry out the PMF simulations, two of the multipolar molecules in
the test system were converted into \ce{q+} and \ce{q-} ions and
constrained to remain at a fixed distance for the duration of the
simulation. The constrained distance was then varied from
$5 - 12 \mathrm{~\r{A}}$. In the PMF calculations, all simulations were
equilibrated for 500 ps in the NVT ensemble and run for 5 ns in the
microcanonical (NVE) ensemble.  Constraint forces were sampled every
20~fs.
      
\section{Results}
\subsection{Dipolar fluid}
\begin{figure}
\includegraphics[width=5in]{dielectricFinal_Dipole.pdf}
\caption{The polarizability ($\alpha_D$), correction factor ($A$), and
  dielectric constant ($\epsilon$) for the test dipolar fluid. The
  left panels were computed using external fields, and those on the
  right are the result of equilibrium fluctuations.  In the GSF and SP
  methods, the corrections are large in with small values of $\alpha$,
  and a optimal damping coefficient is evident around 0.25 \AA$^{-1}$.
  The dashed lines in the upper panel indicate back-calculation of the
  polarizability using the Ewald estimate
  (Ref. \onlinecite{NeumannI83}) for the dielectric constant. NOTE:
  Change $A_{dipole}$ to $A$.}
\label{fig:dielectricDipole}
\end{figure}
The macroscopic polarizability ($\alpha_D$) for the dipolar fluid is
shown in the upper panels in Fig.  \ref{fig:dielectricDipole}.  The
polarizability obtained from the both perturbation and fluctuation
approaches are in excellent agreement with each other.  The data also
show a stong dependence on the damping parameter for both the Shifted
Potential (SP) and Gradient Shifted force (GSF) methods, while Taylor
shifted force (TSF) is largely independent of the damping
parameter.

The calculated correction factors discussed in section
\ref{sec:corrFactor} are shown in the middle panels.  Because the TSF
method has $A = 1$ for all values of the damping parameter, the
computed polarizabilities need no correction for the dielectric
calculation. The value of $A$ varies with the damping parameter in
both the SP and GSF methods, and inclusion of the correction yields
dielectric estimates (shown in the lower panel) that are generally too
large until the damping reaches $\sim 0.25 \mathrm{~\r{A}}^{-1}$.
Above this value, the dielectric constants are generally in good
agreement with previous simulation results.\cite{NeumannI83}

Figure \ref{fig:dielectricDipole} also contains back-calculations of
the polarizability using the reference (Ewald) simulation
results.\cite{NeumannI83} These are indicated with dashed lines in the
upper panels.  It is clear that the expected polarizability for the SP
and GSF methods are quite close to results obtained from the
simulations.  This indicates that the correction formula for the
dipolar fluid (Eq. \ref{correctionFormula}) is quite sensitive when
the value of $\mathrm{A}$ deviates significantly from unity.

These results also suggest an optimal value for the damping parameter
of ($\alpha \sim 0.2-0.3 \mathrm{\r{A}}^{-1}$ when evaluating
dielectric constants of point dipolar fluids using either the
perturbation and fluctuation approaches for the new real-space
methods.

\begin{figure}
\includegraphics[width=4in]{ScreeningFactor_Dipole.pdf}
\caption{The distance-dependent screening factor, $\epsilon(r)$,
  between two ions immersed in the dipolar fluid. The new methods are
  shown in separate panels, and different values of the damping
  parameter ($\alpha$) are indicated with different symbols. All of
  the methods appear to be converging to the bulk dielectric constant
  ($\sim 65$) at large ion separations. CHECK PLOT}
\label{fig:ScreeningFactor_Dipole}
\end{figure}
We have also evaluated the distance-dependent screening factor,
$\epsilon(r)$, between two oppositely charged ions when they are
placed in the dipolar fluid.  These results were computed using
Eq. \ref{eq:pmf} and are shown in Fig.
\ref{fig:ScreeningFactor_Dipole}.

The screening factor is similar to the dielectric constant, but
measures a local property of the ions in the fluid and depends on both
ion-dipole and dipole-dipole interactions. These interactions depend
on the distance between ions as well as the electrostatic interaction
methods utilized in the simulations. The screening should converge to
the dielectric constant when the field due to ions is small. This
occurs when the ions are separated (or when the damping parameter is
large). In Fig. \ref{fig:ScreeningFactor_Dipole} we observe that for
the higher value of damping alpha \textit{i.e.}
$\alpha = 0.2 \AA^{-1}$ and $0.3 \AA^{-1}$ and large separation
between ions, the screening factor does indeed approach the correct
dielectric constant.

REVISIT THESE THREE PARAGRAPHS:

It is also notable that the TSF method again displays smaller
perturbations away from the correct dielectric screening behavior.  We
also observe that for TSF method yields high dielectric screening even
for lower values of $\alpha$. 

At short distances, the presence of the ions creates a strong local
field that acts to align nearby dipoles nearly perfectly in opposition
to the field from the ions.  This has the effect of increasing the
effective screening when the ions are brought close to one another.

This is due to the absence of corrections for dipole-dipole
interactions as discussed above.  However, the TSF remains quite
perturbative for ion-dipole interactions, an effect that is most
apparent at small ion separations.

\subsection{Quadrupolar fluid}
\begin{figure}
\includegraphics[width=4in]{polarizabilityFinal_Quad.pdf}
\caption{The quadrupole polarizability ($\alpha_Q$), correction factor
  ($B$), and susceptibility ($\chi_Q$) for the test quadrupolar
  fluid. The left panels were computed using external field gradients,
  and those on the right are the result of equilibrium fluctuations.
  The GSF and SP methods allow nearly unmodified use of the
  ``conducting boundary'' or polarizability results in place of the
  bulk susceptibility. NOTE: Change $A_{quad}$ to $B$, change $\chi$
  to $\chi_Q$.}
\label{fig:dielectricQuad}
\end{figure}
The polarizability ($\alpha_Q$), correction factor ($B$), and
susceptibility ($\chi_Q$) for the quadrupolar fluid is plotted against
damping parameter Fig.  \ref{fig:dielectricQuad}.  In quadrupolar
fluids, both the polarizability and susceptibility have units of
$\mathrm{length}^2$ Although the susceptibility has dimensionality, it
is the relevant measure of macroscopic quadrupolar
properties.\cite{JeonI03, JeonII03}. The left panel in
Fig. \ref{fig:dielectricQuad} shows results obtained from the applied
field gradient simulations whereas the results from the equilibrium
fluctuation formula are plotted in the right panels.

The susceptibility for the quadrupolar fluid is obtained from
quadrupolarizability and a correction factor using
Eq. (\ref{eq:quadrupolarSusceptiblity}).  The susceptibilities are
shown in the bottom panels of Fig. \ref{fig:dielectricQuad}. All three
methods: (SP, GSF, and TSF) produce similar susceptibilities over the
range of damping parameters.  This shows that susceptibility derived
using the quadrupolarizability and the correction factors are
essentially independent of the electrostatic method utilized in the
simulation.

\begin{figure}
\includegraphics[width=5in]{ScreeningFactor_Quad.pdf}
\caption{\label{fig:screenQuad} The distance-dependent screening
  factor, $\epsilon(r)$, between two ions immersed in the quadrupolar
  fluid. Results from  the perturbation and fluctuation methods are
  shown in right and central panels. Here the susceptibility is
  calculated from the  simulation and the geometrical factor is
  evaluated analytically, using the field and field-gradient produced
  by ions. The right hand panel shows the screening factor obtained
  from the PMF calculations.}
\end{figure}
A more difficult test of the quadrupolar susceptibility is made by
comparing with direct calculation of the electrostatic screening using
the potential of mean force (PMF). Since the effective dielectric
constant for a quadrupolar fluid depends on the geometry of the field
and field gradient, this is not a physical property of the quadrupolar
fluid. 

The geometrical factor for embedded ions changes with the ion
separation distance. It is therefore reasonable to treat the
dielectric constant as a distance-dependent screening factor.  Since
the quadrupolar molecules couple with the gradient of the field, the
distribution of the quadrupoles will be inhomogeneously distributed
around the point charges.  Hence the distribution of quadrupolar
molecules should be taken into account when computing the geometrical
factors in the presence of this perturbation,
\begin{eqnarray}
G &=& \frac{\int_V g(\mathbf{r})  \left|\nabla \mathbf{E}^\circ
  \right|^2 d\mathbf{r}}{\int_V\left|\mathbf{E}^\circ\right|^2
  d\mathbf{r}}  \nonumber \\ \nonumber \\
 &=& \frac{ 2 \pi \int_{-1}^{1}  \int_{0}^{R} r^2  g(r,
     \cos\theta)  \left|\nabla \mathbf{E}^\circ  \right|^2 dr d\cos\theta }{\int_V\left|\mathbf{E}^\circ\right|^2 d\mathbf{r}}
\label{eq:geometricalFactor}
\end{eqnarray}
where $g(r,\cos\theta)$ is a distribution function for the quadrupoles
with respect to an origin at midpoint of a line joining the two probe
charges.

The effective screening factor is plotted against ion separation
distance in Fig. \ref{fig:screenQuad}.  The screening evaluated from
the perturbation and fluctuation methods are shown in right and
central panels. Here the susceptibility is calculated from the
simulation and the geometrical factor is evaluated analytically, using
the field and field-gradient produced by ions. The right hand panel
shows the screening factor obtained from the PMF calculations.

We note that the screening factor obtained from both the perturbation
and fluctuation formula show good agreement with the screening factor
calculated using PMF method. As there are no large differences in
quadrupole-quadruople interactions for various real-space
methods,\cite{PaperI,PaperII} we generally good agreement for the
screening factors using any of the real space methods.


\section{Discussion}
We have used the perturbation and fluctuation formula to evaluate
dielectric properties for dipolar and quadrupolar fluids implementing
SP, GSF, and TSF methods in the simulation. Similarly the correction
factors for the both dipolar and quadrupolar systems have also been
calculated for all SP, GSF, and TSF methods.  We have also derived
screening factor between two ions immersed in the dipolar and
quadrupolar fluids using PMF method. The result from the perturbation
and fluctuation formula are compared with PMF method to test the
accuracy of the simulation result.

For the dipolar fluid, the polarizability evaluated using the perturbation and fluctuation methods show excellent agreement with each other. The dielectric constant evaluated using polarizability and correction factor agree with the previous simulation results \cite{NeumannI83} for the damping parameter 0.25 - 0.3$\AA^{-1}$ and the electrostatic interaction methods, SP and GSF implemented in the simulation. Since the correction factor for TSF, $\mathrm{A_{dipole}} = 1$, it produces dielectric constant more closer to previous simulation results \cite{NeumannI83} for all values damping parameter. This method also produces best dielectric constant for damping parameter 0.15 - 0.25$\AA^{-1}$. We have also found that correction formula (equation ~\ref{correctionFormula}) is very sensitive, when the correction factor ($\mathrm{A_{dipole}}$) much away from the 1 for the SP and GSF methods. Therefore it is better to choose damping parameter from 0.25 - 0.3$\AA^{-1}$ to evaluate dielectric constant using SP and GSF methods. The screening factor between ions is found to be closer to dielectric constant when the distance between the ions are far apart ($\sim 9 \AA$) and damping parameter $0.2$ to $0.3 \AA^{-1}$.

The quadpolarizability evaluated from the both perturbation and fluctuation methods shows excellent agreement with each other for the case of quadrupolar fluid. The susceptibility is calculated from the quadpolarizability and correction factor tends to produce same result for all; SP, GSF, and TSF, methods. Similarly the screening factor calculated using susceptibility and geometrical factor shows excellent agreement with the result obtained from the PMF method.
   

\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}
\section{Gradient of the field due to quadrupolar polarization}
\label{singularQuad}
In this section, we will discuss the gradient of the field produced by
quadrupolar polarization. For this purpose, we consider a distribution
of charge ${\rho}(r)$ which gives rise to an electric field
$\vec{E}(r)$ and gradient of the field $\vec{\nabla} \vec{E}(r)$
throughout space. The total gradient of the electric field over volume
due to the all charges within the sphere of radius $R$ is given by
(cf. Jackson equation 4.14):
\begin{equation}
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r = -\int_{r=R} R^2 \vec{E}\;\hat{n}\; d\Omega
\label{eq:8}
\end{equation}
where $d\Omega$ is the solid angle and $\hat{n}$ is the normal vector
of the surface of the sphere which is equal to
$sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} +
cos[\theta]\hat{z}$
in spherical coordinates.  For the charge density ${\rho}(r')$, the
total gradient of the electric field can be written as (cf. Jackson
equation 4.16),
\begin{equation}
\int_{r<R} \vec{\nabla}\vec{E}\; d^3r=-\int_{r=R} R^2\; \vec{\nabla}\Phi\; \hat{n}\; d\Omega  =-\frac{1}{4\pi\;\epsilon_o}\int_{r=R} R^2\; \vec{\nabla}\;\left(\int \frac{\rho(r')}{|\vec{r}-\vec{r'}|}\;d^3r'\right) \hat{n}\; d\Omega
\label{eq:9}
\end{equation}
The radial function in the equation (\ref{eq:9}) can be expressed in
terms of spherical harmonics as (cf. Jackson equation 3.70),
\begin{equation}
\frac{1}{|\vec{r} - \vec{r'}|} = 4\pi \sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\frac{{r^l_<}}{{r^{l+1}_>}}\;{Y^*}_{lm}(\theta', \phi')\;Y_{lm}(\theta, \phi)
\label{eq:10}
\end{equation}
If the sphere completely encloses the charge density then $ r_< = r'$ and $r_> = R$. Substituting equation (\ref{eq:10}) into (\ref{eq:9}) we get,
\begin{equation}
\begin{split}
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r &=-\frac{R^2}{\epsilon_o}\int_{r=R} \; \vec{\nabla}\;\left(\int \rho(r')\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\frac{{r'^l}}{{R^{l+1}}}\;{Y^*}_{lm}(\theta', \phi')\;Y_{lm}(\theta, \phi)\;d^3r'\right) \hat{n}\; d\Omega \\
 &= -\frac{R^2}{\epsilon_o}\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\int \rho(r')\;{r'^l}\;{Y^*}_{lm}(\theta', \phi')\left(\int_{r=R}\vec{\nabla}\left({R^{-(l+1)}}\;Y_{lm}(\theta, \phi)\right)\hat{n}\; d\Omega \right)d^3r
'
\end{split}
\label{eq:11}
\end{equation} 
The gradient of the product of radial function and spherical harmonics
is given by (cf. Arfken, p.811 eq. 16.94):
\begin{equation}
\begin{split}
\vec{\nabla}\left[ f(r)\;Y_{lm}(\theta, \phi)\right] = &-\left(\frac{l+1}{2l+1}\right)^{1/2}\; \left[\frac{\partial}{\partial r}-\frac{l}{r} \right]f(r)\; Y_{l, l+1, m}(\theta, \phi)\\ &+ \left(\frac{l}{2l+1}\right)^{1/2}\left[\frac
{\partial}{\partial r}+\frac{l}{r} \right]f(r)\; Y_{l, l-1, m}(\theta, \phi).
\end{split}
\label{eq:12}
\end{equation}
Using equation (\ref{eq:12}) we get,
\begin{equation}
\vec{\nabla}\left({R^{-(l+1)}}\;Y_{lm}(\theta, \phi)\right) = [(l+1)(2l+1)]^{1/2}\; Y_{l,l+1,m}(\theta, \phi) \; \frac{1}{R^{l+2}},
\label{eq:13}
\end{equation}
where $ Y_{l,l+1,m}(\theta, \phi)$ is the vector spherical harmonics
which can be expressed in terms of spherical harmonics as shown in
below (cf. Arfkan p.811),
\begin{equation}
Y_{l,l+1,m}(\theta, \phi) = \sum_{m_1, m_2} C(l+1,1,l|m_1,m_2,m)\; {Y_{l+1}}^{m_1}(\theta,\phi)\; \hat{e}_{m_2},
\label{eq:14}
\end{equation}
where $C(l+1,1,l|m_1,m_2,m)$ is a Clebsch-Gordan coefficient and
$\hat{e}_{m_2}$ is a spherical tensor of rank 1 which can be expressed
in terms of Cartesian coordinates,
\begin{equation}
{\hat{e}}_{+1} = - \frac{\hat{x}+i\hat{y}}{\sqrt{2}},\quad {\hat{e}}_{0} = \hat{z},\quad and \quad {\hat{e}}_{-1} = \frac{\hat{x}-i\hat{y}}{\sqrt{2}}
\label{eq:15}
\end{equation} 
The normal vector $\hat{n} $ can be expressed in terms of spherical tensor of rank 1 as shown in below,
\begin{equation}
\hat{n} = \sqrt{\frac{4\pi}{3}}\left(-{Y_1}^{-1}{\hat{e}}_1 -{Y_1}^{1}{\hat{e}}_{-1} + {Y_1}^{0}{\hat{e}}_0 \right)
\label{eq:16}
\end{equation}
The surface integral of the product of $\hat{n}$ and
${Y_{l+1}}^{m_1}(\theta, \phi)$ gives,
\begin{equation}
\begin{split}
\int \hat{n}\;{Y_{l+1}}^{m_1}\;d\Omega &= \int \sqrt{\frac{4\pi}{3}}\left(-{Y_1}^{-1}{\hat{e}}_1 -{Y_1}^{1}{\hat{e}}_{-1} + {Y_1}^{0}{\hat{e}}_0 \right)\;{Y_{l+1}}^{m_1}\; d\Omega \\
&=  \int \sqrt{\frac{4\pi}{3}}\left({{Y_1}^{1}}^* {\hat{e}}_1 +{{Y_1}^{-1}}^* {\hat{e}}_{-1} + {{Y_1}^{0}}^* {\hat{e}}_0 \right)\;{Y_{l+1}}^{m_1}\; d\Omega \\
&=   \sqrt{\frac{4\pi}{3}}\left({\delta}_{l+1, 1}\;{\delta}_{1, m_1}\;{\hat{e}}_1 + {\delta}_{l+1, 1}\;{\delta}_{-1, m_1}\;{\hat{e}}_{-1}+ {\delta}_{l+1, 1}\;{\delta}_{0, m_1} \;{\hat{e}}_0\right),
\end{split}
\label{eq:17}
\end{equation}
where ${Y_{l}}^{-m} = (-1)^m\;{{Y_{l}}^{m}}^* $ and
$ \int {{Y_{l}}^{m}}^*\;{Y_{l'}}^{m'}\;d\Omega =
\delta_{ll'}\delta_{mm'} $.
Non-vanishing values of equation \ref{eq:17} require $l = 0$,
therefore the value of $ m = 0 $. Since the values of $ m_1$ are -1,
1, and 0 then $m_2$ takes the values 1, -1, and 0, respectively
provided that $m = m_1 + m_2$.  Equation \ref{eq:11} can therefore be
modified,
\begin{equation}
\begin{split}
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r = &- \sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\int \rho(r')\;{Y^*}_{00}(\theta', \phi')[ C(1, 1, 0|-1,1,0)\;{\hat{e}_{-1}}{\hat{e}_{1}}\\  &+ C(1, 1, 0|-1,1,0)\;{\hat{e}_{1}}{\hat{e}_{-1}}+C(
1, 1, 0|0,0,0)\;{\hat{e}_{0}}{\hat{e}_{0}} ]\; d^3r'.
\end{split}
\label{eq:18} 
\end{equation}
After substituting ${Y^*}_{00} = \frac{1}{\sqrt{4\pi}} $ and using the
values of the Clebsch-Gorden coefficients:  $ C(1, 1, 0|-1,1,0) =
\frac{1}{\sqrt{3}}, \;   C(1, 1, 0|-1,1,0)= \frac{1}{\sqrt{3}}$ and $
C(1, 1, 0|0,0,0) = -\frac{1}{\sqrt{3}}$ in equation \ref{eq:18} we
obtain,
\begin{equation}
\begin{split}
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r &= -\sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\int \rho(r')\;d^3r'\left({\hat{e}_{-1}}{\hat{e}_{1}}+{\hat{e}_{1}}{\hat{e}_{-1}}-{\hat{e}_{0}}{\hat{e}_{0}}\right)\\
&= - \sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\;C_{total}\;\left({\hat{e}_{-1}}{\hat{e}_{1}}+{\hat{e}_{1}}{\hat{e}_{-1}}-{\hat{e}_{0}}{\hat{e}_{0}}\right).
\end{split}
\label{eq:19} 
\end{equation}
Equation (\ref{eq:19}) gives the total gradient of the field over a
sphere due to the distribution of the charges. For quadrupolar fluids
the total charge within a sphere is zero, therefore
$ \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = 0 $.  Hence the quadrupolar
polarization produces zero net gradient of the field inside the
sphere.

\section{Applied field or field gradient}
\label{Ap:fieldOrGradient}

To satisfy the condition $ \nabla . E = 0 $, within the box of molecules we have taken electrostatic potential in the following form
\begin{equation}
\begin{split}
\phi(x, y, z) =\; &-g_o \left(\frac{1}{2}(a_1\;b_1 - \frac{cos\psi}{3})\;x^2+\frac{1}{2}(a_2\;b_2 - \frac{cos\psi}{3})\;y^2 + \frac{1}{2}(a_3\;b_3 - \frac{cos\psi}{3})\;z^2 \right. \\
& \left. + \frac{(a_1\;b_2 + a_2\;b_1)}{2} x\;y + \frac{(a_1\;b_3 + a_3\;b_1)}{2} x\;z +  \frac{(a_2\;b_3 + a_3\;b_2)}{2} y\;z \right),
\end{split}
\label{eq:appliedPotential}
\end{equation} 
where $a = (a_1, a_2, a_3)$ and $b = (b_1, b_2, b_3)$ are basis vectors  determine coefficients in x, y, and z direction. And $g_o$ and $\psi$ are overall strength of the potential and angle between basis vectors respectively. The electric field derived from the above potential is,
\[\bf{E} 
=\frac{g_o}{2} \left(\begin{array}{ccc}
2(a_1\; b_1 - \frac{cos\psi}{3})\;x \;+  (a_1\; b_2 \;+ a_2\; b_1)\;y + (a_1\; b_3 \;+ a_3\; b_1)\;z \\
 (a_2\; b_1 \;+ a_1\; b_2)\;x + 2(a_2\; b_2 \;- \frac{cos\psi}{3})\;y +  (a_2\; b_3 \;+ a_3\; b_3)\;z \\
(a_3\; b_1 \;+ a_3\; b_2)\;x +  (a_3\; b_2 \;+ a_2\; b_3)y + 2(a_3\; b_3 \;- \frac{cos\psi}{3})\;z 
\end{array} \right).\] 
The gradient of the applied field derived from the potential can be written in the following form,
\[\nabla\bf{E} 
= \frac{g_o}{2}\left(\begin{array}{ccc}
2(a_1\; b_1 - \frac{cos\psi}{3}) &  (a_1\; b_2 \;+ a_2\; b_1) & (a_1\; b_3 \;+ a_3\; b_1)\;z \\
 (a_2\; b_1 \;+ a_1\; b_2) & 2(a_2\; b_2 \;- \frac{cos\psi}{3}) & (a_2\; b_3 \;+ a_3\; b_3)\;z \\
(a_3\; b_1 \;+ a_3\; b_2) & (a_3\; b_2 \;+ a_2\; b_3) & 2(a_3\; b_3 \;- \frac{cos\psi}{3})\;z 
\end{array} \right).\] 
\newpage

\bibliography{dielectric_new}

\end{document}
Revision:	4401
Committed:	Mon Mar 28 20:25:13 2016 UTC (8 years, 3 months ago) by gezelter
Content type:	application/x-tex
File size:	73355 byte(s)
Log Message:	Latest updates