trunk/multipole/dielectric_new.tex

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

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{multirow}
\usepackage{bm}% bold math
\usepackage{natbib}
\usepackage{times}
\usepackage{mathptmx}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}
\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}
  Note: This manuscript is a work in progress.   

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

\maketitle

\section{Introduction}

Over the past several years, there has been increasing interest in
pairwise methods for correcting electrostatic interactions in computer
simulations of condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06}
These techniques were initially developed from the observations and efforts of
Wolf {\it et al.}  and their work towards an $\mathcal{O}(N)$
Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is
able to obtain excellent agreement with Madelung energies in ionic
crystals.\cite{Wolf99} Later, Zahn \textit{et al.} and Fennell and Gezelter extended this method which incorporates Wolf's electrostatic energy and modified it to conserve the total energy in molecular dynamic simulation.\cite{Zahn02, Fennell06} In the previous two papers
we developed three new generalized real space methods: Shifted potential (SP), Gradeint shifted force (GSF), and Taylor shifted force (TSF).\cite{PaperI, PaperII} These methods evaluate electrostatic interactions for higher order multipoles (dipoles and quadrupoles) using finite cutoff sphere with the neutralization of the electrostatic moment within the cutoff sphere. Furthermore, extra terms added to the potential energy so that force and torque vanish smoothly at the cutoff radius. This ensures that the total energy is conserved in a molecular dynamic simulation.

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. Since dielectric properties are macroscopic properties, all interactions between molecules in an entire system are significantly important. But it is computationally infeasible to consider every interactions between molecules in the macroscopically large system. Therefore small molecular system with periodic boundary condition and finite cutoff region of interactions is usually considered in computer simulations. While calculating dielectric properties, the formula should be modified in such a way so that it can accommodate behaviour of electrostatic neutrality and smoothness of energy, force and torque at the cutoff radius. Previously many studies have been conducted to calculate dipolar and quadrupolar dielectric properties using computer simulations. \cite{Kirkwood39, Onsagar36,LoganI81, LoganII82, LoganIII82} But these methods do not specifically take account of the cutoff behavior common in real-space electrosatic methods. In 1983 Neumann proposed a general formula for evaluating dielectric properties for dipolar fluid using real-space cutoff methods. \cite{Neumann83} In the same year Steinhauser and Neumann used  this formula to evaluate the correct dielectric constant for the Stockmayer fluid using two different methods: Ewald-Kornfield (EK) and reaction field (RF) methods. \cite{Neumann-Steinhauser83} This formula contains a correction factor which is equal to $\frac{3}{4 \pi} $ times volume integral of the dipole-dipole interactions for a given electrostatic cutoff method (See equation \ref{dipole-diopleTensor}).\cite{Neumann83} Similarly Zahn \textit{et al.}\cite{Zahn02} also evaluated correction factor for dipole-dipole interaction using damped shifted charge-charge kernel (see equation \ref{dipole-chargeTensor}). This later generalized by Izvekove \textit{et al.}, which is equal to $\frac{3}{4 \pi} $. \cite{Izvekov:2008wo} When the correction factor is equal to  $\frac{3}{4 \pi} $, the expression for the dielectric constant reduces to widely-used \textit{conducting boundary} formula (see equation (\ref{correctionFormula})). Many studies have also been conducted to understand solvation theory using dielectric properties of  quadrupolar fluid.\cite{JeonI03, JeonII03, Chitanvis96}. But these studies do not use correction factor straight forwardly to evaluate correct dielectric properties for quadrupolar fluid.

In this paper we are proposing general consecutive formulas for calculating the dielectric properties for quadrupolar fluid. Furthermore we have also evaluated the correction factor for SP, GSF, and TSF method for both dipolar and quadrupolar fluid considering charge-charge, dipole-dipole or quadrupole-quadrupole interactions. The relation between quadrupolar susceptibility and dielectric constant is not straight forward for quadrupolar fluid  as in the dipolar case. The dielectric constant depends on the geometry of the external field perturbation.\cite{Ernst92} We have also calculated the geometrical factor for two ions immersed quadrupolar system to evaluate dielectric constant from the quadrupolar susceptibility. We have used three different methods: i) external field perturbation, ii) fluctuation formula,  and iii) the potential of mean force, to study dielectric properties of the dipolar and quadrupolar system. In the external field perturbation, 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 should also be modified accordingly. The potential of mean force method calculates dielectric constant from the potential energy between ions before and after dielectric material is introduced. All of these different methods for calculating dielectric properties will be discussed in detail in the following sections: \ref{subsec:perturbation}, \ref{subsec:fluctuation}, and \ref{sec:PMF}.

