NFFT based Ewald summation for electrostatic systems with charges and dipoles

NFFT based Ewald summation for electrostatic systems with charges and

dipoles

Michael Hofmann

Chemnitz University of Technology, Faculty of Computer Science, 09107 Chemnitz, Germany

Franziska Nestler

, Michael Pippig

Chemnitz University of Technology, Faculty of Mathematics, 09107 Chemnitz, Germany

The efficient computation of Coulomb interactions in charged particle systems is of great importance in the field of molecular dynamics simulations. It is widely known that an approximation can be realized based on the Ewald summation approach and the fast Fourier transform (FFT). In the present paper we consider particle systems containing a mixture ofN point charges as well as point dipoles. New cutoff errors in the Ewald summation formulas concerning charge-dipole interactions are derived and, moreover, validated by numerical examples. Furthermore, we present for the first time anO(NlogN) particle mesh algorithm for computing mixed charge-dipole interactions based on the FFT for nonequispaced data (NFFT). We present first numerical results for charge-dipole systems, showing that the introduced method can be tuned to a high precision and verifying the O(NlogN) scaling. In order to calculate the interactions with dipoles efficiently, two new variants of the NFFT, namely the Hessian NFFT as well as the adjoint gradient NFFT, are derived and implemented. In the context of NFFT, these new variants are of great importance on their own. The presented particle mesh method is an extension of the particle-particle NFFT (P²NFFT) framework.

Therefore, all the formerly derived P²NFFT features, which cover for instance the treatment of arbitrary combinations of periodic and non-periodic boundary conditions, the handling of triclinic box shapes and a massively parallel implementation, are now also supported for mixed charge-dipole as well as pure dipole systems. The algorithms are publicly available as a part of the ScaFaCoS library.

Key words and phrases: Ewald summation, nonequispaced fast Fourier transform, par- ticle methods, charged particle systems, dipole-dipole interactions, charge-dipole in- teractions, periodic boundary conditions, NFFT, P2NFFT, P3M, ScaFaCoS

2000 AMS Mathematics Subject Classification : 65T

1 Introduction

In this paper we consider electrostatic systems consisting of point charges as well as point dipoles and show how the interactions between these particles can be computed efficiently via the non-

∗michael.hofmann@informatik.tu-chemnitz.de

†franziska.nestler@mathematik.tu-chemnitz.de

‡michael.pippig.tuc@gmail.com

equispaced FFT (NFFT). We present a particle mesh method capable of calculating electrostatic interactions in systems containing a mixture of charges and dipoles with an arithmetic complexity of O(NlogN), where N denotes the number of present particles. In the following we introduce the considered electrostatic quantities.

Given three linearly independent vectors `1,`2,`3 ∈R³, the particles are assumed to be dis- tributed in a box spanned by these three vectors. If periodic boundary conditions are applied in all three dimensions, the underlying space is the three-dimensional torus

T :=`1T×`2T×`3T:={λ1`1+λ2`2+λ3`3|λ1, λ2, λ3∈T}, whereT:=R/Z'[−¹/2,¹/2).

For N_c chargesq_j ∈ Rat positions x_j ∈ T, j = 1, . . . , N_c, as well as N_d dipoles µ_j ∈ R³ at positionsx_j∈ T,j=N_c+ 1, . . . , N_c+N_d, the overall electrostatic energy can be written as

U = 1 2

N

X

j=1

ξjφ(j) (1.1)

with the potentials

φ(j) = X

n∈Z³

0 N

X

i=1

ξ_i

kxij+Lnk (1.2)

= X

n∈Z³ 0 Nc

X

i=1

q_i

kxij+Lnk − X

n∈Z³

0 N

X

i=N_c+1

µ^>_i (x_ij+Ln) kxij+Lnk³ .

Thereby, we define the difference vectorsxij :=xi−xj, the Euclidean norm is denoted by k · k andN :=N_c+N_drepresents the total number of particles. Furthermore, the prime indicates that forn=0all terms withi=j are omitted and the regular matrixL is given byL:= [`<sub>1</sub>,`₂,`₃].

The operatorsξ_j are defined via ξj:=

