1import re 

2import numpy as np 

3from ase.units import Bohr, Angstrom, Hartree, eV, Debye 

4 

5 

6class OctopusIOError(IOError): 

7 pass # Cannot find output files 

8 

9 

10def read_eigenvalues_file(fd): 

11 unit = None 

12 

13 for line in fd: 

14 m = re.match(r'Eigenvalues\s*\[(.+?)\]', line) 

15 if m is not None: 

16 unit = m.group(1) 

17 break 

18 line = next(fd) 

19 assert line.strip().startswith('#st'), line 

20 

21 # fermilevel = None 

22 kpts = [] 

23 eigs = [] 

24 occs = [] 

25 

26 for line in fd: 

27 m = re.match(r'#k.*?\(\s*(.+?),\s*(.+?),\s*(.+?)\)', line) 

28 if m: 

29 k = m.group(1, 2, 3) 

30 kpts.append(np.array(k, float)) 

31 eigs.append({}) 

32 occs.append({}) 

33 else: 

34 m = re.match(r'\s*\d+\s*(\S+)\s*(\S+)\s*(\S+)', line) 

35 if m is None: 

36 m = re.match(r'Fermi energy\s*=\s*(\S+)\s*', line) 

37 assert m is not None 

38 # We can also return the fermilevel but so far we just read 

39 # it from the static/info instead. 

40 # fermilevel = float(m.group(1)) 

41 else: 

42 spin, eig, occ = m.group(1, 2, 3) 

43 

44 if not eigs: 

45 # Only initialized if kpoint header was written 

46 eigs.append({}) 

47 occs.append({}) 

48 

49 eigs[-1].setdefault(spin, []).append(float(eig)) 

50 occs[-1].setdefault(spin, []).append(float(occ)) 

51 

52 nkpts = len(kpts) 

53 nspins = len(eigs[0]) 

54 nbands = len(eigs[0][spin]) 

55 

56 kptsarr = np.array(kpts, float) 

57 eigsarr = np.empty((nkpts, nspins, nbands)) 

58 occsarr = np.empty((nkpts, nspins, nbands)) 

59 

60 arrs = [eigsarr, occsarr] 

61 

62 for arr in arrs: 

63 arr.fill(np.nan) 

64 

65 for k in range(nkpts): 

66 for arr, lst in [(eigsarr, eigs), (occsarr, occs)]: 

67 arr[k, :, :] = [lst[k][sp] for sp 

68 in (['--'] if nspins == 1 else ['up', 'dn'])] 

69 

70 for arr in arrs: 

71 assert not np.isnan(arr).any() 

72 

73 eigsarr *= {'H': Hartree, 'eV': eV}[unit] 

74 return kptsarr, eigsarr, occsarr 

75 

76 

77def read_static_info_stress(fd): 

78 stress_cv = np.empty((3, 3)) 

79 

80 headers = next(fd) 

81 assert headers.strip().startswith('T_{ij}') 

82 for i in range(3): 

83 line = next(fd) 

84 tokens = line.split() 

85 vec = np.array(tokens[1:4]).astype(float) 

86 stress_cv[i] = vec 

87 return stress_cv 

88 

89 

90def read_static_info_kpoints(fd): 

91 for line in fd: 

92 if line.startswith('List of k-points'): 

93 break 

94 

95 tokens = next(fd).split() 

96 assert tokens == ['ik', 'k_x', 'k_y', 'k_z', 'Weight'] 

97 bar = next(fd) 

98 assert bar.startswith('---') 

99 

100 kpts = [] 

101 weights = [] 

102 

103 for line in fd: 

104 # Format: index kx ky kz weight 

105 m = re.match(r'\s*\d+\s*(\S+)\s*(\S+)\s*(\S+)\s*(\S+)', line) 

106 if m is None: 

107 break 

108 kxyz = m.group(1, 2, 3) 

109 weight = m.group(4) 

110 kpts.append(kxyz) 

111 weights.append(weight) 

112 

113 ibz_k_points = np.array(kpts, float) 

114 k_point_weights = np.array(weights, float) 

115 return dict(ibz_k_points=ibz_k_points, k_point_weights=k_point_weights) 

116 

117 

118def read_static_info_eigenvalues(fd, energy_unit): 

119 

120 values_sknx = {} 

121 

122 nbands = 0 

123 fermilevel = None 

124 for line in fd: 

125 line = line.strip() 

126 if line.startswith('#'): 

127 continue 

128 if not line[:1].isdigit(): 

129 m = re.match(r'Fermi energy\s*=\s*(\S+)', line) 

130 if m is not None: 

131 fermilevel = float(m.group(1)) * energy_unit 

132 break 

133 

134 tokens = line.split() 

135 nbands = max(nbands, int(tokens[0])) 

136 energy = float(tokens[2]) * energy_unit 

137 occupation = float(tokens[3]) 

