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} $$



In [5]:

    
%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:



In [6]:

    
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



In [7]:

    
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



In [8]:

    
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.]")









    Out[8]:





<matplotlib.text.Text at 0x515eb90>



In [ ]: