Coverage for gpaw/test/sphere/test_rshe.py: 100%

1"""Test functionality to expand a function on a spherical grid in real 

2spherical harmonics.""" 

3 

4import pytest 

5import itertools 

6 

7import numpy as np 

8 

9from gpaw.atom.radialgd import EquidistantRadialGridDescriptor 

10from gpaw.sphere.lebedev import Y_nL 

11from gpaw.sphere.rshe import (RealSphericalHarmonicsExpansion, 

12 calculate_reduced_rshe) 

13 

14 

15def generate_combinations(): 

16 """Generate combinations of one and two spherical harmonic indices. 

17 """ 

18 L_L = np.arange(Y_nL.shape[1]) 

19 oneL_i = list(itertools.combinations(L_L, 1)) 

20 twoL_i = list(itertools.combinations(L_L, 2)) 

21 

22 Lcomb_i = oneL_i + twoL_i 

23 

24 return Lcomb_i 

25 

26 

27@pytest.mark.parametrize('Lcomb', generate_combinations()) 

28def test_rshe(Lcomb): 

29 """Test the ability to regenerate a given function based on a 

30 real spherical harmonics expansion.""" 

31 

32 # Build the angular dependence of the function on the Lebedev quadrature 

33 f_n = np.zeros(Y_nL.shape[0]) 

34 for L in Lcomb: 

35 f_n += Y_nL[:, L] 

36 # Build the radial dependence as a Lorentzian 

37 rgd = EquidistantRadialGridDescriptor(h=0.2, # grid spacing 

38 N=11) 

39 f_g = 0.25 / (np.pi * (0.25**2 + rgd.r_g**2.)) 

40 f_ng = f_n[:, np.newaxis] * f_g[np.newaxis] 

41 

42 # Test the real spherical harmonics expansion without basis reduction 

43 rshe = RealSphericalHarmonicsExpansion.from_spherical_grid(rgd, f_ng, Y_nL) 

44 assert rshe.evaluate_on_quadrature() == pytest.approx(f_ng) 

45 

46 # Test the ability to reduce the expansion via minimum weights 

47 rshe, _ = calculate_reduced_rshe(rgd, f_ng, Y_nL, wmin=1e-8) 

48 assert len(rshe.L_M) == len(Lcomb) 

49 assert all([L in rshe.L_M for L in Lcomb]) 

50 assert all([int(np.sqrt(L)) in rshe.l_M for L in Lcomb]) 

51 assert rshe.evaluate_on_quadrature() == pytest.approx(f_ng) 

52 

53 # Test the ability to reduce the expansion by an lmax 

54 Lmax = max(Lcomb) 

55 lmax = int(np.ceil(np.sqrt((Lmax + 1)) - 1)) 

56 if lmax < 4: 

57 rshe, _ = calculate_reduced_rshe(rgd, f_ng, Y_nL, lmax=lmax) 

58 assert len(rshe.L_M) == (lmax + 1)**2 

59 assert rshe.evaluate_on_quadrature() == pytest.approx(f_ng)