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

1# fmt: off 

2 

3import re 

4from typing import List, Tuple, Union 

5 

6import numpy as np 

7 

8from ase.atoms import Atoms 

9from ase.calculators.singlepoint import ( 

10 SinglePointDFTCalculator, 

11 SinglePointKPoint, 

12) 

13 

14 

15def index_startswith(lines: List[str], string: str) -> int: 

16 for i, line in enumerate(lines): 

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

18 return i 

19 raise ValueError 

20 

21 

22def index_pattern(lines: List[str], pattern: str) -> int: 

23 repat = re.compile(pattern) 

24 for i, line in enumerate(lines): 

25 if repat.match(line): 

26 return i 

27 raise ValueError 

28 

29 

30def read_forces(lines: List[str], 

31 ii: int, 

32 atoms: Atoms) -> Tuple[List[Tuple[float, float, float]], int]: 

33 f = [] 

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

35 try: 

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

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

38 except (ValueError, IndexError) as m: 

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

40 return f, i 

41 

42 

43def read_stresses(lines: List[str], 

44 ii: int,) -> Tuple[List[Tuple[float, float, float]], int]: 

45 s = [] 

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

47 try: 

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

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

50 except (ValueError, IndexError) as m: 

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

52 return s, i 

53 

54 

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

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

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

58 

59 # read charge 

60 try: 

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

62 except ValueError: 

63 q = None 

64 else: 

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

66 

67 blocks = [] 

68 i1 = 0 

69 for i2, line in enumerate(lines): 

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

71 if i1 > 0: 

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

73 i1 = i2 

74 blocks.append(lines[i1:]) 

75 

76 images: List[Atoms] = [] 

77 for lines in blocks: 

78 try: 

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

80 except ValueError: 

81 pass 

82 else: 

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

84 del lines[i + 2] # old format 

85 cell: List[Union[float, List[float]]] = [] 

86 pbc = [] 

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

88 words = line.split() 

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

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

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

92 else: # new format with GUC 

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

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

95 

96 symbols = [] 

97 positions = [] 

98 magmoms = [] 

99 for line in lines[1:]: 

100 words = line.split() 

101 if len(words) < 5: 

102 break 

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

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

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

106 if len(words) > 5: 

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

108 magmoms.append(mom) 

109 if len(symbols): 

110 atoms = Atoms(symbols=symbols, positions=positions, 

111 cell=cell, pbc=pbc) 

112 else: 

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

114 

115 try: 

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

117 word = lines[ii].split() 

118 kx = int(word[2]) 

119 ky = int(word[4]) 

120 kz = int(word[6]) 

121 bz_kpts = (kx, ky, kz) 

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

123 except (ValueError, TypeError, IndexError): 

124 bz_kpts = None 

125 ibz_kpts = None 

126 

127 try: 

128 i = index_startswith(lines, 'energy contributions relative to') 

129 except ValueError: 

130 e = energy_contributions = None 

131 else: 

132 energy_contributions = {} 

133 for line in lines[i + 2:i + 13]: 

134 fields = line.split(':') 

135 if len(fields) == 2: 

136 name = fields[0] 

137 energy = float(fields[1]) 

138 energy_contributions[name] = energy 

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

140 e = energy 

141 break 

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]