\section{Boltzmann average for orientational polarization}
The dielectric properties of the system is mainly arise from two different ways: i) the applied field distort the charge distributions so it produces an induced multipolar moment in each molecule; and ii) the applied field tends to line up originally randomly oriented molecular moment towards the direction of the applied field. In this study, we basically focus on the orientational contribution in the dielectric properties. If we consider a system of molecules in the presence of external field perturbation, the perturbation experienced by any molecule will not be only due to external field or field gradient but also due to the field or field gradient produced by the all other molecules in the system. In the following subsections \ref{subsec:boltzAverage-Dipole} and \ref{subsec:boltzAverage-Quad}, we will discuss about the molecular polarization only due to external field perturbation. The contribution of the field or field gradient due to all other molecules will be taken into account while calculating correction factor in the section \ref{sec:corrFactor}.

\subsection{Dipole}
\label{subsec:boltzAverage-Dipole}
Consider a system of molecules with permenent dipole moment $p_o$. In the absense of external field, thermal agitation makes dipole randomly oriented therefore there is no net dipole moment. But external field tends them to line up in the direction of applied field. Here we have considered net field acting due to all other molecules is considered to be zero.  Therefore the total Hamiltonian of the molecule is,\cite{Jackson98}

\begin{equation}
H = H_o - \bf{p_o} .\bf{E},
\end{equation}
where $H_o$ is a function of the internal coordinates of the molecule. Now Boltzmann average of the dipole moment is given by,
\begin{equation}
\braket{p_{mol}} = \frac{\displaystyle\int d\Omega\; p_o\; cos\theta\;  e^{\frac{p_oE\; cos\theta}{k_B T}}}{\displaystyle\int d\Omega\; e^{\frac{p_oE\;cos\theta}{k_B T}}},
\end{equation}
where $\bf{E}$ is selected along z-axis. If we consider applied field is small i.e. $\frac{p_oE\; cos\theta}{k_B T} << 1$ then we get,

\begin{equation}
\braket{p_{mol}}  \approx \frac{1}{3}\frac{{p_o}^2}{k_B T}E,
\end{equation}
where $ \alpha_p = \frac{1}{3}\frac{{p_o}^2}{k_B T}$ is a molecular polarizability. The orientational polarization depends inversely on the temperature and applied field must overcome the thermal agitation. 


