Coverage for gpaw/response/gw_bands.py: 9%

1import pickle 

2 

3import numpy as np 

4from ase.utils import gcd 

5from scipy.interpolate import InterpolatedUnivariateSpline 

6 

7from gpaw import GPAW 

8from gpaw.spinorbit import soc_eigenstates 

9from gpaw.kpt_descriptor import to1bz 

10 

11 

12class GWBands: 

13 """This class defines the GW_bands properties""" 

14 

15 def __init__(self, 

16 calc=None, 

17 comm=None, 

18 gw_file=None, 

19 kpoints=None, 

20 bandrange=None): 

21 

22 self.calc = GPAW(calc, communicator=comm, txt=None) 

23 if gw_file is not None: 

24 self.gw_file = pickle.load(open(gw_file, 'rb'), encoding='bytes') 

25 self.kpoints = kpoints 

26 if bandrange is None: 

27 self.bandrange = np.arange(self.calc.get_number_of_bands()) 

28 else: 

29 self.bandrange = bandrange 

30 

31 self.gd = self.calc.wfs.gd.new_descriptor() 

32 self.kd = self.calc.wfs.kd 

33 

34 self.acell_cv = self.gd.cell_cv 

35 self.bcell_cv = 2 * np.pi * self.gd.icell_cv 

36 

37 def find_k_along_path(self, plot_BZ=True): 

38 """Finds the k-points along the bandpath present in the 

39 original calculation""" 

40 kd = self.kd 

41 acell_cv = self.acell_cv 

42 bcell_cv = self.bcell_cv 

43 kpoints = self.kpoints 

44 

45 if plot_BZ: 

46 """Plotting the points in the Brillouin Zone""" 

47 kp_1bz = to1bz(kd.bzk_kc, acell_cv) 

48 

49 bzk_kcv = np.dot(kd.bzk_kc, bcell_cv) 

50 kp_1bz_v = np.dot(kp_1bz, bcell_cv) 

51 

52 import matplotlib.pyplot as plt 

53 plt.plot(bzk_kcv[:, 0], bzk_kcv[:, 1], 'xg') 

54 plt.plot(kp_1bz_v[:, 0], kp_1bz_v[:, 1], 'ob') 

55 for ik in range(1, len(kpoints)): 

56 kpoint1_v = np.dot(kpoints[ik], bcell_cv) 

57 kpoint2_v = np.dot(kpoints[ik - 1], bcell_cv) 

58 plt.plot([kpoint1_v[0], kpoint2_v[0]], [kpoint1_v[1], 

59 kpoint2_v[1]], '--vr') 

60 

61 """Finding the points along given directions""" 

62 print('Finding the kpoints along the path') 

63 N_c = kd.N_c 

64 wpts_xc = kpoints 

65 

66 x_x = [] 

67 k_xc = [] 

68 k_x = [] 

69 x = 0. 

70 X = [] 

71 for nwpt in range(1, len(wpts_xc)): 

72 X.append(x) 

73 to_c = wpts_xc[nwpt] 

74 from_c = wpts_xc[nwpt - 1] 

75 vec_c = to_c - from_c 

76 print('From ', from_c, ' to ', to_c) 

77 Nv_c = (vec_c * N_c).round().astype(int) 

78 Nv = abs(gcd(gcd(Nv_c[0], Nv_c[1]), Nv_c[2])) 

79 print(Nv, ' points found') 

80 dv_c = vec_c / Nv 

81 dv_v = np.dot(dv_c, bcell_cv) 

82 dx = np.linalg.norm(dv_v) 

83 if nwpt == len(wpts_xc) - 1: 

84 # X.append(Nv * dx) 

85 Nv += 1 

86 for n in range(Nv): 

87 k_c = from_c + n * dv_c 

88 bzk_c = to1bz(np.array([k_c]), acell_cv)[0] 

89 ikpt = kd.where_is_q(bzk_c, kd.bzk_kc) 

90 x_x.append(x) 

91 k_xc.append(k_c) 

92 k_x.append(ikpt) 

93 x += dx 

94 X.append(x_x[-1]) 

95 if plot_BZ is True: 

96 for ik in range(len(k_xc)): 

97 ktemp_xcv = np.dot(k_xc[ik], bcell_cv) 

98 plt.plot(ktemp_xcv[0], ktemp_xcv[1], 'xr', markersize=10) 

99 plt.show() 

100 

101 return x_x, k_xc, k_x, X 

102 

103 def get_dft_eigenvalues(self): 

104 Nk = len(self.calc.get_ibz_k_points()) 

105 bands = np.arange(self.bandrange[0], self.bandrange[-1]) 

106 e_kn = np.array([self.calc.get_eigenvalues(kpt=k)[bands] 

107 for k in range(Nk)]) 

108 return e_kn 

109 

110 def get_vacuum_level(self, plot_pot=False): 

111 """Finds the vacuum level through Hartree potential""" 

