NFFT basierte schnelle Ewald-Summation für gemischt periodische Randbedingungen

(1)

Fast summation Coulomb interactions Fast Ewald Results Conclusion

NFFT basierte schnelle Ewald-Summation für gemischt periodische Randbedingungen

Franziska Nestler

Fakultät für Mathematik

Mecklenburger Workshop

Approximationsmethoden und schnelle Algorithmen Hasenwinkel, 18. März 2014

(2)

The 3d-NFFT (FFT for nonequispaced data)

Notation

• define the torusT:=R/Z'[−¹/2,¹/2)

• forM ∈2NsetI^M :={−^M/2, . . . ,^M/2−1}³⊂Z³

NFFT: f(xj) := X

k∈I_M

fˆke^2πik·x^j xj∈T³, j= 1, . . . , N

FFT: f(j) := X

k∈I_M

fˆke^2πik·^j^/^M j∈ IM

, N:=|IM|=M³

adjoint NFFT: h(k) :=

N

X

j=1

fje^−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

• forM ∈2NsetI^M :={−^M/2, . . . ,^M/2−1}³⊂Z³

NFFT: f(xj) := X

k∈I_M

fˆke^2πik·x^j xj∈T³, j= 1, . . . , N FFT: f(j) := X

k∈I_M

fˆke^2πik·^j^/^M j∈ IM, N:=|IM|=M³

j=1

∈ I

(5)

The 3d-NFFT (FFT for nonequispaced data)

Notation

• forM ∈2NsetI^M :={−^M/2, . . . ,^M/2−1}³⊂Z³

NFFT: f(xj) := X

k∈I_M

fˆke^2πik·x^j xj∈T³, j= 1, . . . , N FFT: f(j) := X

k∈I_M

fˆke^2πik·^j^/^M j∈ IM, N:=|IM|=M³

adjoint NFFT: h(k) :=

N

X

j=1

fje^−2πik·x^j k∈ I^M

(6)

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

• xi−xj∈[−L, L]

• constructKB, claim smoothness inx=±L

• regularized kernelKRis smooth and periodic with period h:= 2(L+δ)

f(xj) =

N

X

i=1

ciKR(xi−xj)

≈

N

X

i=1

ci M/2−1

X

l=−^M/2

ˆb_le^2πil(xⁱ^−x^j^)/h =

NFFT

