Comparison of simulation methods for charged systems of slab geometry
Using the specific model of a system of like charged ions confined between
 two planar like charged surfaces, we compare the predictions for the energy
 and density profile of four simulation methods available to treat the long
 range Coulomb interactions in systems periodic...
Збережено в:
| Опубліковано в: : | Condensed Matter Physics |
|---|---|
| Дата: | 2001 |
| Автори: | , , , |
| Формат: | Стаття |
| Мова: | Англійська |
| Опубліковано: |
Інститут фізики конденсованих систем НАН України
2001
|
| Онлайн доступ: | https://nasplib.isofts.kiev.ua/handle/123456789/120533 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | Comparison of simulation methods for
 charged systems of slab geometry / M. Mazars, J.-M. Caillol, J.-J. Weis, D. Levesque // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 697-713. — Бібліогр.: 33 назв. — англ. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| _version_ | 1860259160540577792 |
|---|---|
| author | Mazars, M. Caillol, J.-M. Weis, J.-J. Levesque, D. |
| author_facet | Mazars, M. Caillol, J.-M. Weis, J.-J. Levesque, D. |
| citation_txt | Comparison of simulation methods for
 charged systems of slab geometry / M. Mazars, J.-M. Caillol, J.-J. Weis, D. Levesque // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 697-713. — Бібліогр.: 33 назв. — англ. |
