Coverage for ase / io / octopus / output.py: 97.67%

1import re 

2 

3import numpy as np 

4 

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

6 

7 

8class OctopusIOError(IOError): 

9 pass # Cannot find output files 

10 

11 

12def read_eigenvalues_file(fd): 

13 unit = None 

14 

15 for line in fd: 

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

17 if m is not None: 

18 unit = m.group(1) 

19 break 

20 line = next(fd) 

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

22 

23 # fermilevel = None 

24 kpts = [] 

25 eigs = [] 

26 occs = [] 

27 

28 for line in fd: 

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

30 if m: 

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

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

33 eigs.append({}) 

34 occs.append({}) 

35 else: 

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

37 if m is None: 

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

39 assert m is not None 

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

41 # it from the static/info instead. 

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

43 else: 

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

45 

46 if not eigs: 

47 # Only initialized if kpoint header was written 

48 eigs.append({}) 

49 occs.append({}) 

50 

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

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

53 

54 nkpts = len(kpts) 

55 nspins = len(eigs[0]) 

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

57 

58 kptsarr = np.array(kpts, float) 

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

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

61 

62 arrs = [eigsarr, occsarr] 

63 

64 for arr in arrs: 

65 arr.fill(np.nan) 

66 

67 for k in range(nkpts): 

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

69 arr[k, :, :] = [ 

70 lst[k][sp] for sp in (['--'] if nspins == 1 else ['up', 'dn']) 

71 ] 

72 

73 for arr in arrs: 

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

75 

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

77 return kptsarr, eigsarr, occsarr 

78 

79 

80def read_static_info_stress(fd): 

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

82 

83 headers = next(fd) 

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

85 for i in range(3): 

86 line = next(fd) 

87 tokens = line.split() 

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

89 stress_cv[i] = vec 

90 return stress_cv 

91 

92 

93def read_static_info_kpoints(fd): 

94 for line in fd: 

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

96 break 

97 

98 tokens = next(fd).split() 

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

100 bar = next(fd) 

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

102 

103 kpts = [] 

104 weights = [] 

105 

106 for line in fd: 

107 # Format: index kx ky kz weight 

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

109 if m is None: 

110 break 

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

112 weight = m.group(4) 

113 kpts.append(kxyz) 

114 weights.append(weight) 

115 

116 ibz_kpoints = np.array(kpts, float) 

117 kpoint_weights = np.array(weights, float) 

118 return dict(ibz_kpoints=ibz_kpoints, kpoint_weights=kpoint_weights) 

119 

120 

121def read_static_info_eigenvalues(fd, energy_unit): 

122 values_sknx = {} 

123 

124 nbands = 0 

125 fermilevel = None 

126 for line in fd: 

127 line = line.strip() 

128 if line.startswith('#'): 

129 continue 

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

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

132 if m is not None: 

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

134 break 

135 

136 tokens = line.split() 

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

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

139 occupation = float(tokens[3]) 

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

141 

142 nspins = len(values_sknx) 

143 if nspins == 1: 

144 val = [values_sknx['--']] 

145 else: 

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

147 val = np.array(val, float) 

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

149 assert remainder == 0 

150 

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

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

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

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

155 assert eps_skn.flags.contiguous 

156 d = dict( 

157 nspins=nspins, 

158 nkpts=nkpts, 

159 nbands=nbands, 

160 eigenvalues=eps_skn, 

161 occupations=occ_skn, 

162 ) 

163 if fermilevel is not None: 

164 d.update(fermi_level=fermilevel) 

165 return d 

166 

167 

168def read_static_info_energy(fd, energy_unit): 

169 def get(name): 

170 for line in fd: 

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

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

173 

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

175 

176 

177def read_static_info(fd): 

178 results = {} 

179 

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

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

182 

183 for line in fd: 

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

185 results.update(read_static_info_kpoints(fd)) 

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

187 unit = get_energy_unit(line) 

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

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

190 unit = get_energy_unit(line) 

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

192 elif line.startswith('Total stress tensor ['): 

193 if '[H/b^3]' in line: 

194 stress = read_static_info_stress(fd) 

195 stress *= Hartree / Bohr**3 

196 results.update(stress=stress) 

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

198 line = next(fd) 

199 values = line.split() 

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

201 line = next(fd) 

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

203 line = next(fd) 

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

205 # Reading Local Magnetic Moments 

206 magmoms = [] 

207 for line in fd: 

208 if line == '\n': 

209 break # there is no more thing to search for 

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

211 values = line.split() 

212 magmoms.append(float(values[-1])) 

213 results['magmoms'] = np.array(magmoms) 

214 elif line.startswith('Dipole'): 

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

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

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

218 elif line.startswith('Forces'): 

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

220 forceunit = {'[eV/A]': eV / Angstrom, '[H/b]': Hartree / Bohr}[ 

221 forceunitspec 

222 ] 

223 forces = [] 

224 line = next(fd) 

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

226 for line in fd: 

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

228 break 

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

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

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

232 

233 return results