Coverage for gpaw/test/test_multipoletest.py: 100%

1from math import sqrt, pi 

2 

3import numpy as np 

4import pytest 

5 

6from gpaw.setup import create_setup 

7from gpaw.grid_descriptor import GridDescriptor 

8from gpaw.lfc import LFC 

9from gpaw.xc import XC 

10 

11n = 60 # 40 /8 * 10 

12a = 10.0 

13 

14 

15@pytest.mark.ci 

16def test_multipole(): 

17 gd = GridDescriptor((n, n, n), (a, a, a)) 

18 c_LL = np.identity(9, float) 

19 a_Lg = gd.zeros(9) 

20 xc = XC('LDA') 

21 for soft in [False]: 

22 s = create_setup('Cu', xc, lmax=2) 

23 ghat_l = s.ghat_l 

24 ghat_Lg = LFC(gd, [ghat_l]) 

25 ghat_Lg.set_positions([(0.54321, 0.5432, 0.543)]) 

26 a_Lg[:] = 0.0 

27 ghat_Lg.add(a_Lg, {0: c_LL} if ghat_Lg.my_atom_indices else {}) 

28 for l in range(3): 

29 for m in range(2 * l + 1): 

30 L = l**2 + m 

31 a_g = a_Lg[L] 

32 Q0 = gd.integrate(a_g) / sqrt(4 * pi) 

33 Q1_m = -gd.calculate_dipole_moment(a_g) / sqrt(4 * pi / 3) 

34 print(Q0) 

35 if l == 0: 

36 Q0 -= 1.0 

37 Q1_m[:] = 0.0 

38 elif l == 1: 

39 Q1_m[(m + 1) % 3] -= 1.0 

40 print(soft, l, m, Q0, Q1_m) 

41 assert abs(Q0) < 2e-6 

42 assert (abs(Q1_m) < 3e-5).all() 

43 b_Lg = np.reshape(a_Lg, (9, -1)) 

44 S_LL = np.inner(b_Lg, b_Lg) 

45 gd.comm.sum(S_LL) 

46 S_LL.ravel()[::10] = 0.0 

47 print(max(abs(S_LL).ravel())) 

48 assert max(abs(S_LL).ravel()) < 3e-4