(qj :j∈ {1, . . . , Nc}, (µ^>_j∇x_j) :j∈ {Nc+ 1, . . . , N}.

The expressions for the energy (1.1) and the potentials (1.2) are the same as for charge-charge systems [14, 28], where the chargesqj have simply been replaced by the operators ξj. Note that the infinite sum in (1.2) is conditionally convergent, provided that we have charge neutrality, i.e.,

N_c

X

j=1

qj= 0. (1.3)

Thus, the value of the energy (1.1) strongly depends on the underlying summation order. Gener- ally, a spherical order of summation is applied, see [28] for more details.

We are also interested in the computation of the fields E(j), energies U(j) and forcesF(j) of the single particles, which are given by

E(j) :=−∇_xjφ(j) = X

n∈Z³ 0

N_c

X

i=1

qixij,n

kxij,nk³ + X

n∈Z³

0 N

X

i=N_c+1

3µ^>_i xij,n

kxij,nk⁵x_ij,n− µ_i kxij,nk³

,

where we setxij,n:=xij+Ln,

U(j) :=ξjφ(j) =

(q_jφ(j) :j∈ {1, . . . , Nc},

−µ^>_jE(j) :j∈ {Nc+ 1, . . . , N}, (1.4)

and

F(j) :=−∇xjU(j) =

(q_jE(j) :j∈ {1, . . . , N_c},

G(j)µ_j :j∈ {Nc+ 1, . . . , N}, (1.5) whereG(j)∈R^3×3 is the negative electric field gradient, i.e.,

G(j) := − ∇_xj∇^>_x

jφ(j)

= X

n∈Z³ 0 Nc

X

i=1

q_i

kxij,nk³I_3,3− 3q_i

kxij,nk⁵x_ij,nx^>_ij,n

+ X

n∈Z³

0 N

X

i=N_c+1

15µ^>_i x_ij,n

kxij,nk⁷ x_ij,nx^>_ij,n− 3

kxij,nk⁵(µ_ix^>_ij,n+x_ij,nµ^>_i )−3µ^>_i x_ij,n kxij,nk⁵I_3,3

. We remark that the computation of sums of a very similar structure is required in order to solve Stokes equations in three dimensions, for which fast algorithms have been proposed recently, see [44, 26].

For the dipole particles, i.e., forj=Nc+ 1, . . . , N, the torques are given by

τ(j) :=µ_j×E(j). (1.6) Note that we can further simplify the expression of the overall energy (1.1) via

U =1 2

Nc

X

j=1

qjφ(j)−1 2

N

X

j=Nc+1

µ^>_jE(j), (1.7) i.e., the overall energy can be obtained easily after having computed the potentials of the charges as well as the fields of the dipoles.

There are already various algorithms for the computation of interactions in charged particle systems, such as the particle-particle particle-mesh (P³M) method [12, 19, 8], the particle-particle NFFT (P²NFFT) [18, 39], further so called particle mesh methods [7, 13], the fast multipole method [17] as well as multigrid-based methods [4], just to mention a few. Note that [2] includes a detailed comparison between different efficient methods for the 3d-periodic Coulomb problem.

The results show that the P³M and P²NFFT solvers rank among the best methods in this field.

Furthermore, for lattice structured systems an approach based on tensor approximation techniques has been proposed recently, see [25].

For systems containing point dipoles the P³M for dipolar systems has already been introduced and applied successfully, see [6, 5]. In this paper we show how the introduced quantities can be approximated efficiently for systems containing charges as well as dipoles. Therefor, the P²NFFT method [39, 33, 34] is extended appropriately. Consequently, we present for the first time an O(NlogN) algorithm for computing electrostatic interactions in charge-dipole systems. We re- mark that the modular structure enables the combination of the extended modules with all former developments, which covers for instance the usage on massively parallel architectures and the application of various types of periodic boundary conditions, i.e.,

φ(j) :=X

n∈S

0 N

X

i=1

ξi

kx_ij+Lnk

withS=Z^p× {0}^3−p and the numberp∈ {0,1,2,3}of dimensions subject to periodic boundary conditions. In the present paper we concentrate on 3d-periodic boundary conditions, i.e.,p= 3, as introduced above. Note that the generalization of the method for the computation of dipole-dipole interactions has already been presented in [31]. We emphasize that the issue of implementing as well as testing the described approach was left open in [31]. In the present work we present

for the first time details on the recently developed implementation as well as numerical results.

We introduce our extension of the method, which enables the treatment of mixed charge-dipole systems, including pure charge and pure dipole systems as special cases.

Our method is based on the Ewald summation approach [14, 28], which we introduce in Section 2.

In Section 3 we discuss the arising root mean square (rms) errors in the forces, which are introduced by truncating the Ewald sums. We revisit the known error estimates for charge-charge [27] as well as for dipole-dipole [46] systems and derive new estimates concerning the interactions between charges and dipoles. The NFFT is briefly described in Section 4. In the charge-charge case we need the NFFT as well as the adjoint NFFT in order to approximate the potentials of the charges and in addition the so called gradient NFFT for the approximation of the acting forces, see [39, 40, 38].

For systems containing both charges and dipoles we need new NFFT modules, which we call Hessian NFFT and adjoint gradient NFFT. These new NFFT modules are described in Section 4, as well. The final algorithm, which is publicly available as a part of the ScaFaCoS library [1], is summarized in Section 5, while its implementation is described in Section 6. In Section 7 we present new numerical results, verifying the correctness of the presented error estimates and showing that the method can be tuned to a high precision. Furthermore, we demonstrate the O(NlogN) scaling with the help of a proper numerical test. A summary is given in Section 8.

2 Ewald summation

Numerous methods in the field of particle simulation are based on the so called Ewald summation technique [14]. Thereby, the Ewald splitting

1

r =erfc(αr)

r +erf(αr) r

is used, where α >0 is called Ewald or splitting parameter, erf(x) = 2π^−1/2Rx

0 e^−t2dt is the well known error function and erfc(x) = 1−erf(x) is the complementary error function.

If the splitting is applied to (1.1) by setting r:=kx_ij+Lnkin each summand, the energy U is split into an absolutely converging short range part

U^short=1 2

N

X

j=1

ξjφ^short(j) with φ^short(j) := X

n∈Z³

0 N

X

i=1

ξi

erfc(αkxij+Lnk)

kxij+Lnk (2.1) as well as a long range partU^longincluding all erf terms.

A transform of the long range part into Fourier space yieldsU^long=U^F−U^self+U⁰ with the following three sums. First, the Fourier sum is given by

U^F=1 2

N

X

j=1

ξ_jφ^F(j) with φ^F(j) := 1 πV

X

k∈Z³

ψ(k)ˆ

N

X

i=1

ξ_ie^{2πik>L−1xi}

!

e^{−2πik>L−1xj}, (2.2)

where V :=|det(L)| denotes the volume of the simulation box and the Fourier coefficients ˆψ(k) are given by

ψ(k) :=ˆ







e^{−π2kL−>kk2/α2}

kL^−>kk² :k6=0,

0 :k=0.

(2.3)

Second, the self energy U^self := α

√π

Nc

X

j=1

q²_j+ 2α³ 3√

π

N

X

j=N_c+1

kµ_jk²=: 1 2

Nc

X

j=1

q_jφ^self(j)−1 2

N

X

j=N_c+1

µ^>_j E^self(j) (2.4)

is subtracted in order to correct for the included terms withkxij+Lnk= 0. Finally, the correction termU⁰is defined as

U⁰=1 2

N

X

j=1

ξ_jφ⁰(j) with φ⁰(j) :=−2π 3V

N

X

i=1

ξ_ikx_ijk². (2.5) In summary we have

U =U^short+U^F−U^self+U⁰. For a detailed derivation we refer to [14, 28].

Remark 2.1. The correction term is the only term representing the order of summation as well as the type of the surrounding medium. In general, the prefactor of the correction energy is

− 2π 3V(1 +),

whereis the dielectric constant of the surrounding medium. For= 0 we apply vacuum boundary conditions and end up with (2.5), whereas= +∞corresponds to metallic boundary conditions.

For systems containing only charges the term U⁰ is generally known as the dipole correction term and is more precisely given by

U⁰=− π 3V

Nc

X

i,j=1

qiqjkxijk²= 2π 3V

Nc

X

i=1

qixi

2

,

where the second identity follows from the charge neutrality condition (1.3). The correction term for mixed systems (2.5) is obtained by replacing the chargesqi by the operatorsξi. Expressed in terms of the potentials we obtain

φ⁰(j) =−2π 3V

N

X

i=1

ξikxijk²= 4π

3V D^>_c xj−1 2

Nc

X

i=1

qikxik²−

N

X

i=N_c+1

µ^>_i xi+D^>_dxj

! , where we define the overall dipole moments

D_c :=

Nc

X

i=1

q_ix_i and D_d:=

N

X

i=N_c+1

µ_i. For the fields we obtain

E⁰(j) =−∇xjφ⁰(j) =−4π

3V (D_c+D_d). In summary we obtain for the correction energy

U⁰=1 2

Nc

X

j=1

qjφ⁰(j)−1 2

N

X

j=Nc+1

µ^>_jE⁰(j) = 2π 3V

kDck²+ 2D^>_cDd+kDdk² .

3 Cutoff errors in the Ewald formulas – an extension to charge-dipole systems

In the field of molecular dynamics simulations the root mean square (rms) error in the forces is commonly considered in order to measure the accuracy of an algorithm. We denote the rms force errors for the charges and for the dipoles by

∆F_c:=

v u u t

1 Nc

Nc

X

j=1

kF(j)−F_≈(j)k² (3.1)

and

∆F_d:=

v u u t

1 N_d

N

X

j=N_c+1

kF(j)−F_≈(j)k², (3.2)

respectively, whereF_≈(j) denote approximations of the exact forcesF(j).

Both, the short range part as well as the Fourier space part, are converging exponentially fast and can thus be truncated. In the following we show how the arising rms errors in the forces caused by the truncation of the Ewald sums can be estimated. Thereby, we revisit the estimates for charge-charge as well as dipole-dipole interactions, as already considered in the literature, and present new estimates concerning interactions between charges and dipoles.

The NFFT based approximation of the truncated Fourier space part, as described above, intro- duce further approximation errors, which are not discussed here.

3.1 Rms force error in the short range part

In the following we denote the kernel function of the short range part shortly by f(r) := erfc(αr)

r . (3.3)

Furthermore, we set

xij,n:=xij+Ln and rij,n:=kxij,nk.

Rms force errors for the charges

For the charges, i.e., for allj = 1, . . . , Nc, the short range part of the force F(j) =−qj∇x_jφ(j) can be written as

F^short(j) =F^short_c.c. (j) +F^short_c.d. (j), where the first term

F^short_c.c. (j) =−qj∇x_j N_c

X

i=1

X

n∈Z³

0qif(rij,n) =qj N_c

X

i=1

X

n∈Z³

0qif⁰(rij,n)xij,n

r_ij,n, includes all charge-charge (c.c.) interactions, while the second term

F^short_c.d. (j) =−qj∇xj

N

X

i=N_c+1

X

n∈Z³

µ^>_i ∇xif(r_ij,n)

=qj N

X

i=Nc+1

X

n∈Z³

f⁰⁰(rij,n)µ^>_i xij,n

x_ij,n

r²_ij,n +f⁰(rij,n) µ_i rij,n

−µ^>_i xij,n

x_ij,n r_ij,n³

!

(3.4) includes all contributions originating from charge-dipole (c.d.) interactions.

From now we denote byrcut>0 the near field cutoff radius and approximate the short range part via considering only distances rij,n smaller than rcut. The resulting approximation error, measured in the Euclidean norm, is denoted by

δF_c.c.^short(j) :=

q_j

Nc

X

i=1

q_i X

n∈Z3 rij,n>rcut

0 f⁰(r_ij,n)xij,n

rij,n

.

In the following we present an estimate of the corresponding rms error

∆F_c.c.^short:=

v u u t

1 Nc

Nc

X

j=1

δF_c.c.^short(j)²,

which has already been considered in [27], before extending the theory to charge-dipole interac- tions. Thereby, we will denote by

Q:=

Nc

X

i=1

q²_i the sum over all squared charge values.

The error is estimated as presented in [27] via δF_c.c.^short(j)²≈ q²_jQ

V Z 2π

0

dφ Z π

0

sinθdθ Z ∞

r_cut

f⁰(r)r

r

2

r²dr, where

f⁰(r) =− 2α

√πre^−α2r2−erfc(αr) r² . Thus, we obtain

δF_c.c.^short(j)²≈4πq_j²Q V

Z ∞ r_cut

2α

√πre^−α2r2+erfc(αr) r²

2

r²dr

≈q²_jQrcut

V α² e^−2α2rcut2 2α

√πrcut

+ 1

√παr_cut^{3 2}

, see [27]. Thereby, the asymptotic expansion formula

Z ∞ A

e^−Bx2f(x)dx≈e^−BA2f(A)

2AB if d

dx f(x)

2Bx f(x)∀x≥A, (3.5) is applied two times. In case of the complementary error function we obtain

erfc(x) = 2

√π Z ∞

x

e^−x2dx≈ e^−x2

√πx, (3.6)

which is used first, followed by approximating the outer integral. Of course, this is only allowed provided thatr_cutis chosen not too small.

The corresponding rms error can thus be approximated by

∆F_c.c.^short:=

v u u t

1 Nc

Nc

X

j=1

δF_c.c.^short(j)²≈ Q√ rcut

√NcV αe^−α2rcut2 2α

√πrcut

+ 1

√παr_cut³

. (3.7)

In the following we derive a corresponding estimate concerning the charge-dipole interaction, which is done in an analog manner. Thereby we denote by θ :=^(µ_i,r) the angle between the dipole momentµ_i and the vectorr. We obtain for the single contributions to the error, cf. (3.4),

f⁰⁰(r)µ^>_i rr

r² +f⁰(r)µ_i

r −µ^>_i rr r³

2

=

f⁰⁰(r)kµ_ikcos(θ)r

r +f⁰(r)µ_i

r − kµ_ikcos(θ)r r²

2

=_f00(r) r −^f0_r^(r)2

²

r²kµ_ik²cos²θ+^f0_r^(r)2²kµ_ik²+ 2_f00(r) r −^f0_r^(r)2

_f0(r)

r cos(θ)kµ_ikµ^>_i r

=kµ_ik²_f0(r)²

r² + cos²θ

f⁰⁰(r)²−^f0_r^(r)2²

.

Thus, we can approximate the error measured in theL2-norm by

δF_c.d.^short(j)²≈q_j² V

N

X

i=Nc+1

Z 2π 0

Z π 0

sinθ Z ∞

rcut

r²

f⁰⁰(r)µ^>_i rr

r²+f(r)⁰µ_i

r −µ^>_i rr r³

2

drdθdφ

=2πq²_j V

N

X

i=N_c+1

kµ_ik² Z π

0

sinθ Z ∞

rcut

f⁰(r)²+ cos²θ r²f⁰⁰(r)²−f⁰(r)² drdθ

=2πq²_j V

N

X

i=N_c+1

kµ_ik² Z ∞

rcut

4

3f⁰(r)²+²₃r²f⁰⁰(r)²dr, (3.8) analogously to the approximation ofδF_c.c.^short(j).

The second derivativef⁰⁰ of the function, as given in (3.3), can be written in the form

f⁰⁰(r) =4α³

√πe^−α2r2+ 4α

√πr²e^−α2r2+ 2erfc(αr)

r³ ≈ e^−α2r2

√παr⁴ 4α⁴r⁴+ 4α²r²+ 2 .

Thereby, the approximation is again obtained via applying the asymptotic expansion formula (3.6) of the complementary error function.

By inserting the obtained approximations for the derivatives off into (3.8) and making use of the asymptotic expansion (3.5) we obtain

δF_c.d.^short(j)²≈ 4πq²_jM

3V ·e^−2α2rcut2 4rcutα² · 1

πα²r_cut⁶ 2B_cut² + (4α⁴r_cut⁴ + 2B_cut)² with

Bcut:= 1 + 2α²r_cut² and

M:=

N

X

i=N_c+1

kµ_ik².

The corresponding rms error can be approximated by

∆F_c.d.^short:=

v u u t

1 Nc

N_c

X

j=1

δF_c.d.^short(j)²≈e^−α2rcut2 α²r_cut³

r 2QM 3NcV rcut

(3B_cut² + 8α⁴r_cut⁴ Bcut+ 8α⁸r_cut⁸ ). (3.9)

Assuming that the contributions originating from interactions with charges and those resulting from interactions with dipoles are independent, we may approximate the rms near field error in the forces of the charges via

∆F_c^short≈ q

(∆F_c.c.^short)²+ ∆F_c.d.^short2

. (3.10)

Rms force errors for the dipole particles

For the dipoles, i.e., for all j = Nc + 1, . . . , N, the short range parts of the forces F(j) =

−∇x_jµ^>_j∇x_jφ(j) are given by

F^short(j) =F^short_d.c. (j) +F^short_d.d. (j),

where

F^short_d.c. (j) =− ∇xjµ^>_j ∇xj

Nc

X

i=1

X

n∈Z³

q_if(r_ij,n)

=−

Nc

X

i=1

qi

X

n∈Z³

f⁰⁰(rij,n)µ^>_jxij,n

x_ij,n

r_ij,n² +f⁰(rij,n) µ_j rij,n

−µ^>_jxij,n

r³_ij,n xij,n

! ,

F^short_d.d. (j) =− ∇xjµ^>_j ∇xj

N

X

i=N_c+1

X

n∈Z³

0µ^>_i ∇xif(rij,n)

=−

N

X

i=Nc+1

X

n∈Z³

0C(rij,n) (µ^>_i xij,n)µ_j+ (µ^>_jxij,n)µ_i+ (µ^>_i µ_j)xij,n

+

N

X

i=Nc+1

X

n∈Z³

0D(rij,n)(µ^>_i xij,n)(µ^>_jxij,n)xij,n.

The functionsC andD are given by C(r) :=f⁰⁰(r)

r² −f⁰(r)

r³ = 4α³

√πr²e^−α2r2+ 6α

√πr⁴e^−α2r2+3erfc(αr) r⁵ , D(r) :=−f⁰⁰⁰(r)

r³ + 3f⁰⁰(r)

r⁴ −3f⁰(r) r² =

8α⁵

√πr² + 20α³

√πr⁴ + 30α

√πr⁶

e^−α2r2+15erfc(αr) r⁷ , see [46].

The same steps, which led to the approximation ofδF_c.d.^short(j) and ∆F_c.d.^short, give δF_d.c.^short(j)²≈ kµ_jk²Qe^−2α2rcut2

3V α⁴r_cut⁷ 2B_cut² + (4α⁴r_cut⁴ + 2B_cut² ) for allj=N_c+ 1, . . . , N as well as

∆F_d.c.^short:=

v u u t

1 N_d

N

X

j=Nc+1

δF_d.c.^short(j)²≈e^−α2rcut2 α²r_cut³

r 2QM

3N_dV r_cut(3B_cut² + 8α⁴r_cut⁴ Bcut+ 8α⁸r_cut⁸ ).

(3.11) In the case of dipole-dipole interactions we refer to the derivation of the corresponding rms errors as presented in [46]. The final result reads as

∆F_d.d.^short:=

v u u t

1 Nd

N

X

j=Nc+1

δF_d.d.^short(j)²

≈e^−α2rcut2 M α²r_cut⁴

s 1 r_cutN_dV

13

6 C_cut² + 2

15D_cut² +13

15C_cutD_cut

, (3.12)

where

C_cut:= 4α⁴r_cut⁴ + 6α²r_cut² + 3,

D_cut:= 8α⁶r_cut⁶ + 20α⁴r_cut⁴ + 30α²r_cut² + 15.

Assuming that the contributions originating from interactions with charges and those resulting from interactions with the other dipoles are independent, we may approximate the rms near field error in the forces of the dipoles via

∆F_d^short≈ q

∆F_d.c.^short2

+ ∆F_d.d.^short2

. (3.13)

3.2 Rms force error in the Fourier space part

For someM ∈2N³we define the index setIM by IM :=

−^M₂¹, . . . ,^M₂¹−1 ×. . .×

−^M₂³, . . . ,^M₂³ −1 .

Since the Fourier coefficients (2.3) decrease exponentially fast for growing k we can replace the infinite lattice by a finite lattice

Z³7→ IM, whereM ∈2N³has to be chosen appropriately.

In the following we denote by

˜`<sub>1</sub>,˜`₂,˜`₃

the dual lattice vectors fulfilling ˜`<sub>i</sub>·`_j =δ_ij. The dual lattice vectors are simply the columns of the matrixL^−>.

We obtain

kL^−>kk²=kk1˜`1+k2˜`2+k3˜`3k². Furthermore, we set

nkl:= ˜`kט`l

and choose the cutoff parametersc₁, c₂, c₃∈R+ such that

c1k˜`1k · |cos^(˜`1,n23)|=c2k˜`2k · |cos^(˜`2,n13)|=c3k˜`3k · |cos^(˜`3,n12)|=β, (3.14) whereβ >0. This means that the three heights of the parallelepiped

