SampTA2017SamplingTheoryandApplications,12thInternationalConferenceJuly3–7,2017,Tallinn,Estonia FranziskaNestlerChemnitzUniversityofTechnologyFacultyofMathematics FastEwaldsummationforelectrostaticsystemswithchargesanddipolesforvarioustypesofperiodicbound

Fast Ewald summation for electrostatic systems with charges and dipoles

Fast Ewald summation for electrostatic systems with charges and dipoles for various types of periodic boundary conditions

Franziska Nestler Chemnitz University of Technology

Faculty of Mathematics

SampTA 2017

Sampling Theory and Applications, 12th International Conference

July 3 – 7, 2017, Tallinn, Estonia

Fast Ewald summation for electrostatic systems with charges and dipoles

1 Introduction: Electrostatic interactions in particle systems

2 Nonequispaced FFT and fast summation

3 P²NFFT for systems with charges and dipoles

4 Numerical results

5 Summary

Introduction: Electrostatic interactions in particle systems

The Coulomb problem

Charged particle system: Let N charges q

at positions

be given.

Are interested in the electrostatic potentials (x

:=

):

φ(j) :=

N

X

i=1 i6=j

q

kxijk

, j = 1, . . . , N.

(N

)?

3d-periodic boundary conditions Choose

:=

,

[

/

,

/

)

and set

φ(j) :=

n∈S N

X

i=1 i6=jifn=0

q

kxij

+ Ln

→

crystals, . . .

→conditional convergence for electrical neutral particle systems

Introduction: Electrostatic interactions in particle systems

The Coulomb problem

Charged particle system: Let N charges q

at positions

be given.

Are interested in the electrostatic potentials (x

:=

):

φ(j) :=

N

X

i=1 i6=j

q

kxijk

, j = 1, . . . , N.

(N

)?

3d-periodic boundary conditions Choose

:=

,

and set

φ(j) :=

n∈S N

X

i=1 i6=jifn=0

q

kxij

+ Ln

→

crystals, . . .

L

→conditional convergence for electrical neutral particle systems

Introduction: Electrostatic interactions in particle systems

The Coulomb problem

Charged particle system: Let N charges q

at positions

be given.

Are interested in the electrostatic potentials (x

:=

):

φ(j) :=

N

X

i=1 i6=j

q

kxijk

, j = 1, . . . , N.

(N

)?

3d-periodic boundary conditions Choose

:=

,

[

/

,

/

)

and set

φ(j) :=

n∈S N

X

i=1 i6=jifn=0

q

kxij

+ Ln

→

crystals, . . .

→conditional convergence for electrical neutral particle systems

Introduction: Electrostatic interactions in particle systems

The Coulomb problem

Charged particle system: Let N charges q

at positions

be given.

Are interested in the electrostatic potentials (x

:=

):

φ(j) :=

N

X

i=1 i6=j

q

kxijk

, j = 1, . . . , N.

(N

)?

3d-periodic boundary conditions Choose

,

[

/

,

/

)

and set

φ(j) :=

n∈S N

X

i=1 i6=jifn=0

q

kxij

+

→

crystals, . . .

→conditional convergence for electrical neutral particle systems

Introduction: Electrostatic interactions in particle systems

Extension to systems with charges and dipoles

Add

•

N

charges q

at positions

, j = 1, . . . , N

,

•

N

dipoles with

at positions

, j = N

+ 1, . . . , N

+ N

. Total number of particles N := N

+ N

.

Replace the charges q

by the operators ξ

:

q

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

Electrostatic potentials:

φ(j) :=

n∈S N

X

i=1 i6=jifn=0

ξi

kxij

+ Ln

, j = 1, . . . , N.

Introduction: Electrostatic interactions in particle systems

Periodic boundary conditions

3d-periodic

S

:=

crystals, . . .

L L

L

1

open boundary conditions

:=

0

1

2d-periodic

S

