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}

\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,Fennell06,ZMM papers}
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 zeroing
out of the net multipole moments inside the cutoff spheres.\cite{ZMM
  papers}

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 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. 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 dielectric
response 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 the
simulation depends on the details of the electrostatic interaction
methods that were utilized 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 modern simulations, bulk materials are usually 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 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 that the expression for the
dielectric constant reduces to widely-used \textit{conducting
  boundary} formula for real-space cutoff methods that have first
derivatives that vanish at the cutoff sphere.

One of the primary topics of this paper will be derivation of these
correction factors for the 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 correction
formula for different real-space methods has not yet been developed.

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

We have also calculated the screening behavior for two ions immersed
in multipolar fluids to estimate the distance dependence of dielectric
screening in both dipolar and quadrupolar fluids.  We have used three
different methods to compare our 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}

Using external perturbations, the net polarization of the system is
observed as a linear response of the applied field perturbation, where
proportionality constant is determined by the electrostatic
interaction between the electrostatic multipoles at a given
temperature. The fluctuation formula observes the time average
fluctuation of the multipolar moment as a function of temperature. The
average fluctuation value of the system is determined by the
multipole-multipole interactions between molecules at a given
temperature. Since the expression of the electrostatic interaction
energy, force, and torque in the real space electrostatic methods are
different from their original definition, both fluctuation and
external field perturbation formula must also be modified
accordingly. 

The potential of mean force allows calculation of an effective
dielectric constant or screening factor from the potential energy
between ions before and after dielectric material is introduced.

\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{p_{mol}} = \epsilon_0 \alpha_p \mathbf{E},
\end{equation}
where $ \alpha_p = \frac{1}{3}\frac{{p_o}^2}{\epsilon_0 k_B T}$ is the
contribution to molecular polarizability due solely to reorientation
dynamics. The orientational polarization depends inversely on the
temperature as the applied field must overcome the thermal motion.

Likewise, a condensed phase system of dipolar molecules will also
polarize along the direction of an applied field. The net dipolar
polarization $\mathbf{P}$ of the system is,
\begin{equation}
\textbf{P} = \epsilon_o \alpha_{D} \mathbf{E}. 
\label{pertDipole}
\end{equation} 
The constant $\alpha_D$ is a macroscopic polarizability, which is an
emergent property of the dipolar fluid in a given density and
temperature and is related to the static dielectric constant.

\subsection{Fluctuation Formula}

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 orientationally-polarized molecules 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{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 molecules, the net dipole moment, $\bf{M}$ at thermal
equilibrium of temperature T in the presence of applied field
$\bf{E}$, the average dipolar polarization can be expressed in terms
of fluctuation of the net dipole moment,
\begin{equation}
\braket{\bf{P}} = \epsilon_o \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E}
\label{flucDipole}
\end{equation}
This is similar to the formula for Boltzmann average of single dipolar
molecule. Here $\braket{\bf{P}}$ is average polarization and
$ \braket{\textbf{M}^2}-{\braket{\textbf{M}}}^2$ is the net dipole
fluctuation at temperature T. For the limiting case
$\textbf{E} \rightarrow 0 $, the ensemble average of both the net
dipole moment $\braket{\textbf{M}}$ and dipolar polarization
$\braket{\bf{P}}$ tends to vanish but $\braket{\bf{M}^2}$ will still
be non-zero. The dipolar macroscopic polarizability can be written as,
\begin{equation}
\alpha_D = \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}
\end{equation}
This is a macroscopic property of dipolar material which is true even
if applied field $ \textbf{E} \rightarrow 0 $. 
\subsection{Correction Factors}
In the presence of an external field $ \textbf{E}$ polarization $\textbf{E}$ will be induced in a dipolar system. The total electrostatic field (or Maxwell electric field) at point $\bf{r}$ in a system is,\cite{NeumannI83}
\begin{equation}
\textbf{E}(\textbf{r}) = \textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int d^3r' \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}.
\end{equation}

We can 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)
\label{dipole-chargeTensor}
\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) 
\label{dipole-diopleTensor}
\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.
 