{x1˜`1+x2˜`2+x3˜`3:x∈ Ic} with

Ic:= [−^c1/2,^c1/2]×[−^c2/2,^c2/2]×[−^c3/2,^c3/2]

are all of length β. Thus, the parallelepiped contains the ball of radius ^β/2. The Fourier space cutoffM ∈2N³ is consequently chosen such that

M ≈c= (c1, c2, c3), (3.15)

i.e., we round to the next even integer in each component.

Remark 3.1. The conditions (3.14) are equivalent to c1

k`1k = c2

k`2k = c3

k`3k =β. (3.16)

As an example, the vectorsn23and`1 are linearly dependent. Thus, we obtain c1k˜`1k · |cos^(˜`1,n23)|=c1k˜`1k · |cos^(˜`1,`1)|=c1

|˜`1·`1| k`₁k = c1

k`₁k.

Example 3.2. We consider a box spanned by the three vectors `1= (1,0,0), `2 = (1,1,0) and

`3= (0,0,1). The dual lattice vectors are given by ˜`1= (1,−1,0), ˜`2= (0,1,0) and ˜`3= (0,0,1).

Forβ:= 8 we obtain

c= (8,8√

2,8) and M := (8,12,8), which ensures

{x∈R³:kxk ≤4} ⊂ {x1˜`1+x2˜`1+x3˜`1,x∈ Ic}.

In contrast, for the lattice vectors`<sub>1</sub>= (1,0,0), `₂ = (0,1,0) and `₃ = (0,0,1) we obtain M = (8,8,8). For a graphical illustration in two dimensions see Figure 3.1.

k2 ∈ I¹² k1∈ I8

r= ^β₂ = 4

1

k2∈ I⁸ k1∈ I8

r= ^β₂ = 4

1

Figure 3.1: The sets{k₁˜`<sub>1</sub>+k<sub>2</sub>˜`₂ :k₁ ∈ I_M1, k₂ ∈ I_M2}for `<sub>1</sub> = (1,0), `₂ = (1,1) (left) and for

`<sub>1</sub>= (1,0),`₂= (0,1) (right).

Rms force errors for the charges

For the charges, i.e., j = 1, . . . , Nc, we decompose the Fourier space part of the force F(j) =

−qj∇x_jφ(j) via

F^F(j) =F^F_c.c.(j) +F^F_c.d.(j), where

F^F_c.c.(j) =−q_j∇_xj πV

X

k∈Z³

ψ(k)ˆ

Nc

X

i=1

q_ie^{2πik>L−1xij} = 2iq_j V

Nc

X

i=1

q_i X

k∈Z³

ψ(k)(Lˆ ^−>k) e^{2πik>L−1xij},

F^F_c.d.(j) =−q_j∇_xj πV

X

k∈Z³

ψ(k)ˆ

N

X

i=N_c+1

µ^>_i ∇_xie^{2πik>L−1xij}

=−4πq_j V

N

X

i=N_c+1

kµ_ik X

k∈Z³

ψ(k)kLˆ ^−>kkcos(γ_i,k)(L^−>k) e^{2πik>L−1xij}, (3.17)

andγi,k :=^(µ_i,L^−>k) denotes the angle between a lattice vector L^−>k and a dipole moment µ_i.

If the infinite sum is replaced by a finite sum over k∈ IM, the Fourier space partF^F_c.c.(j) is approximated via

F^F_c.c.(j)≈F^F,≈_c.c.(j) :=2iq_j V

Nc

X

i=1

q_i X

k∈IM

ψ(k)(Lˆ ^−>k) e^{2πik>L−1xij} resulting in an error

F^F_c.c.(j)−F^F,≈_c.c.(j) =

Nc

X

i=1

qiχ(xij), where

χ(x) := 2i V

X

k∈Z³\IM

ψ(k)(Lˆ ^−>k) e^2πik>L−1x.

Note that the rms force errors in the Fourier space part have already been computed in [27].

However, for the sake of completeness, we sketch the derivation of the error, following the common approach as presented and used in [9, 46], for instance. The result is slightly different from the one given in [27].

At first, we compute the quadratic mean χ²:= 1

V Z

x∈LT

χ(x)^>χ(x) dx¯ = 4 V²

X

k∈Z³\IM

ψ(k)ˆ ²kL^−>kk²= 4 V²

X

k∈Z³\IM

e^{−2π2kL−>kk2/α2} kL^−>kk²