Coverage for gpaw/test/response/test_na_plasmon.py: 97%

1import pytest 

2import numpy as np 

3 

4from ase import Atoms 

5from gpaw import GPAW, PW 

6from gpaw.mpi import world 

7from gpaw.test import findpeak 

8from gpaw.utilities import compiled_with_sl 

9from gpaw.response.df import DielectricFunction 

10from gpaw.response.symmetry import QSymmetryAnalyzer 

11 

12# Comparing the plasmon peaks found in bulk sodium for two different 

13# atomic structures. Testing for identical plasmon peaks. Not using 

14# physical sodium cell. 

15 

16 

17@pytest.mark.dielectricfunction 

18@pytest.mark.response 

19def test_response_na_plasmon(in_tmp_dir): 

20 a = 4.23 / 2.0 

21 a1 = Atoms('Na', 

22 scaled_positions=[[0, 0.1, 0]], 

23 cell=(a, a, a), 

24 pbc=True) 

25 

26 # parallel calculations must have domain = 1 

27 parallel = {'band': 1} 

28 if world.size > 1 and compiled_with_sl(): 

29 parallel.update({'domain': 1}) 

30 

31 # Expanding along x-direction 

32 a2 = Atoms('Na2', 

33 scaled_positions=[[0, 0.1, 0], [0.5, 0.1, 0]], 

34 cell=(2 * a, a, a), 

35 pbc=True) 

36 

37 a1.calc = GPAW(mode=PW(250), 

38 kpts={'size': (4, 4, 4), 'gamma': True}, 

39 parallel=parallel, 

40 # txt='small.txt', 

41 ) 

42 

43 # Kpoint sampling should be halved in the expanded direction. 

44 a2.calc = GPAW(mode=PW(250), 

45 kpts={'size': (2, 4, 4), 'gamma': True}, 

46 parallel=parallel, 

47 # txt='large.txt', 

48 ) 

49 

50 a1.get_potential_energy() 

51 a2.get_potential_energy() 

52 

53 # Use twice as many bands for expanded structure 

54 a1.calc.diagonalize_full_hamiltonian(nbands=20) 

55 a2.calc.diagonalize_full_hamiltonian(nbands=40) 

56 

57 a1.calc.write('gs_Na_small.gpw', 'all') 

58 a2.calc.write('gs_Na_large.gpw', 'all') 

59 

60 # Settings that should yield the same result 

61 settings = [ 

62 {'qsymmetry': QSymmetryAnalyzer(pointgroup, timerev)} 

63 for pointgroup in [False, True] 

64 for timerev in [False, True]] 

65 

66 # Test block parallelization (needs scalapack) 

67 if world.size > 1 and compiled_with_sl(): 

68 settings.append({'qsymmetry': True, 'nblocks': 2}) 

69 

70 # Calculate the dielectric functions 

71 dfs0 = [] # Arrays to check for self-consistency 

72 dfs1 = [] 

73 dfs2 = [] 

74 dfs3 = [] 

75 dfs4 = [] 

76 dfs5 = [] 

77 

78 # list of intensities to compare against. Intensity values matched 

79 # to 10-3 w/ higher tol. Speeding up test degraded the agreement to 10-2 

80 # Added additional intensity difference check test with tol 10-3 

81 I_diffs = {'x': [0.008999, 0.007558, 0.008999, 0.007558, 0.006101], 

82 'y': [0.004063, 0.004063, 0.004063, 0.004063, 0.004063], 

83 'z': [0.004063, 0.005689, 0.004244, 0.005689, 0.005689]} 

84 for idx, kwargs in enumerate(settings): 

85 df1 = DielectricFunction('gs_Na_small.gpw', 

86 ecut=40, 

87 rate=0.001, 

88 **kwargs) 

89 

90 df1NLFCx, df1LFCx = df1.get_dielectric_function(direction='x') 

91 df1NLFCy, df1LFCy = df1.get_dielectric_function(direction='y') 

92 df1NLFCz, df1LFCz = df1.get_dielectric_function(direction='z') 

93 

94 df2 = DielectricFunction('gs_Na_large.gpw', 

95 ecut=40, 

96 rate=0.001, 

97 **kwargs) 

98 

99 df2NLFCx, df2LFCx = df2.get_dielectric_function(direction='x') 

100 df2NLFCy, df2LFCy = df2.get_dielectric_function(direction='y') 

101 df2NLFCz, df2LFCz = df2.get_dielectric_function(direction='z') 

102 

103 dfs0.append(df1NLFCx) 

104 dfs1.append(df1LFCx) 

105 dfs2.append(df1NLFCy) 

106 dfs3.append(df1LFCy) 

107 dfs4.append(df1NLFCz) 

108 dfs5.append(df1LFCz) 

109 

110 # Compare plasmon frequencies and intensities: x, y, z 

111 # x values 

112 w_w = df1.chi0calc.wd.omega_w 

113 w1, I1 = findpeak(w_w, -(1. / df1LFCx).imag) 

114 w2, I2 = findpeak(w_w, -(1. / df2LFCx).imag) 

115 I_diff = abs(I1 - I2) 

116 # test that the frequency for 2 settings are aprx equal 

117 assert w1 == pytest.approx(w2, abs=1e-2) 

118 # test that the intensity difference is within some tol 

119 assert I_diff == pytest.approx(I_diffs['x'][idx], abs=5e-3) 

120 # test that the intensities are aprx equal 

121 assert I1 == pytest.approx(I2, abs=1e-2) 

122 

123 # y values 

124 w1, I1 = findpeak(w_w, -(1. / df1LFCy).imag) 

125 w2, I2 = findpeak(w_w, -(1. / df2LFCy).imag) 

126 I_diff = abs(I1 - I2) 

127 assert w1 == pytest.approx(w2, abs=1e-2) 

128 assert I_diff == pytest.approx(I_diffs['y'][idx], abs=5e-3) 

129 assert I1 == pytest.approx(I2, abs=1e-2) 

130 

131 # z values 

132 w1, I1 = findpeak(w_w, -(1. / df1LFCz).imag) 

133 w2, I2 = findpeak(w_w, -(1. / df2LFCz).imag) 

134 I_diff = abs(I1 - I2) 

135 assert w1 == pytest.approx(w2, abs=1e-2) 

136 assert I_diff == pytest.approx(I_diffs['z'][idx], abs=5e-3) 

137 assert I1 == pytest.approx(I2, abs=1e-2) 

138 

139 # Check for self-consistency 

140 for i, dfs in enumerate([dfs0, dfs1, dfs2, dfs3, dfs4, dfs5]): 

141 while len(dfs): 

142 df = dfs.pop() 

143 for df2 in dfs: 

144 assert np.max(np.abs((df - df2) / df)) < 2e-3