Coverage for gpaw/test/response/test_diamond_absorption.py: 100%

1import pytest 

2import numpy as np 

3from ase.units import Bohr 

4from ase.build import bulk 

5from gpaw import GPAW, FermiDirac 

6from gpaw.response.df import DielectricFunction, read_response_function 

7from gpaw.test import findpeak 

8 

9 

10@pytest.mark.dielectricfunction 

11@pytest.mark.response 

12@pytest.mark.parametrize('eshift', [None, 4]) 

13@pytest.mark.libxc 

14def test_response_diamond_absorption(in_tmp_dir, eshift): 

15 a = 6.75 * Bohr 

16 atoms = bulk('C', 'diamond', a=a) 

17 

18 calc = GPAW(mode='pw', 

19 kpts=(3, 3, 3), 

20 eigensolver='rmm-diis', 

21 occupations=FermiDirac(0.001)) 

22 

23 atoms.calc = calc 

24 atoms.get_potential_energy() 

25 calc.write('C.gpw', 'all') 

26 

27 if eshift is None: 

28 eM1_ = 9.727 

29 eM2_ = 9.548 

30 w0_ = 10.7782 

31 I0_ = 5.47 

32 w_ = 10.7532 

33 I_ = 5.98 

34 else: 

35 eM1_ = 6.993 

36 eM2_ = 6.904 

37 w0_ = 14.784 

38 I0_ = 5.47 

39 w_ = 14.757 

40 I_ = 5.998 

41 

42 # Test the old interface to the dielectric constant 

43 df = DielectricFunction('C.gpw', frequencies=(0.,), eta=0.001, ecut=50, 

44 hilbert=False, eshift=eshift) 

45 eM1, eM2 = df.get_macroscopic_dielectric_constant() 

46 assert eM1 == pytest.approx(eM1_, abs=0.01) 

47 assert eM2 == pytest.approx(eM2_, abs=0.01) 

48 

49 # ----- RPA dielectric function ----- # 

50 dfcalc = DielectricFunction( 

51 'C.gpw', eta=0.25, ecut=50, 

52 frequencies=np.linspace(0, 24., 241), hilbert=False, eshift=eshift) 

53 eps = dfcalc.get_literal_dielectric_function() 

54 

55 # Test the new interface to the dielectric constant 

56 eM1, eM2 = eps.dielectric_constant() 

57 assert eM1 == pytest.approx(eM1_, abs=0.01) 

58 assert eM2 == pytest.approx(eM2_, abs=0.01) 

59 

60 # Test the macroscopic dielectric function 

61 omega_w, eps0M_w, epsM_w = eps.macroscopic_dielectric_function().arrays 

62 w0, I0 = findpeak(omega_w, eps0M_w.imag) 

63 assert w0 == pytest.approx(w0_, abs=0.01) 

64 assert I0 / (4 * np.pi) == pytest.approx(I0_, abs=0.1) 

65 w, I = findpeak(omega_w, epsM_w.imag) 

66 assert w == pytest.approx(w_, abs=0.01) 

67 assert I / (4 * np.pi) == pytest.approx(I_, abs=0.1) 

68 

69 # Test polarizability 

70 omega_w, a0rpa_w, arpa_w = eps.polarizability().arrays 

71 w0, I0 = findpeak(omega_w, a0rpa_w.imag) 

72 assert w0 == pytest.approx(w0_, abs=0.01) 

73 assert I0 == pytest.approx(I0_, abs=0.01) 

74 w, I = findpeak(omega_w, arpa_w.imag) 

75 assert w == pytest.approx(w_, abs=0.01) 

76 assert I == pytest.approx(I_, abs=0.01) 

77 

78 # Test that the macroscopic dielectric function can be calculated also from 

79 # the inverse dielectric function and the bare dielectric function 

80 epsinv = dfcalc.get_inverse_dielectric_function() 

81 _, _, epsM_frominv_w = epsinv.macroscopic_dielectric_function().arrays 

82 assert epsM_frominv_w == pytest.approx(epsM_w, rel=1e-6) 

83 epsbare = dfcalc.get_bare_dielectric_function() 

84 _, _, epsM_frombare_w = epsbare.macroscopic_dielectric_function().arrays 

85 assert epsM_frombare_w == pytest.approx(epsM_w, rel=1e-6) 

86 

87 # ----- TDDFT absorption spectra ----- # 

88 

89 # Absorption spectrum calculation ALDA 

90 if eshift is None: 

91 w_ = 10.7562 

92 I_ = 5.8803 

93 else: 

94 w_ = 14.7615 

95 I_ = 5.7946 

96 

97 epsinv = dfcalc.get_inverse_dielectric_function(xc='ALDA', rshelmax=0) 

98 # Here we base the check on a written results file 

99 epsinv.polarizability().write(filename='ALDA_pol.csv') 

100 dfcalc.context.comm.barrier() 

101 omega_w, a0alda_w, aalda_w = read_response_function('ALDA_pol.csv') 

102 

103 assert a0alda_w == pytest.approx(a0rpa_w, rel=1e-4) 

104 w, I = findpeak(omega_w, aalda_w.imag) 

105 assert w == pytest.approx(w_, abs=0.01) 

106 assert I == pytest.approx(I_, abs=0.1) 

107 

108 # Absorption spectrum calculation long-range kernel 

109 if eshift is None: 

110 w_ = 10.2906 

111 I_ = 5.6955 

112 else: 

113 w_ = 14.2901 

114 I_ = 5.5508 

115 

116 epsinv = dfcalc.get_inverse_dielectric_function(xc='LR0.25') 

117 omega_w, a0lr_w, alr_w = epsinv.polarizability().arrays 

118 

119 assert a0lr_w == pytest.approx(a0rpa_w, rel=1e-4) 

120 w, I = findpeak(omega_w, alr_w.imag) 

121 assert w == pytest.approx(w_, abs=0.01) 

122 assert I == pytest.approx(I_, abs=0.1) 

123 

124 # Absorption spectrum calculation Bootstrap 

125 if eshift is None: 

126 w_ = 10.4600 

127 I_ = 6.0263 

128 else: 

129 w_ = 14.2626 

130 I_ = 5.3896 

131 

132 epsinv = dfcalc.get_inverse_dielectric_function(xc='Bootstrap') 

133 omega_w, a0btsr_w, abtsr_w = epsinv.polarizability().arrays 

134 

135 assert a0btsr_w == pytest.approx(a0rpa_w, rel=1e-4) 

136 w, I = findpeak(omega_w, abtsr_w.imag) 

137 assert w == pytest.approx(w_, abs=0.02) 

138 assert I == pytest.approx(I_, abs=0.2) 

139 

140 # import matplotlib.pyplot as plt 

141 # plt.plot(omega_w, a0rpa_w.imag, label='IP') 

142 # plt.plot(omega_w, arpa_w.imag, label='RPA') 

143 # plt.plot(omega_w, aalda_w.imag, label='ALDA') 

144 # plt.plot(omega_w, alr_w.imag, label='LR0.25') 

145 # plt.plot(omega_w, abtsr_w.imag, label='Bootstrap') 

146 # plt.legend() 

147 # plt.show()