112 vHt_g = self.calc.get_electrostatic_potential() 

113 vHt_z = np.mean(np.mean(vHt_g, axis=0), axis=0) 

114 

115 if plot_pot: 

116 import matplotlib.pyplot as plt 

117 plt.plot(vHt_z) 

118 plt.show() 

119 return vHt_z[0] 

120 

121 def get_spinorbit_corrections(self, return_spin=True, return_wfs=False, 

122 bands=None, dft=False, 

123 eig_file=None): 

124 """Gets the spinorbit corrections to the eigenvalues""" 

125 calc = self.calc 

126 bandrange = self.bandrange 

127 

128 if not dft: 

129 try: 

130 e_kn = self.gw_file['qp'][0] 

131 except KeyError: 

132 e_kn = self.gw_file[b'qp'][0] 

133 else: 

134 if eig_file is not None: 

135 e_kn = pickle.load(open(eig_file))[0] 

136 else: 

137 e_kn = self.get_dft_eigenvalues()[ 

138 :, bandrange[0]:bandrange[-1] + 1] 

139 

140 # this will fail - please write a test! 

141 soc = soc_eigenstates( 

142 calc, 

143 n1=bandrange[0], 

144 n2=bandrange[-1] + 1, 

145 eigenvalues=e_kn[np.newaxis]) 

146 eSO_nk = soc.eigenvalues().T 

147 e_kn = eSO_nk.T 

148 return e_kn 

149 

150 def get_gw_bands(self, nk_Int=50, interpolate=False, SO=False, 

151 dft=False, eig_file=None, vac=False): 

152 """Getting Eigenvalues along the path""" 

153 kd = self.kd 

154 if SO: 

155 e_kn = self.get_spinorbit_corrections(return_wfs=True, 

156 dft=dft, 

157 eig_file=eig_file) 

158 elif eig_file is not None: 

159 e_kn = pickle.load(open(eig_file))[0] 

160 else: 

161 if not dft: 

162 try: 

163 e_kn = self.gw_file['qp'][0] 

164 except KeyError: 

165 e_kn = self.gw_file[b'qp'][0] 

166 else: 

167 e_kn = self.get_dft_eigenvalues() 

168 e_kn = np.sort(e_kn, axis=1) 

169 

170 bandrange = self.bandrange 

171 ef = self.calc.get_fermi_level() 

172 if vac: 

173 evac = self.get_vacuum_level() 

174 else: 

175 evac = 0.0 

176 x_x, k_xc, k_x, X = self.find_k_along_path(plot_BZ=False) 

177 

178 k_ibz_x = np.zeros_like(k_x) 

179 eGW_kn = np.zeros((len(k_x), e_kn.shape[1])) 

180 for n in range(e_kn.shape[1]): 

181 for ik in range(len(k_x)): 

182 ibzkpt = kd.bz2ibz_k[k_x[ik]] 

183 k_ibz_x[ik] = ibzkpt 

184 eGW_kn[ik, n] = e_kn[ibzkpt, n] 

185 

186 N_occ = (eGW_kn[0] < ef).sum() 

187 print(N_occ, bandrange[0]) 

188 # N_occ = int(self.calc.get_number_of_electrons()/2) 

189 print(' ') 

190 if SO: 

191 print('The number of Occupied bands is:', N_occ + 2 * bandrange[0]) 

192 else: 

193 print('The number of Occupied bands is:', N_occ + bandrange[0]) 

194 gap = (eGW_kn[:, N_occ].min() - eGW_kn[:, N_occ - 1].max()) 

195 print('The bandgap is: %f' % gap) 

196 

197 vbm = eGW_kn[:, N_occ - 1].max() - evac 

198 cbm = eGW_kn[:, N_occ].min() - evac 

199 

200 if interpolate: 

201 xfit_k = np.linspace(x_x[0], x_x[-1], nk_Int) 

202 xfit_k = np.append(xfit_k, x_x) 

203 xfit_k = np.sort(xfit_k) 

204 nk_Int = len(xfit_k) 

205 efit_kn = np.zeros((nk_Int, eGW_kn.shape[1])) 

206 for n in range(eGW_kn.shape[1]): 

207 fit_e = InterpolatedUnivariateSpline(x_x, eGW_kn[:, n]) 

208 efit_kn[:, n] = fit_e(xfit_k) 

209 

210 results = {'x_k': xfit_k, 

211 'X': X, 

212 'e_kn': efit_kn - evac, 

213 'ef': ef - evac, 

214 'gap': gap, 

215 'vbm': vbm, 

216 'cbm': cbm} 

217 

218 return results 

219 else: 

220 results = {'x_k': x_x, 

221 'X': X, 

222 'k_ibz_x': k_ibz_x, 

223 'e_kn': eGW_kn - evac, 

224 'ef': ef - evac, 

225 'gap': gap, 

226 'vbm': vbm, 

227 'cbm': cbm} 

228 return results