NFFT based Ewald summation for electrostatic systems with charges and dipoles

(1)

NFFT based Ewald summation

for electrostatic systems with charges and dipoles

NFFT based Ewald summation

for electrostatic systems with charges and dipoles

Franziska Nestler Technische Universität Chemnitz

Faculty of Mathematics joint work with M. Pippig and M. Hofmann

Stuttgart, 11/2015

(2)

NFFT based Ewald summation

for electrostatic systems with charges and dipoles

1 Introduction NFFT

Fast summation The Coulomb problem

2 Fast Ewald summation for charged particle systems The 3d-periodic case

Mixed periodic boundary conditions

3 Extension to systems containing charges and dipoles

(3)

Introduction NFFT

The 1d-NFFT (FFT for nonequispaced data)

Notation

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

• forM∈2NsetIM :={−^M/2, . . . ,^M/2−1}

NFFT: f(xj) := X

k∈IM

fˆke^−2πikx^j xj∈T, j= 1, . . . , N

(inverse) FFT: f(j) := X

k∈IM

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

, N:=M

adjoint NFFT: h(k) :=

N

X

j=1

fje^2πikx^j k∈ IM

Complexity:O(MlogM+N)

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

(4)

Introduction NFFT

The 1d-NFFT (FFT for nonequispaced data)

Notation

• forM∈2NsetIM :={−^M/2, . . . ,^M/2−1}

NFFT: f(xj) := X

k∈IM

fˆke^−2πikx^j xj∈T, j= 1, . . . , N (inverse) FFT: f(j) := X

k∈IM

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

N

X

j=1

[Dutt, Rokhlin 1993] [Beylkin 1995]

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

(5)

Introduction NFFT

The 1d-NFFT (FFT for nonequispaced data)

Notation

• forM∈2NsetIM :={−^M/2, . . . ,^M/2−1}

NFFT: f(xj) := X

k∈IM

fˆke^−2πikx^j xj∈T, j= 1, . . . , N (inverse) FFT: f(j) := X

k∈IM

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

N

X

j=1

[Dutt, Rokhlin 1993] [Beylkin 1995]

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

(6)

Introduction NFFT

NFFT window function

Goal: Evaluation of the trigonometric polynomial

f(x) :=

M/2−1

X

k=−^M/2

fˆke^−2πikx

inxj∈T,j= 1, . . . , N.

Idea: Approximatefby

f(x)≈

σM/2−1

X

l=−^σM/2

glϕ(x˜ −^l/σM)

| {z }

discrete convolution

,

where

• ϕ˜is a1-periodic function with small support[−^m/σM,−^m/σM]⊂T(mσM),

• σ≥1is called oversampling factor,

• glare up to now unknown.

(7)

Introduction NFFT

−¹₂ 0 ¹₂

1 σM

−σM^m,_σM^m ϕ(x)

˜ ϕ(x)

−¹₂ 0 ¹₂

1

σM f(x)≈

σMP/2−1

