ANew Methodfor Computationof Long Ranged
ANew Methodfor Computationof Long Ranged
ANew Methodfor Computationof Long Ranged
net/publication/226221631
CITATIONS READS
30 406
2 authors:
All content following this page was uploaded by Eugene Yakub on 21 May 2014.
1
European Commission, Joint Research Centre, Institute for
Transuranium Elements, P.O. Box 2340, D-76125 Karlsruhe, Germany;
2
Computer Science Dept., Odessa State Economic University,
Preobrazhenskaya 8, 65026, Odessa, Ukraine
E-mail: yakub@unive.odessa.ua
1. Introduction
633
1
Thereinafter simply called main cell.
Computation of Coulomb Forces in Disordered Systems 635
forces acting between the i-th and the j-th particles obey the Coulomb law:
Qi Qj
Fij = 2 .
4πε0 rij
Here rij = |r i −r j | is the distance between i-th and j-th charged particles and
Qi is the value of the i-th point charge of type α : Q(α) = {Q(1) , . . . , Q(M ) }.
We shall here assume that the electro-neutrality condition
N M
Nα Q(α) = 0
X X
Qi = (1)
i=1 α=1
is satisfied, and the standard PBC are imposed as described in Refs. 1,11.
The total Coulomb energy of N charges in the main cell is:
(C) X
UN = Qi ϕ(r i ) ,
1≤i≤N
Keeping in mind that all orientations of the main cell in an isotropic me-
dia should be equivalent, we can average both sides of (3) over all directions
of the vector n at a fixed distance rij . If the brackets
Z +1 Z π
1
h. . .i = d(cos ϑ) dψ . . . ,
4π −1 −π
indicate averaging, whereas ψ, ϑ are the polar and azimuthal angles defining
the direction of the vector n (n · r = nr cos ϑ), we can determine the pre-
averaged (effective) potential as ϕ2 (rij ) ≡ hϕ2 (r ij )i. Integration of (3) over
all orientations of the vector n gives immediately:
Qj rij
ϕ2 (rij ) = erfc δ
4πε0 rij L
" # )
1 X 1 π 2 n2 2π
+ exp − sin nrij . (4)
2π 2 n>0 n3 δ2 L
The pair effective (pre-averaged) potential (7) has the following properties:
The last property entails some inconsistency with the initial main cell con-
figuration. If the range of interaction does not exceed L/2 each particle in
the cell contributes (or not) to the sum of interactions with the selected one
just one time: either as original object (inside the main cell) or as one of its
“images.” This is not the case for potential (7) because rm > L/2.
The subsequent modification of the simulation algorithm was described
in Ref. 7. Two different zones exist within the effective sphere surrounding an
arbitrary charge in the main cell.7 The first zone contains charged particles
(inside the main cell) or their images (outside it) which contribute to the
sum of pair interactions of the chosen charge just one time. The second zone
contains those charged particles which contribute twice to the interaction en-
ergy with the given charge — both as original charges and as their additional
“phantom” images.7
where N is the number of positive ions. It should be noted that both terms
in (8) are size-dependent but not directly proportional to N . The first con-
tribution in (8) at given temperature T may be expressed in terms of plasma
parameter Γ = Q2 N 1/3 (4πε0 rm kT )−1 as −0.9kT Γ N 2/3 .
(OCP)
Table 1. Comparison of predicted Monte Carlo internal energies −UMC /N kT
for one-component plasma.
Another simple model system suitable for testing of the method is the
Restricted Primitive Ionic Model (RPIM), i.e. two-component fluid build
from oppositely charged hard spheres of one size (diameter σ). We per-
formed an additional set of Monte Carlo simulations for RPIM fluid at
Computation of Coulomb Forces in Disordered Systems 639
(MC)
Fig. 1. Thermal contribution to internal energy ∆E (th) /N kT = (UN −
(bcc)
UN )/N kT of one-component plasma at fixed Γ = 100 and different sizes of the
main cell. Comparison of our Monte Carlo simulation results (open circles) with
data of Hansen13 (solid square), and the “best fit” value recommended in Ref. 14
(solid triangle). Cubic spline interpolation of data presented in Table 3.3. was util-
(bcc)
ized to calculate the Madelung bcc-contribution UN .
Perfect crystalline structures are the worst application case for a pre-
averaged effective potential essentially devised for spatially uniform systems.
Nevertheless, the test on the perfect NaCl structure presented in Ref. 7 proves
the ability of the method to predict accurate values of the electrostatic energy
even in ordered condensed systems (see Table 1 in Ref. 7). In this work we
performed additional tests of the method using three other cubic crystalline
structures.
The values of the Madelung constant for bcc-structure of Coulomb crys-
tal formed in OCP, CsCl lattice and fluorite (CaF2 ) structure are computed
from (6)–(7) and presented in Tables 3.3., 3.3. and 3.3. at different numbers
of charged particles in the main cell.
640 E. Yakub and C. Ronchi
Fig. 2. Deviations of Monte Carlo internal energy from data of Larsen:15 ∆E/kT =
(MC) (MC)
(UN − UN )/N kT at fixed temperature T ∗ = 4πε0 kT σQ−2 = 0.595 and four
different packing fractions plotted against the size of the main cell: 1 – η = 0.005023
(squares), 2 – η = 0.02058 (triangles), 3 – η = 0.04842 (circles), 4 – η = 0.09525
(diamonds).
4. CONCLUDING REMARKS
Pre-averaging the Ewald sums over all spatial directions leads to an en-
ergy formulation which eliminates periodicity artifacts imposed by conven-
tional Ewald scheme and provides a fast method for computation of electro-
static contribution to the energy of disordered dense systems in Monte Carlo
or molecular dynamics computer simulations.7 Additional calculations and
Monte Carlo simulation results presented in Section 3. prove the ability of the
method to produce accurate simulation results at low computational cost.
Monte Carlo simulations of one-component plasma (small cell size N = 64)
presented in Table 1 are in good agreement with other simulation data.13,14
At the same time the comparison presented in Fig. 1 shows that it is very
likely coincidental. The convergence to a constant value at large N is non-
monotonous. Predicted thermal contribution to internal energy ∆E (th) /N kT
at Γ = 100 shown in Fig. 1 is in accordance with data of Hansen13 but is
higher than predicted in Ref. 14. This rather small (about 0.5% of the total
energy) difference is well beyond the total estimated error of both simulations.
This discrepancy may be due to periodicity artifacts and requires additional
study because it may affect the location of OCP crystallization point.
Computation of Coulomb Forces in Disordered Systems 641
Table 2. Predicted values of the Madelung constant for OCP Coulomb crystal (bcc
lattice). L/a is the ratio of the cell size to the lattice constant, in the second and
third columns are numbers of charged particles in the cell and in the equivalent
sphere. The fourth column show their difference (the actual net charge of the
equivalent sphere). The last two columns present the predicted Madelung constants
and their deviations from the exact value.
Nr of particles in: Deviation Madelung constant
L/a
Cell Sphere ∆(N ) This work ∆,%
2 16 15 −1 −0.89230 −0.41
3 54 59 5 −0.90361 +0.86
4 128 137 9 −0.89985 +0.44
5 250 259 9 −0.89921 +0.37
6 432 411 −21 −0.89595 +0.00
7 686 701 15 −0.89523 −0.08
8 1024 965 −59 −0.89411 −0.20
9 1458 1459 1 −0.89510 −0.09
10 2000 1989 −11 −0.89591 +0.00
11 2662 2685 23 −0.89656 +0.07
12 3456 3479 23 −0.89648 +0.06
13 4394 4477 83 −0.89633 +0.04
14 5488 5577 89 −0.89621 +0.03
15 6750 6735 −15 −0.89579 −0.02
Exact = −0.89593
Table 3. Predicted values of the Madelung constant for CsCl lattice. Columns are
the same as in the Table 2, except for the deviations in number of ions and net
charge (columns four and five) differ.
Table 4. Predicted values of the Madelung constant for the fluorite structure.
Columns 3 to 6 show the deviations in number of ions in equivalent spheres centered
on Ca2+ and F− ions, and the net charges within that spheres.
REFERENCES
1. R.W. Hockney and J.W. Eastwood, Computer Simulation Using Particles,
McGraw-Hill (1981).
2. R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
3. See, e.g., S.A. Bonev, B. Militzer, and G. Galli, Phys. Rev. B 69, 01 4101
(2004), and references therein.
4. E. Madelung, Phys. Z. 19, 524 (1918).
5. P.P. Ewald, Ann. Phys.(Leipzig) 64, 253 (1921).
6. G.J. Martyna and M.E. Tuckerman, J. Chem. Phys. 110, 2810 (1999);
Z.-H. Duan and R. Krasny, J. Chem. Phys. 113, 3492 (2000).
7. E. Yakub and C. Ronchi, J. Chem. Phys. 119, 11 556 (2003).
8. H.-Q. Ding, N. Karasawa, and W.A. Goddard, J. Chem. Phys. 97, 4309 (1992).
9. F. Gygi, Ab-initio Molecular Dynamics Code GP 1.20.0, Lawrence Livermore
National Laboratory (1999–2003).
10. W. Andreoni and A. Curioni, Parallel Computing 26, 819 (2000).
11. D.C. Rappaport, The Art of Molecular Dynamics Simulation, Cambridge
University Press (1995).
12. E.M. Lifshitz and L.P. Pitaevskii, Electrodynamics of Continuous Media.
Course of Theoretical Physics Vol. 8. 2. ed., Butterworth-Heinemann,
Oxford... (2000).
13. J.P. Hansen, Phys. Rev. A 8, 3110 (1973).
14. G.S. Stringfellow, H.E. DeWitt, and W.L. Slattery, Phys. Rev. A 41,
1105 (1990).
15. B. Larsen, Chem. Phys. Lett. 27, 47 (1974).
16. E.S. Yakub, J. Low Temp. Phys. 122, 559 (2001).