Coverage for gpaw/test/lcao/test_dipole_transition2.py: 100%

1import numpy as np 

2import pytest 

3 

4from ase.parallel import world, parprint 

5from ase.units import Bohr 

6from gpaw import GPAW 

7from gpaw.lcao.dipoletransition import get_dipole_transitions 

8from gpaw.lrtddft.kssingle import KSSingles 

9 

10 

11@pytest.mark.old_gpaw_only 

12def test_dipole_transition(gpw_files, tmp_path_factory): 

13 """Check dipole matrix-elements for Li.""" 

14 calc = GPAW(gpw_files['bcc_li_lcao'], parallel=dict(sl_auto=True)) 

15 # Initialize calculator if necessary 

16 if not hasattr(calc.wfs, 'C_nM'): 

17 calc.wfs.set_positions 

18 calc.initialize_positions(calc.atoms) 

19 from gpaw.kohnsham_layouts import BlacsOrbitalLayouts 

20 isblacs = isinstance(calc.wfs.ksl, BlacsOrbitalLayouts) # XXX 

21 print(calc.wfs.ksl.using_blacs, isblacs) 

22 dip_skvnm = get_dipole_transitions(calc.wfs) 

23 parprint("Dipole moments calculated") 

24 assert dip_skvnm.shape == (1, 4, 3, 4, 4) 

25 dip_kvnm = dip_skvnm[0] * Bohr 

26 

27 print(world.rank, dip_kvnm[:, 0, 0, 3]) 

28 

29 # check symmetry: abs(d[i,j]) == abs(d[j,i]) 

30 for k in range(4): 

31 for v in range(3): 

32 dip_kvnm[k, v].T == pytest.approx(dip_kvnm[k, v]) 

33 

34 # Check numerical value of a few elements - signs might change! 

35 assert 0.0824 == pytest.approx(abs(dip_kvnm[0, 0, 0, 1]), abs=1e-4) 

36 assert abs(-0.0781 + 0.0282j) == pytest.approx(abs(dip_kvnm[1, 2, 0, 3]), 

37 abs=1e-4) 

38 

39 calc = GPAW(gpw_files['bcc_li_fd']) 

40 # compare to lrtddft implementation 

41 kss = KSSingles() 

42 atoms = calc.atoms 

43 atoms.calc = calc 

44 kss.calculate(calc.atoms, 0) 

45 lrrefv = [] 

46 for ex in kss: 

47 print(ex) 

48 lrrefv.append(-1. * ex.muv * Bohr) 

49 lrrefv = np.array(lrrefv) 

50 

51 # some printout for manual inspection, if wanted 

52 parprint(" lrtddft(fd)(v) raman(v)") 

53 f = "{} {:+.4f} {:+.4f}" 

54 # At gamma at three excited states are degenerate, so order is a problem 

55 parprint(f.format('k=0, 0->1 (x)', lrrefv[0, 0], dip_kvnm[0, 0, 0, 1])) 

56 parprint(f.format('k=0, 0->1 (y)', lrrefv[0, 1], dip_kvnm[0, 1, 0, 1])) 

57 parprint(f.format('k=0, 0->1 (z)', lrrefv[0, 2], dip_kvnm[0, 2, 0, 1])) 

58 parprint(f.format('k=0, 0->2 (x)', lrrefv[1, 0], dip_kvnm[0, 0, 0, 2])) 

59 parprint(f.format('k=0, 0->2 (y)', lrrefv[1, 1], dip_kvnm[0, 1, 0, 2])) 

60 parprint(f.format('k=0, 0->2 (z)', lrrefv[1, 2], dip_kvnm[0, 2, 0, 2])) 

61 parprint(f.format('k=0, 0->3 (x)', lrrefv[2, 0], dip_kvnm[0, 0, 0, 3])) 

62 parprint(f.format('k=0, 0->3 (y)', lrrefv[2, 1], dip_kvnm[0, 1, 0, 3])) 

63 parprint(f.format('k=0, 0->3 (z)', lrrefv[2, 2], dip_kvnm[0, 2, 0, 3])) 

64 parprint("") 

65 # At 2nd kpoint 2nd and 3rd excited state almost degenerate in LCAO, 

66 # order might be different from FD mode 

67 parprint(f.format('k=1, 0->1 (x)', lrrefv[3, 0], dip_kvnm[1, 0, 0, 1])) 

68 parprint(f.format('k=1, 0->1 (y)', lrrefv[3, 1], dip_kvnm[1, 1, 0, 1])) 

69 parprint(f.format('k=1, 0->1 (z)', lrrefv[3, 2], dip_kvnm[1, 2, 0, 1])) 

70 

71 parprint(f.format('k=1, 0->2 (x)', lrrefv[4, 0], dip_kvnm[1, 0, 0, 2])) 

72 parprint(f.format('k=1, 0->2 (y)', lrrefv[4, 1], dip_kvnm[1, 1, 0, 2])) 

73 parprint(f.format('k=1, 0->2 (z)', lrrefv[4, 2], dip_kvnm[1, 2, 0, 2])) 

74 parprint(f.format('k=1, 0->3 (x)', lrrefv[5, 0], dip_kvnm[1, 0, 0, 3])) 

75 parprint(f.format('k=1, 0->3 (y)', lrrefv[5, 1], dip_kvnm[1, 1, 0, 3])) 

76 parprint(f.format('k=1, 0->3 (z)', lrrefv[5, 2], dip_kvnm[1, 2, 0, 3])) 

77 parprint("") 

78 # At the third k-point first and second excited state are degenerate 

79 parprint(f.format('k=2, 0->3 (x)', lrrefv[8, 0], dip_kvnm[2, 0, 0, 3])) 

80 parprint(f.format('k=2, 0->3 (y)', lrrefv[8, 1], dip_kvnm[2, 1, 0, 3])) 

81 parprint(f.format('k=2, 0->3 (z)', lrrefv[8, 2], dip_kvnm[2, 2, 0, 3])) 

82 # 4th k-point not included KSSingles, probably all bands above Fermi level