\subsection{Quadrupole}
\label{subsec:boltzAverage-Quad}
Consider a system of molecules with permanent quadrupole moment $q_{\alpha\beta} $. The average quadrupole moment at temperature T in the presence of uniform applied field gradient is given by,\cite{AduGyamfi78, AduGyamfi81}
\begin{equation}
\braket{q_{\alpha\beta}} \;=\; \frac{\displaystyle\int d\Omega\; e^{-\frac{H}{k_B T}}q_{\alpha\beta}}{\displaystyle\int d\Omega\; e^{-\frac{H}{k_B T}}} \;=\; \frac{\displaystyle\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}q_{\alpha\beta}}{\displaystyle\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}},
\label{boltzQuad}
\end{equation}
where $\int d\Omega = \int_0^{2\pi} \int_0^\pi \int_0^{2\pi}
sin\theta\; d\theta\ d\phi\ d\psi$ is the integration over Euler
angles, $ H = H_o -q_{\mu\nu}\;\partial_\nu E_\mu $ is the energy of
a quadrupole in the gradient of the  
applied field and $ H_o$ is a function of internal coordinates of the molecule. The energy and quadrupole moment can be transformed into body frame using following relation, 
\begin{equation}
\begin{split}
&q_{\alpha\beta} = \eta_{\alpha\alpha'}\;\eta_{\beta\beta'}\;{q}^* _{\alpha'\beta'} \\
&H = H_o - q:\vec{\nabla}\vec{E} = H_o - q_{\mu\nu}\;\partial_\nu E_\mu = H_o -\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu.
\end{split}
\label{energyQuad}
\end{equation}
Here the starred tensors are the components in the body fixed
frame. Substituting equation (\ref{energyQuad}) in the equation (\ref{boltzQuad})
and taking linear terms in the expansion we get,
\begin{equation}
\braket{q_{\alpha\beta}} = \frac{ \int d\Omega \left(1 + \frac{\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu }{k_B T}\right)q_{\alpha\beta}}{ \int d\Omega \left(1 + \frac{\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu }{k_B T}\right)},
\end{equation}
where $\eta_{\alpha\alpha'}$ is the inverse of the rotation matrix that transforms
the body fixed co-ordinates to the space co-ordinates,
\[\eta_{\alpha\alpha'} 
= \left(\begin{array}{ccc}
cos\phi\; cos\psi - cos\theta\; sin\phi\; sin\psi & -cos\theta\; cos\psi\; sin\phi - cos\phi\; sin\psi & sin\theta\; sin\phi \\
cos\psi\; sin\phi + cos\theta\; cos\phi \; sin\psi & cos\theta\; cos\phi\; cos\psi - sin\phi\; sin\psi & -cos\phi\; sin\theta \\
sin\theta\; sin\psi & -cos\psi\; sin\theta & cos\theta
\end{array} \right).\]
Integration of 1st and 2nd terms in the denominator gives $8 \pi^2$
and $8 \pi^2 /3\;\vec{\nabla}.\vec{E}\; Tr(q^*) $ respectively. The
second term vanishes for charge free space
(i.e. $\vec{\nabla}.\vec{E} \; = \; 0)$. Similarly integration of the
1st term in the numerator produces
$8 \pi^2 /3\; Tr(q^*)\delta_{\alpha\beta}$ and the 2nd term produces
$8 \pi^2 /15k_B T (3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'} -
{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'})\partial_\alpha E_\beta$,
if $\vec{\nabla}.\vec{E} \; = \; 0$,
$ \partial_\alpha E_\beta = \partial_\beta E_\alpha$ and
${q}^*_{\alpha'\beta'}= {q}^*_{\beta'\alpha'}$. Therefore the
Boltzmann average of a quadrupole moment can be written as,

\begin{equation}
\braket{q_{\alpha\beta}}\; = \; \frac{1}{3} Tr(q^*)\;\delta_{\alpha\beta} + \frac{{\bar{q_o}}^2}{15k_BT}\;\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. 

\section{Macroscopic Polarizability}
\label{sec:MacPolarizablity}

If we consider a system of dipolar or quadrupolar fluid in the external field perturbation, the net polarization of the system will still be proportional to the applied field perturbation.\cite{Chitanvis96, Stern-Feller03, Salvchov14, Salvchov14_2} In simulation the net polarization of the system is determined by the interaction of molecule with all other molecules as well as external field perturbation. Therefore the macroscopic polarizablity obtained from the simulation always varies with nature of real-space electrostatic interaction methods implemented in the simulation. To determine a susceptibility or dielectric constant of the material (which is a actual physical property of the dipolar or quadrupolar fluid) from the macroscopic polarizablity, we need to incorporate the interaction between molecules which has been discussed in detail in section \ref{sec:corrFactor}. In this section we discuss about the two different methods of calculating macroscopic polarizablity for both dipolar and quadrupolar fluid. 

\subsection{External field perturbation}
\label{subsec:perturbation}
In the presence of uniform electric field $\textbf{E}^o$, a system of dipolar  molecules polarizes along the direction of the applied field (or field gradient). Therefore the net dipolar polarization $ \textbf{P}$ of the system is,
\begin{equation}
\textbf{P} = \epsilon_o \alpha_{D}\; \textbf{E}^o. 
\label{pertDipole}
\end{equation} 
The constant $\alpha_D$ is a macroscopic polarizability, which is a property of the dipolar fluid in a given density and temperature.

Similarly, in the presence of external field gradient the system of quadrupolar molecule polarizes along the direction of applied field gradient therefore the net quadrupolar polarization of the system can 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}
\label{subsec:fluctuation}
For a system of molecules with net dipolar moment $\bf{M}$ at thermal equilibrium of temperature T in the presence of applied field $\bf{E}^o$, the average dipolar polarization can be expressed in terms of fluctuation of the net dipole moment as below,\cite{Stern03}
\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}^o
\label{flucDipole}
\end{equation}
This is similar to the formula for boltzmann average of single dipolar molecule in the subsection \ref{subsec:boltzAverage-Dipole}. 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}^o \rightarrow 0 $, ensemble average of both 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}^o \rightarrow 0 $.  

