Coverage for gpaw/test/maths/test_fsbt.py: 100%

1import numpy as np 

2 

3from gpaw.atom.radialgd import fsbt, EquidistantRadialGridDescriptor 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[1], len(r)) 

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[1], len(g)) 

30 r, f = ggd.fft(f * g) 

31 assert abs(np.exp(-alpha * r**2) - f / 8 / np.pi**3).max() < 2e-3