- Email: [email protected]
Ewald summation method with electrostatic layer correction for interactions of point dipoles in slab geometry
Chemical Physics Letters 400 (2004) 62–67 www.elsevier.com/locate/cplett
Ewald summation method with electrostatic layer correction for interactions of point dipoles in slab geometry A. Bro´dka Institute of Physics, University of Silesia, Uniwersytecka 4, 40-007 Katowice, Poland Received 14 September 2004; in final form 30 September 2004 Available online 6 November 2004
Abstract The Ewald method with electrostatic layer correction term for summing dipole–dipole interactions in slab geometry is presented. Theoretical cut-off errors in the electrostatic layer correction term as well as real- and reciprocal-space terms of the total energy are estimated, and their applicability is confirmed in numerical calculations. It is demonstrated that the summation method requires relative small gap of an empty space in the simulation box, and the theoretical errors can be used in determining optimal parameters of the method. 2004 Elsevier B.V. All rights reserved.
1. Introduction Computer simulations are important tools for the study of many-body systems. In many models it is necessary to consider polar systems, and a proper and efficient summation of long-range dipole–dipole interactions is of great importance. For three-dimensionally periodic systems the Ewald method [1,2] is the widely used technique to calculate the interaction sums. For a slab geometry, which occurs frequently in surface and interfacial systems, the conventional Ewald summation cannot be used since there is no periodicity in the one of three directions. The first Ewald type method for such a system was reported by Heyes [3]. The same expression was obtained using the integral representation of the gamma function and Poisson summation formula [4]. The main disadvantage of the method is low computational efficiency because the energy expression contains a sum over two-dimensional reciprocal lattice and all pairs of dipole moments. Using the approach proposed by Kawata and Mikami [5] to the reciprocal-space sum-
E-mail address: [email protected]. 0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.10.086
mation for two-dimensional reciprocal-space vector different from the zero vector, it was possible to obtain the modified method [6]. Additionally, the term for the zero vector in the two-dimensional reciprocal space can be modified [7], and we obtained EW3DC method in which one has three-dimensional reciprocal-space summation and the shape dependent correction (SDC) term for a thin slab proposed by Smith [8]. A disadvantage of those methods is relatively large empty space that must be introduced in the simulation box to obtain accurate results. However, in both methods the reciprocal-space terms contain single sums over dipole moments, and test molecular dynamics simulations for polar particles [6] showed that the calculations are about 25 times faster than those for the two-dimensional Ewald method. Analogous methods were worked out for Coulomb interactions. Particularly, Yeh and Berkowitz [9] adapted expression reported by Smith [8] for a macroscopic piece of ionic crystal in shape of a thin slab and they introduced an empty space along nonperiodicity direction to minimize an artificial influence from the periodic images in that direction. Arnold and co-workers [10] presented the EW3DLC method introducing the electrostatic layer correction (ELC) term,
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
which eliminates the unwanted interactions with the layers replicated along the non-periodicity direction, and gives possibility to diminish an empty space gap. Recently, energy expression in the EW3DLC method [11] was derived from the EW3DC method [7] for Coulomb interactions. In this Letter, the EW3DLC summation method for dipole–dipole interactions in a system with two-dimensional periodicity is presented, and theoretical cut-off errors of the interaction energy are derived. Moreover, the applicability of the energy and error formulas is tested in numerical calculations.
2. EW3DLC method for dipole moments Let us start from the EW3DLC method for a system of N charges which is periodic in the (x,y)-plane. For charges placed in a rectangular box with dimensions Lx · Ly · h, the simulation box has dimensions Lx · Ly · Lz where Lz > h, and the interaction energy is described by [10,11] N X0 erfcða j rkl þ n jÞ 2p X 1 1X þ qk ql 2 k;l¼1 j rkl þ n j V G6¼0 G2 n N G2 2p a X 2 exp 2 j S q ðGÞj þ M 2z pffiffiffi q2 V 4a p k¼1 k h p X 1 S qþ Gq S q Gq A Gq 6¼0 Gq ðeGq Lz 1Þ i þS q Gq S qþ Gq ;
energy, respectively. The third term is the SDC term, which is characteristic for a slab geometry [8]. The last term in Eq. (1) is the electrostatic layer correction (ELC) term, introduced by Arnold et al. [10]. The energy expression for a system of dipole moments may be obtained replacing the charge qk in Eqs. (1), (2) by the differential operator lk Æ $rk, where lk and rk are the kth dipole moment and its position, respectively. Using this substitution the real- and reciprocal-space terms as well as the SDC and ELC ones can be calculated directly. However, to calculate the selfcorrection energy for dipole–dipole interactions one must use the following representation of the fourth term in Eq. (1): Z a2 N N a X a X 2 pffiffiffi q2k ¼ pffiffiffi lim rkl !0 qk ql t1=2 erkl t dt; p k¼1 2 p k¼1 ql !qk 0 ð3Þ where the integral comes from the second term of Eq. (6) in [4] for the basic simulation box. The calculations lead to the following form of the interaction energy for polar system:
Eq ¼
E¼
ð1Þ
where Mz is the z-component of the total dipole moment of the basic simulation box and the structure factors are defined as follows: S q ðG Þ ¼
N X k¼1
qk e
iGrk
;
S q ðGq Þ ¼
N X
63
N X 0 1X fðlk ll ÞBða; j rkl þ n jÞ ½lk ðrkl þ nÞ 2 k;l¼1 n 2p X 1 ½ll ðrkl þ nÞCða; j rkl þ n jÞg þ V G6¼0 G2 N G2 2p 2a3 X 2 2 exp 2 j SðGÞj þ M 2z pffiffiffi j lk j V 4a 3 p k¼1 p X 1 S þ ðGq ÞS ðGq Þ G L q z A Gq 6¼0 Gq ðe 1Þ ð4Þ þ S ðGq ÞS þ ðGq Þ ;
where qk e
iGq qk Gq zk
e
:
ð2Þ
k¼1
In the above expressions rkl = rk rl, where rk is position of the charge qk in the basic simulation box, and lattice vectors n = (nxLx, nyLy, nzLz) with nx, ny, nz being integers describe positions of the simulation box replicas. The prime in the sum over the lattice vector n indicates that for n = 0 the terms with k = l are to be omitted. erfc (m) is the complementary error function [12], and the convergence parameter a > 0 is chosen for computational convenience. G = 2p(kx/Lx, ky/Ly, kz/Lz) with kx, ky, kz being integers is three-dimensional reciprocal-space vector and V = LxLyLz is the simulation box volume with an empty space. qk and Gq denote projections of the vectors rk and G onto the (x,y)-plane, respectively, and A = LxLy is area of the simulation box base. The first, second and fourth terms in Eq. (1) represent the real-space summation, the reciprocal-space one and self-correction
pffiffiffi Bða; rÞ ¼ erfcðarÞ=r3 þ ð2a= pÞ expða2 r2 Þ=r2 ;
ð5Þ
pffiffiffi Cða; rÞ ¼ 3erfcðarÞ=r5 þ ð2a= pÞð2a2 þ 3=r2 Þ expða2 r2 Þ=r2 ; SðGÞ ¼
N X
ð6Þ
ðlk GÞeiGrk ;
ð7Þ
k¼1 N N X X q S k Gq ¼ ilk Gq lzk Gq eiGq qk eGq zk : S Gq ¼ k¼1
k¼1
ð8Þ lqk
lzk
and denote components In the above equations of a dipole moment lk that are parallel and perpendicular to the (x,y)-plane, and the other symbols have the same meaning as previously. Expressions (4)–(8) describe energy in the EW3DLC method for point dipole moments.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
64
3. Error formulas Using the Ewald type method one has errors arising from cutting-off the summation in the real and reciprocal space. Therefore, for the EW3DLC method accuracy of calculation of the total interaction energy results from the first, second and last term of expression (4), and it can be estimated as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi real 2 2 2 dE dE ð9Þ þ ðdErec Þ þ dELC ;
D
E sin Gq qkl cosh Gq zkl sin G0q qkl cosh G0q zkl
i 1h ¼ d Gq ; G0q d Gq ; G0q 2D
E cosh Gq zkl cosh G0q zkl ; ð14Þ
D
E cos Gq qkl sinh Gq zkl cos G0q qkl sinh G0q zkl
i 1h ¼ d Gq ; G0q þ d Gq ; G0q 2D
E sinh Gq zkl sinh G0q zkl ; ð15Þ
where dEreal and dErec are the real- and reciprocal-space errors, respectively, and dELC means error of the ELC term. To estimate the ELC error one may apply the approach used by Kolafa and Perram [13] to calculate the reciprocal-space error for Coulomb interactions. In that case the ELC error consists of the off-diagonal term and diagonal one vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u XN N D E X LC u1 LC LC 2 dE t dE dEkk : þ ð10Þ kl k;l¼1 2 k¼1
D
In the above equation angular brackets mean the configuration average, and
where d (G,G 0 ) is the Kronecker delta. Now the averages in Eq. (10) are following: 2 dELC kl
k6¼l
q 4p 2 X Gq
l sin h k sin hl cos # lk ; Gq G L Gq e q z 1 A Gq >Gc cos # lql ; Gq cos Gq qkl cosh Gq zkl q þ sin hk cos hl cos # lk ; Gq þ cos hk sin hl cos # lql ; Gq sin Gq qkl sinh Gq zkl cos hk cos hl cos Gq qkl cosh Gq zkl ; ð11Þ
E
2 G2q pl2 X ¼ 2 Gq A ðeGq Lz 1Þ Gq >Gc
cos2 # lqk ; Gq cos2 # lql ; Gq þ 1 cosh2 Gq zkl þ cos2 # lqk ; Gq þ cos2 # lql ; Gq sinh2 Gq zkl ; ð16Þ
p X Gq 2 q cos # l ; G ¼ l2 dELC q 1 : k kk G L Gq e q z 1 A Gq >Gc
dELC kl ¼
dELC kk ¼
2p 2 X l Gq A Gq >Gc
Gq 2 sin hk cos2 # lqk ; Gq cos2 hk ; 1
eGq Lz
ð12Þ where l is value of a dipole moment, Gc = 2pKc/Lx is cut-off radius in the reciprocal space and hk is an angle between the kth dipole moment and z-axis of the laboratory frame. To calculate average values of the terms appearing in Eq. (10) one assumes random orientations of dipole moments with respect to the z-axis of the laboratory frame and for uncorrelated dipole moments one can set
1 sin2 hk ¼ cos2 hk ¼ ; 2
hsin hk cos hl i ¼ 0:
ð13Þ
Moreover, assuming isotropic distribution of dipole moments in the (x,y)-plane and lack of correlation between inter-dipole distance components qkl and zkl one has
ð17Þ Further calculations require additional information about distribution of zkl. The maximum ELC error appears when large values of zkl are the most probable values for the largest number of dipole pairs. Such a case is for the dipole moments located in two layers distanced by h, and probability of the inter-dipole distance zkl = 0 or h is equal to 0.5 that gives possibility to estimate the averages appearing in Eq. (16). The two-dimensional sum over Gq, in Eqs. (16) and (17), is estimated by integrals with respect to polar co-ordinates of the Gq vector X Gq Gq >Gc
F ðGq Þ
Z
A ð2pÞ
2
Z
2p
1
d/ 0
dGq Gq F Gq ; / :
ð18Þ
Gc
Assuming that the x-axis of the polar co-ordinates is along the lqk vector, i.e. #ðlqk ; Gq Þ ¼ /, one has # lql ; Gq ¼ cos / cos #ðlqk ; lql Þ þ sin / sin #ðlqk ; lql Þ ð19Þ and for random orientations of the dipole moments the trigonometric functions of #ðlqk ; lql Þ can be replaced by their configuration average values, similarly as those in Eqs. (13). Calculating the integral with respect to the angle / the averages (16) and (17) can be rewritten in the following forms:
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
2 E
¼
pl4 8A
Z
G3q
1
Gc
ðeGq Lz
1Þ
2
11 þ 9cosh 2Gq h dGq ; ð20Þ
dELC kk
l2 ¼ 4
Z
1
Gc
G2q dGq : eGq Lz 1
ð21Þ
Analytical calculation of the integrals in the above equations is impossible, however, they can be estimate as follows: Z 1 D 2 E pl4 1 dELC G3 e2Gq Lz kl 8A ð1 eGc Lz Þ2 Gc q 11 þ 9 cosh 2Gq h dGq ; ð22Þ
dELC kk
l2 1 4 1 eGc Lz
Z
1
G2q e2Gq Lz dGq :
ð23Þ
Gc
Calculating the integrals with respect to Gq, and using Eq. (10) the upper bound of the ELC error has the following form: rffiffiffiffiffiffi N l2 1 1 p 2Gc h LC dE 9e g1 ðGc ;Lz hÞ 4 eGc Lz 1 2 2A o 1=2 þ22g1 ðGc ;Lz Þþ9e2Gc h g1 ðGc ;Lz þhÞ þg2 ðGc ;Lz Þ ; ð24Þ where g1 ðG; xÞ ¼ G3 =x þ 3G2 =2x2 þ 3G=2x3 þ 3=4x4 ; g2 ðG; xÞ ¼ G2 =x þ 2G=x2 þ 2=x3 :
ð25Þ
The real- and reciprocal-space errors can be estimated using approach proposed by Wang and Holm [14] for a system with three-dimensional periodicity. For arbitrary box size the errors have the following forms: 1=2 2 1 1 2 1 1 2 real N l dE ¼ pffiffiffiffi 2 3 pffiffiffiffi Bc Bc C c þ C c exp a2 r2c ; a 4 6 15 r r V c c ð26Þ ! rffiffiffiffiffiffiffiffiffi Gc 2a2 Gc G2 rec 2 dE ¼ N l 4a ð27Þ þ exp c2 ; 30V 3p 4a where rc 6 min(Lx, Ly, Lz)/2 is real-space cut-off radius, and Bc ¼ 2a2 r2c þ 1;
C c ¼ 4a4 r4c þ 6a2 r2c þ 3:
4. EW3DLC method in numerical calculations To check the working of the EW3DLC method and the error formulas, which are presented in Sections 2 and 3, the interaction energy is calculated for 1000 dipole moments in a cubic box with side L, i.e. Lx = Ly = h = L. Ten configurations of the dipole moments with random distribution of their positions and orientations are generated. The simulation energy error is calculated as the root mean square error D E1=2 2 DE ¼ ðE Eexact Þ ; ð29Þ where Eexact means exact energy, and the angular brackets denote average over the system configurations. Calculation of the exact energy using direct summation in real space is very time consuming and inaccurate, because the sum over nq, i.e. simulation box replicas in the (x,y)-plane, is slowly convergent. Therefore, the exact energy of the system is calculated from the Lekner method [15]. To obtain the relative accuracy of the energy better than 1016, the maximum value of the summation variable m in Eqs. (32–34) of [15] equals 6, and truncation parameter pcut in the inequality (35) of [15] equals 1020. In all calculations the real-space cut-off radius rc equals L/2. First of all, one has to find an optimal value of height of the simulation box with an empty space. Therefore, for Lz ranging from L to 2L, the energy errors described by Eqs. (26) and (27) are calculated as functions of the convergence parameter a for different values of the reciprocal-space cut-off radius Kc. Similarly, error of the ELC term arising form Eq. (24) is calculated as a function of Kc. Examples of the theoretical real- and reciprocal-space errors for Lz = 1.3L are presented in Fig. 1, and the corresponding ELC errors are shown in Fig. 2. For a particular value of Kc one may
∆ δ
µ
dELC kl
δ
∆
D
65
ð28Þ
Taking into account an empty space in the simulation box the formula (26) may lead to overestimation of the real-space error. In the case of the reciprocal-space error the diagonal term in (27) is dominated, and one expects similar behaviour of the error as in the bulk system.
α Fig. 1. The theoretical real- and reciprocal-space errors, dEreal and dErec, and total numerical error, DE, as functions of the convergence parameter, a, and reciprocal-space cut-off radius, Kc, for rc = L/2 and Lz = 1.3L.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
66
∆
µ
cal-space errors reproduce the total numerical ones, and the minimum of DE for particular value of Kc appears approximately for a value of a when dEreal @ dErec. In Fig. 2 there are presented the numerical errors of the ELC term assuming that the exact value of this term is given for Kc = 50. This figure shows almost exponential decay of the theoretical error as well as numerical one with the reciprocal-space cut-off radius. From theoretical point of view this behaviour is comprehensible because the term consisting exp(2Gch) is dominating in the error formula (24), and the ELC error decreases as exp (Gc(L h)).
δ ∆
Fig. 2. The theoretical, dELC, and numerical error, DELC, as functions of the reciprocal-space cut-off radius, Kc, for Lz = 1.3L.
find from Fig. 1 the convergence parameter a for which dEreal @ dErec, and for those parameters the theoretical errors as functions of Lz are presented in Fig. 3. It is clearly visible that the real- and reciprocal-space errors depend on Lz weakly, whereas the ELC error is very high for small values of Lz and it decreases with Lz fast. The results presented in Fig. 3 gives possibility to determine an optimal value of the simulation box height, and we choose Lz = 1.3L for which the theoretical ELC error is an order of magnitude smaller than that arising from cutting-off in the real and reciprocal space. In Fig. 3 there are also showed errors of the total energy, DE, obtained from numerical calculations. This comparison indicates that for small values of Lz the energy error results mainly from the ELC term, whereas for larger Lz the real- and reciprocal-space errors are dominated. The optimal value of the simulation box height used in numerical calculation leads to the total energy errors DE depicted in Fig. 1. The theoretical real- and recipro-
((δ
δ
δ
∆
µ
∆
5. Conclusions The Ewald summation method with the ELC term for dipole–dipole interactions in two-dimensional system is derived. In the presented derivation we apply the well-known procedure, where operators involving dipole moments and gradients replace charges in the formula for Coulomb interactions. Moreover, the upper bound of the ELC term error is calculated. To estimate the real- and reciprocal-space cut-off errors we use the expressions for a system with three-dimensional periodicity. Comparison of the theoretical errors with the numerical ones shows that the real- and reciprocal-space error formulas allows us to determine optimal values of the convergence parameter and reciprocal-space cut-off radius, whereas the expression for the ELC term error gives possibility to choose the simulation box height. Such an estimation of the optimal parameters that give the expected accuracy can be treated as a starting point in using EW3DLC method for simulation of polar system in slab geometry. The EW3DLC method for polar system, similarly as for Coulomb interactions, requires relatively small height of the simulation box. The presented theoretical and numerical calculations show that using Lz = 1.3L the interaction error is practically determined by errors arising from the cut-off in the real space and reciprocal one. In comparison with the EW3DC method [7] and the modified Ewald summation method [6,15] the height of the simulation box with an empty space is reduced about five times. That means, that number of the reciprocal-space vectors along the z-axis in the EW3DLC method is five times smaller, and the method should be five times faster then the other methods.
References Fig. 3. The theoretical real- and reciprocal-space error, ((dEreal)2 + (dErec)2)1/2 ELC error, dELC, and total numerical error, DE, as functions of the simulation box height, Lz, for a = 8.3/L, rc=L/2 and Kc = 12.
[1] S.W. De Leeuw, J.W. Perram, E.R. Smith, Proc. R. Soc. Lond., Ser. A 373 (1980) 27. [2] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon, Oxford University Press, 1989.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67 [3] D.M. Heyes, Surf. Sci. 110 (1981) L619. [4] A. Grzybowski, E. Gwo´z´dz´, A. Bro´dka, Phys. Rev. B 61 (2000) 6706. [5] M. Kawata, M. Mikami, Chem. Phys. Lett. 340 (2001) 157. [6] A. Grzybowski, A. Bro´dka, Chem. Phys. Lett. 361 (2002) 329. [7] A. Bro´dka, A. Grzybowski, J. Chem. Phys. 117 (2002) 8208. [8] E.R. Smith, Proc. R. Soc. Lond., Ser. A 375 (1981) 475. [9] I.-C. Yeh, M.L. Berkowitz, J. Chem. Phys. 111 (1999) 3155.
67
[10] A. Arnold, J. de Joannis, C. Holm, J. Chem. Phys. 117 (2002) 2496. [11] A. Bro´dka, J. Chem. Phys. 121 (2004) 7032. [12] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1965. [13] J. Kolafa, J.W. Perram, Mol. Simulat. 9 (1992) 351. [14] Z. Wang, C. Holm, J. Chem. Phys. 115 (2001) 6351. [15] A. Grzybowski, A. Bro´dka, Mol. Phys. 101 (2003) 1079.
Ewald summation method with electrostatic layer correction for interactions of point dipoles in slab geometry A. Bro´dka Institute of Physics, University of Silesia, Uniwersytecka 4, 40-007 Katowice, Poland Received 14 September 2004; in final form 30 September 2004 Available online 6 November 2004
Abstract The Ewald method with electrostatic layer correction term for summing dipole–dipole interactions in slab geometry is presented. Theoretical cut-off errors in the electrostatic layer correction term as well as real- and reciprocal-space terms of the total energy are estimated, and their applicability is confirmed in numerical calculations. It is demonstrated that the summation method requires relative small gap of an empty space in the simulation box, and the theoretical errors can be used in determining optimal parameters of the method. 2004 Elsevier B.V. All rights reserved.
1. Introduction Computer simulations are important tools for the study of many-body systems. In many models it is necessary to consider polar systems, and a proper and efficient summation of long-range dipole–dipole interactions is of great importance. For three-dimensionally periodic systems the Ewald method [1,2] is the widely used technique to calculate the interaction sums. For a slab geometry, which occurs frequently in surface and interfacial systems, the conventional Ewald summation cannot be used since there is no periodicity in the one of three directions. The first Ewald type method for such a system was reported by Heyes [3]. The same expression was obtained using the integral representation of the gamma function and Poisson summation formula [4]. The main disadvantage of the method is low computational efficiency because the energy expression contains a sum over two-dimensional reciprocal lattice and all pairs of dipole moments. Using the approach proposed by Kawata and Mikami [5] to the reciprocal-space sum-
E-mail address: [email protected]. 0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.10.086
mation for two-dimensional reciprocal-space vector different from the zero vector, it was possible to obtain the modified method [6]. Additionally, the term for the zero vector in the two-dimensional reciprocal space can be modified [7], and we obtained EW3DC method in which one has three-dimensional reciprocal-space summation and the shape dependent correction (SDC) term for a thin slab proposed by Smith [8]. A disadvantage of those methods is relatively large empty space that must be introduced in the simulation box to obtain accurate results. However, in both methods the reciprocal-space terms contain single sums over dipole moments, and test molecular dynamics simulations for polar particles [6] showed that the calculations are about 25 times faster than those for the two-dimensional Ewald method. Analogous methods were worked out for Coulomb interactions. Particularly, Yeh and Berkowitz [9] adapted expression reported by Smith [8] for a macroscopic piece of ionic crystal in shape of a thin slab and they introduced an empty space along nonperiodicity direction to minimize an artificial influence from the periodic images in that direction. Arnold and co-workers [10] presented the EW3DLC method introducing the electrostatic layer correction (ELC) term,
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
which eliminates the unwanted interactions with the layers replicated along the non-periodicity direction, and gives possibility to diminish an empty space gap. Recently, energy expression in the EW3DLC method [11] was derived from the EW3DC method [7] for Coulomb interactions. In this Letter, the EW3DLC summation method for dipole–dipole interactions in a system with two-dimensional periodicity is presented, and theoretical cut-off errors of the interaction energy are derived. Moreover, the applicability of the energy and error formulas is tested in numerical calculations.
2. EW3DLC method for dipole moments Let us start from the EW3DLC method for a system of N charges which is periodic in the (x,y)-plane. For charges placed in a rectangular box with dimensions Lx · Ly · h, the simulation box has dimensions Lx · Ly · Lz where Lz > h, and the interaction energy is described by [10,11] N X0 erfcða j rkl þ n jÞ 2p X 1 1X þ qk ql 2 k;l¼1 j rkl þ n j V G6¼0 G2 n N G2 2p a X 2 exp 2 j S q ðGÞj þ M 2z pffiffiffi q2 V 4a p k¼1 k h p X 1 S qþ Gq S q Gq A Gq 6¼0 Gq ðeGq Lz 1Þ i þS q Gq S qþ Gq ;
energy, respectively. The third term is the SDC term, which is characteristic for a slab geometry [8]. The last term in Eq. (1) is the electrostatic layer correction (ELC) term, introduced by Arnold et al. [10]. The energy expression for a system of dipole moments may be obtained replacing the charge qk in Eqs. (1), (2) by the differential operator lk Æ $rk, where lk and rk are the kth dipole moment and its position, respectively. Using this substitution the real- and reciprocal-space terms as well as the SDC and ELC ones can be calculated directly. However, to calculate the selfcorrection energy for dipole–dipole interactions one must use the following representation of the fourth term in Eq. (1): Z a2 N N a X a X 2 pffiffiffi q2k ¼ pffiffiffi lim rkl !0 qk ql t1=2 erkl t dt; p k¼1 2 p k¼1 ql !qk 0 ð3Þ where the integral comes from the second term of Eq. (6) in [4] for the basic simulation box. The calculations lead to the following form of the interaction energy for polar system:
Eq ¼
E¼
ð1Þ
where Mz is the z-component of the total dipole moment of the basic simulation box and the structure factors are defined as follows: S q ðG Þ ¼
N X k¼1
qk e
iGrk
;
S q ðGq Þ ¼
N X
63
N X 0 1X fðlk ll ÞBða; j rkl þ n jÞ ½lk ðrkl þ nÞ 2 k;l¼1 n 2p X 1 ½ll ðrkl þ nÞCða; j rkl þ n jÞg þ V G6¼0 G2 N G2 2p 2a3 X 2 2 exp 2 j SðGÞj þ M 2z pffiffiffi j lk j V 4a 3 p k¼1 p X 1 S þ ðGq ÞS ðGq Þ G L q z A Gq 6¼0 Gq ðe 1Þ ð4Þ þ S ðGq ÞS þ ðGq Þ ;
where qk e
iGq qk Gq zk
e
:
ð2Þ
k¼1
In the above expressions rkl = rk rl, where rk is position of the charge qk in the basic simulation box, and lattice vectors n = (nxLx, nyLy, nzLz) with nx, ny, nz being integers describe positions of the simulation box replicas. The prime in the sum over the lattice vector n indicates that for n = 0 the terms with k = l are to be omitted. erfc (m) is the complementary error function [12], and the convergence parameter a > 0 is chosen for computational convenience. G = 2p(kx/Lx, ky/Ly, kz/Lz) with kx, ky, kz being integers is three-dimensional reciprocal-space vector and V = LxLyLz is the simulation box volume with an empty space. qk and Gq denote projections of the vectors rk and G onto the (x,y)-plane, respectively, and A = LxLy is area of the simulation box base. The first, second and fourth terms in Eq. (1) represent the real-space summation, the reciprocal-space one and self-correction
pffiffiffi Bða; rÞ ¼ erfcðarÞ=r3 þ ð2a= pÞ expða2 r2 Þ=r2 ;
ð5Þ
pffiffiffi Cða; rÞ ¼ 3erfcðarÞ=r5 þ ð2a= pÞð2a2 þ 3=r2 Þ expða2 r2 Þ=r2 ; SðGÞ ¼
N X
ð6Þ
ðlk GÞeiGrk ;
ð7Þ
k¼1 N N X X q S k Gq ¼ ilk Gq lzk Gq eiGq qk eGq zk : S Gq ¼ k¼1
k¼1
ð8Þ lqk
lzk
and denote components In the above equations of a dipole moment lk that are parallel and perpendicular to the (x,y)-plane, and the other symbols have the same meaning as previously. Expressions (4)–(8) describe energy in the EW3DLC method for point dipole moments.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
64
3. Error formulas Using the Ewald type method one has errors arising from cutting-off the summation in the real and reciprocal space. Therefore, for the EW3DLC method accuracy of calculation of the total interaction energy results from the first, second and last term of expression (4), and it can be estimated as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi real 2 2 2 dE dE ð9Þ þ ðdErec Þ þ dELC ;
D
E sin Gq qkl cosh Gq zkl sin G0q qkl cosh G0q zkl
i 1h ¼ d Gq ; G0q d Gq ; G0q 2D
E cosh Gq zkl cosh G0q zkl ; ð14Þ
D
E cos Gq qkl sinh Gq zkl cos G0q qkl sinh G0q zkl
i 1h ¼ d Gq ; G0q þ d Gq ; G0q 2D
E sinh Gq zkl sinh G0q zkl ; ð15Þ
where dEreal and dErec are the real- and reciprocal-space errors, respectively, and dELC means error of the ELC term. To estimate the ELC error one may apply the approach used by Kolafa and Perram [13] to calculate the reciprocal-space error for Coulomb interactions. In that case the ELC error consists of the off-diagonal term and diagonal one vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u XN N D E X LC u1 LC LC 2 dE t dE dEkk : þ ð10Þ kl k;l¼1 2 k¼1
D
In the above equation angular brackets mean the configuration average, and
where d (G,G 0 ) is the Kronecker delta. Now the averages in Eq. (10) are following: 2 dELC kl
k6¼l
q 4p 2 X Gq
l sin h k sin hl cos # lk ; Gq G L Gq e q z 1 A Gq >Gc cos # lql ; Gq cos Gq qkl cosh Gq zkl q þ sin hk cos hl cos # lk ; Gq þ cos hk sin hl cos # lql ; Gq sin Gq qkl sinh Gq zkl cos hk cos hl cos Gq qkl cosh Gq zkl ; ð11Þ
E
2 G2q pl2 X ¼ 2 Gq A ðeGq Lz 1Þ Gq >Gc
cos2 # lqk ; Gq cos2 # lql ; Gq þ 1 cosh2 Gq zkl þ cos2 # lqk ; Gq þ cos2 # lql ; Gq sinh2 Gq zkl ; ð16Þ
p X Gq 2 q cos # l ; G ¼ l2 dELC q 1 : k kk G L Gq e q z 1 A Gq >Gc
dELC kl ¼
dELC kk ¼
2p 2 X l Gq A Gq >Gc
Gq 2 sin hk cos2 # lqk ; Gq cos2 hk ; 1
eGq Lz
ð12Þ where l is value of a dipole moment, Gc = 2pKc/Lx is cut-off radius in the reciprocal space and hk is an angle between the kth dipole moment and z-axis of the laboratory frame. To calculate average values of the terms appearing in Eq. (10) one assumes random orientations of dipole moments with respect to the z-axis of the laboratory frame and for uncorrelated dipole moments one can set
1 sin2 hk ¼ cos2 hk ¼ ; 2
hsin hk cos hl i ¼ 0:
ð13Þ
Moreover, assuming isotropic distribution of dipole moments in the (x,y)-plane and lack of correlation between inter-dipole distance components qkl and zkl one has
ð17Þ Further calculations require additional information about distribution of zkl. The maximum ELC error appears when large values of zkl are the most probable values for the largest number of dipole pairs. Such a case is for the dipole moments located in two layers distanced by h, and probability of the inter-dipole distance zkl = 0 or h is equal to 0.5 that gives possibility to estimate the averages appearing in Eq. (16). The two-dimensional sum over Gq, in Eqs. (16) and (17), is estimated by integrals with respect to polar co-ordinates of the Gq vector X Gq Gq >Gc
F ðGq Þ
Z
A ð2pÞ
2
Z
2p
1
d/ 0
dGq Gq F Gq ; / :
ð18Þ
Gc
Assuming that the x-axis of the polar co-ordinates is along the lqk vector, i.e. #ðlqk ; Gq Þ ¼ /, one has # lql ; Gq ¼ cos / cos #ðlqk ; lql Þ þ sin / sin #ðlqk ; lql Þ ð19Þ and for random orientations of the dipole moments the trigonometric functions of #ðlqk ; lql Þ can be replaced by their configuration average values, similarly as those in Eqs. (13). Calculating the integral with respect to the angle / the averages (16) and (17) can be rewritten in the following forms:
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
2 E
¼
pl4 8A
Z
G3q
1
Gc
ðeGq Lz
1Þ
2
11 þ 9cosh 2Gq h dGq ; ð20Þ
dELC kk
l2 ¼ 4
Z
1
Gc
G2q dGq : eGq Lz 1
ð21Þ
Analytical calculation of the integrals in the above equations is impossible, however, they can be estimate as follows: Z 1 D 2 E pl4 1 dELC G3 e2Gq Lz kl 8A ð1 eGc Lz Þ2 Gc q 11 þ 9 cosh 2Gq h dGq ; ð22Þ
dELC kk
l2 1 4 1 eGc Lz
Z
1
G2q e2Gq Lz dGq :
ð23Þ
Gc
Calculating the integrals with respect to Gq, and using Eq. (10) the upper bound of the ELC error has the following form: rffiffiffiffiffiffi N l2 1 1 p 2Gc h LC dE 9e g1 ðGc ;Lz hÞ 4 eGc Lz 1 2 2A o 1=2 þ22g1 ðGc ;Lz Þþ9e2Gc h g1 ðGc ;Lz þhÞ þg2 ðGc ;Lz Þ ; ð24Þ where g1 ðG; xÞ ¼ G3 =x þ 3G2 =2x2 þ 3G=2x3 þ 3=4x4 ; g2 ðG; xÞ ¼ G2 =x þ 2G=x2 þ 2=x3 :
ð25Þ
The real- and reciprocal-space errors can be estimated using approach proposed by Wang and Holm [14] for a system with three-dimensional periodicity. For arbitrary box size the errors have the following forms: 1=2 2 1 1 2 1 1 2 real N l dE ¼ pffiffiffiffi 2 3 pffiffiffiffi Bc Bc C c þ C c exp a2 r2c ; a 4 6 15 r r V c c ð26Þ ! rffiffiffiffiffiffiffiffiffi Gc 2a2 Gc G2 rec 2 dE ¼ N l 4a ð27Þ þ exp c2 ; 30V 3p 4a where rc 6 min(Lx, Ly, Lz)/2 is real-space cut-off radius, and Bc ¼ 2a2 r2c þ 1;
C c ¼ 4a4 r4c þ 6a2 r2c þ 3:
4. EW3DLC method in numerical calculations To check the working of the EW3DLC method and the error formulas, which are presented in Sections 2 and 3, the interaction energy is calculated for 1000 dipole moments in a cubic box with side L, i.e. Lx = Ly = h = L. Ten configurations of the dipole moments with random distribution of their positions and orientations are generated. The simulation energy error is calculated as the root mean square error D E1=2 2 DE ¼ ðE Eexact Þ ; ð29Þ where Eexact means exact energy, and the angular brackets denote average over the system configurations. Calculation of the exact energy using direct summation in real space is very time consuming and inaccurate, because the sum over nq, i.e. simulation box replicas in the (x,y)-plane, is slowly convergent. Therefore, the exact energy of the system is calculated from the Lekner method [15]. To obtain the relative accuracy of the energy better than 1016, the maximum value of the summation variable m in Eqs. (32–34) of [15] equals 6, and truncation parameter pcut in the inequality (35) of [15] equals 1020. In all calculations the real-space cut-off radius rc equals L/2. First of all, one has to find an optimal value of height of the simulation box with an empty space. Therefore, for Lz ranging from L to 2L, the energy errors described by Eqs. (26) and (27) are calculated as functions of the convergence parameter a for different values of the reciprocal-space cut-off radius Kc. Similarly, error of the ELC term arising form Eq. (24) is calculated as a function of Kc. Examples of the theoretical real- and reciprocal-space errors for Lz = 1.3L are presented in Fig. 1, and the corresponding ELC errors are shown in Fig. 2. For a particular value of Kc one may
∆ δ
µ
dELC kl
δ
∆
D
65
ð28Þ
Taking into account an empty space in the simulation box the formula (26) may lead to overestimation of the real-space error. In the case of the reciprocal-space error the diagonal term in (27) is dominated, and one expects similar behaviour of the error as in the bulk system.
α Fig. 1. The theoretical real- and reciprocal-space errors, dEreal and dErec, and total numerical error, DE, as functions of the convergence parameter, a, and reciprocal-space cut-off radius, Kc, for rc = L/2 and Lz = 1.3L.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67
66
∆
µ
cal-space errors reproduce the total numerical ones, and the minimum of DE for particular value of Kc appears approximately for a value of a when dEreal @ dErec. In Fig. 2 there are presented the numerical errors of the ELC term assuming that the exact value of this term is given for Kc = 50. This figure shows almost exponential decay of the theoretical error as well as numerical one with the reciprocal-space cut-off radius. From theoretical point of view this behaviour is comprehensible because the term consisting exp(2Gch) is dominating in the error formula (24), and the ELC error decreases as exp (Gc(L h)).
δ ∆
Fig. 2. The theoretical, dELC, and numerical error, DELC, as functions of the reciprocal-space cut-off radius, Kc, for Lz = 1.3L.
find from Fig. 1 the convergence parameter a for which dEreal @ dErec, and for those parameters the theoretical errors as functions of Lz are presented in Fig. 3. It is clearly visible that the real- and reciprocal-space errors depend on Lz weakly, whereas the ELC error is very high for small values of Lz and it decreases with Lz fast. The results presented in Fig. 3 gives possibility to determine an optimal value of the simulation box height, and we choose Lz = 1.3L for which the theoretical ELC error is an order of magnitude smaller than that arising from cutting-off in the real and reciprocal space. In Fig. 3 there are also showed errors of the total energy, DE, obtained from numerical calculations. This comparison indicates that for small values of Lz the energy error results mainly from the ELC term, whereas for larger Lz the real- and reciprocal-space errors are dominated. The optimal value of the simulation box height used in numerical calculation leads to the total energy errors DE depicted in Fig. 1. The theoretical real- and recipro-
((δ
δ
δ
∆
µ
∆
5. Conclusions The Ewald summation method with the ELC term for dipole–dipole interactions in two-dimensional system is derived. In the presented derivation we apply the well-known procedure, where operators involving dipole moments and gradients replace charges in the formula for Coulomb interactions. Moreover, the upper bound of the ELC term error is calculated. To estimate the real- and reciprocal-space cut-off errors we use the expressions for a system with three-dimensional periodicity. Comparison of the theoretical errors with the numerical ones shows that the real- and reciprocal-space error formulas allows us to determine optimal values of the convergence parameter and reciprocal-space cut-off radius, whereas the expression for the ELC term error gives possibility to choose the simulation box height. Such an estimation of the optimal parameters that give the expected accuracy can be treated as a starting point in using EW3DLC method for simulation of polar system in slab geometry. The EW3DLC method for polar system, similarly as for Coulomb interactions, requires relatively small height of the simulation box. The presented theoretical and numerical calculations show that using Lz = 1.3L the interaction error is practically determined by errors arising from the cut-off in the real space and reciprocal one. In comparison with the EW3DC method [7] and the modified Ewald summation method [6,15] the height of the simulation box with an empty space is reduced about five times. That means, that number of the reciprocal-space vectors along the z-axis in the EW3DLC method is five times smaller, and the method should be five times faster then the other methods.
References Fig. 3. The theoretical real- and reciprocal-space error, ((dEreal)2 + (dErec)2)1/2 ELC error, dELC, and total numerical error, DE, as functions of the simulation box height, Lz, for a = 8.3/L, rc=L/2 and Kc = 12.
[1] S.W. De Leeuw, J.W. Perram, E.R. Smith, Proc. R. Soc. Lond., Ser. A 373 (1980) 27. [2] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon, Oxford University Press, 1989.
A. Bro´dka / Chemical Physics Letters 400 (2004) 62–67 [3] D.M. Heyes, Surf. Sci. 110 (1981) L619. [4] A. Grzybowski, E. Gwo´z´dz´, A. Bro´dka, Phys. Rev. B 61 (2000) 6706. [5] M. Kawata, M. Mikami, Chem. Phys. Lett. 340 (2001) 157. [6] A. Grzybowski, A. Bro´dka, Chem. Phys. Lett. 361 (2002) 329. [7] A. Bro´dka, A. Grzybowski, J. Chem. Phys. 117 (2002) 8208. [8] E.R. Smith, Proc. R. Soc. Lond., Ser. A 375 (1981) 475. [9] I.-C. Yeh, M.L. Berkowitz, J. Chem. Phys. 111 (1999) 3155.
67
[10] A. Arnold, J. de Joannis, C. Holm, J. Chem. Phys. 117 (2002) 2496. [11] A. Bro´dka, J. Chem. Phys. 121 (2004) 7032. [12] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1965. [13] J. Kolafa, J.W. Perram, Mol. Simulat. 9 (1992) 351. [14] Z. Wang, C. Holm, J. Chem. Phys. 115 (2001) 6351. [15] A. Grzybowski, A. Bro´dka, Mol. Phys. 101 (2003) 1079.