Analogous formula can also be written for a system with 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}


\section{Potential of mean force}
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 as,
\begin{equation}
w(r) = \int_{\inf}^{r}\braket{\frac{\partial f}{\partial r'}}dr',
\end{equation}
where $\braket{\partial f/\partial r'}$ is the mean constraint force. 
Since the ions have a protecting Lennard-Jones (LJ) potential,
\begin{equation}
w(r) = w_{LJ}(r) + \frac{q_iq_j}{4\pi \epsilon_o \epsilon(r)}U_{method}(r).
\label{eq:pmf}
\end{equation}
Here $w_{LJ}$ is the PMF calculated without electrostatic interactions and $U_{method}(r)$ is the radial function for the charge-charge interaction, which is different for various real space truncation methods.

The quadrupole molecule can only couple with the gradient of the electric field and the region between two opposite point charges has both an electric field and a gradient of the electric field present. Therefore, this methodology should be usable to determine the dielectric constant for both the dipolar and quadrupolar fluid.
\label{sec:PMF}

\section{Correction factor}
\label{sec:corrFactor}
Since equations (\ref{pertDipole}, \ref{pertQuad}, \ref{flucDipole}, and \ref{flucQuad}) provide relation between polarization (dipolar or quadrupolar)  and applied field (uniform field or field gradient), $\chi_d$ (or $ \chi_q$) is actually a macroscopic polarizability (or quadrupolarizability), which is different than the dipolar (or quadrupolar) susceptibility of the fluid. Actual constitutive relation should have a relation between polarization and Maxwell field (or field gradient) at different point in the sample. We can obtain susceptibility of the fluid from its macroscopic polarizability using correction factor evaluated below.    
\subsection{Dipolar system}
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{Neumann83}
\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{sec:perturbation} and \ref{sec: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{Neumann83, 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{Adams:1980rt,Adams:1981fr,Neumann83} 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}
\subsection{Quadrupolar system}
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{Neumann83}) 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{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. In case of dipolar system we calculated macroscopic polarzability using first two methods separately and derived the dielectric constant utilizing equation (\ref{correctionFormula}). Similarly we used equation (\ref{eq:pmf}) to calculate dielectric constant 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 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. 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 evaluate ensemble average of the box dipole or quadrupole moment as a response field or field gradient. The  macroscopic polarizability of the system is derived using ratio between system multipolar moment and applied field or field gradient. This method works properly only 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 i.e. $\braket{A^2} - \braket{A}^2 $ is measured at the fixed temperature and density for a given multipolar fluid. Finally the macroscopic polraizability 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 different constraint distances. Finally we calculated dielectric constant 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 run 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. The short range interaction in the simulations is incorporated with 6-12 Lennard Jones interaction method. 

To derive the box multipolar (dipolar or quadrupolar) moment, we added the component each individual molecule in the space frame and taken ensemble average of the snapshots of the whole 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 quadrupole moment. Similarly the first term in the quadrupolar system is derived using relation $ \braket{Q^2} = \braket{3 Q:Q - TrQ^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 are 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 quadrupolar molecules 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 of the quadrupole moment the components of box quadrupole moments are printed every 100 fs. Each simulations are repeated at different values of applied constant gradients from $ 0 - 9 \times 10^{-2}\; V/\r{A}^2 $. To find 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\;\r{A} - 12\; \r{A} $ at the step of $0.1\; \r{A} $ 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}

\section{Conclusion}

\newpage

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

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

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

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

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

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

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

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

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

\end{document}
Revision:	4389
Committed:	Tue Nov 3 16:30:49 2015 UTC (8 years, 7 months ago) by mlamichh
Content type:	application/x-tex
File size:	61064 byte(s)
Log Message:	added Methodology