Using constitutive relation, the polarization density $\textbf{P}(\textbf{r})$ is given by,
\begin{equation}
\textbf{P}(\textbf{r}) = \epsilon_o\; \chi^*_D \left(\textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int d^3r' \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}\right).
\label{constDipole}
\end{equation}
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}(\kappa) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(1-  \frac{3}{4\pi}\;\frac{\chi^*_D}{\chi^*_D + 3}\; \textbf{T}({\kappa})\right)^{-1}\textbf{E}^o({\kappa}).
\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}
where $\epsilon_{CB}$ is dielectric constant obtained from conducting boundary condition. Equation (\ref{correctionFormula}) calculates actual dielectric constant from the dielectric constant obtained from the conducting boundary condition (which can be obtained directly from the simulation) using correction factor. The correction factor is different for different real-space cutoff methods. The expression for correction factor assuming a single point dipole or two closely spaced point charges for SP, GSF, and TSF method is 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 & $A_\mathrm{charges}$ & $A_\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}$ \\
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}$ & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ \\\hline
\end{tabular}%
}
\end{table}


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

In analogy to the dipolar case, a molecule with a permanent
quadrupole, $q_{\alpha\beta}$ will also align in the presence of a
uniform electric field gradient $\vec{\nabla}\mathbf{E}$.  The average
polarization is given by
\begin{equation}
\braket{q_{\alpha\beta}} = \frac{1}{3} Tr(q^*) \delta_{\alpha\beta} + \alpha_q \partial_\alpha E_\beta,
\end{equation}
where $ \alpha_q = \frac{{\bar{q_o}}^2}{15k_BT} $ is a molecular
quadrupolarizablity and
${\bar{q_o}}^2=
3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$
is a square of the net quadrupole moment of a molecule.

Similarly, in the presence of and external field gradient, a system of
quadrupolar molecule polarizes along the direction of applied field
gradient, giving a system quadrupolar polarization,
be given by,
\begin{equation}
\begin{split}
& {Q}_{\alpha\beta} = \frac{1}{3} Tr({Q}) \delta_{\alpha\beta} +  \epsilon_o \alpha_Q  \partial_{\alpha} E^o_{\beta} 
\\ & or \\
& \frac{1}{3}\;\Theta_{\alpha\beta} =  \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta} 
\end{split}
\label{pertQuad}
\end{equation} 
where $Q_{\alpha\beta}$ is a tensor component of the traced quadrupolar moment of the system, $ \alpha_Q$ is a macroscopic quadrupolarizability has a dimension of $length^{-2}$, and $\Theta_{\alpha\beta} = 3Q_{\alpha\beta}-Tr(Q) $ is the traceless component of the quadrupole moment.    
\subsection{Fluctuation Formula}

Fluctuation formulae for quadrupolar fluids were developed by Logan
\textit{et al.}\cite{LoganI81,LoganII82,LoganIII82}

In analogy with the dipolar case, a fluctuation formula can  be
written for a system comprising quadrupolar molecules,
\begin{equation}
\braket{Q_{\alpha\beta}} = \frac{1}{3} Tr(\textbf{Q}) \delta_{\alpha\beta} + \epsilon_o \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T}{\partial_\alpha E^o_\beta}
\label{flucQuad}
\end{equation}
where $Q_{\alpha\beta}$ is a component of system quadrupole moment,  $\bf{Q}$ is net quadrupolar moment which can be expressed as $\textbf{Q}^2 =3Q_{\alpha\beta}Q_{\alpha\beta}-(Tr\textbf{Q})^2 $. The macroscopic quadrupolarizability is given by,
\begin{equation}
\alpha_Q = \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T}
\label{propConstQuad}
\end{equation}

