On the reaction field for interaction site models of polar systems

14 Pages·0.2 MB·English

by  Igor P. Omelyan

On the reaction field for interaction site
9
9
9
1 models of polar systems

n
a
J

9
Igor P. OMELYAN
1

]
h
Institute for Condensed Matter Physics, National Ukrainian Academy of Sciences
p
-
p 1 Svientsitsky St., UA-290011 Lviv, Ukraine ∗
m
o
c
.
s
c
i
s
y
h
p
[
Abstract

1
v
It is rigorously shown that the fluctuation formula, which is used in simulations
2
3
to calculate the dielectric constant of interaction site models, corresponds to the
0
1
reaction field with an individual site cut-off rather than with the usual molecular
0
9
center of mass truncation. Within themolecular cut-off scheme, amodifiedreaction
9
/
s field is proposed. An influence of the truncation effects is discussed and examined
c
i
s by actual Monte Carlo simulations for a MCY water model.
y
h
p
Keywords: Dielectric constant; Reaction field; Interaction site models; Computer
:
v
i simulations
X
r
a PACS numbers: 77.22.-d; 24.60.-k; 61.20.Ja
∗
E-mail: [email protected]
1
1 Introduction
The calculation of dielectric quantities by computer experiment requires an explicit
consideration of effects associated with the truncation of long-range interactions. The
concretesuccessinthisdirectionhasbeenachievedwithinthereactionfield(RF)geometry
[1–5]. As a result, computer adapted dielectric theories have been proposed [6–10]. In the
framework of these theories, a bulk dielectric constant can be determined on the basis of
a fluctuation formula via correlations obtained in simulations for finite samples. However,
main attention in the previous investigations has been focused on polar systems with the
point dipole interaction. As is now well established, the model of point dipoles can not
reproduce adequately features of real polar liquids.
At the same time, attempts to apply the RF geometry for more realistic interaction
site (IS) models have also been made [11–13]. However, acting within a semiphenomeno-
logical approach, it was not understood how to perform the truncation of intermolecular
potentials. Asaconsequence, themolecular cut-offandtheusualpointdipoleRF(PDRF)
havebeenassumed. Obviously, suchanapproachincludeseffectsconnectedwithfiniteness
of the molecule inconsistently. Indeed, the interdipolar potential is replaced by site-site
Coulomb interactions, whereas the RF is remained in its usual form. An additional com-
plication for IS models consists in a spatial distribution of charges and this fact is not
taken into account by the standard PDRF geometry.
In the present paper we propose two alternative approaches to remedy this situation.
The first one follows from the usual fluctuation formula which is constructed, however,
on the microscopic operator of polarization density for IS models. This leads to an ISRF
geometry, where the cut-off radius is applied with respect to individual charges rather
than to the molecule as a whole. Nevertheless, the molecular cut-off scheme can also
be acceptable, but the reaction field together with the fluctuation formula need to be
corrected. In the second approach a molecular RF (MRF) geometry is proposed and a
new quadrupole term is identified. On the basis of a MCY water model we show that
uncertainties of the dielectric quantities can besignificant if the standard PDRF geometry
is used in computer simulations.
2
2 Interaction site reaction field
We consider an isotropic, classical system of N identical molecules enclosed in volume
V. The microscopic electrostatic field created by the molecules at point r V is equal to
∈
N r ra
Eˆ(r) = q − i = L(r r′)Qˆ(r′)dr′ , (1)
a r ra 3 −
Xi=1Xa | − i| VZ
where ra denotes the position for charge q of ith molecule, Qˆ(r) = q δ(r ra) is
i a i,a a − i
the microscopic operator of charge density, L(ρ) = ∇ 1/ρ and the suPmmation extends
−
over all molecules and charged sites. For the investigation of dielectric properties, it is
more convenient to rewrite the electric field (1) in the polarization representation
4π
Eˆ(r) = T(r r′)Pˆ(r′)dr′ = Pˆ(r)+ lim T(r r′)Pˆ(r′)dr′ . (2)
− − 3 ρ→+0 −
Z Z
V V
ρ<|r−r′|
Here T(ρ) = ∇∇ 1/ρ is the dipole-dipole tensor, Pˆ(r) denotes the microscopic operator
of polarization density, defined as ∇·Pˆ(r) = Qˆ(r), and the singularity lim T(ρ) =
ρ→0
−
4π/3 δ(ρ)I has been avoided, where I is the unit tensor of the second rank. The both
−
charge (1) and polarization (2) representations are equivalent and applicable for infinite
(N,V ) systems.
→ ∞
In simulations, which deal with finite samples, the sum (1) can not be calculated
exactly taking into account an infinitely large number of terms. Therefore, we must
restrict ourselves to a finite set of terms in (1) or to a finite range of the integration in (1)
and (2) for which r r′ R, where R is a cut-off radius. Now the following problem
| − | ≤
appears. How to estimate the cut-off field caused by the integration over unaccessible
region r r′ > R? The solution of this problem has been found for the first time
| − |
for systems with point dipoles in the RF geometry. The result for conducting boundary
conditions is [7, 8]
4π I
Eˆ(r) EˆRF(r) = Pˆ(r)+ lim T(r r′)+ Pˆ(r′)dr′ , (3)
≈ R − 3 ρ→+0 − R3
Z (cid:18) (cid:19)
V, tbc
ρ<|r−r′|≤R
where a cubic finite sample and toroidal boundary conditions (TBC) have been used, so
that R √3 V/2. The additional term I/R3 in the right-hand site of (3) describes the
≤
3
RF which is used for an approximation of the real cut-off field. For a pure spherical cut-
off (SC) without the RF correction, we have EˆSC(r) = γ( r r′ )L(r r′)Qˆ(r′)dr′,
R | − | −
where γ(ρ) = 1 if ρ R and γ(ρ) = 0 otherwise. ObvZiously, that lim EˆSC(r) =
≤ R→∞ R
lim EˆRF(r) = Eˆ(r).
R→∞ R
Let us perform the spatial Fourier transform (k) = dre−ik·r (r) for arbitrary
F F
functions . Then one obtains R
F
EˆSC(k) = L(k)Qˆ(k) , EˆRF(k) = 4πPˆ(k)+ T(k)+4πj1(kR)I Pˆ(k) , (4)
R R − 3 kR
(cid:16) (cid:17)
where
ik 4π j (kR)
L(k) = 4π 1 j (kR) , T(k) = 1 3 1 3kˆkˆ I , (5)
− − 0 k2 − 3 − kR −
(cid:16) (cid:17) (cid:16) (cid:17)(cid:16) (cid:17)
Qˆ(k) = q e−ik·rai = ik·Pˆ(k), k = 2πn/√3 V is one of the allowed wavevectors of
i,a a
−
the recipProcal lattice, n designates a vector with integer components, k = k , kˆ = k/k
| |
and j (z) = sin(z)/z, j (z) = cos(z)/z +sin(z)/z2 are the spherical Bessel functions of
0 1 −
zero and first order, respectively. In view of (5), the relations (4) transform into
j (kR)
EˆSC(k) = 4π 1 j (kR) Pˆ (k) , EˆRF(k) = 4π 1 3 1 Pˆ (k) , (6)
R − − 0 L R − − kR L
(cid:16) (cid:17) (cid:16) (cid:17)
where Pˆ (k) = kˆkˆ·Pˆ(k) = ikQˆ(k)/k2 is the longitudinal component of the microscopic
L
operator of polarization density.
It is easy to see from (6) that the both functions EˆSC(k) and EˆRF(k) tend to the same
R R
value Eˆ(k) = 4πPˆ (k) of the infinite system at R (k = 0). However, the results
L
− → ∞ 6
converge as R−1 for the pure SC scheme, while as R−2 in the RF geometry, i.e., more
quickly, becauseamainpartofthetruncationeffectsistakeninto accountbytheRF.This
is very important in our case, where we hope to reproduce features of infinite systems on
the basis of finite samples. That is why the pure truncation, which is standard for simple
fluids with short-range potentials, is generally not recommended for polar systems with
long-range nature of the dipolar interaction. The influence of the TBC and the difference
between micro- and canonical ensembles are of order N−1 R−3 [14] and, therefore,
∼
they can be excluded from our consideration. It is worth mentioning that electrostatic
fields are pure longitudinal. They can be defined via the longitudinal component of the
microscopic operator of polarization density, that is confirmed by Eq. (6).
4
LetusenclosethesysteminanexternalelectrostaticfieldE (r). Thematerialrelation
0
ˆ
between themacroscopicpolarizationP (k) = P (k) intheweak externalfield andto-
L L
D E
talmacroscopicfieldis4πP (k) = ε (k) 1 E (k),whereε (k)denotesthelongitudinal
L L − L L
(cid:16) (cid:17)
wavevector-dependent dielectric constant. Applying the first-order perturbation theory
with respect to E yields for rigid molecules Vk TP (k) = Pˆ (k)·Pˆ ( k) E (k),
0 B L L L − 0 0
D E
where k and T are Boltzmann’s constant and temperature, respectively, and ... is the
B h i0
equilibrium average in the absence of the external field. Then, taking into account that
E (k) = E (k)+ EˆRF(k) and eliminating E (k), we obtain the fluctuation formula
L 0 R 0
(cid:28) (cid:29)
ε (k) 1 9yG (k)
L − = L = 9yg (k) . (7)
ε (k) 1+27yG (k)j (kR)/(kR) L
L L 1
HereG (k) = Pˆ (k)·Pˆ ( k) Nµ2 isthelongitudinalcomponent ofthefinite-system
L L L − 0
D E .
wavevector-dependent Kirkwood factor, y = 4πNµ2 9Vk T and µ = µ = q ra
B | i| | a a i|
.
denotes the permanent magnitude of molecule’s dipole moment. It is necessaryPto note
that we consider rigid IS molecules so that effects associated with molecular and elec-
tronic polarizabilities are not included in our investigation. In the case of R , we
→ ∞
have j (kR)/(kR) 0 and computer adapted formula (7) reduces to the well-known fluc-
1 →
tuation formula for macroscopic systems in terms of the infinite-system Kirkwood factor
g (k) = lim G (k).
L R→∞ L
As was mentioned earlier, the electric field EˆRF in the form (3), (4) as well as the
R
fluctuation formula (7) have been proposed for the first time to investigate polar systems
ofpoint dipoles[8]. However, actingwithinasemiphenomenological framework, itwasnot
understood how to perform the truncation of the intermolecular potential ϕ at attempts
ij
to extend this formula for IS models. As a result, the molecular cut-off r = r r R,
ij i j
| − | ≤
where r is the center of mass for ith molecule, and the usual PDRF have been suggested
i
[11–13]:
q q µ ·µ
ϕ = a b i j , r R . (8)
ij ra rb − R3 ij ≤
Xa,b | i − j|
It is essentially to emphasize that the fluctuation formula (7) takes into account finite-
ness of the system explicitly by the factor j (kR)/(kR). As a result, if the system size is
1
sufficiently large (terms of order R−2 can be neglected), the bulk (N,V ) dielectric
→ ∞
constant can be reproduced via the finite-system Kirkwood factor G (k) which depends
L
5
on R in a characteristic way. However, to achieve this self-consistency in the evaluation
of the bulk dielectric constant, the equilibrium averaging in G (k) must be calculated for
L
systems with the intermolecular potential which leads exactly to the microscopic electric
field EˆRF(r) (3). As we shall below, the intermolecular potential (8) does not obey this
R
condition.
To derive the exact intermolecular potential in the charge representation, we perform
the inverse Fourier transform EˆRF(r) = 1 dkEˆRF(k)eik·r and obtain using (6)
R (2π)3 R
Z
∞
r ra 6 r ra 2
EˆRF(r) = q − i 1 | − i| j (kR)j (k r ra )dk . (9)
R a r ra 3  − π R 1 1 | − i| 
Xi,a | − i| Z0
 
Taking into account that 6 ∞j (kR)j (kρ)dk = ρ/R2 if ρ R and is equal to R/ρ2 if
π 0 1 1 ≤
ρ > R, we have R
r ra r ra 3
EˆRF(r) = q − i 1 | − i| if r ra R (10)
R Xi,a a|r −rai|3 − R3 ! | − i| ≤
and EˆRF(r) = 0 otherwise, where the first term in the right-hand side is the Coulomb
R
field, while the second contribution corresponds to the RF in the IS description.
In order to understand nature of this field, we consider a spherical cavity of radius R
with the center at point r, embedded in an infinite conducting medium. Let us place a
point charge q at point ra in the cavity, so that r ra R. The totalelectric field ea(r)
a i | − i| ≤ i
atpointr consists ofthefieldduetothechargeq andthefieldcreatedbyinduced charges
a
located on the surface of the cavity. According to the method of electrostatic images [5],
this last field can be presented as the field of an imaginary charge q∗ = q R/ r ra
a − a | − i|
which is located at point r∗a = r R2(r ra)/ r ra 2 outside the sphere. Then
i − − i | − i|
ea(r) = q (r ra)/ r ra 3 +q∗(r r∗a)/ r r∗a 3 = q (r ra)(1/ r ra 3 1/R3)
i a − i | − i| a − i | − i| a − i | − i| −
that is completely in line with the term of sum (10).
In the potential representation (EˆRF(r) = ∇Φ(r)), we obtain Φ(r) = φa(r),
R − i,a i
where φa(r) = q (1/ρa + 1ρa2/R3 +C), ρa = r ra and C is, in general, anParbitrary
i a i 2 i i | − i|
constant which for infinite systems is chosen as φa = 0. In our case, according
i|ρai→∞
to the toroidal boundary conventional, φa = 0 whence C = 3/2R−1. Then the
i|ρai=R −
intermolecular potential of interaction is ϕ = q φa(rb) = q φb(ra) = ϕab,
ij a,b b i j a,b a j i a,b ij
P P P
6
where
1 1 ra rb 2 3
q q + | i − j| , ra rb R
ϕab =  a b rai rbj 2 R3 − 2R! | i − j| ≤ (11)
| − |
ij  0 , ra rb > R
| i − j|

and the site-site cut-off is performed.
It is easily seen from (11) that the ISRF part 1 q q ra rb 2/R3 transforms into
2 a,b a b| i − j|
the usual form µ ·µ /R3 of point dipoles for r P R d only, where d = 2max δa
− i j ij ≤ − | i|
is the diameter of the molecule and δa = ra r . In the case if the molecular rather
i i − i
than the site-site cut-off is applied to the potential (11), this transformation is valid for
arbitraryr R. Moreover, inthelastcasetheconstantC = 3/2R−1 iscanceled owing
ij
≤ −
electroneutrality( q = 0)ofthemoleculeandwerecover theresult(8)ofpreviouswork
a a
[11]. However, thePpotential of interaction (11) corresponds completely to the conditions
at which the fluctuation formula (7) is derived. Therefore, this potential, instead of (8),
must be used in simulations to obtain a correct value for the dielectric constant.
3 Molecular reaction field
In the case of point dipoles, where d +0, q provided µ const, both (8)
a
→ → ∞ →
and (11) representations are identical and reduced to the well-known result
µ ·µ
ϕ = µ ·T(r )·µ i j , r R (12)
ij − i ij j − R3 ij ≤
for the interdipolar interaction in the RF geometry. It is easy to see that in the case of IS
models, the intermolecular potential (8) takes into account effects associated with finite-
ness of the molecule inconsistently. For example, the interdipolar potential is replaced by
the real site-site Coulomb ones, whereas the reaction field is remained in its usual form
of point dipoles. From this point of view a natural question of how to improve the RF
within the molecular cut-off scheme arises. The simplest way to solve this problem lies in
the following.
Let us consider the mentioned above spherical cavity, centered now at some fixed point
r , in the infinite conducting medium. We place an ith molecule in such a way that all
0
sites of the molecule would be located in the cavity. This condition is fulfilled providing
7
r r R R d/2. The potential of a molecular reaction field at point r belonging
| i− 0| ≤ d ≡ −
the cavity can be presented, according to the method of electrostatic images, as
q∗ q R/ρa q
ϕRF(r) = a = a i = a , (13)
i ρ ρ∗a − R 2 − ρa R
Xa | − i| Xa ρ− ρa ρai Xa Ri ρ− ρaρai
(cid:12) (cid:18) i (cid:19) (cid:12) (cid:12) i (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
where ρ = r r and ρa = ra r . Differentiating (13) over r at point r yields
− 0 i i − 0 0
∂ϕRF(r) µ ∂2ϕRF(r) qr0 ∂3ϕRF(r) gr0
i = i , i = i , i = i , ... (14)
∂r r −R3 ∂r∂r r −R5 ∂r∂r∂r r −R7
(cid:12) 0 (cid:12) 0 (cid:12) 0
(cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12)
Here µ = q ρa = q δa is the dipole moment of ith molecule, which does not
i a a i a a i
depend on rPowing elecPtroneutrality of the molecule, while qr0 = q (3ρaρa ρa2I)
0 i a a i i − i
and gr0 are the tensors of quadrupole and octupole moments, corPrespondingly, of ith
i
r
molecule with respect to r . The third rank tensor g 0 has the following components
0 i
gr0 = 3 q 5ρa ρa ρa ρa2(ρa δ +ρa δ +ρa δ ) . It is more convenient to
i αβγ a a iα iβ iγ − i iα βγ iβ αγ iγ αβ
(cid:16) (cid:17)
present muPltipoles of higher order with respect to the molecular center of mass. For the
tensor of quadrupole moment we obtain qr0 = q +w , where q = q (3δaδa δa2I)
i i i i a a i i − i
is the tensor of quadrupole moment of ith molecule with respect toPits center of mass,
w = 3(µ ρ +ρ µ ) 2µ ·ρ Iandρ = r r . Itisnecessary tounderlinethattensorq is
i i i i i − i i i i− 0 i
split into dynamical ω = q δaδa and conservative q δa2I parts for rigid molecules.
i a a i i a a i
Putting r = r and aPssuming d R, we obtainPthe energy of jth molecule in the
0 j ≪
MRF of ith molecule
∂ϕRF(r) 1 ∂2ϕRF(r) µ ·µ 1q :qrj
φRF = µ · i + q : i +... = j i j i +... , (15)
ji j ∂r r 6 j ∂r∂r r − R3 − 6 R5
(cid:12) j (cid:12) j
(cid:12) (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12)
where multipoles of higher order have been neglected. Finally, using the RF potential
ϕRF = (φRF +φRF)/2 yields the desired intermolecular potential
ij ij ji
q q µ ·µ q :q 3(q :µ r +q :µ r )
a b i j i j − i j ij j i ji , r R
ra rb − R3 − 6R5 ij ≤ d
ϕij =  Xa,b | i − j| 0 , r > R (16)
ij d

where equality q:I = 0 has been used.
The total reaction field, created by all molecules at point r near r is
0
8
ρi≤Rd ∂ϕRF(r) M(R ) Q(R )+W(R )
E (r) = i = d + d d ρ+... , (17)
RF − ∂r R3 R5
i
X
where M(R ) = ρi≤Rd µ and Q(R ) = ρi≤Rd q denote the total dipole and own
d i i d i i
quadrupolemomenPt, respectively, withinthesPphereofradiusR andW(R ) = ρi≤Rd w .
d d i i
In the case of point dipoles, we have R R, q ,g ,... 0 and the MRF (P17) trans-
d i i
→ →
forms into M(R)/R3 + W(R)ρ/R5. This last formula shows that the reaction field of
finite systems is inhomogeneous even for point dipoles. Only for macroscopic (R )
→ ∞
systems, we reproduce the well-known homogeneous reaction field M(R)/R3 introduced
by Barker and Watts [3]. For finite IS systems, additional higher multipole terms appear.
This brings, for example, into existence of the new quadrupole-dipole and quadrupole-
quadrupole interactions in the intermolecular potential (16). We note that the idea of
using the higher multipole moments in the RF has been proposed for the first time by
Friedman [5].
However, the modified intermolecular potential (16) still needs to be complemented
by a self-consistent fluctuation formula as this has already been done in the preceding
section by the fluctuation formula (7) for the potential of interaction in the site-site cut-
off scheme (11). Unfortunately, it is not a simple matter to construct fluctuation formulas
in the molecular cut-off approach. This problem will be considered in further studying.
The difference in the RF geometry between IS and PD models lies in the distinction
for their microscopic operators of polarization density. For IS models
Pˆ (k) = ik N e−ik·ri q e−ik·δai = Mˆ (k) ik kˆkˆ: N ω e−ik·ri +... , (18)
L k2 A L − 2 i
i=1 a i=1
X X X
where Mˆ (k) = kˆ N kˆ·µ e−ik·ri is the microscopic operator of polarization density for
L i=1 i
point dipoles and aPn expansion over small parameter k·δa has been made [15]. However,
i
putting Pˆ (k) Mˆ (k) in the microscopic electric field EˆRF(k) (6) at the very begin-
L ≡ L R
ning and taking attempts to perform the inverse Fourier transform, we obtain that the
corresponding integral is divergent in k-space when k . This divergence is involved
→ ∞
by the specific nature of point dipoles for which the parameter k·δa becomes indetermi-
i
nate in the limit k because of δa +0 and the expansion (18) fails. Therefore, we
→ ∞ i →
must manipulate with the full operator Pˆ (k) to obtain the interdipolar potential (12)
L
consequently and let δa +0 at the end of the calculation only.
i →
9
Since µ d and q d2, the quadrupole contribution with respect to the dipole term
∼ ∼
is varied in (16) from of order (d/R)2 at r = 0 to d/R at r = R . Therefore, as far
ij ij d
as the usual intermolecular potential (8) is applied in simulations, the dielectric constant
can not be reproduced with the precision better than d/R. It is evident that using
∼
the modified intermolecular potential (16) will lead to the uncertainties of order (d/R)2.
They decrease at increasing the size of the sample as R−2, i.e., with the same rate as those
connected with the truncation of the potential. Effects of the octupole and higher order
multipole contributions into the MRF are of order (d/R)3 and can be ignored.
4 Applying the ISRF to a MCY water model
In the previous investigations [11–13], the standard PDRF geometry (8) has been
applied to actual simulations of the MCY and TIP4P models. As a result, the static,
frequency-dependent [11, 12] and wavevector-dependent [13] dielectric constant has been
determined. For these models d = 1.837˚A and the cut-off radius R = 9.856˚A has been
used in the simulations. From the afore said in the preceding section, it is expected that
the precision of these calculations can not exceed d/R 20%. We shall show now by
∼
actual calculations that this prediction indeed takes place.
As an example we apply the ISRF geometry (11) to the MCY potential [16]. The
calculations have been performed with the help of Monte Carlo (MC) simulations, details
ofwhicharesimilartothosereportedearlier[13],atthedensity ofρ=1.0g/cm3 andatthe
temperature of T = 292 K, i.e., in the same thermodynamic point and yet with the same
number N = 256 of molecules and cut-off radius R = 9.856˚A as considered in [11, 13].
Our result of the calculation (7) for the longitudinal components of the wavevector-
dependent infinite-system Kirkwood factor g (k) and dielectric constant ε (k) obtained
L L
within the ISRF geometry is presented in Figs. 1 and 2, respectively, as the full circles
connected by the solid curves. For the purpose of comparison, analogous calculations
performed previously [13] within the PDRF are also included in these figures (the open
circles connected by the dashed curves). It is obvious that oscillations observing in the
shape of g (k) and ε (k) obtained within the PDRF method are nonphysical and caused
L L
by the finite molecular size which is assumed to be zero in this approach. At the same
10