14.Oktober2020 FranziskaNestler NFFT-basiertePartikelsimulationfürverschiedenartigeperiodischeRandbedingungen

(1)

NFFT based particle simulation for various types of periodic boundary conditions

NFFT-basierte Partikelsimulation für verschiedenartige periodische Randbedingungen

Franziska Nestler

Seminar: Theorie, Modellierung und Simulation, Institut für Physik, TU Chemnitz

14. Oktober 2020

(2)

NFFT based particle simulation for various types of periodic boundary conditions

Overview

1

Introduction

2

The Ewald summation approach

3

Fourier Approximations

4

Nonequispaced FFT and related methods

5

The particle-particle NFFT (P

²

NFFT) algorithm

6

Numerical Results

7

Summary and discussion

(3)

NFFT based particle simulation for various types of periodic boundary conditions Introduction: What are we interested in?

The Coulomb problem

Charged particle system: Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Are interested in the electrostatic potentials (x

ij

:=

xi

−

xj

):

φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k , j = 1, . . . , N. → O(N

²

)?

3d-periodic boundary conditions Choose S :=

Z³

,

xj

∈ [−

^L

/

2

,

^L

/

2

)

³

and set

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

→conditional convergence for electrical neutral particle systems

(4)

The Coulomb problem

Charged particle system: Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Are interested in the electrostatic potentials (x

ij

:=

xi

−

xj

):

φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k , j = 1, . . . , N. → O(N

²

)?

3d-periodic boundary conditions Choose S :=

Z³

,

xj∈[−^L/2,^L/2)³

and set

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

L

(5)

The Coulomb problem

Charged particle system: Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Are interested in the electrostatic potentials (x

ij

:=

xi

−

xj

):

φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k , j = 1, . . . , N. → O(N

²

)?

3d-periodic boundary conditions Choose S :=

Z³

,

xj

∈ [−

^L

/

2

,

^L

/

2

)

³

and set

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

(6)

The Coulomb problem

Charged particle system: Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Are interested in the electrostatic potentials (x

ij

:=

xi

−

xj

):

φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k , j = 1, . . . , N. → O(N

²

)?

3d-periodic boundary conditions Choose

S:=Z³

,

xj

∈ [−

^L

/

2

,

^L

/

2

)

³

and set

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+

Lnk

→ crystals, . . .

(7)

Extension to systems with charges and dipoles

Add

dipole particles. Given

•

N

c

charges q

j

∈

R

at positions

xj

, j ∈ C,

•

N

d

dipoles with

dipole momentsµ_j∈R³

at positions

xj

, j ∈ D.

Total number of particles N := N

_c

+ N

_d

. Replace the charges q

j

by the operators ξ

j

:

q

j

7→

ξj:=