:=

0

thin liquid films, . . .

L L

1

1d-periodic

S

:=

0

nano channels, . . .

L

1

Introduction: Electrostatic interactions in particle systems

Periodic boundary conditions

3d-periodic

S

:=

crystals, . . .

L L

L

1

open boundary conditions

:=

0

1

2d-periodic

S

:=

0

thin liquid films, . . .

L L

1d-periodic

S

:=

0

nano channels, . . .

L

Franziska Nestler, TU Chemnitz, Faculty of Mathematics 5/17

Introduction: Electrostatic interactions in particle systems

Ewald splitting

Idea of Ewald summation

:

Ewald splitting

1

r =

| {z } long ranged, continuous

+

| {z } singular in 0, short ranged

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10

r

•

erf(x) :=

Rx

0

e

dt (error function),

r→0 erf(αr)

r =√^2απ

•

erfc(x) := 1

erf(x) (complementary error function)

•

α > 0 (scaling parameter)

Introduction: Electrostatic interactions in particle systems

3d-periodic boundary conditions

φ(j) = X

n∈Z³ N

X

i=1 i6=jifn=0

ξ_ierfc(αkxij+Lnk) kx_ij+Lnk + X

n∈Z³ N

X

i=1

ξ_ierf(αkxij+Lnk)

kx_ij+Lnk −φ^self(j)

• short range part:direct evaluation (truncation)

• long range part:

Fourier coefficientsare known analytically.

