1import re 

2from typing import List, Tuple, Union 

3 

4import numpy as np 

5from ase.atoms import Atoms 

6from ase.calculators.singlepoint import (SinglePointDFTCalculator, 

7 SinglePointKPoint) 

8 

9 

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

11 for i, line in enumerate(lines): 

12 if line.startswith(string): 

13 return i 

14 raise ValueError 

15 

16 

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

18 repat = re.compile(pattern) 

19 for i, line in enumerate(lines): 

20 if repat.match(line): 

21 return i 

22 raise ValueError 

23 

24 

25def read_forces(lines: List[str], 

26 ii: int, 

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

28 f = [] 

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

30 try: 

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

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

33 except (ValueError, IndexError) as m: 

34 raise IOError('Malformed GPAW log file: %s' % m) 

35 return f, i 

36 

37 

38def read_stresses(lines: List[str], 

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

40 s = [] 

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

42 try: 

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

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

45 except (ValueError, IndexError) as m: 

46 raise IOError(f'Malformed GPAW log file: {m}') 

47 return s, i 

48 

49 

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

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

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

53 

54 blocks = [] 

55 i1 = 0 

56 for i2, line in enumerate(lines): 

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

58 if i1 > 0: 

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

60 i1 = i2 

61 blocks.append(lines[i1:]) 

62 

63 images: List[Atoms] = [] 

64 for lines in blocks: 

65 try: 

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

67 except ValueError: 

68 pass 

69 else: 

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

71 del lines[i + 2] # old format 

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

73 pbc = [] 

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

75 words = line.split() 

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

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

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

79 else: # new format with GUC 

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

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

82 

83 symbols = [] 

84 positions = [] 

85 magmoms = [] 

86 for line in lines[1:]: 

87 words = line.split() 

88 if len(words) < 5: 

89 break 

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

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

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

93 if len(words) > 5: 

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

95 magmoms.append(mom) 

96 if len(symbols): 

97 atoms = Atoms(symbols=symbols, positions=positions, 

98 cell=cell, pbc=pbc) 

99 else: 

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

101 

102 try: 

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

104 word = lines[ii].split() 

105 kx = int(word[2]) 

106 ky = int(word[4]) 

107 kz = int(word[6]) 

108 bz_kpts = (kx, ky, kz) 

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

110 except (ValueError, TypeError, IndexError): 

111 bz_kpts = None 

112 ibz_kpts = None 

113 

114 try: 

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

116 except ValueError: 

117 e = energy_contributions = None 

118 else: 

119 energy_contributions = {} 

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

121 fields = line.split(':') 

122 if len(fields) == 2: 

123 name = fields[0] 

124 energy = float(fields[1]) 

125 energy_contributions[name] = energy 

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

127 e = energy 

128 break 

129 else: # no break 

130 raise ValueError 

131 

132 try: 

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

134 except ValueError: 

135 eFermi = None 

136 else: 

137 fields = lines[ii].split() 

138 try: 

139 def strip(string): 

140 for rubbish in '[],': 

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

142 return string 

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

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

145 except ValueError: 

146 eFermi = float(fields[-1]) 

147 

148 # read Eigenvalues and occupations 

149 ii1 = ii2 = 1e32 

150 try: 

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

152 except ValueError: 

153 pass 

154 try: 

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

156 except ValueError: 

157 pass 

158 ii = min(ii1, ii2) 

159 if ii == 1e32: 

160 kpts = None 

161 else: 

162 ii += 1 

163 words = lines[ii].split() 

164 vals = [] 

165 while len(words) > 2: 

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

167 ii += 1 

168 words = lines[ii].split() 

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

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

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

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

173 if vals.shape[0] > 3: 

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

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

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

177 # read charge 

178 try: 

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

180 except ValueError: 

181 q = None 

182 else: 

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

184 # read dipole moment 

185 try: 

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

187 except ValueError: 

188 dipole = None 

189 else: 

190 line = lines[ii] 

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

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

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

194 

195 try: 

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

197 except ValueError: 

198 pass 

199 else: 

200 magmoms = [] 

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

202 magmom = lines[j].split()[-1].rstrip(')') 

203 magmoms.append(float(magmom)) 

204 

205 try: 

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

207 except ValueError: 

208 f = None 

209 else: 

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

211 

212 try: 

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

214 except ValueError: 

215 stress_tensor = None 

216 else: 

217 stress_tensor, i = read_stresses(lines, ii) 

218 

219 try: 

220 parameters = {} 

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

222 except ValueError: 

223 name = 'gpaw' 

224 else: 

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

226 # save uncorrected values 

227 parameters.update({ 

228 'calculator': 'gpaw', 

229 'uncorrected_energy': e, 

230 }) 

231 line = lines[ii + 1] 

232 assert line.startswith('energy:') 

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

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

235 

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

237 break 

238 

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

240 n = len(atoms) 

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

242 

243 if magmoms: 

244 atoms.set_initial_magnetic_moments(magmoms) 

245 else: 

246 magmoms = None 

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

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

249 dipole=dipole, magmoms=magmoms, 

250 efermi=eFermi, 

251 bzkpts=bz_kpts, ibzkpts=ibz_kpts, 

252 stress=stress_tensor) 

253 calc.name = name 

254 calc.parameters = parameters 

255 if energy_contributions is not None: 

256 calc.energy_contributions = energy_contributions 

257 if kpts is not None: 

258 calc.kpts = kpts 

259 atoms.calc = calc 

260 

261 images.append(atoms) 

262 

263 if len(images) == 0: 

264 raise IOError('Corrupted GPAW-text file!') 

265 

266 return images[index]