(qj :j∈ C, µ^>_j∇x_j :j∈ D.

Electrostatic potentials:

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

ξi

kx

ij

+ Lnk , j = 1, . . . , N.

(8)

Periodic boundary conditions

3d-periodic

S :=

Z³

→ crystals, . . .

L L

L

1

open boundary conditions S := {0}

³

1

2d-periodic

S :=

Z²

× {0} → thin liquid ﬁlms, . . .

L L

1

1d-periodic

S :=

Z

× {0}

²

→ nano channels, . . .

L

1

(9)

Periodic boundary conditions

3d-periodic

S :=

Z³

→ crystals, . . .

L L

L

1

open boundary conditions S := {0}

³

1

2d-periodic

S :=

Z²

× {0} → thin liquid ﬁlms, . . .

L L

1d-periodic

S :=

Z

× {0}

²

→ nano channels, . . .

L

(10)

Triclinic structures

Also more general box shapes may be treated:

Insert kx

ij

+

Lnk

with a 3 × 3 matrix

L

and matrix-vector products

Ln.

`1

`2

`3

1

`1

`2

1

`1

1 1

But, for simplicity, we restrict to the cubic case in this talk!

(11)

Preliminary studies and PhD project

•

Dissertation: Oct. 2012 – April 2018

Efﬁcient Computation of Electrostatic Interactions in Particle Systems Based on Nonequispaced Fast Fourier Transforms

•

Starting Point:

• FastO(NlogN)algorithm (P²NFFT) for 3d-periodic charged particle systems [Pippig, Potts 2011]

• Fast summation techniques for open boundary conditions [Nieslony, Potts, Steidl 2004]

•

Goal: Uniﬁed framework / fast algorithm for

• 2d- and 1d-periodic boundary conditions,

• open boundary conditions,

• dipole-dipole interactions,

• charge-dipole interactions.

• Error control and automated parameter tuning (as far as possible).

•

DFG Project

Fast Fourier-based Coulomb solvers for partially periodic boundary conditions, together with Prof. Dr. Christian Holm (Institute for Computational

Physics, University of Stuttgart)

(12)

Preliminary studies and PhD project

•

Dissertation: Oct. 2012 – April 2018

Efﬁcient Computation of Electrostatic Interactions in Particle Systems Based on Nonequispaced Fast Fourier Transforms

•

Starting Point:

• FastO(NlogN)algorithm (P²NFFT) for 3d-periodic charged particle systems [Pippig, Potts 2011]

• Fast summation techniques for open boundary conditions [Nieslony, Potts, Steidl 2004]

•

Goal: Uniﬁed framework / fast algorithm for

• 2d- and 1d-periodic boundary conditions,

• open boundary conditions,

• dipole-dipole interactions,

• charge-dipole interactions.

• Error control and automated parameter tuning (as far as possible).

•

DFG Project

Fast Fourier-based Coulomb solvers for partially periodic boundary conditions, together with Prof. Dr. Christian Holm (Institute for Computational

Physics, University of Stuttgart)

(13)

NFFT based particle simulation for various types of periodic boundary conditions The Ewald summation approach

Ewald splitting

Idea of Ewald summation

[Ewald 1921]

:

Ewald splitting

1 r =

erf(αr) r

| {z } long ranged, continuous

+

erfc(αr) r

| {z } singular in 0, short ranged

0 0.2 0.4 0.6 0.8 1 0

2 4 6 8 10

•

erf(x) :=

^√²_πRx

0

e

^−t²

dt (error function),

lim

r→0 erf(αr)

r = ^√^2α_π

•

erfc(x) := 1 − erf(x) (complementary error function)

•

α > 0 (scaling parameter)

(14)

Ewald splitting

Idea of Ewald summation

[Ewald 1921]

:

Ewald splitting

1 r =

erf(αr) r

| {z } long ranged, continuous

+

erfc(αr) r

| {z } singular in 0, short ranged

0 0.2 0.4 0.6 0.8 1 0

2 4 6 8 10

•

erf(x) :=

^√²_πRx

0

e

^−t²

dt (error function),

lim

r→0 erf(αr)

r = ^√^2α_π

•

erfc(x) := 1 − erf(x) (complementary error function)

•

α > 0 (scaling parameter)

(15)

3d-periodic constraints

φ(j) =

X

n∈Z³ N

X

i=1 i6=jifn=0

ξ

i

kx

ij

+ Lnk

=

X

n∈Z³ N

X

i=1 i6=jifn=0

ξierfc(αkxij+Lnk) kxij+Lnk

+

X

n∈Z³ N

X

i=1

ξierf(αkxij+Lnk)

kxij+Lnk

− φ

^self

(j)

short range part:

direct evaluation (truncation), only consider kx

ij

+ Lnk ≤ r

cut

long range part:

Fourier coefﬁcients

are known an- alytically.

φ

^self

(j) =

(_2α

√π

q

j

: charges 0 : dipoles

φ

^long

(j) =

X

k∈Z³\{0} N

X

i=1

ξ

i

e^−π²^kkk²^/(α²^L²⁾

πLkkk²

e

^2πik^>^x^ij^/L

→ good choice of truncation parameters: O(N

^3/2

)

(16)

3d-periodic constraints

φ(j) =

X

n∈Z³ N

X

i=1 i6=jifn=0

ξ

i

kx

ij

+ Lnk

=

X

n∈Z³ N

X

i=1 i6=jifn=0

ξierfc(αkxij+Lnk) kxij+Lnk

+

X

n∈Z³ N

X

i=1

ξierf(αkxij+Lnk)

kxij+Lnk

− φ

^self

(j)

short range part:

direct evaluation (truncation), only consider kx

ij

+ Lnk ≤ r

cut

long range part:Fourier coefﬁcients

are known an- alytically.

φ

^self

(j) =

(_2α

√π

q

j

: charges 0 : dipoles

φ

^long

(j) =

X

k∈Z³\{0}

N

X

i=1

ξ

i

e^−π²^kkk²^/(α²^L²⁾

πLkkk²

e

^2πik^>^x^ij^/L

→ good choice of truncation parameters: O(N

^3/2

)

(17)

3d-periodic constraints

→ How to compute the long range part even more efﬁciently? Idea: FFT for nonequispaced data (NFFT), particle mesh approach → O(N log N) A fast algorithm is obtained based on the identity e

^y−z

= e

^y

· e

^−z

:

φ

^long,3d

(j) =

X

k∈Z³ N

X

i=1

ξ

iΘ^3d_k,α

e

^2πik^>^x^ij^/L

=

X

k∈Z³

Θ^3d_k,α

N

X

i=1

ξ

i

e

^2πik^>^xⁱ^/L

!

e

^−2πik^>^x^j^/L

(18)

Fourier space representations 3d/2d/1d/0d

• Periodic directions:

Compute analytical Fourier transform (exponential decrease for kkk → ∞ or |k| → ∞).

• Inﬂuence of nonperiodic directions:

Fourier coefﬁcients depend on particle distances regarding the nonperiodic dimensions.

φ

^long,3d

(j) =

X

k∈Z³ N

X

i=1

ξ

iΘ^3d_k,α

e

^2πik^>^x^ij^/L

φ

^long,2d

(j) =

X

k∈Z² N

X

i=1

ξ

i

Θ

^2d_k,α

(|x

3,ij

|)e

^2πik^>^x^˜^ij^/L

,

x

= (˜

x, x3

)

φ

^long,1d

(j) =

X

k∈Z N

X

i=1

ξ

i

Θ

^1d_k,α

(k¯

xij

k)e

^2πik^>^x^1,ij^/L

,

x

= (x

1

,

x)

¯

φ

^long,0d

(j) =

N

X

i=1

ξ

i

Θ

^0d_α

(kx

ij

k) with Θ

^0d_α

(r) = erf(αr)

r

(19)

Fourier space representations 3d/2d/1d/0d

• Periodic directions:

Compute analytical Fourier transform (exponential decrease for kkk → ∞ or |k| → ∞).

• Inﬂuence ofnonperiodic directions:

Fourier coefﬁcients depend on particle distances regarding the nonperiodic dimensions.

φ

^long,3d

(j) =

X

k∈Z³ N

X

i=1

ξ

i

Θ

^3d_k,α

e

^2πik^>^x^ij^/L

φ

^long,2d

(j) =

X

k∈Z² N

X

i=1

ξ

iΘ^2d_k,α(|x3,ij|)e^2πik^>^x^˜^ij^/L

,

x

= (˜

x, x3

)

φ

^long,1d

(j) =

X

k∈Z N

X

i=1

ξ

iΘ^1d_k,α(k¯xijk)e^2πik^>^x^1,ij^/L

,

x

= (x

1

,

x)

¯

φ

^long,0d

(j) =

N

X

i=1

ξ

iΘ^0d_α(kxijk)

with Θ

^0d_α

(r) = erf(αr)

r

(20)

Goal: fast algorithm

3d-periodic:

φ

^long,3d

(j) =

X

k∈Z³

Θ

^3d_k,α

N

X

i=1

ξ

i

e

^2πik^>^xⁱ^/L

!

e

^−2πik^>^x^j^/L

Idea for 2d/1d/0d:

Obtain the same structure by applying appropriate periodization approaches to the symmetric / radial functions

Θ

^2d_k,α

(| · |), Θ

^1d_k,α

(k

·

k) and Θ

^0d_α

(k

·

k).

Approximate by Fourier sums, e.g.:

Θ

^2d_k,α

(| · |) ≈

X

l

ˆ b

_(k,l),α

e

^2πil·/h

for each

k.

(k, l) ∈

Z³

⇒ φ

^long,2d

(j) ≈

X

(k,l)∈Z³

ˆb(k,l),α N

X

i=1

ξ

i

e

^2πi

k l

>

˜ xi /L xi,3/h

!

e

^−2πi

k l

>

˜ xj /L xi,3/h

(21)

Goal: fast algorithm

3d-periodic:

φ

^long,3d

(j) =

X

k∈Z³

Θ^3d_k,α

N

X

i=1

ξ

i

e

^2πik^>^xⁱ^/L

!

e

^−2πik^>^x^j^/L

Idea for 2d/1d/0d:

Obtain the same structure by applying appropriate periodization approaches to the symmetric / radial functions

Θ

^2d_k,α

(| · |), Θ

^1d_k,α

(k

·

k) and Θ

^0d_α

(k

·

k).

Approximate by Fourier sums, e.g.:

Θ

^2d_k,α

(| · |) ≈

X

l

ˆb_(k,l),α

e

^2πil·/h

for each

k. (k, l)∈Z³

⇒ φ

^long,2d

(j) ≈

X

(k,l)∈Z³

ˆb(k,l),α N

X

i=1

ξ

i

e

^2πi

k l

>

˜ xi /L xi,3/h

!

e

^−2πi

k l

>

˜ xj /L xi,3/h

(22)

Involved functions

−10 0 −5 0 5 10 2

4 6

−10 0 −5 0 5 10 3

6 −10 0 −5 0 5 10 15

30 Θ

^2d_k,α

(·), Θ

^1d_k,α

(·):

k

or k large enough.

Type A:

Very fast decay.

Θ

^2d_k,α

(·), Θ

^1d_k,α

(·):

k

or k small; Θ

^0d_α

(·).

Type B:

No decay at all or not fast enough.

The functionsΘ^2d_k,α(·),Θ^1d_k,α(·)andΘ^0d_α(·)are only interesting in the interval[−L, L].

In this example we haveL= 10.

(23)

NFFT based particle simulation for various types of periodic boundary conditions Fourier Approximations

Periodization approaches (nonperiodic dimensions)

One-dimensional setting:

Approximate a given nonperiodic function over[−L, L]by a trig. polynomial.

Different techniques may be applied for functions of type A and B, respectively.

−L L

−^h2 h

2 3h

2 C^∞

−L L

−^h2 h

2 3h

h−L 2 C^p−1

Type A:The function is sufﬁciently small

Type B:Construct an interpolating polynomial

outside a not too large interval[−^h₂,^h₂].

within[L, h−L]that ﬁts the derivatives off

Use analytical Fourier transform:

up to a certain degree.

Poisson:f(x)≈¹_h X

|l|≤^M/2

fˆ

l h

e^2πilx/h

via FFT: f(x)˜ ≈ X

|l|≤^M/2

ˆble^2πilx/h

Higher dimensions:The generalization in case ofType Afunctions is straight forward. The periodization approachType Bis possible in a similar fashion for radial kernelsf=f(kxk).

(24)

Periodization approaches (nonperiodic dimensions)

−L L

−^h2 h

2 3h

2 C^∞

−L L

−^h2 h

2 3h

h−L 2 C^p−1

Type A:The function is sufﬁciently small Type B:Construct an interpolating polynomial outside a not too large interval[−^h₂,^h₂]. within[L, h−L]that ﬁts the derivatives off Use analytical Fourier transform: up to a certain degree.

|l|≤^M/2

fˆ

l h

e^2πilx/h via FFT: f(x)˜ ≈ X

|l|≤^M/2

ˆble^2πilx/h

Higher dimensions:The generalization in case ofType Afunctions is straight forward. The periodization approachType Bis possible in a similar fashion for radial kernelsf=f(kxk).

(25)

Periodization approaches (nonperiodic dimensions)

−L L

−^h2 h

2 3h

2 C^∞

−L L

−^h2 h

2 3h

h−L 2 C^p−1

Type A:The function is sufﬁciently small Type B:Construct an interpolating polynomial outside a not too large interval[−^h₂,^h₂]. within[L, h−L]that ﬁts the derivatives off Use analytical Fourier transform: up to a certain degree.

|l|≤^M/2

fˆ

l h

e^2πilx/h via FFT: f(x)˜ ≈ X

|l|≤^M/2

ˆble^2πilx/h

Higher dimensions:The generalization in case ofType Afunctions is straight forward.

The periodization approachType Bis possible in a similar fashion for radial kernelsf=f(kxk).

(26)

Periodic dimensions:

Exponential decrease forkkk → ∞or|k| → ∞ ⇒truncateFourier series.

Nonperiodic dimensions:

Apply appropriate periodization approaches to the involved functions.

Result: φ^long(j)≈ X

κ∈M

ˆbκ

XN i=1

ξie^2πiκ^>^˘^xⁱ

!

e^−2πiκ^>^x^˘^j

withM ⊂Z³ﬁniteand scaled nodes˘x:=L⁻¹x(periodic),x˘:=h⁻¹x(nonper. directions).

Indices Fourier coefﬁcients Approach

3dp: κ=k ˆbk= Θ^p3_k,α=δk,0e^−π²^kkk²^/(α²^{L2 )}

πLkkk² analytic FT

2dp: κ= (k, l) Computeˆb_k,lvia periodization kkksmall:Type B(1d-FFT) of each functionΘ^p2_k,α(| · |) kkklarge:Type A(analytic) 1dp: κ= (k,l) Computeˆb_k,lvia periodization |k|small:Type B(2d-FFT)

of each functionΘ^p1_k,α(k · k) |k|large:Type A(analytic) 0dp: κ=l ApproximateΘ^p0α(k · k) = ^erf(αk·k)k·k Type B(3d-FFT)

(27)

Exponential decrease forkkk → ∞or|k| → ∞ ⇒truncate Fourier series.

κ∈M

ˆbκ

XN i=1

!

withM ⊂Z³ﬁnite and scaled nodes˘x:=L⁻¹x(periodic),x˘:=h⁻¹x(nonper. directions).

2dp: κ= (k, l) Computeˆb_k,lvia periodization kkksmall:Type B(1d-FFT) of each functionΘ^p2_k,α(| · |) kkklarge:Type A(analytic)

1dp: κ= (k,l) Computeˆb_k,lvia periodization |k|small:Type B(2d-FFT) of each functionΘ^p1_k,α(k · k) |k|large:Type A(analytic) 0dp: κ=l ApproximateΘ^p0α(k · k) = ^erf(αk·k)k·k Type B(3d-FFT)

(28)

κ∈M

ˆbκ

XN i=1

!

of each functionΘ^p1_k,α(k · k) |k|large:Type A(analytic)

0dp: κ=l ApproximateΘ^p0α(k · k) = ^erf(αk·k)k·k Type B(3d-FFT)

(29)

κ∈M

ˆbκ

XN i=1

!

of each functionΘ^p1_k,α(k · k) |k|large:Type A(analytic) 0dp: κ=l ApproximateΘ^p0α(k · k) = ^erf(αk·k)k·k Type B(3d-FFT)

(30)

NFFT based particle simulation for various types of periodic boundary conditions Tools: How do we get a fast algorithm?

FFT for nonequispaced data – NFFT

Notation For M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

^d

⊂

Z^d

.

Torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^x^j xj

∈

T^d

, j = 1, . . . , N

(inverse) FFT: f (j) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^j^/^M j

∈

IM

, N := |I

M

| = M

^d

adjoint NFFT:

h(k) :=

N

X

j=1

f

j

e

^2πik^>^x^j k

∈ I

M

Complexity: O(|I

M

| log |I

M

| + N )

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

(31)

FFT for nonequispaced data – NFFT

Notation For M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

^d

⊂

Z^d

. Torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^x^j xj

∈

T^d

, j = 1, . . . , N

(inverse) FFT: f (j) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^j^/^M j

∈ I

M

, N := |I

M

| = M

^d

adjoint NFFT:

h(k) :=

N

X

j=1

f

j

e

^2πik^>^x^j k

∈ I

M

Complexity: O(|I

M

| log |I

M

| + N )

[Dutt, Rokhlin 1993] [Beylkin 1995]

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

(32)

FFT for nonequispaced data – NFFT

Notation For M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

^d

⊂

Z^d

. Torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^x^j xj

∈

T^d

, j = 1, . . . , N

(inverse) FFT: f (j) :=

X

k∈I_M

f ˆ

k

e

^−2πik^>^j^/^M j

∈ I

M

, N := |I

M

| = M

^d

adjoint NFFT:

h(k) :=

N

X

j=1

f

j

e

^2πik^>^x^j k

∈ I

M

Complexity: O(|I

M

| log |I

M

| + N )

[Dutt, Rokhlin 1993] [Beylkin 1995]

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

(33)

NFFT: idea and workﬂow

Trace problem back to equidistant case via convolution with a window function.

inverse FFT:

gl

:= 1 M

^d

X

k∈I_M

1 c

k

( ˜ ϕ)

f ˆ

k

e

^2πik^>^l/M

(l ∈ I

M

)

Convolution with window function:

f(x

j

) ≈

X

l∈I_M

gl

ϕ ˜

xj

−

_M^l

−

¹2

1

0

2 1

M

−_M^m,^m_M

˜ ϕ(x)

−

¹2

1

0

2

(34)

NFFT: idea and workﬂow

Trace problem back to equidistant case via convolution with a window function.

inverse FFT:

gl

:= 1 M

^d

X

k∈I_M

1 ck( ˜ϕ)

f ˆ

k

e

^2πik^>^l/M

(l ∈ I

M

)

Convolution with window function:

f(xj)≈ X

l∈I_M

gl

ϕ ˜

xj

−

_M^l

−

¹2

1

0

2 1

M

−_M^m,^m_M

˜ ϕ(x)

−

¹2

1

0

2

(35)

NFFT: various parameters

f(xj)

≈

X

l∈I_M

gl

ϕ ˜

xj

−

_M^l

−

¹2

1

0

2 1

M

−^m_M,_M^m

˜ ϕ(x)

−

¹2

1

0

2

Parameters:

•

window function ϕ (B-spline, a Bessel function, Gaussian)

→ compactly supported or truncated (only a few non-zero summands)

•

support parameter m ∈

N

(width of window)

•

oversampling factor σ ≥ 1

(36)

NFFT: various parameters

f(xj)

≈

X

l∈I_σM

gl

ϕ ˜

xj

−

_σM^l

−

¹2

1

0

2 1

σM

−_σM^m ,_σM^m

˜ ϕ(x)

−

¹2

1

0

2

Parameters:

•

window function ϕ (B-spline, a Bessel function, Gaussian)

→ compactly supported or truncated (only a few non-zero summands)

•

support parameter m ∈

N

(width of window)

•

oversampling factor σ ≥ 1

(37)

NFFT based particle simulation for various types of periodic boundary conditions P²NFFT for systems with charges and dipoles

3d/2d/1d/0d-periodic boundary conditions

φ(j) =

X

n∈S N

X

i=1 i6=jifn=0

ξi

erfc(αkxij+Lnk) kxij+Lnk

+

X

n∈S N

X

i=1

ξi

erf(αkxij+Lnk)

kxij+Lnk

− φ

^self

(j)

short:

Direct via truncation, only consider kx

ij

+ Lnk ≤ r

cut

.

long:

Combine different NFFT variants.

φ

^long

(j) ≈

X

κ∈I_M

ˆ b

k

X

i∈C

q

i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. NFFT

+

X

i∈D

µ^>_i

∇

x_i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. grad. NFFT

!

e

^−2πiκ^>^x^˘^j

| {z }

NFFT

j ∈ C ∪ D

Computation of the forces:

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

[Hofmann, N., Pippig 2017]

(38)

3d/2d/1d/0d-periodic boundary conditions

φ(j) =

X

n∈S N

X

i=1 i6=jifn=0

ξi

+

X

n∈S N

X

i=1

ξi

erf(αkxij+Lnk)

kxij+Lnk

− φ

^self

(j)

short:

Direct via truncation, only consider kx

ij

+ Lnk ≤ r

cut

.

long:

Combine different NFFT variants.

φ

^long

(j) ≈

X

κ∈I_M

ˆ b

k

X

i∈C

q

i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. NFFT

+

X

i∈D

µ^>_i

∇

x_i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. grad. NFFT

!

e

^−2πiκ^>^x^˘^j

| {z }

NFFT

j ∈ C ∪ D

Computation of the forces:

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

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

(39)

3d/2d/1d/0d-periodic boundary conditions

φ(j) =

X

n∈S N

X

i=1 i6=jifn=0

ξi

+

X

n∈S N

X

i=1

ξi

erf(αkxij+Lnk)

kxij+Lnk

− φ

^self

(j)

short:

Direct via truncation, only consider kx

ij

+ Lnk ≤ r

cut

.

long:

Combine different NFFT variants.

φ

^long

(j) ≈

X

κ∈I_M

ˆ b

k

X

i∈C

q

i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. NFFT

+

X

i∈D

µ^>_i

∇

x_i

e

^2πiκ^>^x^˘ⁱ

| {z }

adj. grad. NFFT

!

e

^−2πiκ^>^x^˘^j

| {z }

NFFT

j ∈ C ∪ D

Computation of the forces:

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

(40)

Visualization of the P

²

NFFT workﬂow

Precomputations

Run P²NFFT Precompute Fourier coefficients

Structure factors

Direct computations

Long range interactions 3d-periodic

ˆ c(κ):=ψ^p3(k)

analytic FT

2d-periodic ˆ c(κ):=ψ^p2(k,l)

1d-FFT∀k

1d-periodic ˆ c(κ):=ψ^p1(k,l)

2d-FFT∀k

0d-periodic ˆ c(κ):=ψ^p0(l)

3d-FFT

Adj. NFFT Sc(κ)

Adj.∇-NFFT Sd(κ)

Multiplication:

fˆκ:=c(κ) · [Sˆ c(κ)+Sd(κ)]

Near field module Self interactions

Correction terms NFFT(fˆκ)

→Potentials ∇-NFFT(fˆκ)

→Fields ∇∇^>-NFFT(fˆκ)

→Field gradients

MD-Simulations Compute energies, forces and torques.

Time evolution: simulate motion of particles.

Get new particle positions.

New time step.

Size of simulation box does not change.

New time step.

Simulation box has to be increased.

(41)

NFFT based particle simulation for various types of periodic boundary conditions Numerical results

Error estimates for 3d-periodic boundary conditions

Root mean square (rms) errors in the forces

∆F :=

vu ut1

N XN j=1

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

Two types of errors:

1

Truncation of Ewald sums

^F^short(j)≈F^short_trunc.(j),F^long(j)≈F^long_trunc.(j)

X

charge-charge interactions (3d-p.)

[Kolafa, Perram 1992] X

dipole-dipole interactions (3d-p.)

[Wang, Holm 2001]

X

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

2

NFFT based approximation errors

^F^longtrunc.(j)≈F^long_NFFT(j)

(42)

Error estimates for 3d-periodic boundary conditions

Root mean square (rms) errors in the forces

∆F :=

vu ut1

N XN j=1

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

Two types of errors:

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

X

charge-charge interactions (3d-p.)

[Kolafa, Perram 1992]

X

dipole-dipole interactions (3d-p.)

[Wang, Holm 2001]

X

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

2

NFFT based approximation errors

^F^longtrunc.(j)≈F^long_NFFT(j)

(43)

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⁻²

r_cut

=3.5 r_cut

= 4.0 r_cut

= 4.5 rcut

= 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.

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

(44)

Approximation errors: P

²

NFFT method

• Usage of NFFT algorithms introduces further approximation errors.

• Choose NFFT parameters such that additional errors are sufﬁciently small.

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

• Choice of parameters for mixed periodic and open constraints:

based on error estimates for 3d-periodicity (heuristically).

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=Nd= 100,L= 10.

3d-periodic ,0d-periodic

(45)

Summary

P

²

NFFT (particle-particle NFFT)

X O(NlogN)method, based on NFFT [Pippig, Potts 2011]

X 3d-periodic, open, mixed periodic (same structure) [N., Pippig, Potts 2015]

X Extension to charge-dipole systems [Hofmann, N., Pippig 2017]

X Triclinic box shapes [N. 2018, Dissertation]

X Code publicly available github.com/scafacos

X Massively parallel [Pippig 2015]

X Numerical results for charge-dipole systems [N. 2018, Dissertation]

X Error analysis (3d-periodic) [N. 2018, Dissertation]

X Automated parameter tuning (publ. in progress) X Supported in ESPResSo

Thank you for your attention!

(46)

Summary

P

²

NFFT (particle-particle NFFT)

X O(NlogN)method, based on NFFT [Pippig, Potts 2011]

X 3d-periodic, open, mixed periodic (same structure) [N., Pippig, Potts 2015]

X Extension to charge-dipole systems [Hofmann, N., Pippig 2017]

X Triclinic box shapes [N. 2018, Dissertation]

X Code publicly available github.com/scafacos

X Massively parallel [Pippig 2015]

X Numerical results for charge-dipole systems [N. 2018, Dissertation]

X Error analysis (3d-periodic) [N. 2018, Dissertation]

X Automated parameter tuning (publ. in progress) X Supported in ESPResSo

Thank you for your attention!