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

1import pytest 

2import numpy as np 

3 

4from gpaw import GPAW 

5from gpaw.test import findpeak 

6from gpaw.bztools import find_high_symmetry_monkhorst_pack 

7from gpaw.response.df import DielectricFunction, read_response_function 

8from gpaw.mpi import world 

9 

10 

11@pytest.mark.dielectricfunction 

12@pytest.mark.tetrahedron 

13@pytest.mark.response 

14def test_response_Na_EELS_RPA_tetra_point_comparison(in_tmp_dir, gpw_files): 

15 gpwname = gpw_files['na_chain'] 

16 

17 # Generate grid compatible with tetrahedron integration 

18 kpts = find_high_symmetry_monkhorst_pack(gpwname, 6.0) 

19 

20 # Calculate the wave functions on the new kpts grid 

21 calc = GPAW(gpwname).fixed_density(kpts=kpts, update_fermi_level=True) 

22 calc.write('Na', 'all') 

23 

24 # Excited state calculation 

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

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

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

28 

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

30 df1 = DielectricFunction(calc='Na', frequencies=w_w, eta=0.2, ecut=50, 

31 hilbert=False, rate=0.2) 

32 df1.get_eels_spectrum(xc='RPA', filename='EELS_Na-PI_q0', q_c=q0_c) 

33 df1.get_eels_spectrum(xc='RPA', filename='EELS_Na-PI_q1', q_c=q1_c) 

34 

35 # Calculate the eels spectrum using tetrahedron integration at both q 

36 df2 = DielectricFunction(calc='Na', eta=0.2, ecut=50, 

37 integrationmode='tetrahedron integration', 

38 hilbert=True, rate=0.2) 

39 df2.get_eels_spectrum(xc='RPA', filename='EELS_Na-TI_q0', q_c=q0_c) 

40 df2.get_eels_spectrum(xc='RPA', filename='EELS_Na-TI_q1', q_c=q1_c) 

41 

42 world.barrier() 

43 omegaP0_w, _, eelsP0_w = read_response_function('EELS_Na-PI_q0') 

44 omegaP1_w, _, eelsP1_w = read_response_function('EELS_Na-PI_q1') 

45 omegaT0_w, _, eelsT0_w = read_response_function('EELS_Na-TI_q0') 

46 omegaT1_w, _, eelsT1_w = read_response_function('EELS_Na-TI_q1') 

47 

48 # calculate tetra & point int. peaks and intensities 

49 wpeakP0, IpeakP0 = findpeak(omegaP0_w, eelsP0_w) 

50 wpeakP1, IpeakP1 = findpeak(omegaP1_w, eelsP1_w) 

51 wpeakT0, IpeakT0 = findpeak(omegaT0_w, eelsT0_w) 

52 wpeakT1, IpeakT1 = findpeak(omegaT1_w, eelsT1_w) 

53 

54 # import matplotlib.pyplot as plt 

55 # plt.subplot(1, 2, 1) 

56 # plt.plot(omegaP0_w, eelsP0_w) 

57 # plt.plot(omegaT0_w, eelsT0_w) 

58 # plt.subplot(1, 2, 2) 

59 # plt.plot(omegaP1_w, eelsP1_w) 

60 # plt.plot(omegaT1_w, eelsT1_w) 

61 # plt.show() 

62 

63 # tetra and point integrators should produce similar results: 

64 # confirm this by comparing the 2 integration methods 

65 assert wpeakT0 == pytest.approx(wpeakP0, abs=0.08) 

66 assert wpeakT1 == pytest.approx(wpeakP1, abs=0.12) 

67 

68 # ensure the wpeak for point & tetra integration do not change 

69 assert wpeakP0 == pytest.approx(3.4811, abs=0.02) 

70 assert wpeakP1 == pytest.approx(3.8076, abs=0.02) 

71 assert wpeakT0 == pytest.approx(3.54, abs=0.02) 

72 assert wpeakT1 == pytest.approx(3.79, abs=0.13) 

73 

74 # ensure the Ipeak for point & tetra integration do not change 

75 assert IpeakP0 == pytest.approx(8.6311, abs=1.) 

76 assert IpeakP1 == pytest.approx(7.7766, abs=1.) 

77 assert IpeakT0 == pytest.approx(8.77, abs=1.) 

78 assert IpeakT1 == pytest.approx(7.51, abs=1.)