1 |
\documentclass[aip,jcp,amsmath,amssymb,preprint]{revtex4-1} |
2 |
|
3 |
\usepackage{graphicx}% Include figure files |
4 |
\usepackage{dcolumn}% Align table columns on decimal point |
5 |
\usepackage{multirow} |
6 |
\usepackage{bm} |
7 |
\usepackage{natbib} |
8 |
\usepackage{times} |
9 |
\usepackage{mathptmx} |
10 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
11 |
\usepackage{url} |
12 |
\usepackage{braket} |
13 |
\usepackage{tabularx} |
14 |
|
15 |
\newcolumntype{Y}{>{\centering\arraybackslash}X} |
16 |
|
17 |
\begin{document} |
18 |
|
19 |
\title{Real space electrostatics for multipoles. III. Dielectric Properties} |
20 |
|
21 |
\author{Madan Lamichhane} |
22 |
\affiliation{Department of Physics, University |
23 |
of Notre Dame, Notre Dame, IN 46556} |
24 |
\author{Thomas Parsons} |
25 |
\affiliation{Department of Chemistry and Biochemistry, University |
26 |
of Notre Dame, Notre Dame, IN 46556} |
27 |
\author{Kathie E. Newman} |
28 |
\affiliation{Department of Physics, University |
29 |
of Notre Dame, Notre Dame, IN 46556} |
30 |
\author{J. Daniel Gezelter} |
31 |
\email{gezelter@nd.edu.} |
32 |
\affiliation{Department of Chemistry and Biochemistry, University |
33 |
of Notre Dame, Notre Dame, IN 46556} |
34 |
|
35 |
\date{\today}% It is always \today, today, |
36 |
% but any date may be explicitly specified |
37 |
|
38 |
\begin{abstract} |
39 |
In the first two papers in this series, we developed new shifted |
40 |
potential (SP), gradient shifted force (GSF), and Taylor shifted |
41 |
force (TSF) real-space methods for multipole interactions in |
42 |
condensed phase simulations. Here, we discuss the dielectric |
43 |
properties of fluids that emerge from simulations using these |
44 |
methods. Most electrostatic methods (including the Ewald sum) |
45 |
require correction to the conducting boundary fluctuation formula |
46 |
for the static dielectric constants, and we discuss the derivation |
47 |
of these corrections for the new real space methods. For quadrupolar |
48 |
fluids, the analogous material property is the quadrupolar |
49 |
susceptibility. As in the dipolar case, the fluctuation formula for |
50 |
the quadrupolar susceptibility has corrections that depend on the |
51 |
electrostatic method being utillized. One of the most important |
52 |
effects measured by both the static dielectric and quadrupolar |
53 |
susceptibility is the ability to screen charges embedded in the |
54 |
fluid. We use potentials of mean force between solvated ions to |
55 |
discuss how geometric factors can lead to distance-dependent |
56 |
screening in both quadrupolar and dipolar fluids. |
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 or ``real space'' methods for computing electrostatic |
65 |
interactions in condensed phase |
66 |
simulations.\cite{Wolf99,Zahn02,Kast03,Beckd05,Ma05,Fennell06,ZMM papers} |
67 |
These techniques were initially developed by Wolf {\it et al.} in |
68 |
their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99} |
69 |
Wolf's method of using cutoff neutralization and electrostatic damping |
70 |
is able to obtain excellent agreement with Madelung energies in ionic |
71 |
crystals.\cite{Wolf99} |
72 |
|
73 |
Zahn \textit{et al.}\cite{Zahn02} and Fennell and Gezelter extended |
74 |
this method using shifted force approximations at the cutoff distance |
75 |
in order to conserve total energy in molecular dynamics |
76 |
simulations.\cite{Fennell06} Other recent advances in real-space |
77 |
methods for systems of point charges have included explicit zeroing |
78 |
out of the net multipole moments inside the cutoff spheres.\cite{ZMM |
79 |
papers} |
80 |
|
81 |
In the previous two papers in this series, we developed three |
82 |
generalized real space methods: shifted potential (SP), gradient |
83 |
shifted force (GSF), and Taylor shifted force (TSF).\cite{PaperI, |
84 |
PaperII} These methods evaluate electrostatic interactions for |
85 |
charges and higher order multipoles using a finite-radius cutoff |
86 |
sphere. The neutralization and damping of local moments within the |
87 |
cutoff sphere is a multipolar generalization of Wolf's sum. In the |
88 |
GSF and TSF methods, additional terms are added to the potential |
89 |
energy so that forces and torques also vanish smoothly at the cutoff |
90 |
radius. This ensures that the total energy is conserved in molecular |
91 |
dynamics simulations. |
92 |
|
93 |
One of the most difficult tests of any new electrostatic method is the |
94 |
fidelity with which that method can reproduce the bulk-phase |
95 |
polarizability or equivalently, the dielectric properties of a |
96 |
fluid. Before the advent of computer simulations, Kirkwood and Onsager |
97 |
developed fluctuation formulae for the dielectric properties of |
98 |
dipolar fluids.\cite{Kirkwood39,Onsagar36} Along with projections of |
99 |
the frequency-dependent dielectric to zero frequency, these |
100 |
fluctuation formulae are now widely used to predict dielectric |
101 |
response of simulated materials. |
102 |
|
103 |
If we consider a system of dipolar or quadrupolar molecules under the |
104 |
influence of an external field or field gradient, the net polarization |
105 |
of the system will largely be proportional to the applied |
106 |
perturbation.\cite{Chitanvis96, Stern-Feller03, SalvchovI14, |
107 |
SalvchovII14} In simulations, the net polarization of the system is |
108 |
also determined by the interactions \textit{between} the |
109 |
molecules. Therefore the macroscopic polarizablity obtained from the |
110 |
simulation depends on the details of the electrostatic interaction |
111 |
methods that were utilized in the simulation. To determine the |
112 |
relevant physical properties of the multipolar fluid from the system |
113 |
fluctuations, the interactions between molecules must be incorporated |
114 |
into the formalism for the bulk properties. |
115 |
|
116 |
In modern simulations, bulk materials are usually treated using |
117 |
periodic replicas of small regions, and this level of approximation |
118 |
requires corrections to the fluctuation formulae that were derived for |
119 |
the bulk fluids. In 1983 Neumann proposed a general formula for |
120 |
evaluating dielectric properties of dipolar fluids using real-space |
121 |
cutoff methods.\cite{NeumannI83} Steinhauser and Neumann used this |
122 |
formula to evaluate the corrected dielectric constant for the |
123 |
Stockmayer fluid using two different methods: Ewald-Kornfield (EK) and |
124 |
reaction field (RF) methods.\cite{NeumannII83} |
125 |
|
126 |
Zahn \textit{et al.}\cite{Zahn02} utilized this approach and evaluated |
127 |
the correction factor for using damped shifted charge-charge |
128 |
kernel. This was later generalized by Izvekov \textit{et |
129 |
al.},\cite{Izvekov08} who that the expression for the |
130 |
dielectric constant reduces to widely-used \textit{conducting |
131 |
boundary} formula for real-space cutoff methods that have first |
132 |
derivatives that vanish at the cutoff sphere. |
133 |
|
134 |
One of the primary topics of this paper will be derivation of these |
135 |
correction factors for the new real space methods. We find that the |
136 |
correction formulae for dipolar molecules depends not only on the |
137 |
methodology being used, but also on whether the molecular dipoles are |
138 |
treated using point charges or point dipoles. We derive correction |
139 |
factors for both cases. |
140 |
|
141 |
In quadrupolar fluids, the relationship between quadrupolar |
142 |
susceptibility and the dielectric constant is not as straightforward |
143 |
as in the dipolar case. The effective dielectric constant depends on |
144 |
the geometry of the external (or internal) field |
145 |
perturbation.\cite{Ernst92} Significant efforts have been made to |
146 |
increase our understanding the dielectric properties of these |
147 |
fluids,\cite{JeonI03,JeonII03,Chitanvis96} although a correction |
148 |
formula for different real-space methods has not yet been developed. |
149 |
|
150 |
In this paper we derive general formulae for calculating the |
151 |
dielectric properties of quadrupolar fluids. We also evaluate the |
152 |
correction factor for SP, GSF, and TSF methods for both dipolar and |
153 |
quadrupolar fluids interacting via point charges, point dipoles or |
154 |
through quadrupole-quadrupole interactions. |
155 |
|
156 |
We have also calculated the screening behavior for two ions immersed |
157 |
in multipolar fluids to estimate the distance dependence of dielectric |
158 |
screening in both dipolar and quadrupolar fluids. We have used three |
159 |
different methods to compare our results with computer simulations: |
160 |
\begin{enumerate} |
161 |
\item responses of the fluid to external perturbations, |
162 |
\item fluctuations of system multipole moments, and |
163 |
\item potentials of mean force between solvated ions, |
164 |
\end{enumerate} |
165 |
|
166 |
Using external perturbations, the net polarization of the system is |
167 |
observed as a linear response of the applied field perturbation, where |
168 |
proportionality constant is determined by the electrostatic |
169 |
interaction between the electrostatic multipoles at a given |
170 |
temperature. The fluctuation formula observes the time average |
171 |
fluctuation of the multipolar moment as a function of temperature. The |
172 |
average fluctuation value of the system is determined by the |
173 |
multipole-multipole interactions between molecules at a given |
174 |
temperature. Since the expression of the electrostatic interaction |
175 |
energy, force, and torque in the real space electrostatic methods are |
176 |
different from their original definition, both fluctuation and |
177 |
external field perturbation formula must also be modified |
178 |
accordingly. |
179 |
|
180 |
The potential of mean force allows calculation of an effective |
181 |
dielectric constant or screening factor from the potential energy |
182 |
between ions before and after dielectric material is introduced. |
183 |
|
184 |
\section{Dipolar Fluids and the Dielectric Constant} |
185 |
|
186 |
Dielectric properties of a fluid arise mainly from responses of the |
187 |
fluid to either applied fields or transient fields internal to the |
188 |
fluid. In response to an applied field, the molecules have electronic |
189 |
polarizabilities, changes to internal bond lengths, and reorientations |
190 |
towards the direction of the applied field. There is an added |
191 |
complication that in the presence of external field, the perturbation |
192 |
experienced by any single molecule is not only due to external field |
193 |
but also to the fields produced by the all other molecules in the |
194 |
system. |
195 |
|
196 |
\subsection{Response to External Perturbations} |
197 |
|
198 |
In the presence of uniform electric field $\mathbf{E}$, an individual |
199 |
molecule with a permanent dipole moment $p_o$ will realign along the |
200 |
direction of the field with an average polarization given by |
201 |
\begin{equation} |
202 |
\braket{p_{mol}} = \epsilon_0 \alpha_p \mathbf{E}, |
203 |
\end{equation} |
204 |
where $ \alpha_p = \frac{1}{3}\frac{{p_o}^2}{\epsilon_0 k_B T}$ is the |
205 |
contribution to molecular polarizability due solely to reorientation |
206 |
dynamics. The orientational polarization depends inversely on the |
207 |
temperature as the applied field must overcome the thermal motion. |
208 |
|
209 |
Likewise, a condensed phase system of dipolar molecules will also |
210 |
polarize along the direction of an applied field. The net dipolar |
211 |
polarization $\mathbf{P}$ of the system is, |
212 |
\begin{equation} |
213 |
\textbf{P} = \epsilon_o \alpha_{D} \mathbf{E}. |
214 |
\label{pertDipole} |
215 |
\end{equation} |
216 |
The constant $\alpha_D$ is a macroscopic polarizability, which is an |
217 |
emergent property of the dipolar fluid in a given density and |
218 |
temperature and is related to the static dielectric constant. |
219 |
|
220 |
\subsection{Fluctuation Formula} |
221 |
|
222 |
Of particular interest is the static dielectric constant, $\epsilon$. |
223 |
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can |
224 |
be calculated for systems of orientationally-polarized molecules via |
225 |
\begin{equation} |
226 |
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
227 |
\label{eq:staticDielectric} |
228 |
\end{equation} |
229 |
where $\epsilon_0$ is the permittivity of free space and |
230 |
$\langle M^2\rangle$ is the fluctuation of the system dipole |
231 |
moment.\cite{Allen89} The numerator in the fractional term in |
232 |
equation (\ref{eq:staticDielectric}) is identical to the quantity |
233 |
calculated in the finite-system Kirkwood $g$ factor ($G_k$): |
234 |
\begin{equation} |
235 |
G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
236 |
\label{eq:KirkwoodFactor} |
237 |
\end{equation} |
238 |
where $\mu$ is the dipole moment of a single molecule of the |
239 |
homogeneous system.\cite{NeumannI83,NeumannII83,Neumann84,Neumann85} The |
240 |
fluctuation term in both equation (\ref{eq:staticDielectric}) and |
241 |
(\ref{eq:KirkwoodFactor}) is calculated as follows, |
242 |
\begin{equation} |
243 |
\begin{split} |
244 |
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
245 |
- \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
246 |
&= \langle M_x^2+M_y^2+M_z^2\rangle |
247 |
- (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
248 |
+ \langle M_x\rangle^2). |
249 |
\end{split} |
250 |
\label{eq:fluctBoxDipole} |
251 |
\end{equation} |
252 |
This fluctuation term can be accumulated during a simulation; however, |
253 |
it converges rather slowly, thus requiring multi-nanosecond simulation |
254 |
times.\cite{Horn04} In the case of tin-foil boundary conditions, the |
255 |
dielectric/surface term of the Ewald summation is equal to zero. |
256 |
|
257 |
|
258 |
|
259 |
For a system of molecules, the net dipole moment, $\bf{M}$ at thermal |
260 |
equilibrium of temperature T in the presence of applied field |
261 |
$\bf{E}$, the average dipolar polarization can be expressed in terms |
262 |
of fluctuation of the net dipole moment, |
263 |
\begin{equation} |
264 |
\braket{\bf{P}} = \epsilon_o \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E} |
265 |
\label{flucDipole} |
266 |
\end{equation} |
267 |
This is similar to the formula for Boltzmann average of single dipolar |
268 |
molecule. Here $\braket{\bf{P}}$ is average polarization and |
269 |
$ \braket{\textbf{M}^2}-{\braket{\textbf{M}}}^2$ is the net dipole |
270 |
fluctuation at temperature T. For the limiting case |
271 |
$\textbf{E} \rightarrow 0 $, the ensemble average of both the net |
272 |
dipole moment $\braket{\textbf{M}}$ and dipolar polarization |
273 |
$\braket{\bf{P}}$ tends to vanish but $\braket{\bf{M}^2}$ will still |
274 |
be non-zero. The dipolar macroscopic polarizability can be written as, |
275 |
\begin{equation} |
276 |
\alpha_D = \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T} |
277 |
\end{equation} |
278 |
This is a macroscopic property of dipolar material which is true even |
279 |
if applied field $ \textbf{E} \rightarrow 0 $. |
280 |
\subsection{Correction Factors} |
281 |
In the presence of an external field $ \textbf{E}$ polarization $\textbf{E}$ will be induced in a dipolar system. The total electrostatic field (or Maxwell electric field) at point $\bf{r}$ in a system is,\cite{NeumannI83} |
282 |
\begin{equation} |
283 |
\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}')}. |
284 |
\end{equation} |
285 |
|
286 |
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}). |
287 |
\begin{figure} |
288 |
\includegraphics[width=3in]{DielectricFigure} |
289 |
\caption{With the real-space electrostatic methods, the effective |
290 |
dipole tensor, $\mathbf{T}$, governing interactions between |
291 |
molecular dipoles is not the same for charge-charge interactions as |
292 |
for point dipoles.} |
293 |
\label{fig:stockmayer} |
294 |
\end{figure} |
295 |
In the case where point charges are interacting via an electrostatic |
296 |
kernel, $v(r)$, the effective {\it molecular} dipole tensor, |
297 |
$\mathbf{T}$ is obtained from two successive applications of the |
298 |
gradient operator to the electrostatic kernel, |
299 |
\begin{equation} |
300 |
\mathbf{T}_{\alpha \beta}(r) = \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta} |
301 |
\left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha} |
302 |
r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r} |
303 |
v^{\prime}(r) \right) |
304 |
\label{dipole-chargeTensor} |
305 |
\end{equation} |
306 |
where $v(r)$ may be either the bare kernel ($1/r$) or one of the |
307 |
modified (Wolf or DSF) kernels. This tensor describes the effective |
308 |
interaction between molecular dipoles ($\mathbf{D}$) in Gaussian |
309 |
units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$. |
310 |
|
311 |
When utilizing the new real-space methods for point dipoles, the |
312 |
tensor is explicitly constructed, |
313 |
\begin{equation} |
314 |
\mathbf{T}_{\alpha \beta}(r) = \delta_{\alpha \beta} v_{21}(r) + |
315 |
\frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r) |
316 |
\label{dipole-diopleTensor} |
317 |
\end{equation} |
318 |
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of |
319 |
the approximation. Although the Taylor-shifted (TSF) and |
320 |
gradient-shifted (GSF) models produce to the same $v(r)$ function for |
321 |
point charges, they have distinct forms for the dipole-dipole |
322 |
interactions. |
323 |
|
324 |
Using constitutive relation, the polarization density $\textbf{P}(\textbf{r})$ is given by, |
325 |
\begin{equation} |
326 |
\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). |
327 |
\label{constDipole} |
328 |
\end{equation} |
329 |
Here $\chi^*_D$ is a dipolar susceptibility can be expressed in terms of dielectric constant as $ \chi^*_D = \epsilon - 1$ which different than macroscopic dipolar polarizability $\alpha_D$ in the sections \ref{subsec:perturbation} and \ref{subsec:fluctuation}. We can split integral into two parts: singular part i.e $|\textbf{r}-\textbf{r}'|\rightarrow 0 $ and non-singular part i.e $|\textbf{r}-\textbf{r}'| > 0 $ . The singular part of the integral can be written as,\cite{NeumannI83, Jackson98} |
330 |
\begin{equation} |
331 |
\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} |
332 |
\label{singular} |
333 |
\end{equation} |
334 |
Substituting equation (\ref{singular}) in the equation (\ref{constDipole}) we get, |
335 |
\begin{equation} |
336 |
\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). |
337 |
\end{equation} |
338 |
For both polarization and electric field homogeneous, this can be easily solved using Fourier transformation, |
339 |
\begin{equation} |
340 |
\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}). |
341 |
\end{equation} |
342 |
For homogeneous applied field Fourier component is non-zero only if $\kappa = 0$. Therefore, |
343 |
\begin{equation} |
344 |
\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}). |
345 |
\label{fourierDipole} |
346 |
\end{equation} |
347 |
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, |
348 |
\begin{equation} |
349 |
\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} |
350 |
\end{equation} |
351 |
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, |
352 |
\begin{equation} |
353 |
\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} |
354 |
\label{correctionFormula} |
355 |
\end{equation} |
356 |
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}. |
357 |
\begin{table} |
358 |
\caption{Expressions for the dipolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter |
359 |
($\alpha$) and the cutoff radius ($r_c$). The Ewald-Kornfeld result |
360 |
derived in Refs. \onlinecite{Adams80,Adams81,NeumannI83} is shown for comparison using the Ewald |
361 |
convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). } |
362 |
\label{tab:A} |
363 |
{% |
364 |
\begin{tabular}{l|c|c|c|} |
365 |
|
366 |
Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ \\ |
367 |
\hline |
368 |
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}$ \\ |
369 |
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} $\\ |
370 |
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} $ \\ |
371 |
Taylor-shifted (TSF) & 1 & 1 \\ |
372 |
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 |
373 |
\end{tabular}% |
374 |
} |
375 |
\end{table} |
376 |
|
377 |
|
378 |
\section{Quadrupolar Fluids and the Quadrupolar Susceptibility} |
379 |
\subsection{Response to External Perturbations} |
380 |
|
381 |
In analogy to the dipolar case, a molecule with a permanent |
382 |
quadrupole, $q_{\alpha\beta}$ will also align in the presence of a |
383 |
uniform electric field gradient $\vec{\nabla}\mathbf{E}$. The average |
384 |
polarization is given by |
385 |
\begin{equation} |
386 |
\braket{q_{\alpha\beta}} = \frac{1}{3} Tr(q^*) \delta_{\alpha\beta} + \alpha_q \partial_\alpha E_\beta, |
387 |
\end{equation} |
388 |
where $ \alpha_q = \frac{{\bar{q_o}}^2}{15k_BT} $ is a molecular |
389 |
quadrupolarizablity and |
390 |
${\bar{q_o}}^2= |
391 |
3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$ |
392 |
is a square of the net quadrupole moment of a molecule. |
393 |
|
394 |
Similarly, in the presence of and external field gradient, a system of |
395 |
quadrupolar molecule polarizes along the direction of applied field |
396 |
gradient, giving a system quadrupolar polarization, |
397 |
be given by, |
398 |
\begin{equation} |
399 |
\begin{split} |
400 |
& {Q}_{\alpha\beta} = \frac{1}{3} Tr({Q}) \delta_{\alpha\beta} + \epsilon_o \alpha_Q \partial_{\alpha} E^o_{\beta} |
401 |
\\ & or \\ |
402 |
& \frac{1}{3}\;\Theta_{\alpha\beta} = \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta} |
403 |
\end{split} |
404 |
\label{pertQuad} |
405 |
\end{equation} |
406 |
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. |
407 |
\subsection{Fluctuation Formula} |
408 |
|
409 |
Fluctuation formulae for quadrupolar fluids were developed by Logan |
410 |
\textit{et al.}\cite{LoganI81,LoganII82,LoganIII82} |
411 |
|
412 |
In analogy with the dipolar case, a fluctuation formula can be |
413 |
written for a system comprising quadrupolar molecules, |
414 |
\begin{equation} |
415 |
\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} |
416 |
\label{flucQuad} |
417 |
\end{equation} |
418 |
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, |
419 |
\begin{equation} |
420 |
\alpha_Q = \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T} |
421 |
\label{propConstQuad} |
422 |
\end{equation} |
423 |
|
424 |
\subsection{Correction Factors} |
425 |
In the presence of the field gradient $\partial_\alpha {E}_\beta $, a |
426 |
non-vanishing quadrupolar polarization (quadrupole moment per unit |
427 |
volume) $\bar{Q}_{\alpha\beta}$ will be induced in the entire volume |
428 |
of a sample. The total field at any point $\vec{r}$ in the sample is |
429 |
given by, |
430 |
\begin{equation} |
431 |
\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' |
432 |
\label{gradMaxwell} |
433 |
\end{equation} |
434 |
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)$. |
435 |
|
436 |
\begin{equation} |
437 |
\begin{split} |
438 |
T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \nabla_\gamma \nabla_\delta\;v(r) |
439 |
\\ &= \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) |
440 |
\\ &+ \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) |
441 |
\\ &+ 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), |
442 |
\end{split} |
443 |
\label{quadCharge} |
444 |
\end{equation} |
445 |
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, |
446 |
|
447 |
\begin{equation} |
448 |
\begin{split} |
449 |
T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \;v_{\gamma\delta}(r) |
450 |
\\ &= \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} |
451 |
\\ &+ \delta_{\gamma\delta} r_\alpha r_\beta \left(\frac{v''_{21}(r)}{r^2}-\frac{v'_{21}(r)}{r^3} \right) |
452 |
\\ &+\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) |
453 |
\\ &+ 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), |
454 |
\end{split} |
455 |
\label{quadDip} |
456 |
\end{equation} |
457 |
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, |
458 |
\begin{equation} |
459 |
\begin{split} |
460 |
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), |
461 |
\end{split} |
462 |
\label{quadRadial} |
463 |
\end{equation} |
464 |
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. |
465 |
\begin{figure} |
466 |
\includegraphics[width=3in]{QuadrupoleFigure} |
467 |
\caption{With the real-space electrostatic methods, the effective |
468 |
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.} |
469 |
\label{fig:quadrupolarFluid} |
470 |
\end{figure} |
471 |
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 |
472 |
field gradient due to quadrupolar fluid vanishes at the singularity (see Appendix \ref{singularQuad}), equation (\ref{gradMaxwell}) can be written as, |
473 |
\begin{equation} |
474 |
\partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) + |
475 |
\frac{1}{4\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 } |
476 |
T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{Q}_{\gamma\delta}(\textbf{r}')\; |
477 |
d^3r'. |
478 |
\end{equation} |
479 |
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} |
480 |
\begin{equation} |
481 |
\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', |
482 |
\end{equation} |
483 |
where $\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$ |
484 |
is the traceless quadrupole moment. The total quadrupolar polarization is written as, |
485 |
\begin{equation} |
486 |
{Q}_{\alpha\beta}(\textbf{r}) = \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q})+\epsilon_o {\chi}^*_Q\;\partial_\alpha E_\beta(\textbf{r}), |
487 |
\label{constQaud} |
488 |
\end{equation} |
489 |
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, |
490 |
\begin{equation} |
491 |
\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) |
492 |
\end{equation} |
493 |
For toroidal boundary conditions, both $\partial_\alpha E_\beta$ and |
494 |
${\Theta}_{\alpha\beta}$ are uniform over the entire space. After |
495 |
performing a Fourier transform (see the Appendix in the Neumann's Paper \cite{NeumannI83}) we get, |
496 |
\begin{equation} |
497 |
\frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa})= |
498 |
\epsilon_o {\chi}^*_Q \;\left[{\partial_\alpha |
499 |
{E^o}_\beta}({\kappa})+ \frac{1}{12\pi |
500 |
\epsilon_o}\;{T}_{\alpha\beta\gamma\delta}({\kappa})\; |
501 |
{{\Theta}}_{\gamma\delta}({\kappa})\right] |
502 |
\end{equation} |
503 |
Since the quadrupolar polarization is in the direction of the applied |
504 |
field, we can write |
505 |
${{\Theta}}_{\gamma\delta}({\kappa}) = |
506 |
{{\Theta}}_{\alpha\beta}({\kappa})$ |
507 |
and |
508 |
${T}_{\alpha\beta\gamma\delta}({\kappa}) = |
509 |
{T}_{\alpha\beta\alpha\beta}({\kappa})$. Therefore we can consider each component of the interaction tensor as scalar and perform calculation. |
510 |
\begin{equation} |
511 |
\begin{split} |
512 |
\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] \\ |
513 |
&= \epsilon_o {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\; |
514 |
{T}_{\alpha\beta\alpha\beta}({\kappa})\right)^{-1} |
515 |
{\partial_\alpha E^o_\beta}({\kappa}) |
516 |
\end{split} |
517 |
\label{fourierQuad} |
518 |
\end{equation} |
519 |
If the field gradient is homogeneous over the |
520 |
entire volume, ${\partial_ \alpha E_\beta}({\kappa}) = 0 $ except at |
521 |
$ {\kappa} = 0$, hence it is sufficient to know |
522 |
${T}_{\alpha\beta\alpha\beta}({\kappa})$ at $ {\kappa} = |
523 |
0$. Therefore equation (\ref{fourierQuad}) can be written as, |
524 |
\begin{equation} |
525 |
\begin{split} |
526 |
\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}) |
527 |
\end{split} |
528 |
\label{fourierQuad2} |
529 |
\end{equation} |
530 |
where $ {T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as, |
531 |
\begin{equation} |
532 |
{T}_{\alpha\beta\alpha\beta}({0}) = \int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r |
533 |
\label{realTensorQaud} |
534 |
\end{equation} |
535 |
|
536 |
In terms of traced quadrupole moment equation (\ref{fourierQuad2}) can be written as, |
537 |
\begin{equation} |
538 |
{{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 |
539 |
\label{tracedConstQuad} |
540 |
\end{equation} |
541 |
Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we get, |
542 |
\begin{equation} |
543 |
\begin{split} |
544 |
&\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}, \\ |
545 |
&{\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} |
546 |
\end{split} |
547 |
\end{equation} |
548 |
Finally the quadrupolar susceptibility cab be written in terms of quadrupolar correction factor ($A_{quad}$) as below, |
549 |
\begin{equation} |
550 |
{\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} |
551 |
\label{eq:quadrupolarSusceptiblity} |
552 |
\end{equation} |
553 |
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} |
554 |
|
555 |
\begin{equation} |
556 |
\epsilon = 1 + \chi^*_Q\; G = 1 + G \; \alpha_Q\left(1 + \alpha_Q\; A_{quad}\right)^{-1} |
557 |
\label{eq:dielectricFromQuadrupoles} |
558 |
\end{equation} |
559 |
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, |
560 |
\begin{equation} |
561 |
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} |
562 |
\label{eq:geometricalFactor} |
563 |
\end{equation} |
564 |
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. |
565 |
\begin{table} |
566 |
\caption{Expressions for the quadrupolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter |
567 |
($\alpha$) and the cutoff radius ($r_c$). The dimension of the correction factor is $ length^{-2}$ in case of quadrupolar fluid.} |
568 |
\label{tab:B} |
569 |
\centering |
570 |
\resizebox{\columnwidth}{!}{% |
571 |
|
572 |
\begin{tabular}{l|c|c|c|c|} |
573 |
|
574 |
Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ &$A_\mathrm{quadrupoles}$ \\\hline |
575 |
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}$ \\ |
576 |
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}$ \\ |
577 |
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)$\\ |
578 |
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 |
579 |
\end{tabular}% |
580 |
} |
581 |
\end{table} |
582 |
|
583 |
\section{Screening of Charges by Multipolar Fluids} |
584 |
In this method, we will measure the interaction between a positive and negative charge at varying distances after introducing a dipolar (or quadrupolar) material between them. The potential of mean force (PMF) between two ions in a liquid is obtained by constraining their distance and measuring the mean constraint force required to hold them at a fixed distance $r.$ The PMF is obtained from a sequence of simulations can be expressed as \cite{Wilfred07}, |
585 |
\begin{equation} |
586 |
w(r) = \int_{r_o}^{r}\braket{\frac{\partial f}{\partial r'}}dr' + 2kT \log(r/r_o) + w(r_o), |
587 |
\label{eq:pmf} |
588 |
\end{equation} |
589 |
where $\braket{\partial f/\partial r'}$ is the mean constraint force, $2kT \log(r/r_o)$ is the Fixman factor, and $r_o$ is the reference position. The potential energy between two charges at separation $r$ is given by, |
590 |
\begin{equation} |
591 |
w(r) = \frac{1}{4\pi\epsilon}\frac{q_1q_2}{r}, |
592 |
\end{equation} |
593 |
where $\epsilon$ is the dielectric constant for the fluid. |
594 |
|
595 |
The quadrupole molecule can couple with the gradient of the electric field and the two opposite point charges produces non-zero value of both electric field as well as gradient of the field. Therefore, this simulation set up can be used to determine the dielectric constant for both dipolar and quadrupolar fluid. |
596 |
\label{sec:PMF} |
597 |
\section{Simulation Methodology} |
598 |
We have used three different simulation methods: i) external field perturbation, ii) fluctuation formula, and iii) potential of mean force (PMF), to calculate dielectric properties for dipolar and quadrupolar fluid. For the dipolar system we calculated macroscopic polarzability using first two methods and derived the dielectric constant using polarizability and correction factor (see equation \ref{correctionFormula}). Similarly we used equation (\ref{eq:pmf}) to calculate screening factor from dipolar fluid using PMF method. For quadrupolar fluid, we have calculated quadrupolarizablity using fluctuation formula and external field perturbation and derived quadrupolar susceptibility using quadrupolarizability and correction factor (equation \ref{eq:quadrupolarSusceptiblity}). Since dielectric constant due to quadrupolar fluid depends on the nature of gradient of the field applied in the system, we have used geometrical factor (in equation \ref{eq:geometricalFactor}) and quadrupolar susceptibility to evaluate dielectric constant for two ions dissolved quadrupolar fluid (see equation \ref{eq:dielectricFromQuadrupoles}) . The the dielectric constant evaluated using equation (\ref{eq:dielectricFromQuadrupoles}) has been compared with the result evaluated from PMF method (i.e. equation \ref{eq:pmf}). We have also used three different test systems for both dipolar and quadrupolar fluids to calculate dielectric properties. The parameters used in the test systems are given in table \ref{Tab:C}. |
599 |
|
600 |
\begin{table} |
601 |
\caption{\label{Tab:C}} |
602 |
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
603 |
& \multicolumn{2}{c|}{LJ parameters} & |
604 |
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
605 |
Test system & $\sigma$& $\epsilon$ & $C$ & $D$ & |
606 |
$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass & $I_{xx}$ & $I_{yy}$ & |
607 |
$I_{zz}$ \\ \cline{6-8}\cline{10-12} |
608 |
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
609 |
\AA\textsuperscript{2})} \\ \hline |
610 |
Stockmayer fluid & 3.41 & 0.2381 & - & 1.4026 &-&-&-& 39.948 & 11.613 & 11.613 & 0.0 \\ |
611 |
Quadrupolar fluid & 2.985 & 0.265 & - & - & 0.0 & 0.0 &-2.139 & 18.0153 & 43.0565 & 43.0565 & 0.0 \\ |
612 |
\ce{q+} & 1.0 & 0.1 & +1 & - & - & - & - & 22.98 & - & - & - \\ |
613 |
\ce{q-} & 1.0 & 0.1 & -1 & - & - & - & - & 22.98 & - & - & - \\ \hline |
614 |
\end{tabularx} |
615 |
\end{table} |
616 |
|
617 |
First test system consists of point dipolar or quadrupolar molecules in the presence of constant field or gradient field. Since there is no isolated charge within the system, the divergence of the field should be zero $ i.e. \vec{\nabla} .\vec{E} = 0$. This condition is satisfied by selecting applied potential as described in Appendix \ref{Ap:fieldOrGradient}. When constant electric field or field gradient applied to the system, the molecules align along the direction of the applied field. We have calculated ensemble average of the box dipole or quadrupole moment as a response field or field gradient. Similarly the macroscopic polarizability of the system is derived using ratio between system multipolar moment and applied field or field gradient. This method works properly when the system is at the linear response region of field or field gradient. |
618 |
|
619 |
Second test system consists of box of point dipolar or quadrupolar molecules is simulated for 1 ns in NVE ensemble after equilibration in the absence of any external perturbation. The fluctuation of the ensemble average of the box multipolar moment, $\braket{A^2} - \braket{A}^2 $ where A is box dipolar or quadrupolar moment, is measured at the fixed temperature and density for a given multipolar fluid. Finally the macroscopic polarizability of the system at a particular density is derived using equation (\ref{flucQuad}). |
620 |
|
621 |
Final system consists of dipolar or quadrupolar fluids with two oppositely charged ions immersed in it. These ions are constraint to be at fixed distance throughout the simulation. We run separate simulations for the different constraint distances. We calculated the screening factor using ratio between the force between the two ions in the absence of medium and the average constraint force during the simulation. Since the constraint force is pretty noisy we run each simulation for long time ($\sim 5$ ns) to reduce simulation error. |
622 |
|
623 |
\subsection{Implementation} |
624 |
We have used real-space electrostatic methods implemented in OpenMD \cite{openmd2.3} software to evaluate electrostatic interactions between the molecules. In our simulations we used all three different real-space electrostatic methods: SP, GSF, and TSF developed in the previous paper \cite{PaperI} in the series. The radius of the cutoff sphere is taken to be $12 \r{A}$. Each real space method can be tuned using different values of damping parameter. We have selected ten different values of damping parameter (unit-${\r{A}}^{-1}$); 0.0, 0.05, 0.1, 0.15, 0.175, 0.2, 0.225, 0.25, 0.3, and 0.35 in our simulations for dipolar system. Since quadrupolar interactions are less sensitive to damping parameter as compared to dipole we selected only four values of dampig alpha 0.0, 0.1 0.2, and 0.3 $\AA-1$ in our simulation. The short range interaction in the simulations is incorporated using 6-12 Lennard Jones interactions. |
625 |
|
626 |
To derive the box multipolar (dipolar or quadrupolar) moment, we added the component each individual molecule and taken ensemble average of each snapshot in the entire simulation. The first component of the fluctuation of the dipolar moment is derived by using relation $\braket{M^2} = \braket{{M_x}^2 + {M_y}^2 + {M_z}^2}$, where $M_x$, $M_y$, and $M_z $ are x, y and z components of the box dipolar moment. Similarly the first term in the quadrupolar system is derived using relation $ \braket{Q^2} = \braket{3 Q:Q - \mathrm{Tr}Q^2} $, where $ Q $ is the box quadrupole moment, double dot represent the outer product of the quadrupolar matrices, and $TrQ$ is the trace of the box quadrupolar moment. The second component of the fluctuation formula has been derived using square of the ensemble average of the box dipole moment. The applied constant field or field gradient in the test systems has been taken in the form described in the Appendix \ref{Ap:fieldOrGradient}. |
627 |
\subsection{Model systems} |
628 |
To evaluate dielectric properties for dipolar systems using perturbation and fluctuation formula methods, we have taken system of 2048 Stockmayer molecules with reduced density $ \rho^* = 0.822$, temperature $T^* = 1.15 $, moment of inertia $I^* = 0.025 $, and dipole moment $ \mu^* = \sqrt{3.0} $. Test systems are equilibrated for 500 ps and run for $1\; ns$ and components of box dipole moment are obtained at every femtosecond. The systems have been run in the presence of constant external field from $ 0 - 10\; \times\; 10^{-4}\;V/{\r{A}}$ in the step of $ 10 ^{-4}\; V/\r{A}$ for each simulation. For pmf method, Two dipolar molecules in the above system are converted into $q+$ and $q-$ ions and constrained to remain in fixed distance in simulation. The constrained distance is varied from $5\;\r{A} - 12\; \r{A} $ for different simulations. In pmf method all simulations are equilibrated for 500 ps in NVT ensemble and run for 5 ns in NVE ensemble to print constraint force at an interval of 20 fs. |
629 |
|
630 |
Quadrupolar systems consists 4000 linear point quadrupoles with density $ 2.338\; g/cm^3$ at temperature $ 500\; ^oK $. For both perturbation and fluctuation methods, test systems are equalibrated for 200 ps in NVT ensemble and run for 500 ps in NVE ensemble. To find the ensemble average of the box quadrupole moment and fluctuation, the components of box quadrupole moments are printed every 100 fs. Each simulations repeated at different values of applied gradients from $ 0 - 9 \times 10^{-2}\; V/\AA^2 $. To calculate dielectric constant using pmf method, two ions in the systems are converted into $q+$ and $q-$ ions and constrained to remain at fixed distance in the simulation. These constraint distances are varied from $5\;\AA - 12\; \AA $ at the step of $0.1\; \AA $ for different simulations. For calculating dielectric constant, the test systems are run for 500 ps to equlibrate and run for 5 ns to print constraint force at a time interval of 20 fs. |
631 |
|
632 |
\section{Results} |
633 |
\subsection{Dipolar fluid} |
634 |
The macroscopic polarizability ($\alpha_D$) for Stockmayer fluid evaluated using perturbation (left) and fluctuation (right) methods as shown in figure \ref{fig:dielectricDipole} . The polarizability obtained from the both perturbation and fluctuation methods show excellent agreement with each other. The result shows that polarizability hugely depends on the values of the damping parameter in the case of Shifted Potential (SP) and Gradient Shifted force (GSF) methods. But the results obtained using Taylor Shifted Force (TSF) show slight dependence on the damping alpha. To get correct dielectric properties from the polarizability we used a correction factor as we discussed in section \ref{sec:corrFactor}. The correction factor for dipolar fluids ($\mathrm{A}_\mathrm{dipole}$ ) for the different values of damping parameter is also plotted in figure \ref{fig:dielectricDipole}. For TSF method, $ \mathrm{A}_\mathrm{dipole} = 1 $ for all the values of damping alpha. Therefore we do not need to perform any correction to obtain dielectric constant from the polarizability. The values of $\mathrm{A}_\mathrm{dipole} $ varies with the damping parameter in the case of SP and GSF methods therefore correction factor should be incorporated to calculate the dielectric constant ($\epsilon$) from the polarzability. The dielectric constant obtained from the polarizability and correction factor is also plotted in the figure \ref{fig:dielectricDipole}. The dielectric constant evaluated using SP and GSF method shows good agreement with the previous simulation result (straight dotted line in the lower panels)\cite{NeumannI83} for the damping parameter from $0.25\; \mathrm{to}\; 0.3$ $\AA-1$. For lower value of the damping alpha the dielectric constant is much deviated away from the previous simulation result. To investigate furthermore we performed inverse transformation of expected dielectric constant \cite{NeumannI83}, using correction factor and calculated expected polarizability for the simulation, which is represented by the dotted line in the upper-left and upper-right panels. In the plot we can clearly see that expected ploarizability for the SP and GSF methods are very close to results obtained from the simulations. Sometime even they are within the range of the error-bar but when the dielectric constant is calculated using simulation result and correction formula, it is found to be deviated much away from the expected result for the lower values of the damping alpha. This tell us that the correction formula (equation \ref{correctionFormula}) for the dipolar fluid is very sensitive when the value of $\mathrm{A}_\mathrm{dipole}$ is deviated much away from 1 for the small values of damping parameter. Therefore it is always better to select damping alpha from 0.25 to 0.3 $\AA^{-1}$ to evaluate dielectric constant using perturbation and fluctuation formula while using SP and GSF methods in simulation. |
635 |
|
636 |
\begin{figure} |
637 |
\includegraphics[width=3in]{dielectricFinal_Dipole.pdf} |
638 |
\caption{In the figure, $\alpha_D$, $\mathrm{A}_\mathrm{dipole}$, and $\epsilon$ are polarizability, correction factor, and dielectric constant for Stockmayer fluid. Plots in the left panel show results for (a) perturbation method and right (b) fluctuation method.} |
639 |
\label{fig:dielectricDipole} |
640 |
\end{figure} |
641 |
|
642 |
We have also evaluated the screening factor between two oppositely charged ions when they are placed in the Stockmayer fluid as shown in figure \ref{fig:ScreeningFactor_Dipole}. The screening factor have been evaluated using equation \ref{eq:pmf}. This screening factor is analogous to dielectric screening but there is subtle differences between them. Actually screening factor measures a local property of the ions in the fluid and depends on ion-dipole and dipole-dipole interactions. These interactions depends on the distance between ions and electrostatic interaction methods implemented in the simulations. But we can expect screening to be close to dielectric constant when the field due to ions is very gentle in perturbation. It is possible when distance between ions are far apart and values of damping parameter relatively higher. In the figure \ref{fig:ScreeningFactor_Dipole} we can observe that for the higher value of damping alpha i.e $\alpha = 0.2 \AA^{-1}$ and $0.3 \AA^{-1}$ and large separation between ions, the screening factor approaches to the dielectric constant. We can also observe that for TSF method has got higher screening factor as compared SP and GSF method even for lower damping alpha. It is because of dipole-dipole interactions do not need any correction factor in the case of TSF as we discussed earlier. But presence of ion-dipole interaction causing local effect even in the case of TSF method for small ions separation. |
643 |
|
644 |
\begin{figure} |
645 |
\includegraphics[width=3in]{ScreeningFactor_Dipole.pdf} |
646 |
\caption{Figure shows screening factor between two oppositely charged ions immersed in Stockmayer fluid as a function of ions-separation for the different values damping alpha and electrostatic methods.} |
647 |
\label{fig:ScreeningFactor_Dipole} |
648 |
\end{figure} |
649 |
|
650 |
|
651 |
\subsection{Quadrupolar fluid} |
652 |
The polarizability ($\alpha_Q$), correction factor ($\mathrm{A_{quad}}$), and susceptibility ($\chi$) for the quadrupolar fluid is plotted against damping parameter in the figure \ref{fig:dielectricQuad}. Different than the dipolar fluid, the polarizability and susceptibility for the quadrupolar fluid have a dimension of $\AA ^2$. Although susceptibility has got dimension, it is a measure of macroscopic quadrupolar properties at a given temperature and density \cite{JeonI03, JeonII03}. The left panel in the plot showing results obtained from the perturbation whereas the results from the fluctuation formula are plotted in the right panel. In the figure, we can observe that the polarizablity evaluated from the perturbation (top-left) shows excellent agreement with the result obtained from the fluctuation formula (top-right). The susceptibility for the quadrupolar fluid is obtained from quadrupolarizability and correction factor using equation \ref{eq:quadrupolarSusceptiblity} which is plotted at the bottom-left and bottom-right in the figure \ref{fig:dielectricQuad}. All three methods: SP, GSF, and TSF tends to produce same value of susceptibility for the given value damping alpha. This shows that susceptibility derived using correction factor and polarizablity is tends to be independent of the electrostatic method implemented in the simulation for quadrupolar fluid. |
653 |
\begin{figure} |
654 |
\includegraphics[width=3in]{polarizabilityFinal_Quad.pdf} |
655 |
\caption{In the figure, $\alpha_Q$, $\mathrm{A}_\mathrm{quad}$, and $\chi$ are polarizability, correction factor, and susceptibility for the quadrupolar fluid. Plots in the left panel show results for (a) perturbation method and right (b) fluctuation method.} |
656 |
\label{fig:dielectricQuad} |
657 |
\end{figure} |
658 |
|
659 |
The actual test of the susceptibility can be made by comparing the results obtained from the perturbation and fluctuation with the result calculated from more direct, potential of mean force (PMF), method. For this purpose, we have calculated the expected dielectric constant from the susceptibility and geometrical factor using equation \ref{eq:dielectricFromQuadrupoles}. Since the dielectric constant for quadrupolar fluid is geometry-dependent, it is not a actual measure of the property of the quadrupolar fluid. This geometrical factor varies with the variation in ions separation as well as damping parameter therefore the dielectric constant can be refereed as screening factor. The screening factor plotted against ions separation, evaluated from the perturbation and fluctuation method are shown in right and central panel of the figure \ref{fig:screenQuad}. Here the susceptibility is calculated from the simulation and the geometrical factor is evaluated analytically, using field and field-gradient produced by ions. The right hand panel shows screening factor plotted against ions separation, for the different values of the damping alpha, obtained from the PMF method. The figure clearly shows that the screening factor obtained from the perturbation and fluctuation formula shows good agreement with the screening factor calculated using PMF method. Since there is no huge differences in quadrupole-quadruople interactions for various real-space methods, we do not observe large differences in the screen factors for SP, GSF, and TSF methods. |
660 |
|
661 |
\begin{figure} |
662 |
\includegraphics[width=3in]{ScreeningFactor_Quad.pdf} |
663 |
\caption{} |
664 |
\label{fig:screenQuad} |
665 |
\end{figure} |
666 |
|
667 |
\section{Discussion} |
668 |
We have used the perturbation and fluctuation formula to evaluate dielectric properties for dipolar and quadrupolar fluids implementing SP, GSF, and TSF methods in the simulation. Similarly the correction factors for the both dipolar and quadrupolar systems have also been calculated for all SP, GSF, and TSF methods. We have also derived screening factor between two ions immersed in the dipolar and quadrupolar fluids using PMF method. The result from the perturbation and fluctuation formula are compared with PMF method to test the accuracy of the simulation result. |
669 |
|
670 |
For the dipolar fluid, the polarizability evaluated using the perturbation and fluctuation methods show excellent agreement with each other. The dielectric constant evaluated using polarizability and correction factor agree with the previous simulation results \cite{NeumannI83} for the damping parameter 0.25 - 0.3$\AA^{-1}$ and the electrostatic interaction methods, SP and GSF implemented in the simulation. Since the correction factor for TSF, $\mathrm{A_{dipole}} = 1$, it produces dielectric constant more closer to previous simulation results \cite{NeumannI83} for all values damping parameter. This method also produces best dielectric constant for damping parameter 0.15 - 0.25$\AA^{-1}$. We have also found that correction formula (equation ~\ref{correctionFormula}) is very sensitive, when the correction factor ($\mathrm{A_{dipole}}$) much away from the 1 for the SP and GSF methods. Therefore it is better to choose damping parameter from 0.25 - 0.3$\AA^{-1}$ to evaluate dielectric constant using SP and GSF methods. The screening factor between ions is found to be closer to dielectric constant when the distance between the ions are far apart ($\sim 9 \AA$) and damping parameter $0.2$ to $0.3 \AA^{-1}$. |
671 |
|
672 |
The quadpolarizability evaluated from the both perturbation and fluctuation methods shows excellent agreement with each other for the case of quadrupolar fluid. The susceptibility is calculated from the quadpolarizability and correction factor tends to produce same result for all; SP, GSF, and TSF, methods. Similarly the screening factor calculated using susceptibility and geometrical factor shows excellent agreement with the result obtained from the PMF method. |
673 |
|
674 |
|
675 |
\appendix |
676 |
\section{Point-multipolar interactions with a spatially-varying electric field} |
677 |
|
678 |
We can treat objects $a$, $b$, and $c$ containing embedded collections |
679 |
of charges. When we define the primitive moments, we sum over that |
680 |
collections of charges using a local coordinate system within each |
681 |
object. The point charge, dipole, and quadrupole for object $a$ are |
682 |
given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively. |
683 |
These are the primitive multipoles which can be expressed as a |
684 |
distribution of charges, |
685 |
\begin{align} |
686 |
C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\ |
687 |
D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\ |
688 |
Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in } a} q_k |
689 |
r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} |
690 |
\end{align} |
691 |
Note that the definition of the primitive quadrupole here differs from |
692 |
the standard traceless form, and contains an additional Taylor-series |
693 |
based factor of $1/2$. In Paper 1, we derived the forces and torques |
694 |
each object exerts on the others. |
695 |
|
696 |
Here we must also consider an external electric field that varies in |
697 |
space: $\mathbf E(\mathbf r)$. Each of the local charges $q_k$ in |
698 |
object $a$ will then experience a slightly different field. This |
699 |
electric field can be expanded in a Taylor series around the local |
700 |
origin of each object. A different Taylor series expansion is carried |
701 |
out for each object. |
702 |
|
703 |
For a particular charge $q_k$, the electric field at that site's |
704 |
position is given by: |
705 |
\begin{equation} |
706 |
E_\gamma + \nabla_\delta E_\gamma r_{k \delta} |
707 |
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta} |
708 |
r_{k \varepsilon} + ... |
709 |
\end{equation} |
710 |
Note that the electric field is always evaluated at the origin of the |
711 |
objects, and treating each object using point multipoles simplifies |
712 |
this greatly. |
713 |
|
714 |
To find the force exerted on object $a$ by the electric field, one |
715 |
takes the electric field expression, and multiplies it by $q_k$, and |
716 |
then sum over all charges in $a$: |
717 |
|
718 |
\begin{align} |
719 |
F_\gamma &= \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta} |
720 |
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta} |
721 |
r_{k \varepsilon} + ... \rbrace \\ |
722 |
&= C_a E_\gamma + D_{a \delta} \nabla_\delta E_\gamma |
723 |
+ Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma + |
724 |
... |
725 |
\end{align} |
726 |
|
727 |
Similarly, the torque exerted by the field on $a$ can be expressed as |
728 |
\begin{align} |
729 |
\tau_\alpha &= \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\ |
730 |
& = \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k |
731 |
r_{k\beta} E_\gamma(\mathbf r_k) \\ |
732 |
& = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma |
733 |
+ 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta |
734 |
E_\gamma + ... |
735 |
\end{align} |
736 |
|
737 |
The last term is essentially identical with form derived by Torres del |
738 |
Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation |
739 |
utilized a traceless form of the quadrupole that is different than the |
740 |
primitive definition in use here. We note that the Levi-Civita symbol |
741 |
can be eliminated by utilizing the matrix cross product in an |
742 |
identical form as in Ref. \onlinecite{Smith98}: |
743 |
\begin{equation} |
744 |
\left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta |
745 |
\left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta} |
746 |
-\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta} |
747 |
\right] |
748 |
\label{eq:matrixCross} |
749 |
\end{equation} |
750 |
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of |
751 |
the matrix indices. In table \ref{tab:UFT} we give compact |
752 |
expressions for how the multipole sites interact with an external |
753 |
field that has exhibits spatial variations. |
754 |
|
755 |
\begin{table} |
756 |
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque |
757 |
$(\mathbf{\tau})$ expressions for a multipolar site embedded in an |
758 |
electric field with spatial variations, $\mathbf{E}(\mathbf{r})$. |
759 |
\label{tab:UFT}} |
760 |
\begin{tabular}{r|ccc} |
761 |
& Charge & Dipole & Quadrupole \\ \hline |
762 |
$U$ & $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\ |
763 |
$\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ & $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\ |
764 |
$\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$ |
765 |
\end{tabular} |
766 |
\end{table} |
767 |
\section{Gradient of the field due to quadrupolar polarization} |
768 |
\label{singularQuad} |
769 |
In this section, we will discuss the gradient of the field produced by |
770 |
quadrupolar polarization. For this purpose, we consider a distribution |
771 |
of charge ${\rho}(r)$ which gives rise to an electric field |
772 |
$\vec{E}(r)$ and gradient of the field $\vec{\nabla} \vec{E}(r)$ |
773 |
throughout space. The total gradient of the electric field over volume |
774 |
due to the all charges within the sphere of radius $R$ is given by |
775 |
(cf. Jackson equation 4.14): |
776 |
\begin{equation} |
777 |
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r = -\int_{r=R} R^2 \vec{E}\;\hat{n}\; d\Omega |
778 |
\label{eq:8} |
779 |
\end{equation} |
780 |
where $d\Omega$ is the solid angle and $\hat{n}$ is the normal vector |
781 |
of the surface of the sphere which is equal to |
782 |
$sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} + |
783 |
cos[\theta]\hat{z}$ |
784 |
in spherical coordinates. For the charge density ${\rho}(r')$, the |
785 |
total gradient of the electric field can be written as (cf. Jackson |
786 |
equation 4.16), |
787 |
\begin{equation} |
788 |
\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 |
789 |
\label{eq:9} |
790 |
\end{equation} |
791 |
The radial function in the equation (\ref{eq:9}) can be expressed in |
792 |
terms of spherical harmonics as (cf. Jackson equation 3.70), |
793 |
\begin{equation} |
794 |
\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) |
795 |
\label{eq:10} |
796 |
\end{equation} |
797 |
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, |
798 |
\begin{equation} |
799 |
\begin{split} |
800 |
\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 \\ |
801 |
&= -\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 |
802 |
' |
803 |
\end{split} |
804 |
\label{eq:11} |
805 |
\end{equation} |
806 |
The gradient of the product of radial function and spherical harmonics |
807 |
is given by (cf. Arfken, p.811 eq. 16.94): |
808 |
\begin{equation} |
809 |
\begin{split} |
810 |
\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 |
811 |
{\partial}{\partial r}+\frac{l}{r} \right]f(r)\; Y_{l, l-1, m}(\theta, \phi). |
812 |
\end{split} |
813 |
\label{eq:12} |
814 |
\end{equation} |
815 |
Using equation (\ref{eq:12}) we get, |
816 |
\begin{equation} |
817 |
\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}}, |
818 |
\label{eq:13} |
819 |
\end{equation} |
820 |
where $ Y_{l,l+1,m}(\theta, \phi)$ is the vector spherical harmonics |
821 |
which can be expressed in terms of spherical harmonics as shown in |
822 |
below (cf. Arfkan p.811), |
823 |
\begin{equation} |
824 |
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}, |
825 |
\label{eq:14} |
826 |
\end{equation} |
827 |
where $C(l+1,1,l|m_1,m_2,m)$ is a Clebsch-Gordan coefficient and |
828 |
$\hat{e}_{m_2}$ is a spherical tensor of rank 1 which can be expressed |
829 |
in terms of Cartesian coordinates, |
830 |
\begin{equation} |
831 |
{\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}} |
832 |
\label{eq:15} |
833 |
\end{equation} |
834 |
The normal vector $\hat{n} $ can be expressed in terms of spherical tensor of rank 1 as shown in below, |
835 |
\begin{equation} |
836 |
\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) |
837 |
\label{eq:16} |
838 |
\end{equation} |
839 |
The surface integral of the product of $\hat{n}$ and |
840 |
${Y_{l+1}}^{m_1}(\theta, \phi)$ gives, |
841 |
\begin{equation} |
842 |
\begin{split} |
843 |
\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 \\ |
844 |
&= \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 \\ |
845 |
&= \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), |
846 |
\end{split} |
847 |
\label{eq:17} |
848 |
\end{equation} |
849 |
where ${Y_{l}}^{-m} = (-1)^m\;{{Y_{l}}^{m}}^* $ and |
850 |
$ \int {{Y_{l}}^{m}}^*\;{Y_{l'}}^{m'}\;d\Omega = |
851 |
\delta_{ll'}\delta_{mm'} $. |
852 |
Non-vanishing values of equation \ref{eq:17} require $l = 0$, |
853 |
therefore the value of $ m = 0 $. Since the values of $ m_1$ are -1, |
854 |
1, and 0 then $m_2$ takes the values 1, -1, and 0, respectively |
855 |
provided that $m = m_1 + m_2$. Equation \ref{eq:11} can therefore be |
856 |
modified, |
857 |
\begin{equation} |
858 |
\begin{split} |
859 |
\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( |
860 |
1, 1, 0|0,0,0)\;{\hat{e}_{0}}{\hat{e}_{0}} ]\; d^3r'. |
861 |
\end{split} |
862 |
\label{eq:18} |
863 |
\end{equation} |
864 |
After substituting ${Y^*}_{00} = \frac{1}{\sqrt{4\pi}} $ and using the |
865 |
values of the Clebsch-Gorden coefficients: $ C(1, 1, 0|-1,1,0) = |
866 |
\frac{1}{\sqrt{3}}, \; C(1, 1, 0|-1,1,0)= \frac{1}{\sqrt{3}}$ and $ |
867 |
C(1, 1, 0|0,0,0) = -\frac{1}{\sqrt{3}}$ in equation \ref{eq:18} we |
868 |
obtain, |
869 |
\begin{equation} |
870 |
\begin{split} |
871 |
\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)\\ |
872 |
&= - \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). |
873 |
\end{split} |
874 |
\label{eq:19} |
875 |
\end{equation} |
876 |
Equation (\ref{eq:19}) gives the total gradient of the field over a |
877 |
sphere due to the distribution of the charges. For quadrupolar fluids |
878 |
the total charge within a sphere is zero, therefore |
879 |
$ \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = 0 $. Hence the quadrupolar |
880 |
polarization produces zero net gradient of the field inside the |
881 |
sphere. |
882 |
|
883 |
\section{Applied field or field gradient} |
884 |
\label{Ap:fieldOrGradient} |
885 |
|
886 |
To satisfy the condition $ \nabla . E = 0 $, within the box of molecules we have taken electrostatic potential in the following form |
887 |
\begin{equation} |
888 |
\begin{split} |
889 |
\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. \\ |
890 |
& \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), |
891 |
\end{split} |
892 |
\label{eq:appliedPotential} |
893 |
\end{equation} |
894 |
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, |
895 |
\[\bf{E} |
896 |
=\frac{g_o}{2} \left(\begin{array}{ccc} |
897 |
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 \\ |
898 |
(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 \\ |
899 |
(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 |
900 |
\end{array} \right).\] |
901 |
The gradient of the applied field derived from the potential can be written in the following form, |
902 |
\[\nabla\bf{E} |
903 |
= \frac{g_o}{2}\left(\begin{array}{ccc} |
904 |
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 \\ |
905 |
(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 \\ |
906 |
(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 |
907 |
\end{array} \right).\] |
908 |
\newpage |
909 |
|
910 |
\bibliography{dielectric_new} |
911 |
|
912 |
\end{document} |