We are solving 1D Poisson equation:
$$
{d^2\over dx^2}V(x) = -4\pi n(x)
$$
Using Fast Fourier Transform (FFT). And then calculating the Hartree energy using:
$$
E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2}
$$


```
%pylab inline
from pylab import plot, semilogy
from numpy import loadtxt, linspace
from numpy.fft import fft, fftfreq
from sympy import exp, sin, pi, var, Integral
```

    
    Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
    For more information, type 'help(pylab)'.


Let's setup the right hand side:


```
var("x")
L = 2 # domain [0, L]
rho = exp(sin(2*pi*x/L))
#rho = sin(2*pi*x/L)
integ = Integral(rho, (x, 0, L)).n()
rho -= integ / L
print "rho(x) =", rho
print "Integral of rho =", Integral(rho, (x, 0, L)).n()
```

    rho(x) = exp(sin(pi*x)) - 1.26606587775201
    Integral of rho = -1.41156085310944e-16



```
D = []
for m in range(1, 15):
    #2m+1 terms -Gmax, ..., 0, ..., Gmax
    n = 2*m+1
    xj = linspace(0, L, n+1)[:-1]
    rhoj = array([float(rho.subs(x, val).n(20)) for val in xj])
    nr = fft(rhoj) / n
    # re-order -> -Gmax, ... Gmax
    f = fftfreq(n)
    idx = f.argsort()
    f = f[idx]
    nr = nr[idx]
    nr[m] = 0  # set G=0 term = 0
    G = array([float((2*pi*j/L)) for j in range(-m, m+1)])
    G[m] = 1
    tmp = array([complex(2*pi*abs(nr[j])**2/G[j]**2) for j in range(n)])
    E = sum(tmp).real * L
    print n, E
    D.append((n, E))
D = array(D)
```

    3 0.857623650249
    5 0.825420490259
    7 0.825230159883
    9 0.825229196177
    11 0.825229191862
    13 0.825229191846
    15 0.825229191846
    17 0.825229191846
    19 0.825229191846
    21 0.825229191846
    23 0.825229191846
    25 0.825229191846
    27 0.825229191846
    29 0.825229191846



```
N = D[:, 0]
E = D[:, 1]
E_exact = E[-1]
figure(figsize=(8, 6), dpi=80)
semilogy(N, abs(E-E_exact), "bo-")
grid()
title("Convergence of an FFT Poisson solver (semilogy)")
xlabel("N (number of PW in each direction)")
ylabel("E  -  E_exact  [a.u.]")

figure(figsize=(8, 6), dpi=80)
loglog(N, abs(E-E_exact), "bo-", 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.]")
```


```

```

ae1d688537a2f522294d2cc554865f9a45ae7da1
ipynb
Jupyter Notebook
tutorial_exercises/FFT Introduction.ipynb
certik/scipy-2013-tutorial
26a1cab3a16402afdc20088cedf47acd9bc58483
BSD-3-Clause
2015-02-28T08:53:05.000Z
2021-12-05T05:37:59.000Z
sympy/FFT Introduction.ipynb
certik/scipy-in-13
418c139ab6e1b0c9acd53e7e1a02b8b930005096
BSD-3-Clause
2021-04-17T15:05:46.000Z
2021-04-17T15:05:46.000Z
sympy/FFT Introduction.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