z }| {

M/2−1

X

l=−^M/2

ˆb_l

adj. NFFT

z }| {

N

X

i=1

cie^2πilxⁱ^/h

!

e^−2πilx^j^/h

(7)

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

• extend interval at the boundaries

f(xj) =

N

X

i=1

ciKR(xi−xj)

≈

N

X

i=1

ci M/2−1

X

l=−^M/2

NFFT

z }| {

M/2−1

X

l=−^M/2

ˆb_l

adj. NFFT

z }| {

N

X

i=1

cie^2πilxⁱ^/h

!

e^−2πilx^j^/h

(8)

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

i=1

≈

N

X

i=1

ci M/2−1

X

l=−^M/2

M/2−1

X

l=−^M/2

ˆb_l

adj. NFFT

z }| {

N

X

i=1

cie^2πilxⁱ^/h

!

e^−2πilx^j^/h

(9)

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

f(xj) =

N

X

i=1

ciKR(xi−xj)

≈

N

X

i=1

ci M/2−1

X

l=−^M/2

ˆb_le^2πil(xⁱ^−x^j^)/h

=

NFFT

z }| {

M/2−1

X

l=−^M/2

ˆb_l

adj. NFFT

z }| {

N

X

i=1

cie^2πilxⁱ^/h

!

e^−2πilx^j^/h

(10)

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

f(xj) =

N

X

i=1

ciKR(xi−xj)

≈

N

X

i=1

ci M/2−1

X

l=−^M/2

NFFT

z }| {

M/2−1

X

l=−^M/2

ˆb_l

adj. NFFT

z }| {

N

X

i=1

cie^2πilxⁱ^/h

!

e^−2πilx^j^/h

(11)

Definition of the Coulomb interaction energy

LetN chargesqj∈Rat positionsxj∈R³ be given.

Coulomb interaction energy (xij:=xi−xj):

E:=1 2

N

X

i,j=1 i6=j

qiqj

kxijk=1 2

N

X

j=1

qjφ(xj) with φ(xj) :=

N

X

i=1 i6=j

qi

kxijk.

3d-periodic boundary conditions chooseS:=Z³,xj∈BT³ and set

φ(xj) :=φS(xj) :=X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

→crystals, . . .

Fast NFFT based algorithm: P²NFFT,O(NlogN)[Pippig, Potts 2011]

(12)

Definition of the Coulomb interaction energy

E:=1 2

N

X

i,j=1 i6=j

qiqj

kxijk=1 2

N

X

j=1

N

X

i=1 i6=j

qi

kxijk.

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

→crystals, . . .

B

1

(13)

Definition of the Coulomb interaction energy

E:=1 2

N

X

i,j=1 i6=j

qiqj

kxijk=1 2

N

X

j=1

N

X

i=1 i6=j

qi

kxijk.

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

→crystals, . . .

1

(14)

Definition of the Coulomb interaction energy

E:=1 2

N

X

i,j=1 i6=j

qiqj

kxijk=1 2

N

X

j=1

N

X

i=1 i6=j

qi

kxijk.

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

→crystals, . . .

1

(15)

Coulomb energy s.t. periodic boundary conditions E(S) :=1

2

N

X

j=1

qjφS(xj) with φS(xj) := X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

Mixed periodic boundary conditions:

2d-periodic

xj∈BT²×R,S:=Z²× {0}

→thin liquid films, . . .

B

1

1d-periodic

xj∈BT×R²,S:=Z× {0}²

→nano channels, . . .

B

1

AssumePN

j=1qj= 0⇒conditional convergence

(16)

2

N

X

j=1

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

2d-periodic

xj∈BT²×R,S:=Z²× {0}

→thin liquid films, . . .

B

1

1d-periodic

xj∈BT×R²,S:=Z× {0}²

→nano channels, . . .

B

1

AssumePN

(17)

2

N

X

j=1

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

2d-periodic

xj∈BT²×R,S:=Z²× {0}

→thin liquid films, . . . FMM:O(N) [Greengard, Rokhlin 1987]

MMM2D:O(N^5/3) [Arnold, Holm 2002]

SE2P:O(NlogN) [Lindbo, Tornberg 2011]

1d-periodic

xj∈BT×R²,S:=Z× {0}²

→nano channels, . . . FMM:O(N) [Greengard, Rokhlin 1987]

MMM1D:O(N²) [Arnold, Holm 2005]

AssumePN

(18)

Ewald summation

Idea of the 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)

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

• α >0(scaling parameter)

(19)

Ewald Summation

We apply 1

r =erfc(αr)

r +erf(αr)

r withr:=kxij+Bnkand obtain

X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk

= X

n∈S N

X

i=1 i6=jifn=0

qierfc(αkx_ij+Bnk) kxij+Bnk +

φ^LS(xj) := X

n∈S N

X

i=1

qi

erf(αkxij+Bnk) kxij+Bnk − 2α

√πqj

• short range partcan be obtained via direct evaluation after truncation

• limr→0erf(αr)

r =^√^2α_π ⇒substractself potential

• write thelong range partas a sum in Fourier space

(20)

Ewald Summation

We apply 1

r =erfc(αr)

r +erf(αr)

X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk = X

n∈S N

X

i=1 i6=jifn=0

n∈Si=1 kxij+Bnk − √π

(21)

Ewald Summation

We apply 1

r =erfc(αr)

r +erf(αr)

X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk = X

n∈S N

X

