Coverage for gpaw/test/radial/test_ylexpand.py: 94%

1from math import pi 

2 

3from ase import Atoms 

4from ase.parallel import parprint 

5 

6from gpaw import GPAW 

7import pytest 

8from gpaw.analyse.expandyl import AngularIntegral, ExpandYl 

9 

10 

11def test_radial_ylexpand(in_tmp_dir): 

12 fname = 'H2.gpw' 

13 donot = '' 

14 donot = 'donot' 

15 try: 

16 calc = GPAW(fname + donot, txt=None) 

17 H2 = calc.get_atoms() 

18 calc.converge_wave_functions() 

19 except FileNotFoundError: 

20 R = 0.7 # approx. experimental bond length 

21 a = 2. 

22 c = 3. 

23 H2 = Atoms('HH', 

24 [(a / 2, a / 2, (c - R) / 2), 

25 (a / 2, a / 2, (c + R) / 2)], 

26 cell=(a, a, c), 

27 pbc=True) 

28 calc = GPAW(mode='fd', gpts=(12, 12, 16), nbands=2, kpts=(1, 1, 2), 

29 convergence={'eigenstates': 1.e-6}, 

30 txt=None, 

31 ) 

32 H2.calc = calc 

33 H2.get_potential_energy() 

34 if not donot: 

35 calc.write(fname) 

36 

37 # Check that a / h = 10 is rounded up to 12 as always: 

38 assert (calc.wfs.gd.N_c == (12, 12, 16)).all() 

39 

40 # AngularIntegral 

41 

42 gd = calc.density.gd 

43 ai = AngularIntegral(H2.positions.mean(0), calc.wfs.gd, Rmax=1.5) 

44 unity_g = gd.zeros() + 1. 

45 average_R = ai.average(unity_g) 

46 integral_R = ai.integrate(unity_g) 

47 for V, average, integral, R, Rm in zip(ai.V_R, average_R, integral_R, 

48 ai.radii(), ai.radii('mean')): 

49 if V > 0: 

50 assert average == pytest.approx(1, abs=1.e-9) 

51 assert integral / (4 * pi * Rm**2) == pytest.approx(1, abs=0.61) 

52 assert Rm / R == pytest.approx(1, abs=0.61) 

53 

54 # ExpandYl 

55 

56 yl = ExpandYl(H2.positions.mean(0), calc.wfs.gd, Rmax=1.5) 

57 

58 def max_index(l): 

59 mi = 0 

60 limax = l[0] 

61 for i, li in enumerate(l): 

62 if limax < li: 

63 limax = li 

64 mi = i 

65 return mi 

66 

67 # check numbers 

68 for n in [0, 1]: 

69 # gl, w = yl.expand(calc.get_pseudo_wave_function(band=n)) 

70 gl, w = yl.expand(calc.wfs.kpt_u[0].psit_nG[n]) 

71 parprint('max_index(gl), n=', max_index(gl), n) 

72 assert max_index(gl) == n 

73 

74 # io 

75 fname = 'expandyl.dat' 

76 yl.to_file(calc, fname)