| collection | DSpace DC |
| container_title | Condensed Matter Physics |
| description | Using the specific model of a system of like charged ions confined between
two planar like charged surfaces, we compare the predictions for the energy
and density profile of four simulation methods available to treat the long
range Coulomb interactions in systems periodic in two directions but bound
in the third one. Monte Carlo simulations demonstrate unambiguously complete
agreement between the results obtained with these methods where
the potential between charges is solution of Poisson’s equation in the simulation
cell with adequate boundary conditions. The practical advantages
of the different methods are assessed.
Користуючись специфічною моделлю систем однойменно заряджених іонів обмежених двома плоскими однойменно зарядженими поверхнями проведено порівняння енергії і профілю густини чотирьох
можливих методів комп’ютерного моделювання для врахування далекосяжних кулонівських взаємодій в системах періодичних за двома напрямками і обмежених за третім напрямком. Методи Монте-Карло показують повну однозначну узгодженість результатів отриманих цими методами, якщо потенціал між зарядами задовільняє
рівнянню Пуасона для симуляційної комірки з відповідними граничними умовами. Оцінено практичні переваги різних методів.
|
| first_indexed | 2025-12-07T18:53:03Z |
| format | Article |
| fulltext |
Condensed Matter Physics, 2001, Vol. 4, No. 4(28), pp. 697–713
Comparison of simulation methods for
charged systems of slab geometry∗
M.Mazars, J.-M.Caillol, J.-J.Weis, D.Levesque
Laboratoire de Physique Théorique,
UMR 8627, Bât. 210, Université de Paris-Sud
91405 Orsay Cedex, France
Received October 18, 2001
Using the specific model of a system of like charged ions confined between
two planar like charged surfaces, we compare the predictions for the ener-
gy and density profile of four simulation methods available to treat the long
range Coulomb interactions in systems periodic in two directions but bound
in the third one. Monte Carlo simulations demonstrate unambiguously com-
plete agreement between the results obtained with these methods where
the potential between charges is solution of Poisson’s equation in the sim-
ulation cell with adequate boundary conditions. The practical advantages
of the different methods are assessed.
Key words: simulations, Monte Carlo, Coulomb interaction, colloids
PACS: 61.20.Ja, 68.15.+e, 82.70.-Dd
1. Introduction
The standard way to treat long range interactions (Coulomb, dipolar, Yukawa
near the Coulomb limit) in simulations of bulk systems is to replicate the basic simu-
lation cell periodically in all three directions of space and apply an Ewald summation
technique (EW3D) [1,2]. However, in many situations of electrochemical, biological
or technological interest, as, for instance, electrolyte solutions between charged sur-
faces, charged lipid bilayers in water, suspensions of colloids between glass plates,
clays [3–5], magnetic thin films [6], Wigner crystals [7] etc., the system is finite in
one direction thus necessitating a modification of the conventional Ewald method.
For systems periodic in two dimensions but bounded along the third one, an Ewald
summation method (EW2D) for electrostatic interactions has first been given by
Parry [8,9] and was later rederived by various authors [10–14]. Unfortunately, the
practical use of the EW2D sum is hampered by the occurrence, in the reciprocal
space term, of a double sum over different particles which, due to the complicated
∗We dedicate the article to our colleague J.-P.Badiali.
c© M.Mazars, J.-M.Caillol, J.-J.Weis, D.Levesque 697
M.Mazars et al.
way the bounded coordinates enter the expression, cannot be reduced, as for the
EW3D case, to a sum of order N, a circumstance which renders the method com-
putationally expensive. Not surprisingly, only few calculations using EW2D have
been reported to date mainly to test the validity of more approximate approaches
[15,16]. These calculations involve tabulation and interpolation of the potentials on
a three-dimensional grid.
Various “time saving” schemes have been proposed to bypass the computational
burden of the EW2D method. The purpose of this paper is to compare some of
these methods for the specific case of a system of charged particles confined between
two charged planar surfaces separated by a distance of the order of several particle
diameters. The methods that have been chosen are those which satisfy exactly the
laws of electrostatics (perfect screening, Green’s theorem etc..) and do not present
numerical instabilities difficult to control. Thus the charged sheet method proposed
by Torrie and Valleau [17,18] and its modification by Boda et al [19] have been
discarded. Indeed, the distribution of the point charges of the ions outside the sim-
ulation cell, located in the image cells of the latter, is represented approximately
by a set of uniformly charged planar sheets. Another method, proposed by Lekner
[20,21] is based on rewriting the sum of forces acting on an ion by a second ion and
its periodic images as an absolutely converging infinite sum. The energy is obtained
by integration of the force. A shortcoming of the method is the need to estimate the
sum with a given precision, to retain a number of terms depending strongly on the
relative position of pairs of particles. Because of this technical reason, emphasized
in more detail in [22], the method was not included in this investigation.
The four methods we have considered, the first three of which have already
appeared in the literature, are EW3D [15,16,23,24], a method developed by Hautman
and Klein [25] (HK), the method of hyperspheres (HSG) [26] and that of concentric
spheres (CS) described below.
A way to treat the long-range Coulomb interactions is to place the slab at the
centre of a parallelepipedic simulation box which has dimension perpendicular to
the slab much larger than the width of the latter and apply the EW3D method.
In this way one effectively simulates an infinite number of parallel slabs. Provided
the region of empty space separating the slabs is sufficiently large, one expects the
influence of the periodic images on the behaviour of the system to become negligible.
If the distance between the surfaces confining the particles is small compared
to the extension of the system along the directions parallel to the surfaces (narrow
slab), the lateral distance s between pairs of particles will be much larger than the
distance z normal to the slab surfaces and it will be appropriate to expand the
Coulomb pair interaction in powers of z/s and apply an Ewald sum to the in-plane
component of the interaction. Such an approach was followed by Hautman and Klein
(HK) [25].
An attractive way to investigate the behaviour of Coulomb particles between
charged surfaces which circumvents the cumbersome Ewald sums is to use a closed
space, e.g. a hypersphere in four-dimensional Euclidean space [26–28]. Not only is the
electrostatic potential solution of Poisson’s equation on the hypersphere obtainable
698
Simulation of confined charged systems
in simple closed form [27], thus avoiding approximation of the interaction potential
(as, for instance, in EW3D due to truncation of the direct and reciprocal space sums),
but also the number of operations necessary to calculate the distances between
particles is reduced with respect to Euclidean space with concomitant speed up of
the simulation procedure. This geometry has been applied previously to the study of
the attraction between two like charged surfaces neutralized by solvated counterions
in an endeavour to comprehend the stability of charged lamellar materials such as
clays and cement [27,29] and to the study of the effective interaction of charged
colloidal particles confined between like charged plates [30].
A system of charged particles occupying the region between the outer and the
inner surfaces of two concentric spheres is also a suitable arrangement to describe,
in the limit of sufficiently large radii of the spheres, the system of charged particles
confined between two planar surfaces. This geometry has the advantage that the
interaction between charges and between charges and surfaces is the usual Coulomb
potential; however, it may require the use of a large number of particles to render
the effect of curvature of the confining surfaces small.
The purpose of this paper is to check, by means of Monte Carlo simulation, the
agreement between the results of the four aforementioned methods for the energy
and density profile of the system described above consisting of N like-charged ions
between parallel like-charged surfaces separated by a distance h. The ions, modelled
as hard spheres of diameter d bear a charge q at their centre and each surface, of
area S, a uniform charge density σ. Such a system can be viewed as a crude mod-
el for lamellar liquid crystals formed by ionic amphiphiles [31] or charged lamellar
materials like clays or cement [27]. The convergence of the results to their thermo-
dynamic limit is estimated by performing simulations with an increasing number of
particles keeping the distance between surfaces fixed. The speed of this convergence
is obviously an important criterion for the appraisal of the relative practical interest
of the four methods.
Expressions for the energy calculated with the different methods will be given in
section 2 together with details of the Monte Carlo simulations. Results for the energy
and density profiles obtained with the four methods are compared in section 3 and
conclusions drawn in section 4.
2. Energy expressions
2.1. EW3D
In this method a square slab of ions of thickness h is placed at the centre of a
simulation box having dimension Lz normal to the surfaces much larger than the
lateral dimension L and the system is extended periodically in the three directions
of space. The slab surfaces are located at z = ±h/2 perpendicular to the z-axis.
The closest approaches of the ions to the impenetrable surfaces are therefore z =
±1
2
(h− d).
To evaluate the total energy of the system which is the sum of the contributions
699
M.Mazars et al.
from the ion-ion Uii, ion-surface Uis and surface-surface Uss interactions it will be
convenient to associate with the ions a uniform neutralizing background (filling the
whole simulation cell) of density −Nq/V and with the charged walls a uniform
background of density −2σS/V = −2σ/Lz where V = LzL
2 = LzS is the volume
of the simulation cell. Because of electroneutrality of the system Nq+2σS = 0, the
backgrounds cancel out exactly and will not contribute to the total energy.
The contribution of the ions and their background to the total energy is given
by
Uii =
q2
2
{
N
∑
i,j
∑
ν
′ erfc(α|rij + ν|)
|rij + ν| +
1
V
N
∑
i,j
j 6=i
∑
G6=0
4π2
G
2 e
−G
2/4α2
eiG·(rj−ri)
− 1
V
N
∑
i,j
j 6=i
π
α2
}
− 1
2
Cw . (1)
In these equations rij=rj−ri where ri and rj are the positions of particles i and j.
The term − 1
2
Cw is the self-energy given by
−1
2
Cw = −1
2
q2N
{
−
∑
ν 6=0
erfc(α|ν|)
|ν| − 1
V
∑
G6=0
4π2
G
2 e
−G
2/4α2
+
2α√
π
+
π
α2
}
. (2)
In equation (1) ν is a vector of components (Lnx,Lny,Lznz) (nx, ny, nz integers)
and the prime in the sum over ν means that the terms i = j must be omitted when
ν = 0. The wave-vectors G which enter the reciprocal space contributions to the
energy have components 2πnx/L, 2πny/L and 2πnz/Lz .
In our calculations the sum on lattice vectors extends over all G subject to
|nx| 6 6, |ny| 6 6, |nz| 6 12 and |n| 6 nmax = 12. The α parameter which
governs the rate of convergence of the real- and reciprocal-space contributions was
taken sufficiently large so that only the terms with ν = 0 had to be retained in
equations (1) and (2).
The electrostatic energy between the ions and the two charged walls including
their compensating backgrounds is given by
Uis =
N
∑
i=1
qVs(zi)−
Nq
V
∫
V
drVs(z), (3)
where Vs(z) is the electric potential associated with the electric field E s generated
by the two charged surfaces and their backgrounds. Due to the symmetry of the
simulation cell, this electric field is normal to the charged surfaces and depends
only on the z-coordinate. Furthermore, as a consequence of the periodic boundary
conditions, it must satisfy Es(−Lz/2) = Es(Lz/2) = 0. Applying Gauss’s theorem
700
Simulation of confined charged systems
✁❆
✟❍
z
y,x
r
0
−
L
z
2
−
h
2
h
2
L
z
2
Figure 1. Simulation setup for EW3D. The slab of width h is placed at the
center of a simulation cell having dimension Lz in the direction normal to the slab
larger than the lateral dimensions L. Periodic boundary conditions are applied
in all three directions. To evaluate the electric field Es(z) created by the charged
surfaces and their backgrounds, Gauss’s theorem is applied to the dotted volume.
(cf. figure 1) one finds
Es(z) =
−8πσ
Lz
z − 4πσ, −Lz/2 6 z 6 −h/2
−8πσ
Lz
z, −h/2 6 z 6 h/2
−8πσ
Lz
z + 4πσ, h/2 6 z 6 Lz/2
(4)
and, by integration, for the potential, choosing Vs(0) = 0 (any additional constant
would leave Uis unchanged),
Vs(z) =
−4πσ
Lz
z2 + 4πσz + 2πσh, −Lz/2 6 z 6 −h/2
−4πσ
Lz
z2, −h/2 6 z 6 h/2
−4πσ
Lz
z2 − 4πσz + 2πσh, h/2 6 z 6 Lz/2.
(5)
In the periodic system the potential created by the charged surfaces is thus parabol-
ic. The potential Vs(z) can also be obtained by a direct integration of the Ewald
potential over the two charged planes resulting in a one-dimensional Fourier series
which can be summed explicitly to yield equation (5). Integration of V s over the
volume of the simulation box leads to
Uis = −4πσq
Lz
N
∑
i=1
z2i −
Nqπσ
3Lz
(3h2 − 6Lzh+ 2L2
z). (6)
701
M.Mazars et al.
Finally, the interaction between the two charged surfaces (including their back-
ground) is
Uss =
1
8π
∫
V
drE2
s (z) =
2πSσ2
3L2
z
[h3 + (Lz − h)3]. (7)
2.2. Hautmann-Klein method
In this method the slab, having the characteristics described above, has periodic
boundary conditions only in the x- and y- directions. The origin of coordinates being
at the centre of the parallelepipedic cell, the energy of the system is given by
U =
q2
2
N
∑
i,j=1
∑
ν
′ 1
|rij + ν| + σq
N
∑
i=1
2
∑
α=1
∫
S
dxαdyα
∑
ν
1
|ri − rα + ν|
+
1
2
σ2
2
∑
α,β=1
∫
S
dx′
αdy
′
α
∫
S
dxβdyβ
∑
ν
1
|r′
α − rβ + ν| , (8)
where rij = rj−ri, ri,rj coordinates of the particles i and j, and the vectors rα and
rβ are the position vectors of points on the surfaces Sα (α = 1, 2) having constant
z coordinate z = h/2 if α or β = 1 and z = −h/2 if α or β = 2. The sum over ν
runs over all vectors of components (nxL, nyL) (nx, ny integers) perpendicular to the
z-direction. The three sums in the expression for U correspond to the interactions
between the particles, the interaction between the particles and the surfaces S1 and
S2 and to that between the two surfaces S1 and S2, respectively. To each of these
sums, due to the infinite number of replicated cells, is associated a divergent con-
tribution. However, if electroneutrality is taken into account these divergent terms
will cancel.
The evaluation of the ion-ion interaction starts with the identity [25]
1
r
=
(
1
r
−
M
∑
n=0
an
z2n
s2n+1
)
+
M
∑
n=0
an
z2n
s2n+1
(9)
with
an =
(−1)n(2n)!
22n(n!)2
.
The added and subtracted term represents the first M + 1 terms in the binomial
expansion of 1/r in powers of z/s where s is the component of r in the plane of the
surface. By introducing a convergence function hn(s;α) for each term 1/s2n+1 the
energy can be divided in a short range part
Us
ii =
1
2
N
∑
i=1
N
∑
j=1
qiqj
(
∑
ν
′ 1
rij,ν
−
M
∑
n=0
anz
2n
ij
hn(sij,ν ;α)
s2n+1
ij,ν
)
(10)
702
Simulation of confined charged systems
and a long range part
U l
ii =
1
2
N
∑
i=1
N
∑
j=1
qiqj
M
∑
n=0
anz
2n
ij
(
∑
ν
′ hn(sij,ν ;α)
s2n+1
ij,ν
)
(11)
which can be evaluated in reciprocal space. In these equations sij,ν = |sj − si + ν|.
With the choice
h0(s;α) = erf(αs) (12)
and
hn(s;α)
s2n+1
=
1
an(2n)!
∇2n
(
h0(s;α)
s
)
(13)
as proposed by Hautman and Klein the short and long ranged parts of the electro-
static energy take the form [25]
Us
ii =
1
2
N
∑
i=1
N
∑
j 6=i
qiqj
(
1
rij
− erf(αsij)
sij
−
M
∑
n=1
1
(2n)!
z2nij ∇2n(
erf(αsij)
sij
)
)
, (14)
where only the ν = 0 term has been retained (see below) and
U l
ii =
π
A
N
∑
i=1
N
∑
j=1
qiqj
M
∑
n=0
1
(2n)!
z2nij
∑
G6=0
G2n−1erfc(G/2α)eiG.sij
− α√
π
N
∑
i=1
q2i −
√
π
αA
(
N
∑
i=1
qi)
2
(15)
with G a two-dimensional vector of components 2π(nx/L, ny/L). In our simulations
we kept terms up to M = 3 which allowed the short range part U s
ii to be limited to
its ν = 0 term even for relatively large values of h. Expressions for hn(s;α) up to
M = 3 are
h1(s;α) = erf(αs)− 2αs√
π
e−α2s2(1 + 2α2s2),
h2(s;α) = erf(αs)− 2αs√
π
e−α2s2(1 +
2
3
α2s2 − 4
9
α4s4 +
8
9
α6s6),
h3(s;α) = erf(αs)− 2αs√
π
e−α2s2(1 +
2
3
α2s2 +
4
15
α4s4 +
8
25
α6s6
− 112
225
α8s8 +
32
225
α10s10).
(16)
The attractive feature of the method is that, by writing the z2n
ij explicitly as
polynomials in zi and zj , U
l
ii becomes a sum of terms each of which involves products
of two functions having general form
Fp(G) =
N
∑
i
qiz
p
i e
iG.si. (17)
703
M.Mazars et al.
In the evaluation of U l
ii the sum on the pairs of particles is replaced by sums on
particles which, obviously, makes the computation faster. The contributions to the
energy of the interactions between the ions and the surfaces and between the surfaces
are easily obtained using the method of de Leeuw and Perram [11] and the identity
∫
S
dxαdyα
∑
ν
f(|ra − r + ν|) ≡
∫ +∞
−∞
dxα
∫ +∞
−∞
dyα f(|ra − r|). (18)
The divergent terms of these two contributions having been eliminated through the
use of the electroneutrality condition, the ion-surface and surface-surface interactions
contribute to the energy U by a constant term equal to 2πσ2V .The expression for
the total energy is therefore given by
U = Us
ii + U l
ii + 2πσ2V. (19)
2.3. Hyperspherical geometry
In this method the Monte Carlo simulations are performed on a hypersphere
S3 in four-dimensional Euclidian space. Two surfaces of angular colatitudes θN and
θS = π − θN separated by a distance h = R(π − 2θN) are located symmetrically on
the opposite sides of the equator (see figure 3 of [27]). N neutralizing ions of charge
q are confined between the two surfaces which bear each a charge density σ. For
given h and σ the colatitude θN is obtained by solving the equation
h sin θN
π − 2θN
=
(
− qN
8πσ
)1/2
. (20)
Similar to the EW3D case described above, neutralizing backgrounds are as-
sociated to both the ions and the charged surfaces. As a detailed derivation of
the different contributions to the potential energy, ion-ion Uii, ion-surface Uis and
surface-surface Uss can be found in [27] only the final expressions are given here
Uii =
q2
2πR
N
∑
i=1
N
∑
j 6=i
[(π − θij)cotθij − 0.5]
+
Nq2
2πR
[(π − 2α)cotα + 0.5− π/sinα], (21)
Uis =
qqs
πR
[(π − 2θi)cotθi − 1] +
4q2s
πR
(θNcotθN − 1), (22)
Uss =
q2s
πR
[1 + (π − 4θN)cotθN ]. (23)
Here the angle α is equal to the ion radius d/2 divided by the hypersphere radius
R, θij is the angular separation between particles i and j and qs = σS (S area of
each surface).
704
Simulation of confined charged systems
2.4. Concentric spheres
In the last method we have explored, the ions occupied the region of volume V
between two concentric spheres of radii rl and rm = h+ rl. The external surface of
the inner sphere, of area Sl = 4πr2l , and the internal surface of the outer sphere, of
area Sm = 4πr2m, bear a charge density σ. The two surfaces being impenetrable, the
distances of the closest approach of an ion to the surfaces are rl+ d/2 and rm−d/2.
The electroneutrality conditions read
σ(Sl + Sm) +Nq = 0. (24)
The total energy of the system is readily obtained as
U =
1
2
N
∑
i=1
N
∑
j 6=i
q2
|ri − rj |
+
N
∑
i=1
Slσq
ri
+ qN
σSm
rm
+
σ2SmSl
rm
+
1
2
σ2S2
m
rm
+
1
2
σ2S2
l
rl
, (25)
where the three first terms in the r.h.s. represent the ion-ion energy, the energy
of the ions with the charged surface Sl and with the surface Sm which creates in
the volume V a constant potential. The three last terms correspond to the energies
associated with the surface charges of Sl and Sm.
3. Results
In our comparisons we did not consider realistic values of q and σ as for instance
envisaged in [27,29]. Our aim being methodological, we chose values of q and σ
allowing for a compelling test of the efficiency of the four methods to predict reliably
the properties of the confined charged system. In the following charges q and σ
will be expressed by using reduced units of charge and distance (βq2/d)1/2 and
d, respectively. In these units the surface charge σ has been fixed to −1 in all our
simulations and the value of q varied from 0.5 to 5. The distance between the surfaces
limiting the system has been taken to be h = 5; however, in order to test the limit of
validity of the HK method for large h a few simulations were performed with h = 8
and h = 12. For given values of h and σ electroneutrality entails that the density of
the system increases for decreasing values of q if the volume of the simulation cell
remains unchanged.
The simulations were carried out in the canonical ensemble. Of the order of
105 MC steps per particle were performed to calculate the energy and the density
profiles when N 6 3000. In the CS geometry, for N ≈ 15000, 104 trial configurations
per particle were generated to calculate these quantities. The statistical error on the
energies is of the order of ±0.02 in units of kT . For the two planar geometries, EW3D
and HK, the density profiles are estimated with an error of ∼ 0.5%; the statistical
error, except for systematic error due to curvature effects, is of the same order for
the HSG method. In the case of the CS geometry, for N ≈ 15000, the error on the
energies is of the order of ±0.03 and on the density profiles of the order of 2− 3%.
In the calculations using the EW3D method, the simulation volume has a square
section of side L and an elongation along the z-axis of Lz varying between 60 and
705
M.Mazars et al.
Table 1. Variation with ionic charge q of the total energy of a system of like ions
confined between planar surfaces bearing a charge density σ = −1 and separated
by a distance h = 5, for the four different simulation methods considered. The
number N of ions is given in brackets, β = 1/kT , T temperature.
q βUEW3D/N βUHK/N βUHSG/N βUCS/N
5.0 –4.30 (320) –4.31 (1024) –4.26 (1024) –4.16 (1024)
–4.30 (980) –4.29 (3000) –4.26 (3544)
–4.28 (14036)
4.0 –1.30 (320) –1.28 (1024) –1.26 (1024)
–1.28 (500) –1.27 (3000)
3.0 0.930 (500) 0.943 (1024) 0.948 (1024)
0.946 (3000)
2.0 2.32 (500) 2.35 (1024) 2.36 (1024) 2.33 (1024)
2.35 (980) 2.35 (3000) 2.35 (3544)
2.36 (14036)
1.0 3.32 (980) 3.32 (1024) 3.23 (1024) 3.51 (1024)
3.29 (3000) 3.36 (3544)
3.30 (14036)
0.75 3.13 (500) 3.15 (1024) 3.02 (1024) 3.57 (1024)
3.17 (980) 3.10 (3000) 3.26 (3544)
3.16 (1620) 3.19 (14036)
90. For h = 5 and L between 15 and 30, the influence of the value of Lz on the
simulation results turned out to be always negligible. Also, within this geometry, if
the system confined between the surfaces has a net dipole moment, a correction to
the energy proportional to the square of the dipole moment normal to the surfaces
can be taken into account to remove the interaction between net dipoles in the
periodically repeated slabs [16]. As the dipole moment in our system was small such
a correction was not considered.
A systematic comparison of the energies for 0.75 6 q 6 5 is made in table 1.
It shows without ambiguity the convergence of the results of the four methods. It
appears that the thermodynamic limit for U is obtained for q > 1 as soon as the
number of particles is of the order of 500 for EW3D and HK and of the order of 2000–
3000 for the HSG and CS method. However, at q ≈ 0.75, i. e. at densities larger than
0.4 at fixed h and σ, at least 10000 particles are necessary for the energy calculated
with CS to agree within 1% with the other methods. This result is obviously related
to the important curvature effects expected for radii rl and rm 6 10. Figure 2 shows,
for the four methods, the variation of U when the ionic charge q varies from 0.5 to 5.
A maximum is observed for q ≈ 0.8. This maximum occurs as a consequence of the
706
Simulation of confined charged systems
0.0 1.0 2.0 3.0 4.0 5.0
q
−6.0
−4.0
−2.0
0.0
2.0
4.0
β
U
/N
EW3D
HK
HSG
CS
Figure 2. Variation with ionic charge q of the total energy of a system of like ions
confined between planar surfaces bearing a charge density σ = −1 and separated
by a distance h = 5, for the four different simulation methods considered.
−2.0 −1.0 0.0 1.0 2.0
z/d
0.0
1.0
2.0
3.0
4.0
5.0
6.0
ρ
(z
)
q = 5.0
q = 4.0
q = 3.0
q = 2.5
q = 2.0
q = 1.5
q = 1.0
q = 0.75
q = 0.625
q = 0.5
Figure 3. Variation with ionic charge q of the density profile ρ(z) of a system of
like ions confined between planar surfaces bearing a charge density σ = −1 and
separated by a distance h = 5. The results shown correspond to the method of
HK and are indistinguishable from those obtained for EW3D.
707
M.Mazars et al.
1.00 1.25 1.50 1.75 2.00
z/d
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
ρ
(z
)
EW3D
HSG
HSG
CS
HK
Figure 4. Ion density profile for a slab of width h = 5, surface charge density
σ = −1 and ionic charge q = 2. Comparison between all four simulation methods.
The profile is shown only for z/d > 1.
lowering of q and saturation of the particles in contact with the confining surfaces
due to steric effects. Indeed layers of particles parallel to the surfaces form at high
density as evidenced by figure 3 which presents the density profiles for the different
values of q. For large q, layers of highly localized particles are in contact with the
surfaces. When the density increases (i.e. q decreases) layers separated by a distance
of order d appear in the volume. The profiles shown in figure 3 are obtained by the
method of HK and are indistinguishable from those given by EW3D.
Figures 4 and 5 allow to compare the density profiles obtained with the methods
EW3D and HK with those evaluated with the methods of hypersphere and concentric
spheres. For q = 2 excellent agreement is found for the four methods taking into
account the statistical errors, especially those for the CS method. At q = 0.75
differences are manifest and demonstrate the persistence of curvature effects; they
are small in the hypersphere geometry but remain important in the CS geometry in
spite of the large value N = 14036. This result at q = 0.75 is in agreement with the
708
Simulation of confined charged systems
−0.5 0.0 0.5 1.0 1.5 2.0
z/d
0.0
0.5
1.0
1.5
2.0
ρ
(z
)
EW3D
HSG
HSG
CS
HK
Figure 5. Ion density profile for a slab of width h = 5, surface charge density
σ = −1 and ionic charge q = 0.75. Comparison between all four simulation
methods. The profile is shown only for z/d > −0.5.
−6.0 −4.0 −2.0 0.0 2.0 4.0 6.0
z/d
0.0
1.0
2.0
3.0
4.0
5.0
6.0
ρ
(z
)
HK
EW3D
Figure 6. Comparison of the density profiles obtained with the HK and EW3D
methods for (from top to bottom) h = 5, σ = −1; h = 8, σ = −1 and h = 12,
σ = −2. All results are for q = 1. For clarity the profiles have been shifted by 0.5
with respect to each other.
709
M.Mazars et al.
Table 2. Comparison of energy U and contact value ρc of the density profile for
the two methods HK and EW3D. The surface charge density is σ = −1 except for
the case h = 12 where σ = −2, β = 1/kT , T temperature. The average density
is ρ = N/V = −2σ/qh.
Ewald 3d Hautman-Klein
h q ρ βU/N ρc βU/N ρc
5.0 5.0 0.08 –4.30 6.30 –4.31 6.38
4.0 0.1 –1.28 6.30 –1.28 6.32
3.0 2/15 0.93 6.30 0.94 6.31
2.5 0.16 1.70 6.30 1.71 6.29
2.0 0.2 2.35 6.30 2.35 6.29
1.5 4/15 2.97 6.35 2.98 6.33
1.0 0.4 3.32 6.53 3.32 6.54
0.75 8/15 3.17 7.22 3.15 7.20
0.625 0.64 2.90 8.36 2.90 8.38
0.5 0.8 2.49 11.6 2.49 11.6
8.0 1.0 0.25 3.46 6.40 3.53 6.39
12.0a 1.0 1/3 9.07 27.4 8.92 29.3
aσ = −2
slow convergence, in the latter geometry, of the energy to its thermodynamic limit
when q < 1.
Extrapolations towards z = ±2 give the values ρc of the profiles when the par-
ticles are in contact with the surfaces. It has been shown that for planar surfaces
separated by a distance much larger than the diameter d of the particles the pressure
is given by [32,33]
P/kT = ρc − E2
c/8πkT, (26)
where Ec, the value of the electric field at the surface, equals 4πσ for planar surfaces.
For the confined systems considered in this work this expression is not strictly valid
anymore. Nonetheless, considering it as an approximation and using the values of
ρc given in table 2, an estimate of P can be obtained. The value of P is small and
positive for q > 1 and rises to a value of 5 at q = 0.5 where hard core interactions are
important. The nearly zero value of the pressure for q > 1 is clearly in accord with
the absence of particles in the central part of the simulation cell (cf. figure 2) and
complete screening of the charged surfaces by the particles in their vicinity implying
710
Simulation of confined charged systems
an almost vanishing electric field in this same region of simulation volume.
Table 2 and figure 6 comparing simulation results for the energy and density
profiles for h = 8 and 12 obtained with the EW3D and HK methods, show that the
latter method remains accurate even for h/L ∼ 0.5. However, reliable results are
obtained only if, in equations 14 and 15 for the energy, M is chosen to be 3.
4. Conclusion
The convergence of the results obtained with the four simulation methods con-
sidered is clearly the main conclusion of this study. It demonstrates the strict equiv-
alence of the various possible routes to take into account the long range of the
Coulomb interactions: the use of periodic boundary conditions, a closed system
without boundary or simply a large system with boundaries. The equivalence is
made realized only by the use of potentials preserving the laws of electrostatics and
thus solution of Poisson’s equation associated with the simulation cell and its bound-
ary conditions. One can remark that the system under investigation is particularly
challenging to establish this equivalence because particles all bear a charge of the
same sign and screening effects therefore result only from interactions between the
particles and the surfaces.
Depending on the type of the system which is simulated, the advantage of using
one method or the other may differ. This work clearly shows that the EW3D and
HK methods apply more favourably to the study of systems confined by planar
surfaces, the thermodynamic limit being reached for N ∼ 500. However, with regard
to computing time, the method of hypersphere is the most efficient despite the
necessity to use systems with larger values of N to make curvature effects negligible.
Obviously, the CS method is unfavourable to study planar interfaces as curvature
effects get small only for N > 10000 − 20000. It has been used here mainly to
establish without ambiguity the equivalence between the Ewald and hypersphere
“Coulomb” potentials and the usual Coulomb potential.
The calculations carried out with the CS method nevertheless show that the ionic
density profiles are affected in the vicinity of a charged interface by the curvature of
the latter. Such modifications of the density profiles near planar charged interfaces
may be present in solutions of inverted micelles or suspensions of charged colloids
of small diameter.
The determination of the relative efficiency of the simulation methods for the
calculation of the properties of confined systems will be complete only when pres-
sure, free energy and surface tension will have been evaluated. Calculation of these
quantities is in progress.
References
1. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. Oxford, Clarendon, 1987.
2. Salin G., Caillol J.-M. // J. Chem. Phys., 2000, vol. 113, p. 10459.
3. Israelachvili J. Intermolecular & Surface Forces. San Diego, Academic, 1992.
711
M.Mazars et al.
4. Evans D.F., Wennerström H. The Colloidal Domain. New York, VCH, 1994.
5. von Grünberg H.H., Mbamala E.C. // J. Phys.: Condens. Matter, 2000, vol. 12,
p. 10349.
6. De’Bell K., MacIsaac A.B., Whitehead J.P. // Rev. Mod. Phys., 2000, vol. 72, p. 225.
7. Weis J.-J., Levesque D., Jorge S. // Phys. Rev. B, 2001, vol. 63, p. 045308.
8. Parry D.E. // Surf. Sci., 1975, vol. 49, p. 433.
9. Parry D.E. // Surf. Sci., 1975, vol. 54, p. 195.
10. Heyes D.M., Barber M., Clarke J.H.R. // J. Chem. Soc. Faraday Trans. II, 1977,
vol. 73, p. 1485.
11. de Leeuw S.W., Perram J.W. // Mol. Phys., 1979, vol. 37, p. 1313.
12. Cichocki S.W., Felderhof B.U. // Mol. Phys., 1989, vol. 67, p. 1373.
13. Harris F.E. // Int. J. Quantum Chem., 1998, vol. 68, p. 385.
14. Grzybowski A., Gwóźdź, Bródka A. // Phys. Rev. B, 2000, vol. 61, p. 6706.
15. Spohr E. // J. Chem. Phys., 1997, vol. 107, p. 6342.
16. Yeh I.-Ch, Berkowitz M.L. // J. Chem. Phys., 1999, vol. 111, p. 3155.
17. Torrie G.M., Valleau J.P. // J. Chem. Phys., 1980, vol. 73, p. 5807.
18. Valleau J.P., Ivkov R., Torrie G.M. // J. Chem. Phys., 1991, vol. 95, p. 520.
19. Boda D., Chan K.-Y., Henderson D. // J. Chem. Phys., 1998, vol. 109, p. 7362.
20. Lekner J. // Physica A, 1989, vol. 157, p. 826.
21. Lekner J. // Physica A, 1991, vol. 176, p. 485.
22. Mazars M. // J. Chem. Phys., 2001, vol. 115, p. 2955.
23. Shelley J.C., Patey G.N. // Mol. Phys., 1996, vol. 88, p. 385.
24. Crozier P.S., Rowley R.L., Spohr E., Henderson D. // J. Chem. Phys., 2000, vol. 112,
p. 9253.
25. Hautman J., Klein M.L. // Mol. Phys., 1992, vol. 75, p. 379.
26. Caillol J.-M., Levesque D. // J. Chem. Phys., 1991, vol. 94, p. 597.
27. Pellenq R.J.-M., Caillol J.-M., Delville A. // J. Phys. Chem. B, 1997, vol. 101, p. 8584.
28. Caillol J.-M. // J. Chem. Phys., 1999, vol. 111, p. 6528.
29. Delville A., Pellenq R.J.-M., Caillol J.-M. // J. Chem. Phys., 1997, vol. 106, p. 7275.
30. Allahyarov E., D’Amico I., Löwen H. // Phys. Rev. E, 1999, vol. 60, p. 3199.
31. Guldbrand P.G., Jönsson B., Wennerström H., Linse P. // J. Chem. Phys., 1984,
vol. 80, p. 2221.
32. Henderson D., Blum L., Lebowitz J.L. // J. Electroanal. Chem., 1979, vol. 102, p. 315.
33. Carnie S.L., Chan D.Y.C. // J. Chem. Phys., 1981, vol. 74, p. 1293.
712
Simulation of confined charged systems
Порівнення методів комп’ютерного моделювання
для заряджених систем, обмежених плоскими
поверхнями
М.Мазар, Ж.-М. Каіло, Ж.-Ж.Вейс, Д.Левек
Лабораторія теоретичної фізики, Університет Парі-Сюд,
91405 Орсей, Франція
Отримано 18 жовтня 2001 р.
Користуючись специфічною моделлю систем однойменно зарядже-
них іонів обмежених двома плоскими однойменно зарядженими по-
верхнями проведено порівняння енергії і профілю густини чотирьох
можливих методів комп’ютерного моделювання для врахування да-
лекосяжних кулонівських взаємодій в системах періодичних за дво-
ма напрямками і обмежених за третім напрямком. Методи Монте-
Карло показують повну однозначну узгодженість результатів отри-
маних цими методами, якщо потенціал між зарядами задовільняє
рівнянню Пуасона для симуляційної комірки з відповідними гранич-
ними умовами. Оцінено практичні переваги різних методів.
Ключові слова: моделювання, Монте Карло, кулонівські взаємодія,
колоїди
PACS: 61.20.Ja, 68.15.+e, 82.70.-Dd
713
714
|
| id | nasplib_isofts_kiev_ua-123456789-120533 |
| institution | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| issn | 1607-324X |
| language | English |
| last_indexed | 2025-12-07T18:53:03Z |
| publishDate | 2001 |
| publisher | Інститут фізики конденсованих систем НАН України |
| record_format | dspace |
| spelling | Mazars, M. Caillol, J.-M. Weis, J.-J. Levesque, D. 2017-06-12T10:08:47Z 2017-06-12T10:08:47Z 2001 Comparison of simulation methods for
 charged systems of slab geometry / M. Mazars, J.-M. Caillol, J.-J. Weis, D. Levesque // Condensed Matter Physics. — 2001. — Т. 4, № 4(28). — С. 697-713. — Бібліогр.: 33 назв. — англ. 1607-324X PACS: 61.20.Ja, 68.15.+e, 82.70.-Dd DOI:10.5488/CMP.4.4.697 https://nasplib.isofts.kiev.ua/handle/123456789/120533 Using the specific model of a system of like charged ions confined between
 two planar like charged surfaces, we compare the predictions for the energy
 and density profile of four simulation methods available to treat the long
 range Coulomb interactions in systems periodic in two directions but bound
 in the third one. Monte Carlo simulations demonstrate unambiguously complete
 agreement between the results obtained with these methods where
 the potential between charges is solution of Poisson’s equation in the simulation
 cell with adequate boundary conditions. The practical advantages
 of the different methods are assessed. Користуючись специфічною моделлю систем однойменно заряджених іонів обмежених двома плоскими однойменно зарядженими поверхнями проведено порівняння енергії і профілю густини чотирьох
 можливих методів комп’ютерного моделювання для врахування далекосяжних кулонівських взаємодій в системах періодичних за двома напрямками і обмежених за третім напрямком. Методи Монте-Карло показують повну однозначну узгодженість результатів отриманих цими методами, якщо потенціал між зарядами задовільняє
 рівнянню Пуасона для симуляційної комірки з відповідними граничними умовами. Оцінено практичні переваги різних методів. en Інститут фізики конденсованих систем НАН України Condensed Matter Physics Comparison of simulation methods for charged systems of slab geometry Порівнення методів комп’ютерного моделювання для заряджених систем, обмежених плоскими поверхнями Article published earlier |
