Particle simulation based on the nonequispaced FFT

(1)

Particle simulation based on the nonequispaced FFT

Franziska Nestler Technische Universität Chemnitz

Faculty of Mathematics joint work with M. Pippig

Strobl16

Time-Frequency Analysis and Related Topics

Strobl, 06/2016

(2)

Particle simulation based on the nonequispaced FFT

1

Introduction: NFFT and fast summation

2

Coulomb interactions and fast Ewald summation

3

Numerical results

4

Extension to systems with dipoles

5

Summary

(3)

Introduction: NFFT and fast summation

The 3d-NFFT (FFT for nonequispaced data)

Notation

•

define the torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

•

for M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

³

⊂

Z³

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik·x^j xj

∈

T³

, j = 1, . . . , N

FFT: f(j) :=

X

k∈I_M

f ˆ

k

e

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

∈ I

M

, N := |I

M

| = M

³

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]

(4)

The 3d-NFFT (FFT for nonequispaced data)

Notation

•

define the torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

•

for M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

³

⊂

Z³

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik·x^j xj

∈

T³

, j = 1, . . . , N

FFT: f(j) :=

X

k∈I_M

f ˆ

k

e

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

∈ I

M

, N := |I

M

| = M

³

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]

(5)

The 3d-NFFT (FFT for nonequispaced data)

Notation

•

define the torus

T

:=

R

/

Z

' [−

¹

/

2

,

¹

/

2

)

•

for M ∈ 2

N

set I

M

:= {−

^M

/

2

, . . . ,

^M

/

2

− 1}

³

⊂

Z³

NFFT:

f(x

j

) :=

X

k∈I_M

f ˆ

k

e

^−2πik·x^j xj

∈

T³

, j = 1, . . . , N

FFT: f(j) :=

X

k∈I_M

f ˆ

k

e

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

∈ I

M

, N := |I

M

| = M

³

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]

(6)

NFFT workflow

1

Deconvolution in Fourier space: O(|I

M

|)

ˆ g

k

:=