\subsection{Correction Factors}
In the presence of the field gradient $\partial_\alpha {E}_\beta $, a
non-vanishing quadrupolar polarization (quadrupole moment per unit
volume) $\bar{Q}_{\alpha\beta}$ will be induced in the entire volume
of a sample. The total field at any point $\vec{r}$ in the sample is
given by,
\begin{equation}
\partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{4\pi \epsilon_o}\int T_{\alpha\beta\gamma\delta}(|{\textbf{r}-\textbf{r}'}|)\;{Q}_{\gamma\delta}(\textbf{r}')\; d^3r'
\label{gradMaxwell}
\end{equation}
where $\partial_\alpha {E^o}_\beta$ is the applied field gradient and $ T_{\alpha\beta\gamma\delta}$ is the quadrupole-quadrupole interaction tensor. We can represent quadrupole as a group of four closely spaced charges, two closely spaced point dipoles or single point quadrupole (see Fig. \ref{fig:quadrupolarFluid}). The quadrupole-quadrupole interaction tensor from the charge representation can obtained from the application of the four successive gradient operator to the electrostatic kernel $v(r)$. 

\begin{equation}
\begin{split}
T_{\alpha\beta\gamma\delta}(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'(r)}{r^3} + \frac{v''(r)}{r^2}\right) 
\\ &+ \left(\delta_{\alpha\beta} r_\gamma r_\delta + 5 \; permutations \right) \left(\frac{3v'(r)}{r^5}-\frac{3v''(r)}{r^4} + \frac{v'''(r)}{r^3}\right) 
\\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(-\frac{15v'(r)}{r^7}+\frac{15v''(r)}{r^6}-\frac{6v'''(r)}{r^5} + \frac{v''''(r)}{r^4}\right), 
\end{split}
\label{quadCharge}
\end{equation}
where $v(r)$ can either be electrostatic kernel for spherical truncation or one of the modified (Wolf or DSF) method. Similarly in point dipole representation the qaudrupole-quadrupole interaction tensor can be obtained from the applications of the two successive gradient in the dipole-dipole interaction tensor,

\begin{equation}
\begin{split}
T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \;v_{\gamma\delta}(r) 
\\ &= \delta_{\alpha\beta}\delta_{\gamma\delta} \frac{v'_{21}(r)}{r} + \left(\delta_{\alpha\gamma}\delta_{\beta\delta}+ \delta_{\alpha\delta}\delta_{\beta\gamma}\right)\frac{v_{22}(r)}{r^2} 
\\ &+ \delta_{\gamma\delta} r_\alpha r_\beta \left(\frac{v''_{21}(r)}{r^2}-\frac{v'_{21}(r)}{r^3} \right)
\\ &+\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'_{22}(r)}{r^3}-\frac{2v_{22}(r)}{r^4}\right) 
\\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(\frac{v''_{22}(r)}{r^4}-\frac{5v'_{22}(r)}{r^5}+\frac{8v_{22}(r)}{r^6}\right), 
\end{split}
\label{quadDip}
\end{equation}
where $v_{\gamma\delta}(r)$ is the electrostatic dipole-dipole interaction tensor, which is different for different electrostatic cut off methods. Similarly $v_{21}(r) \;and\; v_{22}(r)$ are the radial function for different real space cutoff methods defined in Paper I of the series.\cite{PaperI} Using point quadrupole representation the quadrupole-quadrupole interaction can be constructed as,  
\begin{equation}
\begin{split}
T_{\alpha\beta\gamma\delta}(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} \\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(\frac{v_{43}(r)}{r^4}\right), 
\end{split}
\label{quadRadial}
\end{equation}
where $v_{41}(r),\; v_{42}(r), \; \text{and} \; v_{43}(r)$ are defined in Paper I  of the series. \cite{PaperI} They have  different functional forms for different electrostatic cutoff methods.
\begin{figure}
\includegraphics[width=3in]{QuadrupoleFigure}
\caption{With the real-space electrostatic methods, the effective
  quadrupolar tensor, $\mathbf{T}_{\alpha\beta\gamma\delta}(r)$, governing interactions between molecular quadrupoles can be represented by interaction of charges, point dipoles or single point quadrupoles.}
\label{fig:quadrupolarFluid}
\end{figure}
The integral in equation (\ref{gradMaxwell}) can be divided into two parts, $|\textbf{r}-\textbf{r}'|\rightarrow 0 $ and $|\textbf{r}-\textbf{r}'|> 0$. Since the total
field gradient due to quadrupolar fluid vanishes at the singularity (see Appendix \ref{singularQuad}), equation (\ref{gradMaxwell}) can be written as,
\begin{equation}
\partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) +
  \frac{1}{4\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 }
  T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{Q}_{\gamma\delta}(\textbf{r}')\;
  d^3r'. 
\end{equation}
If $\textbf{r} = \textbf{r}'$ is excluded from the integration, the gradient of the electric can be expressed in terms of traceless quadrupole moment as below, \cite{LoganI81}
\begin{equation}
\partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 } T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{\Theta}_{\gamma\delta}(\textbf{r}')\; d^3r',
\end{equation}
where $\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$
is the traceless quadrupole moment. The total quadrupolar polarization is written as,
\begin{equation}
{Q}_{\alpha\beta}(\textbf{r}) =  \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q})+\epsilon_o {\chi}^*_Q\;\partial_\alpha E_\beta(\textbf{r}),
\label{constQaud}
\end{equation}
In the equation (\ref{constQaud}), $\partial_{\alpha}E_{\beta}$ is Maxwell field gradient and ${\chi}^*_Q$ is the actual quadrupolar susceptibility of the fluid which is different than the proportionality constant $\chi_q $ in the equation (\ref{propConstQuad}). In terms of traceless quadrupole moment, equation (\ref{constQaud}) can be written as,
\begin{equation}
\frac{1}{3}{\Theta}_{\alpha\beta}(\textbf{r}) =  \epsilon_o {\chi}^*_Q \; \partial_\alpha E_\beta (\textbf{r})=  \epsilon_o {\chi}^*_Q \left(\partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 } T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{\Theta}_{\gamma\delta}(\textbf{r}')\; d^3r'\right)
\end{equation}
For toroidal boundary conditions, both $\partial_\alpha E_\beta$ and
${\Theta}_{\alpha\beta}$ are uniform over the entire space. After
performing a Fourier transform (see the Appendix in the Neumann's Paper  \cite{NeumannI83}) we get,
\begin{equation}
\frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa})=
\epsilon_o {\chi}^*_Q \;\left[{\partial_\alpha
    {E^o}_\beta}({\kappa})+ \frac{1}{12\pi
    \epsilon_o}\;{T}_{\alpha\beta\gamma\delta}({\kappa})\;
 {{\Theta}}_{\gamma\delta}({\kappa})\right] 
\end{equation} 
Since the quadrupolar polarization is in the direction of the applied
field, we can write
${{\Theta}}_{\gamma\delta}({\kappa}) =
{{\Theta}}_{\alpha\beta}({\kappa})$
and
${T}_{\alpha\beta\gamma\delta}({\kappa}) =
{T}_{\alpha\beta\alpha\beta}({\kappa})$. Therefore we can consider each component of the interaction tensor as scalar and perform calculation.
\begin{equation}
\begin{split}
\frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa}) &= \epsilon_o {\chi}^*_Q \left[{\partial_\alpha E^o_\beta}({\kappa})+ \frac{1}{12\pi \epsilon_o}{T}_{\alpha\beta\alpha\beta}({\kappa})\;{{\Theta}}_{\alpha\beta}({\kappa})\right] \\
&= \epsilon_o {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\;
 {T}_{\alpha\beta\alpha\beta}({\kappa})\right)^{-1}
{\partial_\alpha E^o_\beta}({\kappa}) 
\end{split}
\label{fourierQuad}
\end{equation} 
If the field gradient is homogeneous over the
entire volume, ${\partial_ \alpha E_\beta}({\kappa}) = 0 $ except at
$ {\kappa} = 0$, hence it is sufficient to know
${T}_{\alpha\beta\alpha\beta}({\kappa})$ at $ {\kappa} =
0$. Therefore equation (\ref{fourierQuad}) can be written as,
\begin{equation}
\begin{split}
\frac{1}{3}{{\Theta}}_{\alpha\beta}({0}) &= \epsilon_o {\chi}^*_Q\; \left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1} \partial_\alpha E^o_\beta({0})
\end{split}
\label{fourierQuad2}
\end{equation}
where $ {T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as,
\begin{equation}
{T}_{\alpha\beta\alpha\beta}({0}) = \int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r
\label{realTensorQaud}
\end{equation}

In terms of traced quadrupole moment equation (\ref{fourierQuad2}) can be written as,
\begin{equation}
{{Q}}_{\alpha\beta} = \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q}) + \epsilon_o\; {\chi}^*_Q\left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}\; \partial_\alpha E^o_\beta
\label{tracedConstQuad}
\end{equation}
Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we get,
\begin{equation}
\begin{split}
&\frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\; =\; {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}, \\ 
&{\chi}^*_Q \;=\; \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{1}{4\pi} \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}
\end{split}
\end{equation}
Finally the quadrupolar susceptibility cab be written in terms of quadrupolar correction factor ($A_{quad}$) as below, 
\begin{equation}
{\chi}^*_Q \;=\; \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\; A_{quad}\right)^{-1} = \alpha_Q\left(1 + \alpha_Q\; A_{quad}\right)^{-1}
\label{eq:quadrupolarSusceptiblity}
\end{equation}
where $A_{quad} = \frac{1}{4\pi}\int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r $ has dimension of the $length^{-2}$ is different for different cutoff methods which is listed in Table \ref{tab:B}. The dielectric constant associated with the quadrupolar susceptibility is defined as,\cite{Ernst92}

\begin{equation}
\epsilon = 1 + \chi^*_Q\; G = 1 + G \; \alpha_Q\left(1 + \alpha_Q\;  A_{quad}\right)^{-1}
\label{eq:dielectricFromQuadrupoles}
\end{equation}
where $G = \frac{\displaystyle\int_V |\partial_\alpha E^o_\beta|^2 d^3r}{\displaystyle\int_V{|E^o|}^2 d^3r}$ is a geometrical factor depends on the nature of the external field perturbation. This is true when the quadrupolar fluid is homogeneous over the sample. Since quadrupolar molecule couple with the gradient of the field, the distribution of the quadrupoles is inhomogeneous for varying field gradient. Hence the distribution function should also be taken into account to calculate actual geometrical factor in the presence of non-uniform gradient field. Therefore,
\begin{equation}
G = \frac{\displaystyle\int_V\; g(r, \theta, \phi)\; |\partial_\alpha E^o_\beta|^2 d^3r}{\displaystyle\int_V{|E^o|}^2 d^3r}
\label{eq:geometricalFactor}
\end{equation}
where $g(r,\theta, \phi)$ is a distribution function of the quadrupoles in with respect to origin at the center of line joining two probe charges.  
\begin{table}
  \caption{Expressions for the quadrupolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter
    ($\alpha$) and the cutoff radius ($r_c$). The dimension of the correction factor is $ length^{-2}$ in case of quadrupolar fluid.}
\label{tab:B}
\centering
\resizebox{\columnwidth}{!}{%

\begin{tabular}{l|c|c|c|c|}
       
Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ &$A_\mathrm{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)$ \\\hline
\end{tabular}%
}
\end{table}

\section{Screening of Charges by Multipolar Fluids}
In this method, we will measure the interaction between a positive and negative charge at varying distances after introducing a dipolar (or quadrupolar) material between them. The potential of mean force (PMF) between two ions in a liquid is obtained by constraining their distance and measuring the mean constraint force required to hold them at a fixed distance $r.$ The PMF  is obtained from a sequence of simulations can be expressed as \cite{Wilfred07},
\begin{equation}
w(r) = \int_{r_o}^{r}\braket{\frac{\partial f}{\partial r'}}dr' + 2kT \log(r/r_o) + w(r_o),
\label{eq:pmf}
\end{equation}
where $\braket{\partial f/\partial r'}$ is the mean constraint force, $2kT \log(r/r_o)$ is the Fixman factor, and $r_o$ is the reference position. The potential energy between two charges at separation $r$ is given by,
\begin{equation}
w(r) = \frac{1}{4\pi\epsilon}\frac{q_1q_2}{r},
\end{equation}
where $\epsilon$ is the dielectric constant for the fluid.

The quadrupole molecule can couple with the gradient of the electric field and the two opposite point charges produces non-zero value of both electric field as well as gradient of the field. Therefore, this simulation set up can be used to determine the dielectric constant for both dipolar and quadrupolar fluid.
\label{sec:PMF}
\section{Simulation Methodology}
We have used three different simulation methods: i) external field perturbation, ii) fluctuation formula, and iii) potential of mean force (PMF), to calculate dielectric properties for dipolar and quadrupolar fluid. For the dipolar system we calculated macroscopic polarzability using first two methods and derived the dielectric constant using polarizability and correction factor (see equation \ref{correctionFormula}). Similarly we used equation (\ref{eq:pmf}) to calculate screening factor from dipolar fluid using PMF method. For quadrupolar fluid, we have calculated quadrupolarizablity using fluctuation formula and external field perturbation and derived quadrupolar susceptibility using quadrupolarizability and correction factor (equation \ref{eq:quadrupolarSusceptiblity}). Since dielectric constant due to quadrupolar fluid depends on the nature of gradient of the field applied in the system, we have used geometrical factor (in equation \ref{eq:geometricalFactor}) and quadrupolar susceptibility to evaluate dielectric constant for two ions dissolved quadrupolar fluid (see equation \ref{eq:dielectricFromQuadrupoles}) . The the dielectric constant evaluated using equation (\ref{eq:dielectricFromQuadrupoles}) has been compared with the result evaluated from PMF method (i.e. equation \ref{eq:pmf}). We have also used three different test systems for both dipolar and quadrupolar fluids to calculate dielectric properties. The parameters used in the test systems are given in table \ref{Tab:C}.

\begin{table}
  \caption{\label{Tab:C}}
\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
    Stockmayer 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}

First test system consists of point dipolar or quadrupolar molecules in the presence of constant field or gradient field. Since there is no isolated charge within the system, the divergence of the field should be zero $ i.e. \vec{\nabla} .\vec{E} = 0$. This condition is satisfied by selecting applied potential as described in Appendix \ref{Ap:fieldOrGradient}.  When constant electric field or field gradient applied to the system, the molecules align along the direction of the applied field. We have calculated ensemble average of the box dipole or quadrupole moment as a response field or field gradient. Similarly the  macroscopic polarizability of the system is derived using ratio between system multipolar moment and applied field or field gradient. This method works properly when the system is at the linear response region of field or field gradient.    

Second test system consists of box of point dipolar or quadrupolar molecules is simulated for 1 ns in NVE ensemble after equilibration in the absence of any external perturbation. The fluctuation of the ensemble average of the box multipolar moment, $\braket{A^2} - \braket{A}^2 $ where A is box dipolar or quadrupolar moment,  is measured at the fixed temperature and density for a given multipolar fluid. Finally the macroscopic polarizability of the system at a particular density is derived using equation (\ref{flucQuad}).

Final system consists of dipolar or quadrupolar fluids with two oppositely charged ions immersed in it. These ions are constraint to be at fixed distance throughout the simulation. We run separate simulations for the different constraint distances. We calculated the screening factor using ratio between the force between the two ions in the absence of medium and the average constraint force during the simulation. Since the constraint force is pretty noisy we run each simulation for long time ($\sim 5$ ns) to reduce simulation error.

\subsection{Implementation}
We have used real-space electrostatic methods implemented in OpenMD \cite{openmd2.3} software to evaluate electrostatic interactions between the molecules. In our simulations we used all three different real-space electrostatic methods: SP, GSF, and TSF developed in the previous paper \cite{PaperI} in the series. The radius of the cutoff sphere is taken to be $12 \r{A}$. Each real space method can be tuned using different values of damping parameter. We have selected ten different values of damping parameter (unit-${\r{A}}^{-1}$); 0.0, 0.05, 0.1, 0.15, 0.175, 0.2, 0.225, 0.25, 0.3, and 0.35 in our simulations for dipolar system. Since quadrupolar interactions are less sensitive to damping parameter as compared to dipole we selected only four values of dampig alpha 0.0, 0.1 0.2, and 0.3 $\AA-1$ in our simulation. The short range interaction in the simulations is incorporated using 6-12 Lennard Jones interactions. 

To derive the box multipolar (dipolar or quadrupolar) moment, we added the component each individual molecule and taken ensemble average of each snapshot in the entire simulation. The first component of the fluctuation of the dipolar moment is derived by using relation $\braket{M^2} = \braket{{M_x}^2 + {M_y}^2 + {M_z}^2}$, where $M_x$, $M_y$, and $M_z $ are x, y and z components of the box dipolar moment. Similarly the first term in the quadrupolar system is derived using relation $ \braket{Q^2} = \braket{3 Q:Q - \mathrm{Tr}Q^2} $, where $ Q $ is the box quadrupole moment, double dot represent the outer product of the quadrupolar matrices, and $TrQ$ is the trace of the box quadrupolar moment. The second component of the fluctuation formula has been derived using square of the ensemble average of the box dipole moment. The applied constant field or field gradient in the test systems has been taken in the form described in the Appendix \ref{Ap:fieldOrGradient}.
\subsection{Model systems}
To evaluate dielectric properties for dipolar systems using perturbation and fluctuation formula methods, we have taken system of 2048 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} $. Test systems are equilibrated for 500 ps and run for $1\; ns$ and components of box dipole moment are obtained at every femtosecond. The systems have been run in the presence of constant external field from $ 0 - 10\; \times\; 10^{-4}\;V/{\r{A}}$ in the step of $ 10 ^{-4}\; V/\r{A}$ for each simulation. For pmf method, Two dipolar molecules in the above system are converted into $q+$ and $q-$ ions and constrained to remain in fixed distance in simulation. The constrained distance is varied from $5\;\r{A} - 12\; \r{A} $ for different simulations. In pmf method all simulations are equilibrated for 500 ps in NVT ensemble and run for 5 ns in NVE ensemble to print constraint force at an interval of 20 fs.   

Quadrupolar systems consists 4000 linear point quadrupoles with density $ 2.338\; g/cm^3$ at temperature $ 500\; ^oK $. For both perturbation and fluctuation methods, test systems are equalibrated for 200 ps in NVT ensemble and run for 500 ps in NVE ensemble. To find the ensemble average of the box quadrupole moment and fluctuation, the components of box quadrupole moments are printed every 100 fs. Each simulations repeated at different values of applied gradients from $ 0 - 9 \times 10^{-2}\; V/\AA^2 $. To calculate dielectric constant using pmf method, two ions in the systems are converted into $q+$ and $q-$ ions and constrained to remain at fixed distance in the simulation. These constraint distances are varied from $5\;\AA - 12\; \AA $ at the step of $0.1\; \AA $ for different simulations. For calculating dielectric constant, the test systems are run for 500 ps to equlibrate and run for 5 ns to print constraint force at a time interval of 20 fs.
      
\section{Results}
\subsection{Dipolar fluid}
The macroscopic polarizability ($\alpha_D$) for Stockmayer fluid evaluated using perturbation (left) and fluctuation (right) methods as shown in figure \ref{fig:dielectricDipole} . The polarizability obtained from the both perturbation and fluctuation methods show excellent agreement with each other. The result shows that polarizability hugely depends on the values of the damping parameter in the case of Shifted Potential (SP) and Gradient Shifted force (GSF) methods. But the results obtained using Taylor Shifted Force (TSF) show slight dependence on the damping alpha.  To get correct dielectric properties from the polarizability we used a correction factor as we discussed in section \ref{sec:corrFactor}. The correction factor for dipolar fluids ($\mathrm{A}_\mathrm{dipole}$ ) for the different values of damping parameter is also plotted in figure \ref{fig:dielectricDipole}. For TSF method, $ \mathrm{A}_\mathrm{dipole} = 1 $ for all the values of damping alpha. Therefore we do not need to perform any correction to obtain dielectric constant from the polarizability. The values of $\mathrm{A}_\mathrm{dipole} $ varies with the damping parameter in the case of SP and GSF methods therefore correction factor should be incorporated to calculate the dielectric constant ($\epsilon$) from the polarzability. The dielectric constant obtained from the polarizability and correction factor is also plotted in the figure \ref{fig:dielectricDipole}. The dielectric constant evaluated using SP and GSF method shows good agreement with the previous simulation result (straight dotted line in the lower panels)\cite{NeumannI83} for the damping parameter from $0.25\; \mathrm{to}\; 0.3$ $\AA-1$. For lower value of the damping alpha the dielectric constant is much deviated away from the previous simulation result. To investigate furthermore we performed inverse transformation of expected dielectric constant \cite{NeumannI83}, using correction factor and calculated expected polarizability for the simulation, which is represented by the dotted line in the upper-left and upper-right panels. In the plot we can clearly see that expected ploarizability for the SP and GSF methods are very close to results obtained from the simulations. Sometime even they are within the range of the error-bar but when the dielectric constant is calculated using simulation result and correction formula, it is found to be deviated much away from the expected result for the lower values of the damping alpha. This tell us that the correction formula (equation \ref{correctionFormula}) for the dipolar fluid is very sensitive when the value of $\mathrm{A}_\mathrm{dipole}$ is deviated much away from 1 for the small values of damping parameter. Therefore it is always better to select damping alpha from 0.25 to 0.3 $\AA^{-1}$ to evaluate dielectric constant using perturbation and fluctuation formula while using SP and GSF methods in simulation.    

\begin{figure}
\includegraphics[width=3in]{dielectricFinal_Dipole.pdf}
\caption{In the figure, $\alpha_D$, $\mathrm{A}_\mathrm{dipole}$, and $\epsilon$  are polarizability, correction factor, and dielectric constant for Stockmayer fluid. Plots in the left panel show results for (a) perturbation method and right (b) fluctuation method.}
\label{fig:dielectricDipole}
\end{figure}

We have also evaluated the screening factor between two oppositely charged ions when they are placed in the Stockmayer fluid as shown in figure \ref{fig:ScreeningFactor_Dipole}. The screening factor have been evaluated using equation \ref{eq:pmf}. This screening factor is analogous to dielectric screening but there is subtle differences between them. Actually screening factor measures a local property of the ions in the fluid and depends on ion-dipole and dipole-dipole interactions. These interactions depends on the distance between ions and electrostatic interaction methods implemented in the simulations. But we can expect screening to be close to dielectric constant when the field due to ions is very gentle in perturbation. It is possible when distance between ions are far apart and values of damping parameter relatively higher. In the figure \ref{fig:ScreeningFactor_Dipole} we can observe that for the higher value of damping alpha i.e $\alpha = 0.2 \AA^{-1}$ and $0.3 \AA^{-1}$ and large separation between ions, the screening factor approaches to the dielectric constant. We can also  observe that for TSF method has got higher screening factor as compared SP and GSF method even for lower damping alpha. It is because of dipole-dipole interactions do not need any correction factor in the case of TSF as we discussed earlier. But presence of ion-dipole interaction causing local effect even in the case of TSF method for small ions separation.

\begin{figure}
\includegraphics[width=3in]{ScreeningFactor_Dipole.pdf}
\caption{Figure shows screening factor between two oppositely charged ions immersed in Stockmayer fluid as a function of ions-separation for the different values damping alpha and electrostatic methods.}
\label{fig:ScreeningFactor_Dipole}
\end{figure}


\subsection{Quadrupolar fluid}
The polarizability ($\alpha_Q$), correction factor ($\mathrm{A_{quad}}$), and susceptibility ($\chi$) for the quadrupolar fluid is plotted against damping parameter in the figure \ref{fig:dielectricQuad}. Different than the dipolar fluid, the polarizability and susceptibility for the quadrupolar fluid have a dimension of $\AA ^2$. Although susceptibility has got dimension, it is a measure of macroscopic quadrupolar properties at a given temperature and density \cite{JeonI03, JeonII03}. The left panel in the plot showing results obtained from  the perturbation whereas the results from the fluctuation formula are plotted in the right panel. In the figure, we can observe that the polarizablity evaluated from the perturbation (top-left) shows excellent agreement with the result obtained from the fluctuation formula (top-right). The susceptibility for the quadrupolar fluid is obtained from quadrupolarizability and correction factor using equation \ref{eq:quadrupolarSusceptiblity} which is plotted at the bottom-left and bottom-right in the figure \ref{fig:dielectricQuad}. All three methods: SP, GSF, and TSF tends to produce same value of susceptibility for the given value damping alpha. This shows that susceptibility derived using correction factor and polarizablity is tends to be independent of the electrostatic method implemented in the simulation for quadrupolar fluid.
\begin{figure}
\includegraphics[width=3in]{polarizabilityFinal_Quad.pdf}
\caption{In the figure, $\alpha_Q$, $\mathrm{A}_\mathrm{quad}$, and $\chi$  are polarizability, correction factor, and susceptibility for the quadrupolar fluid. Plots in the left panel show results for (a) perturbation method and right (b) fluctuation method.}
\label{fig:dielectricQuad}
\end{figure}

The actual test of the susceptibility can be made by comparing the results obtained from the perturbation and fluctuation with the result calculated from more direct, potential of mean force (PMF), method. For this purpose, we have calculated the expected dielectric constant from the susceptibility and geometrical factor using equation \ref{eq:dielectricFromQuadrupoles}. Since the dielectric constant for quadrupolar fluid is geometry-dependent, it is not a actual measure of the property of the quadrupolar fluid. This geometrical factor varies with the variation in ions separation as well as damping parameter therefore the dielectric constant can be refereed as screening factor. The screening factor plotted against ions separation, evaluated from the perturbation and fluctuation method are shown in right and central panel of the figure \ref{fig:screenQuad}. Here the susceptibility is calculated from the simulation and the geometrical factor is evaluated analytically, using field and field-gradient produced by ions. The right hand panel shows screening factor plotted against ions separation, for the different values of the damping alpha, obtained from the PMF method. The figure clearly shows that the screening factor obtained from the perturbation and fluctuation formula shows good agreement with the screening factor calculated using PMF method. Since there is no huge differences in quadrupole-quadruople interactions for various real-space methods, we do not observe large differences in the screen factors for SP, GSF, and TSF methods.        

\begin{figure}
\includegraphics[width=3in]{ScreeningFactor_Quad.pdf}
\caption{}
\label{fig:screenQuad}
\end{figure}

\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:	4399
Committed:	Thu Mar 24 20:05:30 2016 UTC (8 years, 3 months ago) by gezelter
Content type:	application/x-tex
File size:	69031 byte(s)
Log Message:	latest version