Coverage for gpaw/test/lrtddft/test_3.py: 100%

1import re 

2 

3import pytest 

4import numpy as np 

5from ase.build import molecule 

6from ase.units import Hartree 

7 

8from gpaw import GPAW 

9from gpaw.mpi import world 

10from gpaw.gauss import Gauss 

11from gpaw.lrtddft import LrTDDFT, photoabsorption_spectrum 

12from gpaw.lrtddft.kssingle import KSSingles 

13 

14 

15@pytest.mark.lrtddft 

16@pytest.mark.slow 

17@pytest.mark.skipif(world.size > 1, reason='test is serial') 

18def test_lrtddft_3(in_tmp_dir): 

19 from io import StringIO 

20 

21 txt = None 

22 N2 = molecule('N2', vacuum=2.0) 

23 

24 calc = GPAW(mode='fd', 

25 h=0.25, 

26 nbands=-2, 

27 spinpol=True, 

28 xc='LDA', 

29 txt=txt) 

30 

31 N2.calc = calc 

32 _ = N2.get_potential_energy() 

33 calc.write('N2_wfs.gpw', 'all') 

34 

35 # selections 

36 for obj in [KSSingles, LrTDDFT]: 

37 # selection using state numbers 

38 el = obj(restrict={'istart': 3, 'jend': 6}, txt=txt) 

39 el.calculate(N2) 

40 if hasattr(obj, 'diagonalize'): 

41 el.diagonalize() 

42 

43 assert len(el) == 8 

44 # selection using an energy range 

45 el = obj(restrict={'energy_range': 8}, txt=txt) 

46 el.calculate(N2) 

47 if hasattr(obj, 'diagonalize'): 

48 el.diagonalize() 

49 

50 assert len(el) == 4 

51 el = obj(restrict={'energy_range': 11.5}, txt=txt) 

52 el.calculate(N2) 

53 

54 if hasattr(obj, 'diagonalize'): 

55 el.diagonalize() 

56 

57 # Used to be == 18, but we lowered vacuum so test runs fast 

58 # and now it's: 

59 assert len(el) == 16 

60 

61 if hasattr(obj, 'diagonalize'): 

62 el.diagonalize(restrict={'energy_range': 8}) 

63 assert len(el) == 4 

64 

65 lr = LrTDDFT(calc, nspins=2) 

66 lr.write('lrtddft3.dat.gz') 

67 lr.diagonalize() 

68 

69 # This is done to test if writing and reading again yields the same result 

70 lr2 = LrTDDFT.read('lrtddft3.dat.gz') 

71 lr2.diagonalize() 

72 

73 # Post processing 

74 

75 Epeak = 19.5 # The peak we want to investigate (this is alone) 

76 Elist = np.asarray([lrsingle.get_energy() * Hartree 

77 for lrsingle in lr]) 

78 n = np.argmin(np.abs(Elist - Epeak)) # index of the peak 

79 

80 E = lr[n].get_energy() * Hartree 

81 osz = lr[n].get_oscillator_strength() 

82 print('Original object :', E, osz[0]) 

83 

84 # Test the output of analyse 

85 sio = StringIO() 

86 lr.analyse(n, out=sio) 

87 s = sio.getvalue() 

88 

89 match = re.findall( 

90 r'%i: E=([0-9]*\.[0-9]*) eV, f=([0-9]*\.[0-9]*)*' % n, s) 

91 Eanalyse = float(match[0][0]) 

92 oszanalyse = float(match[0][1]) 

93 print('From analyse :', Eanalyse, oszanalyse) 

94 assert E == pytest.approx(Eanalyse, 

95 abs=1e-3) # Written precision in analyse 

96 assert osz[0] == pytest.approx(oszanalyse, abs=1e-3) 

97 

98 E2 = lr2[n].get_energy() * Hartree 

99 osz2 = lr2[n].get_oscillator_strength() 

100 print('Written and read object:', E2, osz2[0]) 

101 

102 # Compare values of original and written/read objects 

103 assert E == pytest.approx(E2, abs=1e-4) 

104 for i in range(len(osz)): 

105 assert osz[i] == pytest.approx(osz2[i], abs=1.7e-4) 

106 

107 width = 0.05 

108 photoabsorption_spectrum(lr, 

109 spectrum_file='lrtddft3-spectrum.dat', 

110 width=width) 

111 # We need to be able to check the heights in the spectrum 

112 weight = Gauss(width).get(0) 

113 

114 spectrum = np.loadtxt('lrtddft3-spectrum.dat', usecols=(0, 1)) 

115 idx = (spectrum[:, 0] >= E - 0.1) & (spectrum[:, 0] <= E + 0.1) 

116 peak = np.argmax(spectrum[idx, 1]) + np.nonzero(idx)[0][0] 

117 Espec = spectrum[peak, 0] 

118 oszspec = spectrum[peak, 1] / weight 

119 

120 print('Values from spectrum :', Espec, oszspec) 

121 # Compare calculated values with values written to file 

122 assert E == pytest.approx(Espec, 

123 abs=1e-2) # The spectrum has a low sampling 

124 assert osz[0] == pytest.approx(oszspec, abs=1e-2)