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

1import time 

2import pytest 

3import numpy as np 

4 

5from ase.build import bulk 

6from ase.parallel import parprint 

7 

8from gpaw import GPAW, PW 

9from gpaw.test import findpeak 

10from gpaw.bztools import find_high_symmetry_monkhorst_pack 

11from gpaw.response.df import DielectricFunction, read_response_function 

12from gpaw.mpi import size, world 

13 

14 

15# Affected by https://gitlab.com/gpaw/gpaw/-/issues/840 

16# We are disabling assertions below as necessary, should be reenabled 

17# after fixing 840. 

18@pytest.mark.dielectricfunction 

19@pytest.mark.tetrahedron 

20@pytest.mark.response 

21def test_response_aluminum_EELS_RPA(in_tmp_dir): 

22 assert size <= 4**3 

23 

24 # Ground state calculation 

25 

26 t1 = time.time() 

27 

28 a = 4.043 

29 atoms = bulk('Al', 'fcc', a=a) 

30 atoms.center() 

31 calc = GPAW(mode=PW(200), 

32 nbands=4, 

33 kpts=(4, 4, 4), 

34 parallel={'band': 1}, 

35 xc='LDA') 

36 

37 atoms.calc = calc 

38 atoms.get_potential_energy() 

39 calc.write('Al_gs.gpw') 

40 

41 # Generate grid compatible with tetrahedron integration 

42 kpts = find_high_symmetry_monkhorst_pack('Al_gs.gpw', 2.0) 

43 

44 # Calculate the wave functions on the new kpts grid 

45 calc = GPAW('Al_gs.gpw').fixed_density(kpts=kpts, update_fermi_level=True) 

46 calc.write('Al.gpw', 'all') 

47 

48 t2 = time.time() 

49 

50 # Excited state calculation 

51 q0_c = np.array([0., 0., 0.]) 

52 q1_c = np.array([1 / 4., 0., 0.]) 

53 w_w = np.linspace(0, 24, 241) 

54 

55 # Calculate the eels spectrum using point integration at both q-points 

56 df1 = DielectricFunction(calc='Al.gpw', frequencies=w_w, eta=0.2, ecut=50, 

57 integrationmode='point integration', 

58 hilbert=False, rate=0.2) 

59 df1.get_eels_spectrum(xc='RPA', filename='EELS_Al-PI_q0', q_c=q0_c) 

60 df1.get_eels_spectrum(xc='RPA', filename='EELS_Al-PI_q1', q_c=q1_c) 

61 

62 t3 = time.time() 

63 

64 # Calculate the eels spectrum using tetrahedron integration at q=0 

65 # NB: We skip the finite q-point, because the underlying symmetry 

66 # exploration runs excruciatingly slowly at finite q... 

67 df2 = DielectricFunction(calc='Al.gpw', eta=0.2, ecut=50, 

68 integrationmode='tetrahedron integration', 

69 hilbert=True, rate=0.2) 

70 df2.get_eels_spectrum(xc='RPA', filename='EELS_Al-TI_q0', q_c=q0_c) 

71 

72 t4 = time.time() 

73 

74 parprint('') 

75 parprint('For ground state calc, it took', (t2 - t1) / 60, 'minutes') 

76 parprint('For PI excited state calc, it took', (t3 - t2) / 60, 'minutes') 

77 parprint('For TI excited state calc, it took', (t4 - t3) / 60, 'minutes') 

78 

79 world.barrier() 

80 omegaP0_w, eels0P0_w, eelsP0_w = read_response_function('EELS_Al-PI_q0') 

81 omegaP1_w, eels0P1_w, eelsP1_w = read_response_function('EELS_Al-PI_q1') 

82 omegaT0_w, eels0T0_w, eelsT0_w = read_response_function('EELS_Al-TI_q0') 

83 

84 # calculate the 1 & 2 wpeak and ipeak values for tetra and point int. 

85 wpeak1P0, Ipeak1P0 = findpeak(omegaP0_w, eels0P0_w) 

86 wpeak2P0, Ipeak2P0 = findpeak(omegaP0_w, eelsP0_w) 

87 wpeak1P1, Ipeak1P1 = findpeak(omegaP1_w, eels0P1_w) 

88 wpeak2P1, Ipeak2P1 = findpeak(omegaP1_w, eelsP1_w) 

89 wpeak1T0, Ipeak1T0 = findpeak(omegaT0_w, eels0T0_w) 

90 wpeak2T0, Ipeak2T0 = findpeak(omegaT0_w, eelsT0_w) 

91 

92 # import matplotlib.pyplot as plt 

93 # plt.plot(omegaP0_w, eelsP0_w) 

94 # plt.plot(omegaT0_w, eelsT0_w) 

95 # plt.show() 

96 

97 # tetra and point integrators should produce similar results; however, 

98 # Al converges very slowly w.r.t. kpts so we just make sure the 

99 # values don't change and tests consistency elsewhere 

100 assert wpeak1P0 == pytest.approx(15.7111, abs=0.02) 

101 assert wpeak2P0 == pytest.approx(15.7096, abs=0.02) 

102 assert wpeak1P1 == pytest.approx(15.8402, abs=0.02) 

103 assert wpeak2P1 == pytest.approx(15.8645, abs=0.02) 

104 # assert wpeak1T0 == pytest.approx(20.2119, abs=0.02) # XXX #840 

105 # assert wpeak2T0 == pytest.approx(20.2179, abs=0.02) # XXX #840 

106 

107 assert Ipeak1P0 == pytest.approx(29.40, abs=1.) 

108 assert Ipeak2P0 == pytest.approx(27.70, abs=1.) 

109 assert Ipeak1P1 == pytest.approx(28.39, abs=1.) 

110 assert Ipeak2P1 == pytest.approx(26.89, abs=1.) 

111 # assert Ipeak1T0 == pytest.approx(46.24, abs=1.) # XXX #840 

112 # assert Ipeak2T0 == pytest.approx(44.27, abs=1.) # XXX #840