Coverage for ase / io / gpaw_out.py: 83.73%

1# fmt: off 

2 

3import re 

4 

5import numpy as np 

6 

7from ase.atoms import Atoms 

8from ase.calculators.singlepoint import ( 

9 SinglePointDFTCalculator, 

10 SinglePointKPoint, 

11) 

12 

13 

14def index_startswith(lines: list[str], string: str) -> int: 

15 for i, line in enumerate(lines): 

16 if line.strip().startswith(string): 

17 return i 

18 raise ValueError 

19 

20 

21def index_pattern(lines: list[str], pattern: str) -> int: 

22 repat = re.compile(pattern) 

23 for i, line in enumerate(lines): 

24 if repat.match(line): 

25 return i 

26 raise ValueError 

27 

28 

29def read_forces(lines: list[str], 

30 ii: int, 

31 atoms: Atoms) -> tuple[list[tuple[float, float, float]], int]: 

32 f = [] 

33 for i in range(ii + 1, ii + 1 + len(atoms)): 

34 try: 

35 x, y, z = lines[i].split()[-3:] 

36 f.append((float(x), float(y), float(z))) 

37 except (ValueError, IndexError) as m: 

38 raise OSError(f'Malformed GPAW log file: {m}') 

39 return f, i 

40 

41 

42def read_stresses(lines: list[str], 

43 ii: int,) -> tuple[list[tuple[float, float, float]], int]: 

44 s = [] 

45 for i in range(ii + 1, ii + 4): 

46 try: 

47 x, y, z = lines[i].split()[-3:] 

48 s.append((float(x), float(y), float(z))) 

49 except (ValueError, IndexError) as m: 

50 raise OSError(f'Malformed GPAW log file: {m}') 

51 return s, i 

52 

53 

54def read_gpaw_out(fileobj, index): # -> Union[Atoms, List[Atoms]]: 

55 """Read text output from GPAW calculation.""" 

56 lines = [line.lower() for line in fileobj.readlines()] 

57 

58 # read charge 

59 try: 

60 ii = index_startswith(lines, 'total charge:') 

61 except ValueError: 

62 q = None 

63 else: 

64 q = float(lines[ii].split()[2]) 

65 

66 blocks = [] 

67 i1 = 0 

68 for i2, line in enumerate(lines): 

69 if line == 'positions:\n': 

70 if i1 > 0: 

71 blocks.append(lines[i1:i2]) 

72 i1 = i2 

73 blocks.append(lines[i1:]) 

74 

75 images: list[Atoms] = [] 

76 for lines in blocks: 

77 try: 

78 i = lines.index('unit cell:\n') 

79 except ValueError: 

80 pass 

81 else: 

82 if lines[i + 2].startswith(' -'): 

83 del lines[i + 2] # old format 

84 cell: list[float | list[float]] = [] 

85 pbc = [] 

86 for line in lines[i + 2:i + 5]: 

87 words = line.split() 

88 if len(words) == 5: # old format 

89 cell.append(float(words[2])) 

90 pbc.append(words[1] == 'yes') 

91 else: # new format with GUC 

92 cell.append([float(word) for word in words[3:6]]) 

93 pbc.append(words[2] == 'yes') 

94 

95 symbols = [] 

96 positions = [] 

97 magmoms = [] 

98 for line in lines[1:]: 

99 words = line.split() 

100 if len(words) < 5: 

101 break 

102 n, symbol, x, y, z = words[:5] 

103 symbols.append(symbol.split('.')[0].title()) 

104 positions.append([float(x), float(y), float(z)]) 

105 if len(words) > 5: 

106 mom = float(words[-1].rstrip(')')) 

107 magmoms.append(mom) 

108 if len(symbols): 

109 atoms = Atoms(symbols=symbols, positions=positions, 

110 cell=cell, pbc=pbc) 

111 else: 

112 atoms = Atoms(cell=cell, pbc=pbc) 

113 

114 try: 

115 ii = index_pattern(lines, '\\d+ k-point') 

116 word = lines[ii].split() 

117 kx = int(word[2]) 

118 ky = int(word[4]) 

119 kz = int(word[6]) 

120 bz_kpts = (kx, ky, kz) 

121 ibz_kpts = int(lines[ii + 1].split()[0]) 

122 except (ValueError, TypeError, IndexError): 

123 bz_kpts = None 

124 ibz_kpts = None 

125 

126 try: 

127 i = index_startswith(lines, 'energy contributions relative to') + 2 

128 except ValueError: 

129 e = energy_contributions = None 

130 else: 

131 energy_contributions = {} 

132 while i < len(lines): 

133 fields = lines[i].split(':') 

134 if len(fields) == 2: 

135 name = fields[0] 

136 energy = float(fields[1]) 