`=−^σM/2

g_`ϕ(x˜ −^`/σM)

ind >1dimensions:tensor product approachϕ(x)7→ϕ(x1)·. . .·ϕ(xd).

(8)

Introduction NFFT

NFFT workflow

1 Deconvolution in Fourier space:O(M)

ˆ gk:=

( _f_ˆ

k

ck( ˜ϕ) :k∈ {−^M/2, . . . ,^M/2−1},

0 :else.

2 Inverse FFT:O(σMlog(σM))

gl:= 1 σM

σM/2−1

X

k=−^σM/2

ˆ

gke−2πikl/(σM), l=−^σM/2, . . . ,^σM/2−1.

3 Convolution in spatial domain:O((2m+ 1)N)

f(xj)≈

σM/2−1

X

l=−^σM/2

glϕ(x˜ j−^l/σM), j= 1, . . . , N.

(9)

Introduction NFFT

NFFT workflow

ˆ gk:=

( _f_ˆ

k

ck( ˜ϕ) :k∈ {−^M/2, . . . ,^M/2−1},

0 :else.

gl:= 1 σM

σM/2−1

X

k=−^σM/2

ˆ

gke−2πikl/(σM), l=−^σM/2, . . . ,^σM/2−1.

f(xj)≈

σM/2−1

X

l=−^σM/2

glϕ(x˜ j−^l/σM), j= 1, . . . , N.

(10)

Introduction NFFT

NFFT workflow

ˆ gk:=

( _f_ˆ

k

ck( ˜ϕ) :k∈ {−^M/2, . . . ,^M/2−1},

0 :else.

gl:= 1 σM

σM/2−1

X

k=−^σM/2

ˆ

gke−2πikl/(σM), l=−^σM/2, . . . ,^σM/2−1.

f(xj)≈

σM/2−1

X

l=−^σM/2

glϕ(x˜ j−^l/σM), j= 1, . . . , N.

(11)

Introduction NFFT

Deconvolution:

ˆ

gk:= ˆdk·fˆk. Standard approach:

dˆk:=





 1

ck( ˜ϕ) :k=−^M/2, . . . ,^M/2−1,

0 :else.

Least square optimized NFFT:

dˆk:=







ck( ˜ϕ) P

r∈Z

c²_k+rσM( ˜ϕ) :k=−^M/2, . . . ,^M/2−1,

0 :else.

(12)

Introduction NFFT

Matrix-vector form

1 Deconvolution in Fourier space:gˆk:= ˆdkfˆk

ˆ

g:=D·f,ˆ D∈R^σM×M : diag. matrix with entriesdˆk 2 Inverse FFT:gˆ7→g

g=F ·ˆg, F ∈C^σM×σM :Fourier matrix

3 Convolution in spatial domain:f(xj)≈P^σM/2−1

l=−^σM/2glϕ(x˜ j−^l/σM) f˜=B·g, B∈R^N×σM, bjl= ˜ϕ(xj−^l/σM)

(13)

Introduction NFFT

Adjoint NFFT

NFFT: Compute for allj= 1, . . . , N fj:=

M/2−1

X

k=−^M/2

fˆke^−2πikx^j ⇐⇒ f:=A·fˆ,

whereA= e^−2πikx^j

j,k∈R^N×M. Approximation:f≈BF D·fˆ.

adjoint NFFT: Compute for allk=−^M/2, . . . ,^M/2−1 hk:=

N

X

j=1

fje^2πikx^j ⇐⇒ h=A^∗·f,

whereA^∗= e^2πikx^j

k,j∈R^M×N. Approximation:h≈DF^∗B^T·f.

(14)

Introduction NFFT

Adjoint NFFT

NFFT: Compute for allj= 1, . . . , N fj:=

M/2−1

X

k=−^M/2

fˆke^−2πikx^j ⇐⇒ f:=A·fˆ,

whereA= e^−2πikx^j

j,k∈R^N×M. Approximation:f≈BF D·fˆ.

adjoint NFFT: Compute for allk=−^M/2, . . . ,^M/2−1 hk:=

N

X

j=1

fje^2πikx^j ⇐⇒ h=A^∗·f,

whereA^∗= e^2πikx^j

k,j∈R^M×N. Approximation:h≈DF^∗B^T·f.

(15)

Introduction Fast summation

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.

(16)

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

(17)

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

(18)

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(xij)

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(xij)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`x^ij^/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`.

(19)

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(xij)

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(xij)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`x^ij^/h

= X

`∈IM

ˆb` N

X

i=1

cie^2πi`xⁱ^/h

!

| {z }

adj. NFFT

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

| {z }

NFFT

(20)

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(xij)

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(xij)≈

N

X

i=1

ci

X

`∈IM

ˆb`e^2πi`x^ij^/h= X

`∈IM

ˆb` N

X

i=1

cie^2πi`xⁱ^/h

!

| {z }

adj. NFFT

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

| {z }

NFFT

(21)

Fast decaying kernels

Computef(xj) :=

N

X

i=1

Assume that

• K(·)is sufficiently small outside of[−^h/2,^h/2]⊇[−L, L].

• we know the analytical Fourier transformK(·).ˆ

Use analytical Fourier coefficients (Poisson summation, no FFT-precomputation).

−^h/2 h/2

K(x)≈ X∞ n=−∞

K(x+nh) =1 h

X∞

`=−∞

K(ˆ^`/^h)e^2πi`x/h K(x)

h-periodization

N 1

X

i=1

ciK(xij)≈ 1 h

X

`∈IM

K(ˆ ^`/h)

N

X

i=1

cie^2πi`xⁱ^/h

!

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

(22)

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)

1

• 1d regularization on[−^h/2,^h/2],^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

(23)

Introduction The Coulomb problem

Coulomb interactions

LetNchargesqj∈Rat positionsxj∈[−^L/2,^L/2]³be given.

Coulomb interaction energy s.t. periodic boundary conditions:

U :=1 2

N

X

j=1

qjφ(j) with φ(j) :=X

n∈S N

X

i=1 0 qi

kxij+Lnk, whereS ⊂Z³.

• 0d-periodic:S:={0}³

• 3d-periodic:S:=Z³

• 2d-periodic:S:=Z²× {0}

• 1d-periodic:S:=Z× {0}²

(24)

Coulomb interactions

U :=1 2

N

X

j=1

n∈S N

X

i=1 0 qi

Franziska Nestler 1 TU Chemnitz, Faculty of Mathematics

(25)

Coulomb interactions

U :=1 2

N

X

j=1

n∈S N

X

i=1 0 qi

L L

L

Franziska Nestler 1 TU Chemnitz, Faculty of Mathematics

(26)

Coulomb interactions

U :=1 2

N

X

j=1

n∈S N

X

i=1 0 qi

L

1

(27)

Coulomb interactions

U :=1 2

N

X

j=1

n∈S N

X

i=1 0 qi

• 1d-periodic:S:=Z× {0}² ^L

1

(28)

Fast Ewald summation for charged particle systems

Ewald summation

Apply the well known Ewald splitting 1

r = erfc(αr)

r +erf(αr) r withr:=kxij+Lnkand obtain

φ(j) =X

n∈S N

X

i=1 0 qi

kxij+Lnk

=X

n∈S N

X

i=1

0qierfc(αkxij+Lnk) kxij+Lnk +

φ^L(j) :=X

n∈S N

X

i=1

qi

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

√πqj

• short range partcan be obtained via direct evaluation after truncation

• lim

r→0 erf(αr)

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

• write thelong range partas a sum in Fourier space

(29)

Ewald summation

r = erfc(αr)

φ(j) =X

n∈S N

X

i=1 0 qi

kxij+Lnk =X

n∈S N

X

i=1

φ^L(j) :=X

n∈S N

X

i=1

qi

√πqj

• lim

r→0 erf(αr)

(30)

Ewald summation

r = erfc(αr)

φ(j) =X

n∈S N

X

i=1 0 qi

kxij+Lnk =X

n∈S N

X

i=1

φ^L(j) :=

X

n∈S N

X

i=1

qi

√πqj

• lim

r→0 erf(αr)

(31)

Ewald summation

r = erfc(αr)

φ(j) =X

n∈S N

X

i=1 0 qi

kxij+Lnk =X

n∈S N

X

i=1

φ^L(j) :=X

n∈S N

X

i=1

qi

erf(αkxij+Lnk) kxij+Lnk

− 2α

√πqj

• lim

r→0 erf(αr)

(32)

Fast Ewald summation for charged particle systems The 3d-periodic case

3d-periodic Ewald summation

Applying the splitting we obtain

φ(j) =φ^short(j) +φ^F(j) +φ^self(j) +φ^correct(j), where

φ^short(j) = X

n∈Z³ N

X

i=1

0qierfc(αkxij+Lnk) kxij+Lnk

φ^F(j) = 1 πL

X

k∈Z³\{0}

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

N

X

i=1

qie^2πik^T^xⁱ^/L

!

e^−2πik^T^x^j^/L

φ^self(j) =−2α

√πqj

and the correction termφ^correct(j)represents the applied order of summation or rather the surrounding medium.

(33)

Order of summation

The correction term reads as

φ^correct(j) = 4π V(1 + 2)

N

X

i=1

qixi

!^T

xj− 2π V(1 + 2)

N

X

i=1

qikxik²,

wheredenotes the dielectric constant of the surrounding medium.

Corresponding summation orders:

vacuum b.c.,= 1 X

i,j

X

n∈Z³

0:= lim

R→∞

N

X

i,j=1

X

kLnk≤R 0

spherical shells, whole boxes

metallic b.c.,=∞ 1

(ScaFaCoS)

X

i,j

X

n∈Z³

0:= lim

R→∞

N

X

i,j=1

X

kx_ij+Lnk≤R 0

increasing sharp cutoff

Franziska Nestler TU Chemnitz, Faculty of Mathematics

(34)

NFFT based computation of the long range part

Define the Fourier coefficients ψ(k) :=ˆ

( 1 πL

e^−π²^kkk²^/(α²^L2)

kkk² :k6=0,

0 :k=0.

Truncate the infinite sum and compute the structure factors

S(k) :=

N

X

i=1

qie^2πik^T^xⁱ^/L≈S˜(k),

wherek∈ IM, via theadjoint NFFT.

Finally, the long range parts of the potentials are approximated via theNFFT.

φ^L(j)≈ X

k∈I_M

ψ(k) ˜ˆ S(k) e^−2πik^T^x^j^/L.

(35)

Fields and forces

• Potentials:φ(j) =φ^short(j) +φ^F(j) +φ^self(j) +φ^correct(j)

• Energies:U(j) =qjφ(j)

• Overall energy:U =¹₂PN j=1qjφ(j)

• Fields:E(j) =−∇x_jφ(j) =E^short(j) +E^F(j)+E^correct(j)

• Forces:F(j) =qjE(j)

Fourier space parts of the fields:

E^F(j) =−∇xj

X

k∈I_M

ψ(k)Sˆ (k) e^−2πik^T^x^j^/L

(36)

Gradient NFFT

∇f(xj) =∇x_j

X

k∈I_M

fˆke^−2πik^T^x^j.

• ik-differentiation:

∇f(xj) =−2πi X

k∈I_M

fˆkke^−2πik^T^x^j.

We need an FFT in each of the 3 dimensions.

{fˆk}7→ {ˆ^D gk} 7→ {−2πiˆgkk}^3xFFT7→ {g_l} 7→X

l

g_lϕ(x˜ j−^l/σM)

• Analytical differentiation:

{fˆk}7→ {ˆ^D gk}^1xFFT7→ {gl} 7→X

l

gl∇ϕ(x˜ j−^l/σM)

Apply the differentiation operator to the NFFT window function. Only one FFT is needed.

(37)

Gradient NFFT

∇f(xj) =∇x_j

X

k∈I_M

fˆke^−2πik^T^x^j.

• ik-differentiation:

∇f(xj) =−2πi X

k∈I_M

fˆkke^−2πik^T^x^j.

We need an FFT in each of the 3 dimensions.

{fˆk}7→ {ˆ^D gk} 7→ {−2πiˆgkk}^3xFFT7→ {g_l} 7→X

l

g_lϕ(x˜ j−^l/σM)

• Analytical differentiation:

{fˆk}7→ {ˆ^D gk}^1xFFT7→ {gl} 7→X

l

gl∇ϕ(x˜ j−^l/σM)

Apply the differentiation operator to the NFFT window function. Only one FFT is needed.

(38)

Fast Ewald summation for charged particle systems Mixed periodic boundary conditions

Ewald summation with 2d-periodic b.c.

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

We writex_ij= (˜x_ij, xij,3)and obtain φ^L(j) = X

k∈Z²\{0}

N

X

i=1

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

N

X

i=1

qiΘ^p2(0, xij,3)

2d Fourier series k=0part

Θ^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

e2πkr/L 2Lk erfc_πk

αL+αr

+^e^−2πkr/L_2Lk erfc_πk αL−αr

Θ^p2(0, r) =²

√π

αL2e^−α²^r²+^2rπ_L₂erf(αr)

(39)

Ewald summation with 2d-periodic b.c.

k∈Z²\{0}

N

X

i=1

N

X

i=1

qiΘ^p2(0, 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

e2πkr/L 2Lk erfc

πk αL+αr

+^e^−2πkr/L_2Lk erfc πk αL−αr

Θ^p2(0, r) =²

√π

(40)

Ewald summation with 2d-periodic b.c.

k∈Z²\{0}

N

X

i=1

N

X

i=1

qiΘ^p2(0, 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

e2πkr/L 2Lk erfc

πk αL+αr

+^e^−2πkr/L_2Lk erfc πk αL−αr

Θ^p2(0, r) =²

√π

(41)

Apply 1d fast summation

k= 0and smallk6= 0 k6= 0large enough

0

−L L

−^h/2 h/2

R(k,·) Θ^p2(k,·)

polynomial polynomial

1

−^h/2 h/2

Θ^p2(k, r)≈ X∞ n=−∞

Θ^p2(k, r+nh) Θ^p2(k,·) h-periodization

1

Choose an appropriate periodh >2L:

• embedΘ^p2(k,·)into a smooth function R(k,·)(regularization)

• interpolate the derivatives at the boundary up to orderp∈N

• approximate Fourier coefficients byFFT

• approximate the function by itsh-periodic version

• useanalytical Fourier coefficients (Poisson summation formula)

Result: Θ^p2(k, r)≈

M3/2−1

P

l=−^M3/2

ˆbk,le^2πilr/h , r∈[−L, L].

(42)

NFFT based algorithm

φ^L(j) = X

k∈Z²

e^2πik·^x^˜^ij^/L

N

X

i=1

qiΘ^p2(kkk, xij,3)

≈ X

k∈I_M

e^2πik·^x^˜^ij^/L

N

X

i=1

qiΘ^p2(kkk, xij,3)

≈ X

k∈I_M

X

l∈IM3

ˆbkkk,l







N

X

i=1

qie

2πi



 k1 k2 l



·



 xi1/L xi2/L xi3/h











| {z }

3d adj. NFFT

e

−2πi



 k1 k2 l



·



 xj1/L xj2/L xj3/h





| {z }

3d NFFT,j=1,...,N

(43)

1d and 0d-periodic b.c.

• 2d-periodic:

φ^F(j) = X

k∈Z² N

X

i=1

qie^2πik^T^x^˜^ij^/LΘ^p2(kkk,|xij,3|) Approximate the functionsΘ^p2(kkk,·)by trigonometric polynomials.

• 1d-periodic:

φ^F(j) =X

k∈Z N

X

i=1

qie^2πikx^ij,1^/LΘ^p1(|k|,k˜xijk)

Approximate the functionsΘ^p1(|k|,k · k)by bivariate trigonometric polynomials.

• 0d-periodic (open):

φ^F(j) =

N

X

i=1

qi

erf(αkxijk) kxijk

Approximateerf(αk · k)/k · kby a trivariate trigonometric polynomial.

(44)

1d and 0d-periodic b.c.

• 2d-periodic:

φ^F(j) = X

k∈Z² N

X

i=1

qie^2πik^T^x^˜^ij^/LΘ^p2(kkk,|xij,3|)

Approximate the functionsΘ^p2(kkk,·)by trigonometric polynomials.

• 1d-periodic:

φ^F(j) =X

k∈Z N

X

i=1

qie^2πikx^ij,1^/LΘ^p1(|k|,k˜x_ijk)

φ^F(j) =

N

X

i=1

qi

erf(αkxijk) kxijk

(45)

1d and 0d-periodic b.c.

• 2d-periodic:

φ^F(j) = X

k∈Z² N

X

i=1

qie^2πik^T^x^˜^ij^/LΘ^p2(kkk,|xij,3|)

Approximate the functionsΘ^p2(kkk,·)by trigonometric polynomials.

• 1d-periodic:

φ^F(j) =X

k∈Z N

X

i=1

qie^2πikx^ij,1^/LΘ^p1(|k|,k˜x_ijk)

φ^F(j) =

N

X

i=1

qi

erf(αkxijk) kxijk

(46)

P

²

NFFT, computation of the Fourier space part

3d-periodic: φ^F(j) = X

k∈IM

ψ(k)S(k)ˆ e

−2πikT





 xj,1/L xj,2/L xj,3/L







(k,l)∈IM

ˆbk,lS(k, l)e

−2πi(k,l)T





 xj,1/L xj,2/L xj,3/h







(k,l)∈IM

ˆbk,lS(k,l)e

−2πi(k,l)T





 xj,1/L xj,2/h xj,3/h







l∈IM

ˆblS(l)e

−2πilT





 xj,1/h xj,2/h xj,3/h







Steps:1. adjoint NFFT, 2. multiply with the Fourier coefficients, 3. NFFT:

φ^F≈BF D·D·DF^∗B^T·q