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

1from ase import Atoms 

2from scipy.ndimage import gaussian_filter1d 

3import numpy as np 

4import pytest 

5 

6from gpaw import GPAW, FermiDirac 

7from gpaw import PW 

8from gpaw.bztools import find_high_symmetry_monkhorst_pack 

9from gpaw.response.df import DielectricFunction 

10from gpaw.test import findpeak 

11 

12 

13@pytest.mark.tetrahedron 

14@pytest.mark.dielectricfunction 

15@pytest.mark.response 

16def test_point_tetra_match(in_tmp_dir): 

17 

18 gs_file = 'gs.gpw' 

19 response_file = 'gsresponse.gpw' 

20 

21 # Create graphene lattice 

22 a = 2.5 

23 atoms = Atoms( 

24 symbols='C2', positions=[(0.5 * a, -np.sqrt(3) / 6 * a, 0.0), 

25 (0.5 * a, np.sqrt(3) / 6 * a, 0.0)], 

26 cell=[(0.5 * a, -0.5 * 3**0.5 * a, 0), 

27 (0.5 * a, 0.5 * 3**0.5 * a, 0), 

28 (0.0, 0.0, 3.2)], 

29 pbc=[True, True, False]) 

30 

31 atoms.center(axis=2) 

32 

33 calc = GPAW(mode=PW(400), 

34 kpts={'density': 10.0, 'gamma': True}, 

35 occupations=FermiDirac(0.1)) 

36 

37 atoms.calc = calc 

38 atoms.get_potential_energy() 

39 calc.write(gs_file) 

40 

41 density = 15 

42 kpts = find_high_symmetry_monkhorst_pack(gs_file, density=density) 

43 responseGS = GPAW(gs_file).fixed_density( 

44 kpts=kpts, 

45 parallel={'band': 1}, 

46 nbands=20, 

47 occupations=FermiDirac(0.001), 

48 convergence={'bands': 10}) 

49 responseGS.write(response_file, 'all') 

50 

51 # DF with tetra integration 

52 df = DielectricFunction(response_file, 

53 eta=25e-3, 

54 rate='eta', 

55 frequencies={'type': 'nonlinear', 

56 'domega0': 0.1}, 

57 integrationmode='tetrahedron integration') 

58 df1_tetra, df2_tetra = df.get_dielectric_function(q_c=[0, 0, 0]) 

59 

60 # DF with point integration 

61 df = DielectricFunction(response_file, 

62 frequencies={'type': 'nonlinear', 

63 'domega0': 0.1}, 

64 eta=25e-3, 

65 rate='eta') 

66 df1_point, df2_point = df.get_dielectric_function(q_c=[0, 0, 0]) 

67 

68 omega = df.get_frequencies() 

69 df2_tetra = np.imag(df2_tetra) 

70 df2_point = np.imag(df2_point) 

71 

72 # Do not use frequencies near the w=0 singularity 

73 slicer = [(freq >= 1.5) and (freq <= 20) for freq in omega] 

74 # Convolve with Gaussian to smoothen the curve 

75 sigma = 1.9 

76 df2_gauss = gaussian_filter1d(df2_point[slicer], sigma) 

77 

78 rms_diff_tetra_point = np.sqrt(np.sum((df2_tetra[slicer] 

79 - df2_point[slicer])**2) / np.sum(slicer)) 

80 rms_diff_tetra_gauss = np.sqrt(np.sum((df2_tetra[slicer] 

81 - df2_gauss)**2) / np.sum(slicer)) 

82 

83 assert rms_diff_tetra_point < 1.35 

84 assert rms_diff_tetra_gauss < 1.10 

85 assert rms_diff_tetra_point * 0.80 > rms_diff_tetra_gauss 

86 

87 freq1, amp1 = findpeak(omega[slicer], df2_tetra[slicer]) 

88 freq2, amp2 = findpeak(omega[slicer], df2_gauss) 

89 

90 assert freq1 - freq2 < 0.1 

91 

92 # # Plot the figures 

93 # import matplotlib.pyplot as plt 

94 # plt.figure(figsize=(6, 6)) 

95 # plt.plot(omega[slicer], df2_tetra[slicer], label='Tetrahedron') 

96 # plt.plot(omega[slicer], df2_point[slicer], label='Point sampling') 

97 # plt.plot(omega[slicer], df2_gauss, label='Gaussian convolution') 

98 # plt.xlabel('Frequency (eV)') 

99 # plt.ylabel('$\\mathrm{Im}\\varepsilon$') 

100 # plt.xlim(0, 6) 

101 # plt.ylim(0, 20) 

102 # plt.legend() 

103 # plt.tight_layout() 

104 # plt.savefig('graphene_eps.png', dpi=300)