( _f_ˆ

k

c_k( ˜ϕ)

:

k

∈ I

M

,

0 : else.

2

Inverse FFT: O(|I

σM

| log |I

σM

|) g

l

:= 1

|I

σM

|

X

k∈I_σM

ˆ

g

k

e

−2πik·l/(σM)

,

l

∈ I

σM

.

3

Convolution in spatial domain: O(N )

Approximate f by a sum of translates of a 1-periodic

window-function.

f(x

j

) ≈

X

l∈I_σM

g

l

ϕ(x ˜

j

−

^l

/

σM

), j = 1, . . . , N.

Oversampling factor σ ≥ 1: we can use more coefficients than given Fourier coefficients.

(7)

NFFT workflow

1

Deconvolution in Fourier space: O(|I

M

|)

ˆ g

k

:=

( _f_ˆ

k

c_k( ˜ϕ)

:

k

∈ I

M

,

0 : else.

2

Inverse FFT: O(|I

σM

| log |I

σM

|) g

l

:= 1

|I

σM

|

X

k∈I_σM

ˆ

g

k

e

−2πik·l/(σM)

,

l

∈ I

σM

.

3

Convolution in spatial domain: O(N )

Approximate f by a sum of translates of a 1-periodic

f(x

j

) ≈

X

l∈I_σM

g

l

ϕ(x ˜

j

−

^l

/

σM

), j = 1, . . . , N.

Oversampling factor σ ≥ 1: we can use more coefficients than given Fourier coefficients.

(8)

NFFT workflow

1

Deconvolution in Fourier space: O(|I

M

|)

ˆ g

k

:=

( _f_ˆ

k

c_k( ˜ϕ)

:

k

∈ I

M

,

0 : else.

2

Inverse FFT: O(|I

σM

| log |I

σM

|) g

l

:= 1

|I

σM

|

X

k∈I_σM

ˆ

g

k

e

−2πik·l/(σM)

,

l

∈ I

σM

.

3

Convolution in spatial domain: O(N )

Approximate f by a sum of translates of a 1-periodic

f(x

j

) ≈

X

l∈I_σM

g

l

ϕ(x ˜

j

−

^l

/

σM

), j = 1, . . . , N.

Oversampling factor σ ≥ 1: we can use more coefficients than given Fourier coefficients.

(9)

The 3d-NFFT (FFT for nonequispaced data)

Further implemented variants:

•

Gradient NFFT: Approximate

∇f(x

j

) = −

X

k∈I_M

2πik f ˆ

k

e

^−2πik·x^j

∀ j = 1, . . . , N.

•

Hessian NFFT: Approximate Hf(x

j

) = −

X

k∈I_M

4π

²kk^>

f ˆ

k

e

^−2πik·x^j

∀ j = 1, . . . , N.

•

Adjoint gradient NFFT (f

_j

∈

C³

): Approximate

N

X

j=1

f^>_j

∇

x

e

^2πik·x _x=x

j

=

N

X

j=1

2πif

^>_jk

e

^2πik·x^j

∀

k

∈ I

M

.

Complexity: O(|I

M

| log |I

M

| + N )

X [M. Pippig, PNFFT library, https://github.com/mpip/pfft]

(10)

ik and analytic differentiation

...for the gradient NFFT:

1

Differentiation in Fourier space (ik differentiation):

∇f(x

j

) = −

X

k∈I_M

2πik f ˆ

k

e

^−2πik·x^j

Compute three 3D-FFTs in the second step (one in each dimension).

2

Analytic differentiation: apply gradient to the window function. Modify the last step of the NFFT as follows:

∇f(x

j

) ≈

X

l∈I_σM

g

l

∇ ϕ(x ˜

j

−

^l

/

σM

)

→ only one 3D-FFT has to be computed in the second step.

(11)

ik and analytic differentiation

...for the gradient NFFT:

1

Differentiation in Fourier space (ik differentiation):

∇f(x

j

) = −

X

k∈I_M

2πik f ˆ

k

e

^−2πik·x^j

Compute three 3D-FFTs in the second step (one in each dimension).

2

Analytic differentiation: apply gradient to the window function. Modify the last step of the NFFT as follows:

∇f(x

j

) ≈

X

l∈I_σM

g

l

∇ ϕ(x ˜

j

−

^l

/

σM

)

→ only one 3D-FFT has to be computed in the second step.

(12)

Fast summation based on NFFT in 1d

[Potts, Steidl 2003]

Computef(xj) :=

N

X

i=1

ciK(xi−xj)forxj∈[−^L/2,^L/2], j= 1, . . . , N.

0

−L L

K(x)

1

• Extend the interval at both ends.

• Construct asmooth transition(match derivatives up to some orderp∈N).

→Two-point-Taylor interpolation

• Theregularized kernel functionis smooth and periodic.

(13)

Fast summation based on NFFT in 1d

Computef(xj) :=

N

X

i=1

0

−L L

−(L+δ) L+δ

∂j

∂xjKB(L) =_∂x^∂^jjK(L)

K(x)

KB(x) KB(x)

1

(14)

Fast summation based on NFFT in 1d

Computef(xj) :=

N

X

i=1

0

−L L

−(L+δ) L+δ

K(x)

KR(x)

KB(x) KB(x)

1

(15)

Fast summation based on NFFT in 1d

Computef(xj) :=

N

X

i=1

K(x) =KR(x) ∀x∈[−L, L]

f(xj) =

N

X

i=1

ciK_R(xi−xj)

0

−L L

−(L+δ) L+δ

K(x)

KR(x)

KB(x) KB(x)

1

Far field computations (use FFT and NFFT): K_Ris periodic with period2(L+δ) =:h. f(xj) =

N

X

i=1

ciKR(xi−xj)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`(xⁱ^−x^j^)/h

=X

`∈IM

ˆb` N

X

i=1

cie^2πi`xⁱ^/h

!

| {z }

adj. NFFT

e^−2πi`x^j^/h

| {z }

NFFT

*Use the FFT to approximate the FKˆb`.

(16)

Fast summation based on NFFT in 1d

Computef(xj) :=

N

X

i=1

K(x) =KR(x) ∀x∈[−L, L]

f(xj) =

N

X

i=1

ciK_R(xi−xj)

0

−L L

−(L+δ) L+δ

K(x)

KR(x)

KB(x) KB(x)

1

Far field computations (use FFT and NFFT):

KRis periodic with period2(L+δ) =:h.

f(xj) =

N

X

i=1

ciKR(xi−xj)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`(xⁱ^−x^j^)/h

=X

`∈IM

ˆb` N

X

i=1

cie^2πi`xⁱ^/h

!

| {z }

adj. NFFT

e^−2πi`x^j^/h

| {z }

NFFT

(17)

Fast summation based on NFFT in 1d

Computef(xj) :=

N

X

i=1

K(x) =KR(x) ∀x∈[−L, L]

f(xj) =

N

X

i=1

ciK_R(xi−xj)

0

−L L

−(L+δ) L+δ

K(x)

KR(x)

KB(x) KB(x)

1

Far field computations (use FFT and NFFT):

KRis periodic with period2(L+δ) =:h.

f(xj) =

N

X

i=1

ciKR(xi−xj)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`(xⁱ^−x^j^)/h=X

`∈IM

ˆb` N

X

i=1

cie^2πi`xⁱ^/h

!

| {z }

adj. NFFT

e^−2πi`x^j^/h

| {z }

NFFT

(18)

For d ≥ 2 dimensions, radial kernels

Computef(xj) :=

N

X

i=1

ciK(kxi−xjk)forj= 1, . . . , N. Assumekxi−xjk ≤L.

d^j

dr^jKB(L+δ) = 0 d^j

dr^jKB(L) =_dr^d^jjK(L)

L

−L L+δ

−(L+δ)

K(x)

• 1d regularization on[−^h/2,^h/2]with h:= 2(L+δ)

• claim vanishing derivatives at the boundary

• rotate and extend the function to the torushT² (constant value over the striped area)

• result: periodically smooth function in 2 variables

• approximate by bivariate trigonometric polynomials (2d FFT)

Franziska Nestler TU Chemnitz, Faculty of Mathematics

(19)

Coulomb interactions and fast Ewald summation

Definition of the Coulomb interaction energy

Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Coulomb interaction energy (x

ij

:=

xi

−

xj

):

U := 1 2

N

X

i,j=1 i6=j

q

i

q

j

kx

ij

k = 1 2

N

X

j=1

q

j

φ(j) with φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k .

3d-periodic boundary conditions choose S :=

Z³

,

xj

∈ L

T³

and set

φ(j) := φ

S

(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

Fast NFFT based algorithm: P

²

NFFT, O(N log N)

[Pippig, Potts 2011]

(20)

Definition of the Coulomb interaction energy

Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Coulomb interaction energy (x

ij

:=

xi

−

xj

):

U := 1 2

N

X

i,j=1 i6=j

q

i

q

j

kx

ij

k = 1 2

N

X

j=1

q

j

φ(j) with φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k .

3d-periodic boundary conditions choose S :=

Z³

,

xj∈LT³

and set

φ(j) := φ

S

(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+

Lnk

→ crystals, . . .

L

Fast NFFT based algorithm: P

²

NFFT, O(N log N)

(21)

Definition of the Coulomb interaction energy

Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Coulomb interaction energy (x

ij

:=

xi

−

xj

):

U := 1 2

N

X

i,j=1 i6=j

q

i

q

j

kx

ij

k = 1 2

N

X

j=1

q

j

φ(j) with φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k .

3d-periodic boundary conditions choose

S:=Z³

,

xj

∈ L

T³

and set

φ(j) := φ

S

(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

Fast NFFT based algorithm: P

²

NFFT, O(N log N)

(22)

Definition of the Coulomb interaction energy

Let N charges q

j

∈

R

at positions

xj

∈

R³

be given.

Coulomb interaction energy (x

ij

:=

xi

−

xj

):

U := 1 2

N

X

i,j=1 i6=j

q

i

q

j

kx

ij

k = 1 2

N

X

j=1

q

j

φ(j) with φ(j) :=

N

X

i=1 i6=j

q

i

kx

ij

k .

3d-periodic boundary conditions choose S :=

Z³

,

xj

∈ L

T³

and set

φ(j) := φ

S

(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

→ crystals, . . .

Fast NFFT based algorithm: P

²

NFFT, O(N log N)

(23)

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

r

•

erf(x) :=

^√²_πRx

0

e

^−t²

dt (error function), lim

r→0erf(αr) r

=

^√^2α_π

•

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

•

α > 0 (scaling parameter)

(24)

Open boundary conditions

φ(j) =

N

X

i=1 i6=j

q

i

kx

ij

k =

N

X

i=1 i6=j

qi

erfc(αkxijk) kxijk

+

N

X

i=1

qi

erf(αkxijk) kxijk

− 2α

√ π q

j

• short range part:

direct evaluation (truncation)

• long range part:

use NFFT based fast summation for radial kernels (d = 3) φ

^long

(j) ≈

X

k∈I_M

ˆ b

k N

X

i=1

q

i

e

^2πik·xⁱ^/h

!

| {z }

adj. NFFT

e

^−2πik·x^j^/h

| {z }

NFFT

Computation of the forces:

F

(j) := −q

j

∇φ(j) = −q

j

∇φ

^short

(j)

| {z }

direct

−q

j

∇φ

^long

(j)

| {z }

via gradient NFFT

(25)

Open boundary conditions

φ(j) =

N

X

i=1 i6=j

q

i

kx

ij

k =

N

X

i=1 i6=j

qi

erfc(αkxijk) kxijk

+

N

X

i=1

qi

erf(αkxijk) kxijk

− 2α

√ π q

j

direct evaluation (truncation)

use NFFT based fast summation for radial kernels (d = 3) φ

^long

(j) ≈

X

k∈I_M

ˆ b

k N

X

i=1

q

i

e

^2πik·xⁱ^/h

!

| {z }

adj. NFFT

e

^−2πik·x^j^/h

| {z }

NFFT

Computation of the forces:

F

(j) := −q

j

∇φ(j) = −q

j

∇φ

^short

(j)

| {z }

direct

−q

j

∇φ

^long

(j)

| {z }

via gradient NFFT

(26)

3d-periodic boundary conditions

φ(j) =

X

n∈Z³ N

X

i=1 i6=jifn=0

qierfc(αkxij+Lnk) kxij+Lnk

+

X

n∈Z³ N

X

i=1

qierf(αkxij+Lnk) kxij+Lnk

− 2α

√ π q

j

direct evaluation (truncation)

Fourier coefficients are known analytically.

φ

^long

(j) =

X

k∈Z³

e

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

kkk

²

N

X

i=1

q

i

e

^2πik·xⁱ^/L

!

| {z }

adj. NFFT

e

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

| {z }

NFFT

• coefficients tend to zero exponentially fast

• truncate:Z³7→ IM

(27)

Mixed periodicity

φ(j) :=

X

n∈S N

X

i=1 i6=jifn=0

q

i

kx

ij

+ Lnk

2d-periodic

xj

∈ L

T²

×

R

, S :=

Z²

× {0}

→ thin liquid films, . . .

L

1

•

Compute analytic solution with respect to the 2 periodic dimensions.

•

Use fast summation in 1d for the non-periodic dimension.

1d-periodic

xj

∈ L

T

×

R²

, S :=

Z

× {0}

²

→ nano channels, . . .

L

1

•

Compute analytic solution with respect to the 1 periodic dimension.

•

Use fast summation in 2d for the

non-periodic dimensions.

(28)

Numerical results

Parameter choice

Parameter tuning: quite well understood for the 3d-periodic case.

Choose the NFFT mesh size in the mixed periodic case as follows:

2d-periodc 3d-periodic

periodic dims. non periodic dim. periodic dims.

box length L

h:= 2L+ 2δ L

# grid points M M

3

:= 2M + P M

0

−L L

−(L+δ) L+δ

regularization regularization

periodic

non periodic

P/2 M M P/2

grid points

Franziska Nestler TU Chemnitz, Faculty of Mathematics

(29)

Numerical results

Parameter choice

If the parameters are chosen appropriately, the achieved rms errors are comparable:

0.5 1 1.5 2

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

splitting parameterα

∆FZ2×{0}

M= 16, P= 8, h= 25 M= 32, P= 16, h= 25 M= 64, P= 30, h≈24.67 M= 128, P= 44, h≈23.44 M= 256, P= 76, h≈22.97

0.5 1 1.5 2

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

∆FZ3

Figure:Achieved rms force errors for the 2d-periodic (left) compared to the 3d-periodic (right) case. (L= 10,N= 300)

Large particle systems:

• Tune parameters for a small system.

• Use tuned parameters also for lager systems (same particle density), chooseM³∼N.

• Complexity:O(NlogN).

10³ 10⁴ 10⁵ 10⁶

10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³

#charges

time[s]

∼N

∼NlogN

Figure:Comparison of the computation times (2dp:*, 3dp:4).

(30)

Numerical results

Parameter choice

If the parameters are chosen appropriately, the achieved rms errors are comparable:

0.5 1 1.5 2

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

∆FZ2×{0}

M= 16, P= 8, h= 25 M= 32, P= 16, h= 25 M= 64, P= 30, h≈24.67 M= 128, P= 44, h≈23.44 M= 256, P= 76, h≈22.97

0.5 1 1.5 2

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

∆FZ3

Figure:Achieved rms force errors for the 2d-periodic (left) compared to the 3d-periodic (right) case. (L= 10,N= 300)

Large particle systems:

• Tune parameters for a small system.

• Use tuned parameters also for lager systems (same particle density), chooseM³∼N.

• Complexity:O(NlogN).

10³ 10⁴ 10⁵ 10⁶

10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³

#charges

time[s]

∼N

∼NlogN

Figure:Comparison of the computation times (2dp:*, 3dp:4).

(31)

Extension to systems with dipoles

Systems with charges and dipoles

Given

•

N

1

charges q

j

∈

R

at positions

xj

(j = 1, . . . , N

1

)

•

N

2

dipoles with dipole moments

µ_j

∈

R³

at positions

xj

(j = N

1

+ 1, . . . , N

1

+ N

2

) Replace the charges q

j

by the operators ξ

j

:

q

j

7→ ξ

j

:=

(

q

j

: j ∈ {1, . . . , N

1

},

µ^>_j

∇

x_j

: j ∈ {N

1

+ 1, . . . , N

1

+ N

2

}.

Electrostatic energy and potentials:

U := 1 2

N1+N2

X

j=1

ξ

j

φ(j) with φ(j) :=

X

n∈Z³ N1+N2

X

i=1 i6=jifn=0

ξ

i

kx

ij

+ Lnk .

(32)

Extension to systems with dipoles

NFFT based computation of the long range part

Potentials (for charges and dipoles):

φ

^long

(j) ≈

X

k∈I_M

ˆ b

k







N1

X

i=1

q

i

e

^2πik·xⁱ^/L

| {z }

adj. NFFT

+

N1+N2

X

i=N1+1

µ^>_i

∇

x_i

e

^2πik·xⁱ^/L

| {z }

adj. gradient NFFT







e

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

| {z }

NFFT

Forces of the charges:

F

(j) = −q

j

∇

x_j

φ(j),

gradient NFFT

instead of NFFT.

Forces of the dipoles:

F

(j) = −[∇

x_j

∇

^>x_j

φ(j)] ·

µ_j

F^long

(j) ≈ −



∇xj

∇

^>xj

X

k∈I_M

ˆ b

k

(S

c

(k) + S

d

(k)) e

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





| {z }

Hessian NFFT

·µ

_j

(33)

Summary

P

²

NFFT for charged particle systems

Long range parts of the potentials: same structure for all types of periodic boundary conditions.

[N., Pippig, Potts, J. Comput. Phys. 2015], [N., Pippig, Potts, Comput. Trends in Solvation and Transport in Liquids 2015]

• 3d-periodic: coefficientsˆbkare known analytically (Ewald summation)

• 0d-periodic: NFFT based fast summation (3D)

• 1d-periodic: Ewald summation (1D) and fast summation (2D)

Long range parts of the forces: Use gradient NFFT.

Extended method for computation of interactions with dipoles.[N., Appl. Numer. Math. 2016]

New NFFT modules: Hessian NFFT, adjoint gradient NFFT.[M. Pippig, PNFFT library, https://github.com/mpip/pfft] P²NFFT is publicly available at www.scafacos.de

Numerical results for systems with charges and dipoles: going to be published.

Thank you for your attention!

(34)

Summary

P

²

NFFT for charged particle systems

• 2d-periodic: Ewald summation (2D) and fast summation (1D) Long range parts of the forces: Use gradient NFFT.

New NFFT modules: Hessian NFFT, adjoint gradient NFFT.[M. Pippig, PNFFT library, https://github.com/mpip/pfft] P²NFFT is publicly available at www.scafacos.de

Thank you for your attention!

(35)

Summary

P

²

NFFT for charged particle systems

New NFFT modules: Hessian NFFT, adjoint gradient NFFT.[M. Pippig, PNFFT library, https://github.com/mpip/pfft]

P²NFFT is publicly available at www.scafacos.de

Thank you for your attention!

(36)

Summary

P

²

NFFT for charged particle systems

Thank you for your attention!

(37)

Summary

P

²

NFFT for charged particle systems

Thank you for your attention!