Coverage for gpaw/test/utilities/test_ldos.py: 100%

1import numpy as np 

2from ase import Atom, Atoms 

3 

4import gpaw.mpi as mpi 

5from gpaw import GPAW, FermiDirac, PoissonSolver 

6import pytest 

7from gpaw.utilities.dos import RawLDOS, raw_orbital_LDOS, raw_wignerseitz_LDOS 

8 

9 

10def test_utilities_ldos(in_tmp_dir, gpaw_new): 

11 comms = [mpi.world.new_communicator(np.array([r])) 

12 for r in range(mpi.size)] 

13 comm = comms[mpi.rank] 

14 

15 Hnospin = Atoms([Atom('H')], cell=[5, 5, 5], pbc=False) 

16 Hspin = Atoms([Atom('H', magmom=1)], cell=[5, 5, 5], pbc=False) 

17 LiH = Atoms([Atom('Li', [.0, .0, .41]), 

18 Atom('H', [.0, .0, -1.23])]) 

19 Hnospin.center() 

20 Hspin.center() 

21 LiH.center(vacuum=3.0) 

22 

23 # This is needed for the Wigner-Seitz test to give 

24 # architecture-independent results: 

25 LiH.translate(0.003234) 

26 

27 calc = GPAW(mode='fd', gpts=(24, 24, 24), communicator=comm) 

28 Hnospin.calc = calc 

29 e_Hnospin = Hnospin.get_potential_energy() 

30 energies, sweight = raw_orbital_LDOS(calc, a=0, spin=0, angular='s') 

31 energies, pdfweight = raw_orbital_LDOS(calc, a=0, spin=0, angular='pdf') 

32 

33 calc = GPAW(mode='fd', 

34 gpts=(24, 24, 24), 

35 occupations=FermiDirac(width=0, fixmagmom=True), 

36 poissonsolver=PoissonSolver('fd'), 

37 hund=True, 

38 communicator=comm) 

39 Hspin.calc = calc 

40 e_Hspin = Hspin.get_potential_energy() 

41 energies, sweight_spin = raw_orbital_LDOS(calc, a=0, spin=0, angular='s') 

42 

43 calc = GPAW(mode='fd', gpts=(32, 32, 40), nbands=2, 

44 poissonsolver=PoissonSolver('fd'), 

45 communicator=comm) 

46 LiH.calc = calc 

47 e_LiH = LiH.get_potential_energy() 

48 energies, Li_orbitalweight = raw_orbital_LDOS(calc, a=0, spin=0, 

49 angular=None) 

50 energies, H_orbitalweight = raw_orbital_LDOS(calc, a=1, spin=0, 

51 angular=None) 

52 energies, Li_wzweight = raw_wignerseitz_LDOS(calc, a=0, spin=0) 

53 energies, H_wzweight = raw_wignerseitz_LDOS(calc, a=1, spin=0) 

54 

55 if not gpaw_new: 

56 n_a = calc.get_wigner_seitz_densities(spin=0) 

57 print(n_a) 

58 assert n_a.sum() == pytest.approx(0.0, abs=1e-5) 

59 assert n_a[1] == pytest.approx(0.737, abs=0.001) 

60 

61 print(sweight, pdfweight) 

62 print(sweight_spin) 

63 print(Li_wzweight) 

64 print(H_wzweight) 

65 

66 assert sweight[0] == pytest.approx(1.0, abs=0.06) 

67 assert pdfweight[0] == pytest.approx(0.0, abs=0.0001) 

68 assert sweight_spin[0] == pytest.approx(1.14, abs=0.06) 

69 assert ((Li_wzweight - [.13, 0.93]).round(2) == 0).all() 

70 assert ((H_wzweight - [0.87, 0.07]).round(2) == 0).all() 

71 assert ((Li_wzweight + H_wzweight).round(5) == 1).all() 

72 

73 print(Li_orbitalweight) 

74 print(H_orbitalweight) 

75 # HOMO s py pz px *s 

76 Li_orbitalweight[0] -= [.5, .0, .6, .0, .0] 

77 H_orbitalweight[0] -= [.7, .0, .0, .0, .0] 

78 

79 # LUMO s py pz px *s 

80 Li_orbitalweight[1] -= [1.0, .0, 0.9, .0, .0] 

81 H_orbitalweight[1] -= [0.1, .0, 0.0, .0, .0] 

82 

83 assert not Li_orbitalweight.round(1).any() 

84 assert not H_orbitalweight.round(1).any() 

85 

86 ldos = RawLDOS(calc) 

87 fname = 'ldbe' 

88 ldos.by_element_to_file(fname + '.dat', shift=False) 

89 ldos.by_element_to_file(fname + '_2.0.dat', 2.0, shift=False) 

90 ldos.by_element_to_file(fname + '_indx0.dat', indices=[0]) 

91 # the hydrogen entries are missing for index 0 only 

92 if mpi.world.rank == 0: 

93 assert (np.loadtxt(fname + '_indx0.dat').shape[1] + 3 == 

94 np.loadtxt(fname + '.dat').shape[1]) 

95 

96 energy_tolerance = 0.001 

97 assert e_Hnospin == pytest.approx(0.153991, abs=energy_tolerance) 

98 assert e_Hspin == pytest.approx(-0.782309, abs=energy_tolerance) 

99 assert e_LiH == pytest.approx(-3.74582, abs=energy_tolerance)