137 energy_contributions[name] = energy 

138 if name in ['zero kelvin', 'extrapolated']: 

139 e = energy 

140 break 

141 i += 1 

142 else: # no break 

143 raise ValueError 

144 

145 try: 

146 ii = index_pattern(lines, '(fixed )?fermi level(s)?:') 

147 except ValueError: 

148 eFermi = None 

149 else: 

150 fields = lines[ii].split() 

151 try: 

152 def strip(string): 

153 for rubbish in '[],': 

154 string = string.replace(rubbish, '') 

155 return string 

156 eFermi = [float(strip(fields[-2])), 

157 float(strip(fields[-1]))] 

158 except ValueError: 

159 eFermi = float(fields[-1]) 

160 

161 # read Eigenvalues and occupations 

162 ii1 = ii2 = 1e32 

163 try: 

164 ii1 = index_startswith(lines, 'band eigenvalues occupancy') 

165 except ValueError: 

166 pass 

167 try: 

168 ii2 = index_startswith(lines, 'band eigenvalues occupancy') 

169 except ValueError: 

170 pass 

171 ii = min(ii1, ii2) 

172 if ii == 1e32: 

173 kpts = None 

174 else: 

175 ii += 1 

176 words = lines[ii].split() 

177 vals = [] 

178 while len(words) > 2: 

179 vals.append([float(w) for w in words]) 

180 ii += 1 

181 words = lines[ii].split() 

182 vals = np.array(vals).transpose() 

183 kpts = [SinglePointKPoint(1, 0, 0)] 

184 kpts[0].eps_n = vals[1] 

185 kpts[0].f_n = vals[2] 

186 if vals.shape[0] > 3: 

187 kpts.append(SinglePointKPoint(1, 1, 0)) 

188 kpts[1].eps_n = vals[3] 

189 kpts[1].f_n = vals[4] 

190 # read dipole moment 

191 try: 

192 ii = index_startswith(lines, 'dipole moment:') 

193 except ValueError: 

194 dipole = None 

195 else: 

196 line = lines[ii] 

197 for x in '()[],': 

198 line = line.replace(x, '') 

199 dipole = np.array([float(c) for c in line.split()[2:5]]) 

200 

201 try: 

202 ii = index_startswith(lines, 'local magnetic moments') 

203 except ValueError: 

204 pass 

205 else: 

206 magmoms = [] 

207 for j in range(ii + 1, ii + 1 + len(atoms)): 

208 line = lines[j] 

209 if '#' in line: # new GPAW format 

210 magmom = line.split()[-4].split(']')[0] 

211 else: 

212 magmom = line.split()[-1].rstrip(')') 

213 magmoms.append(float(magmom)) 

214 

215 try: 

216 ii = lines.index('forces in ev/ang:\n') 

217 except ValueError: 

218 f = None 

219 else: 

220 f, i = read_forces(lines, ii, atoms) 

221 

222 try: 

223 ii = lines.index('stress tensor:\n') 

224 except ValueError: 

225 stress_tensor = None 

226 else: 

227 stress_tensor, i = read_stresses(lines, ii) 

228 

229 try: 

230 parameters = {} 

231 ii = index_startswith(lines, 'vdw correction:') 

232 except ValueError: 

233 name = 'gpaw' 

234 else: 

235 name = lines[ii - 1].strip() 

236 # save uncorrected values 

237 parameters.update({ 

238 'calculator': 'gpaw', 

239 'uncorrected_energy': e, 

240 }) 

241 line = lines[ii + 1] 

242 assert line.startswith('energy:') 

243 e = float(line.split()[-1]) 

244 f, i = read_forces(lines, ii + 3, atoms) 

245 

246 if len(images) > 0 and e is None: 

247 break 

248 

249 if q is not None and len(atoms) > 0: 

250 n = len(atoms) 

251 atoms.set_initial_charges([q / n] * n) 

252 

253 if magmoms: 

254 atoms.set_initial_magnetic_moments(magmoms) 

255 else: 

256 magmoms = None 

257 if e is not None or f is not None: 

258 calc = SinglePointDFTCalculator(atoms, energy=e, forces=f, 

259 dipole=dipole, magmoms=magmoms, 

260 efermi=eFermi, 

261 bzkpts=bz_kpts, ibzkpts=ibz_kpts, 

262 stress=stress_tensor) 

263 calc.name = name 

264 calc.parameters = parameters 

265 if energy_contributions is not None: 

266 calc.energy_contributions = energy_contributions 

267 if kpts is not None: 

268 calc.kpts = kpts 

269 atoms.calc = calc 

270 

271 images.append(atoms) 

272 

273 if len(images) == 0: 

274 raise OSError('Corrupted GPAW-text file!') 

275 

276 return images[index]