1import numpy as np
2
3from gpaw.atom.radialgd import fsbt, RadialGridDescriptor as RGD
4
5
6def test_maths_fsbt():
7    N = 1024
8    L = 50.0
9    h = L / N
10    alpha = 5.0
11    r = np.linspace(0, L, N, endpoint=False)
12    G = np.linspace(0, np.pi / h, N // 2 + 1)
13    n = np.exp(-alpha * r**2)
14
15    for l in range(7):
16        f = fsbt(l, n * r**l, r, G)
17        f0 = (np.pi**0.5 / alpha**(l + 1.5) / 2**l * G**l / 4 *
18              np.exp(-G**2 / (4 * alpha)))
19        tol = 3e-6 * 10**(-7 + l)
20        print(l, abs(f - f0).max(), 'tol=', tol)
21        assert abs(f - f0).max() < tol
22
23    rgd = RGD(r, r * 0 + r[1])
24    g, f = rgd.fft(n * r)
25    f0 = 4 * np.pi**1.5 / alpha**1.5 / 4 * np.exp(-g**2 / 4 / alpha)
26    assert abs(f - f0).max() < 1e-6
27
28    # This is how to do the inverse FFT:
29    ggd = RGD(g, g * 0 + g[1])
30    r, f = ggd.fft(f * g)
31    assert abs(np.exp(-alpha * r**2) - f / 8 / np.pi**3).max() < 2e-3
32