```
%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import sqrt, sin, pi, var, integrate, Symbol, S, Integral
var("x y z")
n = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
Vh = -((x-1)**4 + (y-1)**4 + (z-1)**4) + 2*(x**2+y**2+z**2)
def laplace(f):
    return f.diff(x, 2) + f.diff(y, 2) + f.diff(z, 2)
print "n ="
n
```


```
print "Vh ="
Vh
```

Let's check charge neutrality:


```
integrate(n, (x, 0, 2), (y, 0, 2), (z, 0, 2))
```

$V_H$ is produced by the charge density $n$:


```
laplace(Vh) / (-4*pi)
```


```
_.simplify()
```

The Hartree energy $E_H$ is: 


```
integrate(Vh*n, (x, 0, 2), (y, 0, 2), (z, 0, 2))/2
```

We are solving a Poisson equation:
$$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$
with
$$n(x, y, z) = {3\over\pi} ((x-1)^2 + (y-1)^2 + (z-1)^2 - 1)$$
on the domain $[0, 2] \times [0, 2] \times [0, 2]$ with periodic boundary condition. This charge density is net neutral. The solution is given by:
$$
V_H({\bf x}) = \int {n({\bf y})\over |{\bf x}-{\bf y}|} d^3 y
$$
where ${\bf x}=(x, y, z)$. In reciprocal space
$$
V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}
$$
where ${\bf G}$ are the reciprocal space vectors. I am interested in the Hartree energy:
$$
E_H = {1\over 2} \int {n({\bf x}) n({\bf y})\over |{\bf x}-{\bf y}|} d^3 x d^3 y
= {1\over 2} \int V_H({\bf x}) n({\bf x}) d^3 x
$$
In reciprocal space this becomes (after discretization):
$$
E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2}
$$
The ${\bf G}=0$ term is omitted, which effectively makes the charge density net neutral (and since it is already neutral, then everything is consistent).

For the test problem above, this can be evaluated analytically and one gets:
$$
E_H = {128\over 35\pi} = 1.16410...
$$
Now we implement the FFT solver in NumPy below and plot the convergence graphs with respect to this exact energy.


This problem has also been discussed at:

http://scicomp.stackexchange.com/questions/7097/convergence-rate-of-fft-poisson-solver/7169

[Computational Science StackExchange](http://scicomp.stackexchange.com/) is a great site to ask these questions!


```
from numpy import empty, pi, meshgrid, linspace, sum
from numpy.fft import fftn, fftfreq
E_exact = 128/(35*pi)
print "Hartree Energy (exact):      %.15f" % E_exact
f = open("conv.txt", "w")
for N in range(3, 200, 10):
    print "N =", N
    L = 2.
    x1d = linspace(0, L, N+1)[:-1]
    x, y, z = meshgrid(x1d, x1d, x1d)

    nr = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
    ng = fftn(nr) / N**3

    G1d = N * fftfreq(N) * 2*pi/L
    kx, ky, kz = meshgrid(G1d, G1d, G1d)
    G2 = kx**2+ky**2+kz**2
    G2[0, 0, 0] = 1  # omit the G=0 term

    tmp = 2*pi*abs(ng)**2 / G2
    tmp[0, 0, 0] = 0  # omit the G=0 term
    E = sum(tmp) * L**3
    print "Hartree Energy (calculated): %.15f" % E
    f.write("%d %.15f\n" % (N, E))
f.close()
```

    Hartree Energy (exact):      1.164104726615006
    N = 3
    Hartree Energy (calculated): 2.446338611821947
    N = 13
    Hartree Energy (calculated): 1.213786470434117
    N = 23
    Hartree Energy (calculated): 1.179663010471846
    N = 33
    Hartree Energy (calculated): 1.171624324575022
    N = 43
    Hartree Energy (calculated): 1.168524684513086
    N = 53
    Hartree Energy (calculated): 1.167011254992898
    N = 63
    Hartree Energy (calculated): 1.166160629019379
    N = 73
    Hartree Energy (calculated): 1.165635416812430
    N = 83
    Hartree Energy (calculated): 1.165288523120677
    N = 93
    Hartree Energy (calculated): 1.165047479421731
    N = 103
    Hartree Energy (calculated): 1.164873217880501
    N = 113
    Hartree Energy (calculated): 1.164743164819577
    N = 123
    Hartree Energy (calculated): 1.164643537882899
    N = 133
    Hartree Energy (calculated): 1.164565535654283
    N = 143
    Hartree Energy (calculated): 1.164503323719686
    N = 153
    Hartree Energy (calculated): 1.164452910534660
    N = 163
    Hartree Energy (calculated): 1.164411490567212
    N = 173
    Hartree Energy (calculated): 1.164377045217771
    N = 183
    Hartree Energy (calculated): 1.164348092059736
    N = 193
    Hartree Energy (calculated): 1.164323522587376



```
from pylab import plot, semilogy
from numpy import loadtxt, pi
D = loadtxt("conv.txt")
N = D[:, 0]
E = D[:, 1]
E_exact = 128/(35*pi)
```


```
figure(figsize=(8, 6), dpi=80)
semilogy(N, E-E_exact, "k-")
grid()
title("Convergence of an FFT Poisson solver (semilogy)")
xlabel("N (number of PW in each direction)")
ylabel("E  -  E_exact  [a.u.]")
savefig("fft_convergence_semilogy.png")
```


```
figure(figsize=(8, 6), dpi=80)
calculated = E-E_exact
loglog(N, E-E_exact, "ko", label="calculated")
predicted = 1/N**2
predicted = predicted * (calculated[-1] /  predicted[-1])
loglog(N, predicted, "g-", label="$1/N^2$")
grid()
title("Convergence of an FFT Poisson solver (loglog)")
xlabel("N (number of PW in each direction)")
ylabel("E  -  E_exact  [a.u.]")
legend()
savefig("fft_convergence_loglog.png")
```

# Exercise

Rerun the above FFT calculation with the following potential:


```
nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4
```

Plot convergence study. Why is the convergence rate different?

## Solution


```
from numpy import empty, pi, meshgrid, linspace, sum, exp, sin
from numpy.fft import fftn, fftfreq
f = open("conv2.txt", "w")
for N in range(3, 30, 2):
    print "N =", N
    L = 2.
    x1d = linspace(0, L, N+1)[:-1]
    x, y, z = meshgrid(x1d, x1d, x1d)

    nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4
    ng = fftn(nr) / N**3

    G1d = N * fftfreq(N) * 2*pi/L
    kx, ky, kz = meshgrid(G1d, G1d, G1d)
    G2 = kx**2+ky**2+kz**2
    G2[0, 0, 0] = 1  # omit the G=0 term

    tmp = 2*pi*abs(ng)**2 / G2
    tmp[0, 0, 0] = 0  # omit the G=0 term
    E = sum(tmp) * L**3
    print "Hartree Energy (calculated): %.15f" % E
    f.write("%d %.15f\n" % (N, E))
f.close()
```

    N = 3
    Hartree Energy (calculated): 1.580426011908313
    N = 5
    Hartree Energy (calculated): 1.414548280420217
    N = 7
    Hartree Energy (calculated): 1.413994578897476
    N = 9
    Hartree Energy (calculated): 1.413991759887669
    N = 11
    Hartree Energy (calculated): 1.413991748461275
    N = 13
    Hartree Energy (calculated): 1.413991748422897
    N = 15
    Hartree Energy (calculated): 1.413991748422800
    N = 17
    Hartree Energy (calculated): 1.413991748422796
    N = 19
    Hartree Energy (calculated): 1.413991748422795
    N = 21
    Hartree Energy (calculated): 1.413991748422802
    N = 23
    Hartree Energy (calculated): 1.413991748422796
    N = 25
    Hartree Energy (calculated): 1.413991748422800
    N = 27
    Hartree Energy (calculated): 1.413991748422803
    N = 29
    Hartree Energy (calculated): 1.413991748422796



```
from pylab import plot, semilogy
from numpy import loadtxt, pi
D = loadtxt("conv2.txt")
N = D[:, 0]
E = D[:, 1]
E_exact = E[-1]

figure(figsize=(8, 6), dpi=80)
calculated = E-E_exact
semilogy(N, E-E_exact, "ko-", label="calculated")
grid()
title("Convergence of an FFT Poisson solver (loglog)")
xlabel("N (number of PW in each direction)")
ylabel("E  -  E_exact  [a.u.]")
legend()
savefig("fft_convergence_loglog.png")
```


```

```

da2f8f853147f3ab1b72632e903396053e2bd6d0
ipynb
Jupyter Notebook
tutorial_exercises/FFT convergence.ipynb
certik/scipy-2013-tutorial
26a1cab3a16402afdc20088cedf47acd9bc58483
BSD-3-Clause
2015-02-28T08:53:05.000Z
2021-12-05T05:37:59.000Z
sympy/FFT convergence.ipynb
certik/scipy-in-13
418c139ab6e1b0c9acd53e7e1a02b8b930005096
BSD-3-Clause
2021-04-17T15:05:46.000Z
2021-04-17T15:05:46.000Z
sympy/FFT convergence.ipynb
certik/scipy-in-13
418c139ab6e1b0c9acd53e7e1a02b8b930005096
BSD-3-Clause
2015-03-11T00:25:21.000Z
2021-08-25T14:52:40.000Z
Qwen/Qwen-72B
1. YES
2. YES