ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipole/dielectric_new.tex
Revision: 4363
Committed: Wed Oct 21 14:30:18 2015 UTC (8 years, 8 months ago) by mlamichh
Content type: application/x-tex
File size: 53607 byte(s)
Log Message:
-m "added methodoloy"
 

File Contents

# Content
1 \documentclass[%
2 aip,
3 jcp,
4 amsmath,amssymb,
5 preprint,%
6 % reprint,%
7 %author-year,%
8 %author-numerical,%
9 ]{revtex4-1}
10
11 \usepackage{graphicx}% Include figure files
12 \usepackage{dcolumn}% Align table columns on decimal point
13 \usepackage{multirow}
14 \usepackage{bm}% bold math
15 \usepackage{natbib}
16 \usepackage{times}
17 \usepackage{mathptmx}
18 \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
19 \usepackage{url}
20 \usepackage{braket}
21
22 \begin{document}
23
24 \title{Real space electrostatics for multipoles. III. Dielectric Properties}
25
26 \author{Madan Lamichhane}
27 \affiliation{Department of Physics, University
28 of Notre Dame, Notre Dame, IN 46556}
29 \author{Thomas Parsons}
30 \affiliation{Department of Chemistry and Biochemistry, University
31 of Notre Dame, Notre Dame, IN 46556}
32 \author{Kathie E. Newman}
33 \affiliation{Department of Physics, University
34 of Notre Dame, Notre Dame, IN 46556}
35 \author{J. Daniel Gezelter}
36 \email{gezelter@nd.edu.}
37 \affiliation{Department of Chemistry and Biochemistry, University
38 of Notre Dame, Notre Dame, IN 46556}
39
40 \date{\today}% It is always \today, today,
41 % but any date may be explicitly specified
42
43 \begin{abstract}
44 Note: This manuscript is a work in progress.
45
46 We report on the dielectric properties of the shifted potential
47 (SP), gradient shifted force (GSF), and Taylor shifted force (TSF)
48 real-space methods for multipole interactions that were developed in
49 the first two papers in this series. We find that some subtlety is
50 required for computing dielectric properties with the real-space
51 methods, particularly when using the common fluctuation formulae.
52 Three distinct methods for computing the dielectric constant are
53 investigated, including the standard fluctuation formulae,
54 potentials of mean force between solvated ions, and direct
55 measurement of linear solvent polarization in response to applied
56 fields and field gradients.
57 \end{abstract}
58
59 \maketitle
60
61 \section{Introduction}
62
63 Over the past several years, there has been increasing interest in
64 pairwise methods for correcting electrostatic interactions in computer
65 simulations of condensed molecular
66 systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06}
67 These techniques were initially developed from the observations and efforts of
68 Wolf {\it et al.} and their work towards an $\mathcal{O}(N)$
69 Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is
70 able to obtain excellent agreement with Madelung energies in ionic
71 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
72 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.
73
74 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.
75
76 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}.
77
78 \section{Boltzmann average for orientational polarization}
79 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}.
80
81 \subsection{Dipole}
82 \label{subsec:boltzAverage-Dipole}
83 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}
84
85 \begin{equation}
86 H = H_o - \bf{p_o} .\bf{E},
87 \end{equation}
88 where $H_o$ is a function of the internal coordinates of the molecule. Now Boltzmann average of the dipole moment is given by,
89 \begin{equation}
90 \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}}},
91 \end{equation}
92 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,
93
94 \begin{equation}
95 \braket{p_{mol}} \approx \frac{1}{3}\frac{{p_o}^2}{k_B T}E,
96 \end{equation}
97 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.
98
99
100 \subsection{Quadrupole}
101 \label{subsec:boltzAverage-Quad}
102 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}
103 \begin{equation}
104 \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}}},
105 \label{boltzQuad}
106 \end{equation}
107 where $\int d\Omega = \int_0^{2\pi} \int_0^\pi \int_0^{2\pi}
108 sin\theta\; d\theta\ d\phi\ d\psi$ is the integration over Euler
109 angles, $ H = H_o -q_{\mu\nu}\;\partial_\nu E_\mu $ is the energy of
110 a quadrupole in the gradient of the
111 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,
112 \begin{equation}
113 \begin{split}
114 &q_{\alpha\beta} = \eta_{\alpha\alpha'}\;\eta_{\beta\beta'}\;{q}^* _{\alpha'\beta'} \\
115 &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.
116 \end{split}
117 \label{energyQuad}
118 \end{equation}
119 Here the starred tensors are the components in the body fixed
120 frame. Substituting equation (\ref{energyQuad}) in the equation (\ref{boltzQuad})
121 and taking linear terms in the expansion we get,
122 \begin{equation}
123 \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)},
124 \end{equation}
125 where $\eta_{\alpha\alpha'}$ is the inverse of the rotation matrix that transforms
126 the body fixed co-ordinates to the space co-ordinates,
127 \[\eta_{\alpha\alpha'}
128 = \left(\begin{array}{ccc}
129 cos\phi\; cos\psi - cos\theta\; sin\phi\; sin\psi & -cos\theta\; cos\psi\; sin\phi - cos\phi\; sin\psi & sin\theta\; sin\phi \\
130 cos\psi\; sin\phi + cos\theta\; cos\phi \; sin\psi & cos\theta\; cos\phi\; cos\psi - sin\phi\; sin\psi & -cos\phi\; sin\theta \\
131 sin\theta\; sin\psi & -cos\psi\; sin\theta & cos\theta
132 \end{array} \right).\]
133 Integration of 1st and 2nd terms in the denominator gives $8 \pi^2$
134 and $8 \pi^2 /3\;\vec{\nabla}.\vec{E}\; Tr(q^*) $ respectively. The
135 second term vanishes for charge free space
136 (i.e. $\vec{\nabla}.\vec{E} \; = \; 0)$. Similarly integration of the
137 1st term in the numerator produces
138 $8 \pi^2 /3\; Tr(q^*)\delta_{\alpha\beta}$ and the 2nd term produces
139 $8 \pi^2 /15k_B T (3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'} -
140 {q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'})\partial_\alpha E_\beta$,
141 if $\vec{\nabla}.\vec{E} \; = \; 0$,
142 $ \partial_\alpha E_\beta = \partial_\beta E_\alpha$ and
143 ${q}^*_{\alpha'\beta'}= {q}^*_{\beta'\alpha'}$. Therefore the
144 Boltzmann average of a quadrupole moment can be written as,
145
146 \begin{equation}
147 \braket{q_{\alpha\beta}}\; = \; \frac{1}{3} Tr(q^*)\;\delta_{\alpha\beta} + \frac{{\bar{q_o}}^2}{15k_BT}\;\partial_\alpha E_\beta,
148 \end{equation}
149 where $ \alpha_q = \frac{{\bar{q_o}}^2}{15k_BT} $ is a molecular quadrupolarizablity and ${\bar{q_o}}^2=
150 3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$ is a square of the net quadrupole moment of a molecule.
151
152 \section{Macroscopic Polarizability}
153 \label{sec:MacPolarizablity}
154
155 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.
156
157 \subsection{External field perturbation}
158 \label{subsec:perturbation}
159 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,
160 \begin{equation}
161 \textbf{P} = \epsilon_o \alpha_{D}\; \textbf{E}^o.
162 \label{pertDipole}
163 \end{equation}
164 The constant $\alpha_D$ is a macroscopic polarizability, which is a property of the dipolar fluid in a given density and temperature.
165
166 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,
167 \begin{equation}
168 \begin{split}
169 & {Q}_{\alpha\beta} = \frac{1}{3}\; Tr({Q})\; \delta_{\alpha\beta} + \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta}
170 \\ & or \\
171 & \frac{1}{3}\;\Theta_{\alpha\beta} = \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta}
172 \end{split}
173 \label{pertQuad}
174 \end{equation}
175 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.
176
177
178 \subsection{Fluctuation formula}
179 \label{subsec:fluctuation}
180 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}
181 \begin{equation}
182 \braket{\bf{P}} = \epsilon_o \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E}^o
183 \label{flucDipole}
184 \end{equation}
185 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,
186 \begin{equation}
187 \alpha_D = \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}
188 \end{equation}
189 This is a macroscopic property of dipolar material which is true even if applied field $ \textbf{E}^o \rightarrow 0 $.
190
191 Analogous formula can also be written for a system with quadrupolar molecules,
192 \begin{equation}
193 \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}
194 \label{flucQuad}
195 \end{equation}
196 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,
197 \begin{equation}
198 \alpha_Q = \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T}
199 \label{propConstQuad}
200 \end{equation}
201
202
203 \section{Potential of mean force}
204 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,
205 \begin{equation}
206 w(r) = \int_{\inf}^{r}\braket{\frac{\partial f}{\partial r'}}dr',
207 \end{equation}
208 where $\braket{\partial f/\partial r'}$ is the mean constraint force.
209 Since the ions have a protecting Lennard-Jones (LJ) potential,
210 \begin{equation}
211 w(r) = w_{LJ}(r) + \frac{q_iq_j}{4\pi \epsilon_o \epsilon(r)}U_{method}(r).
212 \label{eq:pmf}
213 \end{equation}
214 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.
215
216 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.
217 \label{sec:PMF}
218
219 \section{Correction factor}
220 \label{sec:corrFactor}
221 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.
222 \subsection{Dipolar system}
223 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}
224 \begin{equation}
225 \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}')}.
226 \end{equation}
227
228 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}).
229 \begin{figure}
230 \includegraphics[width=3in]{DielectricFigure}
231 \caption{With the real-space electrostatic methods, the effective
232 dipole tensor, $\mathbf{T}$, governing interactions between
233 molecular dipoles is not the same for charge-charge interactions as
234 for point dipoles.}
235 \label{fig:stockmayer}
236 \end{figure}
237 In the case where point charges are interacting via an electrostatic
238 kernel, $v(r)$, the effective {\it molecular} dipole tensor,
239 $\mathbf{T}$ is obtained from two successive applications of the
240 gradient operator to the electrostatic kernel,
241 \begin{equation}
242 \mathbf{T}_{\alpha \beta}(r) = \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta}
243 \left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha}
244 r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r}
245 v^{\prime}(r) \right)
246 \label{dipole-chargeTensor}
247 \end{equation}
248 where $v(r)$ may be either the bare kernel ($1/r$) or one of the
249 modified (Wolf or DSF) kernels. This tensor describes the effective
250 interaction between molecular dipoles ($\mathbf{D}$) in Gaussian
251 units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$.
252
253 When utilizing the new real-space methods for point dipoles, the
254 tensor is explicitly constructed,
255 \begin{equation}
256 \mathbf{T}_{\alpha \beta}(r) = \delta_{\alpha \beta} v_{21}(r) +
257 \frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r)
258 \label{dipole-diopleTensor}
259 \end{equation}
260 where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of
261 the approximation. Although the Taylor-shifted (TSF) and
262 gradient-shifted (GSF) models produce to the same $v(r)$ function for
263 point charges, they have distinct forms for the dipole-dipole
264 interactions.
265
266 Using constitutive relation, the polarization density $\textbf{P}(\textbf{r})$ is given by,
267 \begin{equation}
268 \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).
269 \label{constDipole}
270 \end{equation}
271 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}
272 \begin{equation}
273 \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}
274 \label{singular}
275 \end{equation}
276 Substituting equation (\ref{singular}) in the equation (\ref{constDipole}) we get,
277 \begin{equation}
278 \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).
279 \end{equation}
280 For both polarization and electric field homogeneous, this can be easily solved using Fourier transformation,
281 \begin{equation}
282 \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}).
283 \end{equation}
284 For homogeneous applied field Fourier component is non-zero only if $\kappa = 0$. Therefore,
285 \begin{equation}
286 \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}).
287 \label{fourierDipole}
288 \end{equation}
289 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,
290 \begin{equation}
291 \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}
292 \end{equation}
293 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,
294 \begin{equation}
295 \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}
296 \label{correctionFormula}
297 \end{equation}
298 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}.
299 \begin{table}
300 \caption{Expressions for the dipolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter
301 ($\alpha$) and the cutoff radius ($r_c$). The Ewald-Kornfeld result
302 derived in Refs. \onlinecite{Adams:1980rt,Adams:1981fr,Neumann83} is shown for comparison using the Ewald
303 convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). }
304 \label{tab:A}
305 {%
306 \begin{tabular}{l|c|c|c|}
307
308 Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ \\
309 \hline
310 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}$ \\
311 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} $\\
312 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} $ \\
313 Taylor-shifted (TSF) & 1 & 1 \\
314 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
315 \end{tabular}%
316 }
317 \end{table}
318 \subsection{Quadrupolar system}
319 In the presence of the field gradient $\partial_\alpha {E}_\beta $, a
320 non-vanishing quadrupolar polarization (quadrupole moment per unit
321 volume) $\bar{Q}_{\alpha\beta}$ will be induced in the entire volume
322 of a sample. The total field at any point $\vec{r}$ in the sample is
323 given by,
324 \begin{equation}
325 \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'
326 \label{gradMaxwell}
327 \end{equation}
328 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)$.
329
330 \begin{equation}
331 \begin{split}
332 T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \nabla_\gamma \nabla_\delta\;v(r)
333 \\ &= \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)
334 \\ &+ \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)
335 \\ &+ 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),
336 \end{split}
337 \label{quadCharge}
338 \end{equation}
339 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,
340
341 \begin{equation}
342 \begin{split}
343 T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \;v_{\gamma\delta}(r)
344 \\ &= \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}
345 \\ &+ \delta_{\gamma\delta} r_\alpha r_\beta \left(\frac{v''_{21}(r)}{r^2}-\frac{v'_{21}(r)}{r^3} \right)
346 \\ &+\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)
347 \\ &+ 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),
348 \end{split}
349 \label{quadDip}
350 \end{equation}
351 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,
352 \begin{equation}
353 \begin{split}
354 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),
355 \end{split}
356 \label{quadRadial}
357 \end{equation}
358 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.
359 \begin{figure}
360 \includegraphics[width=3in]{QuadrupoleFigure}
361 \caption{With the real-space electrostatic methods, the effective
362 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.}
363 \label{fig:quadrupolarFluid}
364 \end{figure}
365 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
366 field gradient due to quadrupolar fluid vanishes at the singularity (see Appendix \ref{singularQuad}), equation (\ref{gradMaxwell}) can be written as,
367 \begin{equation}
368 \partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) +
369 \frac{1}{4\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 }
370 T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{Q}_{\gamma\delta}(\textbf{r}')\;
371 d^3r'.
372 \end{equation}
373 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}
374 \begin{equation}
375 \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',
376 \end{equation}
377 where $\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$
378 is the traceless quadrupole moment. The total quadrupolar polarization is written as,
379 \begin{equation}
380 {Q}_{\alpha\beta}(\textbf{r}) = \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q})+\epsilon_o {\chi}^*_Q\;\partial_\alpha E_\beta(\textbf{r}),
381 \label{constQaud}
382 \end{equation}
383 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,
384 \begin{equation}
385 \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)
386 \end{equation}
387 For toroidal boundary conditions, both $\partial_\alpha E_\beta$ and
388 ${\Theta}_{\alpha\beta}$ are uniform over the entire space. After
389 performing a Fourier transform (see the Appendix in the Neumann's Paper \cite{Neumann83}) we get,
390 \begin{equation}
391 \frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa})=
392 \epsilon_o {\chi}^*_Q \;\left[{\partial_\alpha
393 {E^o}_\beta}({\kappa})+ \frac{1}{12\pi
394 \epsilon_o}\;{T}_{\alpha\beta\gamma\delta}({\kappa})\;
395 {{\Theta}}_{\gamma\delta}({\kappa})\right]
396 \end{equation}
397 Since the quadrupolar polarization is in the direction of the applied
398 field, we can write
399 ${{\Theta}}_{\gamma\delta}({\kappa}) =
400 {{\Theta}}_{\alpha\beta}({\kappa})$
401 and
402 ${T}_{\alpha\beta\gamma\delta}({\kappa}) =
403 {T}_{\alpha\beta\alpha\beta}({\kappa})$. Therefore we can consider each component of the interaction tensor as scalar and perform calculation.
404 \begin{equation}
405 \begin{split}
406 \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] \\
407 &= \epsilon_o {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\;
408 {T}_{\alpha\beta\alpha\beta}({\kappa})\right)^{-1}
409 {\partial_\alpha E^o_\beta}({\kappa})
410 \end{split}
411 \label{fourierQuad}
412 \end{equation}
413 If the field gradient is homogeneous over the
414 entire volume, ${\partial_ \alpha E_\beta}({\kappa}) = 0 $ except at
415 $ {\kappa} = 0$, hence it is sufficient to know
416 ${T}_{\alpha\beta\alpha\beta}({\kappa})$ at $ {\kappa} =
417 0$. Therefore equation (\ref{fourierQuad}) can be written as,
418 \begin{equation}
419 \begin{split}
420 \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})
421 \end{split}
422 \label{fourierQuad2}
423 \end{equation}
424 where $ {T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as,
425 \begin{equation}
426 {T}_{\alpha\beta\alpha\beta}({0}) = \int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r
427 \label{realTensorQaud}
428 \end{equation}
429
430 In terms of traced quadrupole moment equation (\ref{fourierQuad2}) can be written as,
431 \begin{equation}
432 {{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
433 \label{tracedConstQuad}
434 \end{equation}
435 Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we get,
436 \begin{equation}
437 \begin{split}
438 &\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}, \\
439 &{\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}
440 \end{split}
441 \end{equation}
442 Finally the quadrupolar susceptibility cab be written in terms of quadrupolar correction factor ($A_{quad}$) as below,
443 \begin{equation}
444 {\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}
445 \label{eq:quadrupolarSusceptiblity}
446 \end{equation}
447 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}
448
449 \begin{equation}
450 \epsilon = 1 + \chi^*_Q\; G = 1 + G \; \alpha_Q\left(1 + \alpha_Q\; A_{quad}\right)^{-1}
451 \label{eq:dielectricFromQuadrupoles}
452 \end{equation}
453 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,
454 \begin{equation}
455 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}
456 \label{eq:geometricalFactor}
457 \end{equation}
458 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.
459 \begin{table}
460 \caption{Expressions for the quadrupolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter
461 ($\alpha$) and the cutoff radius ($r_c$). The dimension of the correction factor is $ length^{-2}$ in case of quadrupolar fluid.}
462 \label{tab:B}
463 \centering
464 \resizebox{\columnwidth}{!}{%
465
466 \begin{tabular}{l|c|c|c|c|}
467
468 Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ &$A_\mathrm{quadrupoles}$ \\\hline
469 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}$ \\
470 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}$ \\
471 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)$\\
472 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
473 \end{tabular}%
474 }
475 \end{table}
476 \section{Methodology}
477 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}). We use equation (\ref{eq:pmf}) to calculated 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 derive dielectric constant for two ions dissolved quadrupolar fluid using 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}).
478
479 This research consists of three different types of system for the both dipolar and quadrupolar fluids. First system consist of point dipolar or quadrupolar molecules in the presence of constant electric 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$. The following potential satisfy zero divergence field stated above,
480 \begin{equation}
481 \begin{split}
482 \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. \\
483 & \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),
484 \end{split}
485 \label{eq:appliedPotential}
486 \end{equation}
487 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.
488
489
490 \section{Result}
491
492 \section{Conclusion}
493
494 \newpage
495
496 \appendix
497 \section{Point-multipolar interactions with a spatially-varying electric field}
498
499 We can treat objects $a$, $b$, and $c$ containing embedded collections
500 of charges. When we define the primitive moments, we sum over that
501 collections of charges using a local coordinate system within each
502 object. The point charge, dipole, and quadrupole for object $a$ are
503 given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively.
504 These are the primitive multipoles which can be expressed as a
505 distribution of charges,
506 \begin{align}
507 C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\
508 D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\
509 Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in } a} q_k
510 r_{k\alpha} r_{k\beta} . \label{eq:quadrupole}
511 \end{align}
512 Note that the definition of the primitive quadrupole here differs from
513 the standard traceless form, and contains an additional Taylor-series
514 based factor of $1/2$. In Paper 1, we derived the forces and torques
515 each object exerts on the others.
516
517 Here we must also consider an external electric field that varies in
518 space: $\mathbf E(\mathbf r)$. Each of the local charges $q_k$ in
519 object $a$ will then experience a slightly different field. This
520 electric field can be expanded in a Taylor series around the local
521 origin of each object. A different Taylor series expansion is carried
522 out for each object.
523
524 For a particular charge $q_k$, the electric field at that site's
525 position is given by:
526 \begin{equation}
527 E_\gamma + \nabla_\delta E_\gamma r_{k \delta}
528 + \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
529 r_{k \varepsilon} + ...
530 \end{equation}
531 Note that the electric field is always evaluated at the origin of the
532 objects, and treating each object using point multipoles simplifies
533 this greatly.
534
535 To find the force exerted on object $a$ by the electric field, one
536 takes the electric field expression, and multiplies it by $q_k$, and
537 then sum over all charges in $a$:
538
539 \begin{align}
540 F_\gamma &= \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta}
541 + \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
542 r_{k \varepsilon} + ... \rbrace \\
543 &= C_a E_\gamma + D_{a \delta} \nabla_\delta E_\gamma
544 + Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma +
545 ...
546 \end{align}
547
548 Similarly, the torque exerted by the field on $a$ can be expressed as
549 \begin{align}
550 \tau_\alpha &= \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\
551 & = \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k
552 r_{k\beta} E_\gamma(\mathbf r_k) \\
553 & = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma
554 + 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta
555 E_\gamma + ...
556 \end{align}
557
558 The last term is essentially identical with form derived by Torres del
559 Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation
560 utilized a traceless form of the quadrupole that is different than the
561 primitive definition in use here. We note that the Levi-Civita symbol
562 can be eliminated by utilizing the matrix cross product in an
563 identical form as in Ref. \onlinecite{Smith98}:
564 \begin{equation}
565 \left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta
566 \left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta}
567 -\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta}
568 \right]
569 \label{eq:matrixCross}
570 \end{equation}
571 where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of
572 the matrix indices. In table \ref{tab:UFT} we give compact
573 expressions for how the multipole sites interact with an external
574 field that has exhibits spatial variations.
575
576 \begin{table}
577 \caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
578 $(\mathbf{\tau})$ expressions for a multipolar site embedded in an
579 electric field with spatial variations, $\mathbf{E}(\mathbf{r})$.
580 \label{tab:UFT}}
581 \begin{tabular}{r|ccc}
582 & Charge & Dipole & Quadrupole \\ \hline
583 $U$ & $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\
584 $\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ & $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\
585 $\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$
586 \end{tabular}
587 \end{table}
588 \section{Gradient of the field due to quadrupolar polarization}
589 \label{singularQuad}
590 In this section, we will discuss the gradient of the field produced by
591 quadrupolar polarization. For this purpose, we consider a distribution
592 of charge ${\rho}(r)$ which gives rise to an electric field
593 $\vec{E}(r)$ and gradient of the field $\vec{\nabla} \vec{E}(r)$
594 throughout space. The total gradient of the electric field over volume
595 due to the all charges within the sphere of radius $R$ is given by
596 (cf. Jackson equation 4.14):
597 \begin{equation}
598 \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = -\int_{r=R} R^2 \vec{E}\;\hat{n}\; d\Omega
599 \label{eq:8}
600 \end{equation}
601 where $d\Omega$ is the solid angle and $\hat{n}$ is the normal vector
602 of the surface of the sphere which is equal to
603 $sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} +
604 cos[\theta]\hat{z}$
605 in spherical coordinates. For the charge density ${\rho}(r')$, the
606 total gradient of the electric field can be written as (cf. Jackson
607 equation 4.16),
608 \begin{equation}
609 \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
610 \label{eq:9}
611 \end{equation}
612 The radial function in the equation (\ref{eq:9}) can be expressed in
613 terms of spherical harmonics as (cf. Jackson equation 3.70),
614 \begin{equation}
615 \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)
616 \label{eq:10}
617 \end{equation}
618 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,
619 \begin{equation}
620 \begin{split}
621 \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 \\
622 &= -\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
623 '
624 \end{split}
625 \label{eq:11}
626 \end{equation}
627 The gradient of the product of radial function and spherical harmonics
628 is given by (cf. Arfken, p.811 eq. 16.94):
629 \begin{equation}
630 \begin{split}
631 \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
632 {\partial}{\partial r}+\frac{l}{r} \right]f(r)\; Y_{l, l-1, m}(\theta, \phi).
633 \end{split}
634 \label{eq:12}
635 \end{equation}
636 Using equation (\ref{eq:12}) we get,
637 \begin{equation}
638 \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}},
639 \label{eq:13}
640 \end{equation}
641 where $ Y_{l,l+1,m}(\theta, \phi)$ is the vector spherical harmonics
642 which can be expressed in terms of spherical harmonics as shown in
643 below (cf. Arfkan p.811),
644 \begin{equation}
645 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},
646 \label{eq:14}
647 \end{equation}
648 where $C(l+1,1,l|m_1,m_2,m)$ is a Clebsch-Gordan coefficient and
649 $\hat{e}_{m_2}$ is a spherical tensor of rank 1 which can be expressed
650 in terms of Cartesian coordinates,
651 \begin{equation}
652 {\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}}
653 \label{eq:15}
654 \end{equation}
655 The normal vector $\hat{n} $ can be expressed in terms of spherical tensor of rank 1 as shown in below,
656 \begin{equation}
657 \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)
658 \label{eq:16}
659 \end{equation}
660 The surface integral of the product of $\hat{n}$ and
661 ${Y_{l+1}}^{m_1}(\theta, \phi)$ gives,
662 \begin{equation}
663 \begin{split}
664 \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 \\
665 &= \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 \\
666 &= \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),
667 \end{split}
668 \label{eq:17}
669 \end{equation}
670 where ${Y_{l}}^{-m} = (-1)^m\;{{Y_{l}}^{m}}^* $ and
671 $ \int {{Y_{l}}^{m}}^*\;{Y_{l'}}^{m'}\;d\Omega =
672 \delta_{ll'}\delta_{mm'} $.
673 Non-vanishing values of equation \ref{eq:17} require $l = 0$,
674 therefore the value of $ m = 0 $. Since the values of $ m_1$ are -1,
675 1, and 0 then $m_2$ takes the values 1, -1, and 0, respectively
676 provided that $m = m_1 + m_2$. Equation \ref{eq:11} can therefore be
677 modified,
678 \begin{equation}
679 \begin{split}
680 \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(
681 1, 1, 0|0,0,0)\;{\hat{e}_{0}}{\hat{e}_{0}} ]\; d^3r'.
682 \end{split}
683 \label{eq:18}
684 \end{equation}
685 After substituting ${Y^*}_{00} = \frac{1}{\sqrt{4\pi}} $ and using the
686 values of the Clebsch-Gorden coefficients: $ C(1, 1, 0|-1,1,0) =
687 \frac{1}{\sqrt{3}}, \; C(1, 1, 0|-1,1,0)= \frac{1}{\sqrt{3}}$ and $
688 C(1, 1, 0|0,0,0) = -\frac{1}{\sqrt{3}}$ in equation \ref{eq:18} we
689 obtain,
690 \begin{equation}
691 \begin{split}
692 \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)\\
693 &= - \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).
694 \end{split}
695 \label{eq:19}
696 \end{equation}
697 Equation (\ref{eq:19}) gives the total gradient of the field over a
698 sphere due to the distribution of the charges. For quadrupolar fluids
699 the total charge within a sphere is zero, therefore
700 $ \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = 0 $. Hence the quadrupolar
701 polarization produces zero net gradient of the field inside the
702 sphere.
703
704
705
706 \newpage
707
708 \bibliography{multipole}
709
710 \end{document}