| spellingShingle | Comparison of simulation methods for charged systems of slab geometry Mazars, M. Caillol, J.-M. Weis, J.-J. Levesque, D. |
| title | Comparison of simulation methods for charged systems of slab geometry |
| title_alt | Порівнення методів комп’ютерного моделювання для заряджених систем, обмежених плоскими поверхнями |
| title_full | Comparison of simulation methods for charged systems of slab geometry |
| title_fullStr | Comparison of simulation methods for charged systems of slab geometry |
| title_full_unstemmed | Comparison of simulation methods for charged systems of slab geometry |
| title_short | Comparison of simulation methods for charged systems of slab geometry |
| title_sort | comparison of simulation methods for charged systems of slab geometry |
| url | https://nasplib.isofts.kiev.ua/handle/123456789/120533 |
| work_keys_str_mv | AT mazarsm comparisonofsimulationmethodsforchargedsystemsofslabgeometry AT cailloljm comparisonofsimulationmethodsforchargedsystemsofslabgeometry AT weisjj comparisonofsimulationmethodsforchargedsystemsofslabgeometry AT levesqued comparisonofsimulationmethodsforchargedsystemsofslabgeometry AT mazarsm porívnennâmetodívkompûternogomodelûvannâdlâzarâdženihsistemobmeženihploskimipoverhnâmi AT cailloljm porívnennâmetodívkompûternogomodelûvannâdlâzarâdženihsistemobmeženihploskimipoverhnâmi AT weisjj porívnennâmetodívkompûternogomodelûvannâdlâzarâdženihsistemobmeženihploskimipoverhnâmi AT levesqued porívnennâmetodívkompûternogomodelûvannâdlâzarâdženihsistemobmeženihploskimipoverhnâmi |