φ^self(j) = (_2α

√πqj :charges

0 :dipoles

φ^long(j) = X

k∈Z³\{0} N

X

i=1

ξie^{−π2kkk2/(α2L2)}

πLkkk² e^2πik>xij/L

→good choice of truncation parameters:O(N^3/2)

How to compute the long range part even more efficiently?

(N log N)

•

3d-periodic constraints: FFT for nonequispaced data –

•

open and mixed periodic b.c.:

Introduction: Electrostatic interactions in particle systems

3d-periodic boundary conditions

φ(j) = X

n∈Z³ N

X

i=1 i6=jifn=0

ξ_ierfc(αkxij+Lnk) kx_ij+Lnk + X

n∈Z³ N

X

i=1

ξ_ierf(αkxij+Lnk)

kx_ij+Lnk −φ^self(j)

• short range part:direct evaluation (truncation)

• long range part:Fourier coefficientsare known analytically. φ^self(j) = (_2α

√πqj :charges

0 :dipoles

φ^long(j) = X

k∈Z³\{0} N

X

i=1

ξie^{−π2kkk2/(α2L2)}

πLkkk² e^2πik>xij/L

→good choice of truncation parameters:O(N^3/2)

How to compute the long range part even more efficiently?

(N log N)

•

3d-periodic constraints: FFT for nonequispaced data –

•

open and mixed periodic b.c.:

Introduction: Electrostatic interactions in particle systems

3d-periodic boundary conditions

φ(j) = X

n∈Z³ N

X

i=1 i6=jifn=0

ξ_ierfc(αkxij+Lnk) kx_ij+Lnk + X

n∈Z³ N

X

i=1

ξ_ierf(αkxij+Lnk)

kx_ij+Lnk −φ^self(j)

• short range part:direct evaluation (truncation)

• long range part:Fourier coefficientsare known analytically. φ^self(j) = (_2α

√πqj :charges

0 :dipoles

φ^long(j) = X

k∈Z³\{0} N

X

i=1

ξie^{−π2kkk2/(α2L2)}

πLkkk² e^2πik>xij/L

→good choice of truncation parameters:O(N^3/2)

How to compute the long range part even more efficiently?

(N log N )

•

3d-periodic constraints: FFT for nonequispaced data –

•

open and mixed periodic b.c.:

Nonequispaced FFT and fast summation

FFT for nonequispaced data – NFFT

Notation For M

2

set

:=

/

, . . . ,

/

1

.

Torus

:=

/

[

/

,

/

)

NFFT:

f(x

) :=

k∈IM

f ˆ

e

, j = 1, . . . , N

(inverse) FFT: f(j) :=

k∈IM

f ˆ

e

, N :=

= M

adjoint NFFT:

h(k) :=

N

X

j=1

f

e

Complexity:

(

log

+ N )

[Dutt, Rokhlin 1993] [Beylkin 1995] [Potts, Steidl, Tasche 2001] [Greengard, Lee 2004]

Nonequispaced FFT and fast summation

FFT for nonequispaced data – NFFT

Notation For M

2

set

:=

/

, . . . ,

/

1

. Torus

:=

/

[

/

,

/

)

NFFT:

f(x

) :=

k∈IM

f ˆ

e

, j = 1, . . . , N

(inverse) FFT: f(j) :=

k∈IM

f ˆ

e

, N :=

= M

adjoint NFFT:

h(k) :=

N

X

j=1

f

e

Complexity:

(

log

+ N )

[Dutt, Rokhlin 1993] [Beylkin 1995]

[Potts, Steidl, Tasche 2001] [Greengard, Lee 2004]

Nonequispaced FFT and fast summation

FFT for nonequispaced data – NFFT

Notation For M

2

set

:=

/

, . . . ,

/

1

. Torus

:=

/

[

/

,

/

)

NFFT:

f(x

) :=

k∈IM

f ˆ

e

, j = 1, . . . , N

(inverse) FFT: f(j) :=

k∈IM

f ˆ

e

, N :=

= M

adjoint NFFT:

h(k) :=

N

X

j=1

f

e

Complexity:

(

log

+ N )

[Dutt, Rokhlin 1993] [Beylkin 1995]

[Potts, Steidl, Tasche 2001] [Greengard, Lee 2004]

Nonequispaced FFT and fast summation

Further implemented variants (d = 3)

•

Gradient NFFT: Approximate

∇

f(x

) =

k∈IM

2πik f ˆ

e

j = 1, . . . , N.

•

Hessian NFFT: Approximate

f(x

) =

k∈IM

4π

f ˆ

e

j = 1, . . . , N.

•

Adjoint gradient NFFT: For given

approximate

N

X

j=1

f^>_j∇^x

e

j

=

N

X

j=1

2πif

e

.

Complexity:

(

log

+ N )

Nonequispaced FFT and fast summation

Further implemented variants (d = 3)

•

Gradient NFFT: Approximate

∇

f(x

) =

k∈IM

2πik f ˆ

e

j = 1, . . . , N.

•

Hessian NFFT: Approximate

f(x

) =

k∈IM

4π

f ˆ

e

j = 1, . . . , N.

•

Adjoint gradient NFFT: For given

approximate

N

X

j=1

f^>_j∇^x

e

j

=

N

X

j=1

2πif

e

.

Complexity:

(

log

+ N )

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−L L

K(x)

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes. Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−L L

−^h₂ ^h₂ ^h−L K(x)

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes. Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−L L

−^h₂ ^h₂ ^h−L K(x) KB(x)

C^p−1

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes. Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−^h₂ ^h₂

KR(x) C^p−1

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes.

Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−^h₂ ^h₂

KR(x) C^p−1

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes.

Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>

= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

NFFT based fast summation in 1d

Computef(xj) :=

N

X

i=1

ciK(xi−xj) withxj∈[−^L/2,^L/2] ⇒xi−xj∈[−L, L].

EmbedK(·)into a periodic function:

−^h₂ ^h₂

KR(x) C^p−1

→Two-point-Taylor interpolation

periodh >2L. . .

KR(x)≈ X

`∈IM

ˆb`e<sup>2πi`x/h</sup>

Sampleperiodic func.at equispaced nodes.

Use the FFT to approximate the FKˆb`.

f(xj)≈

N

X

i=1

ci

X

`∈IM

ˆb<sub>`</sub>e<sup>2πi`(xi−xj)/h</sup>= X

`∈IM

ˆb_`

N

X

i=1

cie^2πi`xi/h

!

| {z }

adj. NFFT

e^−2πi`xj/h

| {z }

NFFT

Nonequispaced FFT and fast summation

For d ≥ 2 dimensions, radial kernels

Computef(xj) :=

N

X

i=1

ciK(kxi−xjk), assumekxi−xjk ≤L.

−^h2 −L L h

2 h−L ^3h

2 C^p−1

K(x)

KB(x)

→regularization in 1d

→two-point-Taylor interpolation

→modification: claim vanishing derivatives in^h/2

K(kxk) KB(kxk)

−L L

−^h2 h

2 3h

h−L 2

−L L

−^h₂ h

2 Rotation inR^d:

KR(x) :=







K(kxk) :kxk ≤L KB(kxk) :L <kxk ≤^h/2

KB(^h/2) :else

→h-periodization with respect to all dimensions is a smooth function

→approximateKR(·)by ad-variate trigonometric polynomial (FFT) and proceed analogously to 1d

Nonequispaced FFT and fast summation

For d ≥ 2 dimensions, radial kernels

Computef(xj) :=

N

X

i=1

ciK(kxi−xjk), assumekxi−xjk ≤L.

−^h2 −L L h

2 h−L ^3h

2 C^p−1

K(x)

KB(x)

→regularization in 1d

→two-point-Taylor interpolation

→modification: claim vanishing derivatives in^h/2

K(kxk) KB(kxk)

−L L

−^h2 h

2 3h

h−L 2

−L L

−^h2 h

2 Rotation inR^d:

KR(x) :=







K(kxk) :kxk ≤L KB(kxk) :L <kxk ≤^h/2

KB(^h/2) :else

→h-periodization with respect to all dimensions is a smooth function

→approximateKR(·)by ad-variate trigonometric polynomial (FFT) and proceed analogously to 1d

P²NFFT for systems with charges and dipoles

3d-periodic boundary conditions

φ(j) =

n∈Z³ N

X

i=1 i6=jifn=0

ξierfc(αkx_ij+Lnk) kxij+Lnk

+

n∈Z³ N

X

i=1

ξierf(αkx_ij+Lnk)

kxij+Lnk −

φ

(j)

direct evaluation (truncation)

long:

Fourier coefficients are known analytically. ˆ b

:=

,

=

φ

(j)

k∈IM

ˆ b

X

i=1

q

e

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇^xi

e

| {z }

adj. grad. NFFT

!

e

| {z }

NFFT

Computation of the forces:

charges: F(j) =−∇^xjφ(j)·qj →gradient NFFT dipoles: F(j) =−∇^xj∇^>x_jφ(j)·µ_j →Hessian NFFT

[Hofmann, N., Pippig 2016 (preprint)]

P²NFFT for systems with charges and dipoles

3d-periodic boundary conditions

φ(j) =

n∈Z³ N

X

i=1 i6=jifn=0

ξierfc(αkx_ij+Lnk) kxij+Lnk

+

n∈Z³ N

X

i=1

ξierf(αkx_ij+Lnk)

kxij+Lnk −

φ

(j)

direct evaluation (truncation)

long:

Fourier coefficients are known analytically. ˆ b

:=

,

=

φ

(j)

k∈IM

ˆ b

X

i=1

q

e

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i ∇^xi

e

| {z }

adj. grad. NFFT

!

e

| {z }

NFFT

Computation of the forces:

charges: F(j) =−∇^xjφ(j)·qj →gradient NFFT dipoles: F(j) =−∇^xj∇^>x_jφ(j)·µ_j →Hessian NFFT

Fast algorithm: particle-particle NFFT,O(NlogN)[Pippig, Potts 2011] [Hofmann, N., Pippig 2016 (preprint)]

P²NFFT for systems with charges and dipoles

3d-periodic boundary conditions

φ(j) =

n∈Z³ N

X

i=1 i6=jifn=0

ξierfc(αkx_ij+Lnk) kxij+Lnk

+

n∈Z³ N

X

i=1

ξierf(αkx_ij+Lnk)

kxij+Lnk −

φ

(j)

direct evaluation (truncation)

long:

Fourier coefficients are known analytically. ˆ b

:=

,

=

φ

(j)

k∈IM

ˆ b

X

i=1

q

e

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i ∇^xi

e

| {z }

adj. grad. NFFT

!

e

| {z }

NFFT

Computation of the forces:

charges: F(j) =−∇^xjφ(j)·qj →gradient NFFT dipoles: F(j) =−∇^xj∇^>x_jφ(j)·µ_j →Hessian NFFT

P²NFFT for systems with charges and dipoles

Open and mixed periodic boundary conditions

φ^long(j) regularization / fast summation

0dp =

N

X

i=1

ξi

erf(αkxijk) kxijk

erf(αkxijk) kx_ijk ≈ X

`∈IM

ˆb`e<sup>2πi`>xij/h</sup> d= 3

1dp

≈ X

k∈IM N

X

i=1

ξie^2πikxij,1/LΘ^p1_k,α(kx˜ijk) Θ^p1_k,α(kx˜ijk)≈ X

`∈IM

ˆbk,`e<sup>2πi`>x˜ij/h</sup> d= 2

2dp

≈ X

k∈IM N

X

i=1

ξ_ie^{2πik>˜xij/L}Θ^p2_k,α(|x_ij,3|) Θ^p2_k,α(|x_ij,3|)≈ X

`∈IM

ˆb<sub>k,`</sub>e<sup>2πi`xij,3/h</sup> d= 1

→Ewald summation for 1d- and 2d-periodic constraints [M. Porto 2000] [Grzybowski, Gwó´zd´z, Bródka 2000]

φ^long(j)≈ X

(∗,∗)∈IM

ˆb(∗,∗) Nc

X

i=1

qie^{2πi(∗,∗)>x˘i}

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇x_ie^{2πi(∗,∗)>x˘i}

| {z }

adj. grad. NFFT

!

e^{−2πi(∗,∗)>x˘j}

| {z }

NFFT

→same structure as for 3d-periodic boundary conditions

* periodic dimensions * non periodic dimensions

P²NFFT for systems with charges and dipoles

Open and mixed periodic boundary conditions

φ^long(j) regularization / fast summation

0dp =

N

X

i=1

ξi

erf(αkxijk) kxijk

erf(αkxijk) kxijk ≈ X

`∈IM

ˆb`e<sup>2πi`>xij/h</sup> d= 3

1dp

≈ X

k∈IM N

X

i=1

ξie^2πikxij,1/LΘ^p1_k,α(kx˜ijk) Θ^p1_k,α(kx˜ijk)≈ X

`∈IM

ˆbk,`e<sup>2πi`>x˜ij/h</sup> d= 2

2dp

≈ X

k∈IM N

X

i=1

ξ_ie^{2πik>˜xij/L}Θ^p2_k,α(|x_ij,3|) Θ^p2_k,α(|x_ij,3|)≈ X

`∈IM

ˆb<sub>k,`</sub>e<sup>2πi`xij,3/h</sup> d= 1

→Ewald summation for 1d- and 2d-periodic constraints [M. Porto 2000] [Grzybowski, Gwó´zd´z, Bródka 2000]

φ^long(j)≈ X

(∗,∗)∈IM

ˆb(∗,∗) Nc

X

i=1

qie^{2πi(∗,∗)>x˘i}

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇x_ie^{2πi(∗,∗)>x˘i}

| {z }

adj. grad. NFFT

!

e^{−2πi(∗,∗)>x˘j}

| {z }

NFFT

→same structure as for 3d-periodic boundary conditions

* periodic dimensions * non periodic dimensions

P²NFFT for systems with charges and dipoles

Open and mixed periodic boundary conditions

φ^long(j) regularization / fast summation

0dp =

N

X

i=1

ξi

erf(αkxijk) kxijk

erf(αkxijk) kxijk ≈ X

`∈IM

ˆb`e<sup>2πi`>xij/h</sup> d= 3

1dp ≈ X

k∈IM N

X

i=1

ξie^2πikxij,1/LΘ^p1_k,α(kx˜ijk)

Θ^p1_k,α(kx˜ijk)≈ X

`∈IM

ˆbk,`e<sup>2πi`>x˜ij/h</sup> d= 2

2dp ≈ X

k∈IM N

X

i=1

ξ_ie^{2πik>˜xij/L}Θ^p2_k,α(|x_ij,3|)

Θ^p2_k,α(|x_ij,3|)≈ X

`∈IM

ˆb<sub>k,`</sub>e<sup>2πi`xij,3/h</sup> d= 1

→Ewald summation for 1d- and 2d-periodic constraints [M. Porto 2000] [Grzybowski, Gwó´zd´z, Bródka 2000]

φ^long(j)≈ X

(∗,∗)∈IM

ˆb(∗,∗) Nc

X

i=1

qie^{2πi(∗,∗)>x˘i}

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇x_ie^{2πi(∗,∗)>x˘i}

| {z }

adj. grad. NFFT

!

e^{−2πi(∗,∗)>x˘j}

| {z }

NFFT

→same structure as for 3d-periodic boundary conditions

* periodic dimensions

* non periodic dimensions

P²NFFT for systems with charges and dipoles

Open and mixed periodic boundary conditions

φ^long(j) regularization / fast summation

0dp =

N

X

i=1

ξi

erf(αkxijk) kxijk

erf(αkxijk) kxijk ≈ X

`∈IM

ˆb`e<sup>2πi`>xij/h</sup> d= 3

1dp ≈ X

k∈IM N

X

i=1

ξie^2πikxij,1/LΘ^p1_k,α(kx˜ijk) Θ^p1_k,α(kx˜ijk)≈ X

`∈IM

ˆbk,`e<sup>2πi`>x˜ij/h</sup> d= 2

2dp ≈ X

k∈IM N

X

i=1

ξ_ie^{2πik>˜xij/L}Θ^p2_k,α(|x_ij,3|) Θ^p2_k,α(|x_ij,3|)≈ X

`∈IM

ˆb<sub>k,`</sub>e<sup>2πi`xij,3/h</sup> d= 1

→Ewald summation for 1d- and 2d-periodic constraints [M. Porto 2000] [Grzybowski, Gwó´zd´z, Bródka 2000]

φ^long(j)≈ X

(∗,∗)∈IM

ˆb(∗,∗) Nc

X

i=1

qie^{2πi(∗,∗)>x˘i}

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇x_ie^{2πi(∗,∗)>x˘i}

| {z }

adj. grad. NFFT

!

e^{−2πi(∗,∗)>x˘j}

| {z }

NFFT

→same structure as for 3d-periodic boundary conditions

P²NFFT for systems with charges and dipoles

Open and mixed periodic boundary conditions

φ^long(j) regularization / fast summation

0dp =

N

X

i=1

ξi

erf(αkxijk) kxijk

erf(αkxijk) kxijk ≈ X

`∈IM

ˆb`e<sup>2πi`>xij/h</sup> d= 3

1dp ≈ X

k∈IM N

X

i=1

ξie^2πikxij,1/LΘ^p1_k,α(kx˜ijk) Θ^p1_k,α(kx˜ijk)≈ X

`∈IM

ˆbk,`e<sup>2πi`>x˜ij/h</sup> d= 2

2dp ≈ X

k∈IM N

X

i=1

ξ_ie^{2πik>˜xij/L}Θ^p2_k,α(|x_ij,3|) Θ^p2_k,α(|x_ij,3|)≈ X

`∈IM

ˆb<sub>k,`</sub>e<sup>2πi`xij,3/h</sup> d= 1

→Ewald summation for 1d- and 2d-periodic constraints [M. Porto 2000] [Grzybowski, Gwó´zd´z, Bródka 2000]

φ^long(j)≈ X

(∗,∗)∈IM

ˆb_(∗,∗)

Nc

X

i=1

qie^{2πi(∗,∗)>x˘i}

| {z }

adj. NFFT

+

N

X

i=Nc+1

µ^>_i∇x_ie^{2πi(∗,∗)>˘xi}

| {z }

adj. grad. NFFT

!

e^{−2πi(∗,∗)>x˘j}

| {z }

NFFT

→same structure as for 3d-periodic boundary conditions

* periodic dimensions * non periodic dimensions

Numerical results

Error estimates for 3d-periodic boundary conditions

Root mean square (rms) errors in the forces

∆F:=

v u u t 1 N

N

X

j=1

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

Two types of errors:

1

Truncation of Ewald sums

X

charge-charge interactions (3d-p.)

dipole-dipole interactions (3d-p.)

X

extension to charge-dipole systems (3d-p.)

2

NFFT based approximation errors

Numerical results

Error estimates for 3d-periodic boundary conditions

Root mean square (rms) errors in the forces

∆F:=

v u u t 1 N

N

X

j=1

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

Two types of errors:

1 Truncation of Ewald sums ^Fshort(j)≈F^short_trunc.(j),F^long(j)≈F^long_trunc.(j)

X

charge-charge interactions (3d-p.)

X

dipole-dipole interactions (3d-p.)

X

extension to charge-dipole systems (3d-p.)

2

NFFT based approximation errors

Numerical results

Truncation errors in the Ewald sums (3d-periodic)

Based on existing estimates: Ewald truncation errors are predicted accurately.

charges:

0.5 1 1.5 2

10⁻¹⁴ 10⁻¹¹ 10⁻⁸ 10⁻⁵ 10⁻²

r_cut

=3.5 rcut

=4.0 r_cut

=4.5 rcut

= 5.0 M= 16

M=24 M=32

α

rmsforceerror

dipoles:

0.5 1 1.5 2

10⁻¹⁴ 10⁻¹¹ 10⁻⁸ 10⁻⁵ 10⁻²

rcut

=3.5 r_cut

=4.0 rcut

=4.5 r_cut

= 5.0 M= 16

M=24 M=32

α

rmsforceerror

Figure:

• Solid: Measured rms force errors in a small particle system.Nc=Nd= 50, L= 10.

• Dotted: Predicted rms force errors in the short range part.

• Dasheded: Predicted rms force errors in the Fourier space part.

Numerical results

Approximation errors: P

NFFT method

• Usage of NFFT algorithms introduces further approximation errors.

• Choose NFFT parameters such that additional errors are sufficiently small.

→Error analysis: Quite well understood for 3d-periodic boundary conditions.

• Use the same parameters for open and mixed periodic boundary conditions.

• Mesh size: apply the same resolution in per. and non per. dimensions→M^per L

=! M^np h =:β.

0.5 1 1.5 2

10⁻⁷ 10⁻⁶ 10⁻⁵ 10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰

rcut= 3 rcut= 4 rcut= 5

β= 2

β= 3

β= 4

α

rmsforceerror

charges

0.5 1 1.5 2

10⁻⁷ 10⁻⁶ 10⁻⁵ 10⁻⁴ 10⁻³ 10⁻² 10⁻¹

10⁰ rcut= 3

rcut= 4 rcut= 5

β= 2

β= 3 β= 4

α

dipoles

system:Nc=N_d= 100,L= 10.

3d-periodic

,0d-periodic