i=1 i6=jifn=0

φ^LS(xj) :=

X

n∈S N

X

i=1

qi

erf(αkxij+Bnk) kxij+Bnk − 2α

√πqj

r = ^√^2α_π ⇒substractself potential

(22)

Ewald Summation

r r r

X

n∈S N

X

i=1 i6=jifn=0

qi

kxij+Bnk = X

n∈S N

X

i=1 i6=jifn=0

φ^LS(xj) := X

n∈S N

X

i=1

qi

erf(αkxij+Bnk) kxij+Bnk

− 2α

√πqj

(23)

Ewald summation s.t. 2d-periodic boundary conditions

Long range part under 2d-periodic b.c. [Grzybowski, Gwóźdź, Bródka 2000]

We writexij= (˜xij, xij,3)and obtain forφ^L_Z2×{0}(xj) 1

2B X

k∈Z²\{0}

N

X

i=1

qie^2πik·˜^x^ij^/BΘ^p2(kkk, xij,3)−2√π B²

N

X

i=1

qiΘ^p2₀ (xij,3)

2d Fourier series k=0part

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5

Θ^p2(k, r) k= 1 k=√

2 k=√

3 . . .

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

e^2πkr/Berfc(αB^πk+αr)+e^−2πkr/Berfc(αB^πk−αr)

k

Θ^p2₀ (r) =_α¹e^−α²^r²+r√ πerf(αr)

(24)

Ewald summation s.t. 2d-periodic boundary conditions

2B X

k∈Z²\{0}

N

X

i=1

N

X

i=1

qiΘ^p2₀ (xij,3)

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5

Θ^p2(k, r) k= 1 k=√

2 k=√

3 . . .

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 1.2 1.4

e^2πkr/Berfc(_αB^πk+αr)+e^−2πkr/Berfc(_αB^πk−αr)

k

(25)

Ewald summation s.t. 2d-periodic boundary conditions

2B X

k∈Z²\{0}

N

X

i=1

N

X

i=1

qiΘ^p2₀ (xij,3)

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5

Θ^p2(k, r) k= 1 k=√

2 k=√

3 . . .

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

e^2πkr/Berfc(_αB^πk+αr)+e^−2πkr/Berfc(_αB^πk−αr)

k

(26)

A fast algorithm I

The 2d Fourier sum

1 2B

P

k∈Z²\{0}e^2πik·˜^x^ij^/BPN

i=1qiΘ^p2(kkk, xij,3)

∀k: Θ^p2(k,·)∈C^∞(R) Regularization:

• xij,3 are absolutely bounded

• we are just interested in Θ^p2(k,·)over a finite interval (white area)

• ε∈(0,¹/2)determines the size of the areas at the boundaries

0

−^h/2+hε ^h/2−hε

−^h/2 h/2

∂^j

∂rjKB(k,h/2−hε) = 1

2B∂^j

∂r^jΘ^p2(k,h/2−hε) 1

2BΘ^p2(k,·)

KB(k,·) KB(k,·)

1

.

| {z }

KR(k,·) :hT→R

Franziska Nestler TU Chemnitz

(27)

A fast algorithm I

The 2d Fourier sum

1 2B

P

∀k: Θ^p2(k,·)∈C^∞(R) Regularization:

• xij,3 are absolutely bounded

• we are just interested in Θ^p2(k,·)over a finite interval (white area)

• construct a smooth and periodic functionKR(k,·)

• ε∈(0,¹/2)determines the size of the areas at the boundaries

0

−^h/2+hε ^h/2−hε

−^h/2 h/2

∂^j

∂rjKB(k,h/2−hε) = 1

2B∂^j

∂r^jΘ^p2(k,h/2−hε) 1

2BΘ^p2(k,·)

KB(k,·) KB(k,·)

1

.

| {z }

KR(k,·) :hT→R

Franziska Nestler TU Chemnitz

(28)