138 values_sknx.setdefault(tokens[1], []).append((energy, occupation)) 

139 

140 nspins = len(values_sknx) 

141 if nspins == 1: 

142 val = [values_sknx['--']] 

143 else: 

144 val = [values_sknx['up'], values_sknx['dn']] 

145 val = np.array(val, float) 

146 nkpts, remainder = divmod(len(val[0]), nbands) 

147 assert remainder == 0 

148 

149 eps_skn = val[:, :, 0].reshape(nspins, nkpts, nbands) 

150 occ_skn = val[:, :, 1].reshape(nspins, nkpts, nbands) 

151 eps_skn = eps_skn.transpose(1, 0, 2).copy() 

152 occ_skn = occ_skn.transpose(1, 0, 2).copy() 

153 assert eps_skn.flags.contiguous 

154 d = dict(nspins=nspins, 

155 nkpts=nkpts, 

156 nbands=nbands, 

157 eigenvalues=eps_skn, 

158 occupations=occ_skn) 

159 if fermilevel is not None: 

160 d.update(efermi=fermilevel) 

161 return d 

162 

163 

164def read_static_info_energy(fd, energy_unit): 

165 def get(name): 

166 for line in fd: 

167 if line.strip().startswith(name): 

168 return float(line.split('=')[-1].strip()) * energy_unit 

169 return dict(energy=get('Total'), free_energy=get('Free')) 

170 

171 

172def read_static_info(fd): 

173 results = {} 

174 

175 def get_energy_unit(line): # Convert "title [unit]": ---> unit 

176 return {'[eV]': eV, '[H]': Hartree}[line.split()[1].rstrip(':')] 

177 

178 for line in fd: 

179 if line.strip('*').strip().startswith('Brillouin zone'): 

180 results.update(read_static_info_kpoints(fd)) 

181 elif line.startswith('Eigenvalues ['): 

182 unit = get_energy_unit(line) 

183 results.update(read_static_info_eigenvalues(fd, unit)) 

184 elif line.startswith('Energy ['): 

185 unit = get_energy_unit(line) 

186 results.update(read_static_info_energy(fd, unit)) 

187 elif line.startswith('Stress tensor'): 

188 assert line.split()[-1] == '[H/b^3]' 

189 stress = read_static_info_stress(fd) 

190 stress *= Hartree / Bohr**3 

191 results.update(stress=stress) 

192 elif line.startswith('Total Magnetic Moment'): 

193 if 0: 

194 line = next(fd) 

195 values = line.split() 

196 results['magmom'] = float(values[-1]) 

197 

198 line = next(fd) 

199 assert line.startswith('Local Magnetic Moments') 

200 line = next(fd) 

201 assert line.split() == ['Ion', 'mz'] 

202 # Reading Local Magnetic Moments 

203 mag_moment = [] 

204 for line in fd: 

205 if line == '\n': 

206 break # there is no more thing to search for 

207 line = line.replace('\n', ' ') 

208 values = line.split() 

209 mag_moment.append(float(values[-1])) 

210 

211 results['magmoms'] = np.array(mag_moment) 

212 elif line.startswith('Dipole'): 

213 assert line.split()[-1] == '[Debye]' 

214 dipole = [float(next(fd).split()[-1]) for i in range(3)] 

215 results['dipole'] = np.array(dipole) * Debye 

216 elif line.startswith('Forces'): 

217 forceunitspec = line.split()[-1] 

218 forceunit = {'[eV/A]': eV / Angstrom, 

219 '[H/b]': Hartree / Bohr}[forceunitspec] 

220 forces = [] 

221 line = next(fd) 

222 assert line.strip().startswith('Ion') 

223 for line in fd: 

224 if line.strip().startswith('---'): 

225 break 

226 tokens = line.split()[-3:] 

227 forces.append([float(f) for f in tokens]) 

228 results['forces'] = np.array(forces) * forceunit 

229 elif line.startswith('Fermi'): 

230 tokens = line.split() 

231 unit = {'eV': eV, 'H': Hartree}[tokens[-1]] 

232 eFermi = float(tokens[-2]) * unit 

233 results['efermi'] = eFermi 

234 

235 if 'ibz_k_points' not in results: 

236 results['ibz_k_points'] = np.zeros((1, 3)) 

237 results['k_point_weights'] = np.ones(1) 

238 if 0: # 'efermi' not in results: 

239 # Find HOMO level. Note: This could be a very bad 

240 # implementation with fractional occupations if the Fermi 

241 # level was not found otherwise. 

242 all_energies = results['eigenvalues'].ravel() 

243 all_occupations = results['occupations'].ravel() 

244 args = np.argsort(all_energies) 

245 for arg in args[::-1]: 

246 if all_occupations[arg] > 0.1: 

247 break 

248 eFermi = all_energies[arg] 

249 results['efermi'] = eFermi 

250 

251 return results