A fast algorithm II

The 2d Fourier sum(∗)

1 2B

P

• truncate the infinite sum,Z²→ IM^˜ :={−^M^˜/2. . . ,^M^˜/2−1}²

1

2BΘ^p2(k, xij,3) =KR(k, xij,3)≈

M/2−1

X

l=−^M/2

ˆbk,le^2πilx^ij,3^/h

• we obtain withxî:= ^x_Bî,1,^x_Bî,2,^xî,3_h

the approximation

• analog for thek=0term

(29)

A fast algorithm II

1 2B

P

• forM ∈2Nlarge enough we find for eachk=kkk 6= 0 (k∈ IM^˜)

1

M/2−1

X

l=−^M/2

the approximation

(∗)≈ X

k∈I_M_˜\{0} M/2−1

X

l=−^M/2

ˆb_kkk,l







N

X

i=1

qie

2πi







k1

k2

l





·ˆxi





 e

−2πi







k1

k2

l





·xˆj

(30)

A fast algorithm II

1 2B

P

1

M/2−1

X

l=−^M/2

the approximation

(∗)≈ X

k∈I_M_˜\{0}

M/2−1

X

l=−^M/2

ˆb_kkk,l







N

X

i=1

qie

2πi







k1

k2

l





·ˆxi





 e

−2πi







k1

k2

l





·xˆj

(31)

A fast algorithm II

1 2B

P

1

M/2−1

X

l=−^M/2

the approximation

(∗)≈ X

k∈I_M_˜\{0}

M/2−1

X

l=−^M/2

ˆb_kkk,l







N

X

i=1

qie

2πi







k1

k2

l





·ˆxi







| {z }

3d adj. NFFT

e

−2πi







k1

k2

l





·xˆj

| {z }

3d NFFT

(32)

A fast algorithm II

1 2B

P

1

M/2−1

X

l=−^M/2

the approximation

(∗)≈ X

k∈I_M_˜\{0}

M/2−1

X

l=−^M/2

ˆb_kkk,l







N

X

i=1

qie

2πi







k1

k2

l





·ˆxi







| {z }

3d adj. NFFT

e

−2πi







k1

k2

l





·xˆj

| {z }

3d NFFT

(33)

Ewald summation s.t. 1d periodic boundary conditions

Long range part under 1d-periodic b.c. [Porto 2000]

We writexij= (xij,1,x˜ij)and obtain forφ^L_Z_×{0}2(xj) 2

B X

k∈Z\{0}

N

X

i=1

qie^2πikx^ij,1^/BΘ^p1(|k|,kx˜ijk)− 1 B

N

X

i=1

qiΘ^p1₀ (kx˜ijk)

1d Fourier series k= 0part

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Θ^p1(k, r) k= 1 k= 2 k= 3 . . .

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Θ^p1(k, r) = Zα

0 1 texp

−^π_B²2^kt²2−r²t²

dt Θ^p1₀ (r) =γ+ Γ(0, α²r²) + ln(α²r²)

(34)

Ewald summation s.t. 1d periodic boundary conditions

Long range part under 1d-periodic b.c. [Porto 2000]

We writexij= (xij,1,x˜ij)and obtain forφ^L_Z_×{0}2(xj) 2

B X

k∈Z\{0}

N

X

i=1

qie^2πikx^ij,1^/BΘ^p1(|k|,kx˜ijk)− 1 B

N

X

i=1

qiΘ^p1₀ (kx˜ijk)

1d Fourier series k= 0part

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Θ^p1(k, r) k= 1 k= 2 k= 3 . . .

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5 2 2.5 3 3.5

Θ^p1(k, r) = Zα

0 1 texp

−^π_B²2^kt²2 −r²t² dt

Θ^p1₀ (r) =γ+ Γ(0, α²r²) + ln(α²r²)

NFFT basierte schnelle Ewald-Summation für gemischt periodische Randbedingungen