Coverage for /builds/debichem-team/python-ase/ase/io/gpaw_out.py: 83.25%

1import re 

2from typing import List, Tuple, Union 

3 

4import numpy as np 

5 

6from ase.atoms import Atoms 

7from ase.calculators.singlepoint import ( 

8 SinglePointDFTCalculator, 

9 SinglePointKPoint, 

10) 

11 

12 

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

14 for i, line in enumerate(lines): 

15 if line.startswith(string): 

16 return i 

17 raise ValueError 

18 

19 

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

21 repat = re.compile(pattern) 

22 for i, line in enumerate(lines): 

23 if repat.match(line): 

24 return i 

25 raise ValueError 

26 

27 

28def read_forces(lines: List[str], 

29 ii: int, 

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

31 f = [] 

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

33 try: 

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

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

36 except (ValueError, IndexError) as m: 

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

38 return f, i 

39 

40 

41def read_stresses(lines: List[str], 

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

43 s = [] 

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

45 try: 

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

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

48 except (ValueError, IndexError) as m: 

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

50 return s, i 

51 

52 

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

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

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

56 

57 blocks = [] 

58 i1 = 0 

59 for i2, line in enumerate(lines): 

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

61 if i1 > 0: 

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

63 i1 = i2 

64 blocks.append(lines[i1:]) 

65 

66 images: List[Atoms] = [] 

67 for lines in blocks: 

68 try: 

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

70 except ValueError: 

71 pass 

72 else: 

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

74 del lines[i + 2] # old format 

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

76 pbc = [] 

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

78 words = line.split() 

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

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

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

82 else: # new format with GUC 

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

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

85 

86 symbols = [] 

87 positions = [] 

88 magmoms = [] 

89 for line in lines[1:]: 

90 words = line.split() 

91 if len(words) < 5: 

92 break 

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

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

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

96 if len(words) > 5: 

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

98 magmoms.append(mom) 

99 if len(symbols): 

100 atoms = Atoms(symbols=symbols, positions=positions, 

101 cell=cell, pbc=pbc) 

102 else: 

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

104 

105 try: 

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

107 word = lines[ii].split() 

108 kx = int(word[2]) 

109 ky = int(word[4]) 

110 kz = int(word[6]) 

111 bz_kpts = (kx, ky, kz) 

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

113 except (ValueError, TypeError, IndexError): 

114 bz_kpts = None 

115 ibz_kpts = None 

116 

117 try: 

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

119 except ValueError: 

120 e = energy_contributions = None 

121 else: 

122 energy_contributions = {} 

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

124 fields = line.split(':') 

125 if len(fields) == 2: 

126 name = fields[0] 

127 energy = float(fields[1]) 

128 energy_contributions[name] = energy 

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

130 e = energy 

131 break 

132 else: # no break 

133 raise ValueError 

134 

135 try: 

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

137 except ValueError: 

138 eFermi = None 

139 else: 

140 fields = lines[ii].split() 

141 try: 

142 def strip(string): 

143 for rubbish in '[],': 

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

145 return string 

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

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

148 except ValueError: 

149 eFermi = float(fields[-1]) 

150 

151 # read Eigenvalues and occupations 

152 ii1 = ii2 = 1e32 

153 try: 

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

155 except ValueError: 

156 pass 

157 try: 

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

159 except ValueError: 

160 pass 

161 ii = min(ii1, ii2) 

162 if ii == 1e32: 

163 kpts = None 

164 else: 

165 ii += 1 

166 words = lines[ii].split() 

167 vals = [] 

168 while len(words) > 2: 

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

170 ii += 1 

171 words = lines[ii].split() 

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

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

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

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

176 if vals.shape[0] > 3: 

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

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

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

180 # read charge 

181 try: 

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

183 except ValueError: 

184 q = None 

185 else: 

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

187 # read dipole moment 

188 try: 

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

190 except ValueError: 

191 dipole = None 

192 else: 

193 line = lines[ii] 

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

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

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

197 

198 try: 

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

200 except ValueError: 

201 pass 

202 else: 

203 magmoms = [] 

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

205 line = lines[j] 

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

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

208 else: 

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

210 magmoms.append(float(magmom)) 

211 

212 try: 

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

214 except ValueError: 

215 f = None 

216 else: 

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

218 

219 try: 

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

221 except ValueError: 

222 stress_tensor = None 

223 else: 

224 stress_tensor, i = read_stresses(lines, ii) 

225 

226 try: 

227 parameters = {} 

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

229 except ValueError: 

230 name = 'gpaw' 

231 else: 

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

233 # save uncorrected values 

234 parameters.update({ 

235 'calculator': 'gpaw', 

236 'uncorrected_energy': e, 

237 }) 

238 line = lines[ii + 1] 

239 assert line.startswith('energy:') 

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

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

242 

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

244 break 

245 

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

247 n = len(atoms) 

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

249 

250 if magmoms: 

251 atoms.set_initial_magnetic_moments(magmoms) 

252 else: 

253 magmoms = None 

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

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

256 dipole=dipole, magmoms=magmoms, 

257 efermi=eFermi, 

258 bzkpts=bz_kpts, ibzkpts=ibz_kpts, 

259 stress=stress_tensor) 

260 calc.name = name 

261 calc.parameters = parameters 

262 if energy_contributions is not None: 

263 calc.energy_contributions = energy_contributions 

264 if kpts is not None: 

265 calc.kpts = kpts 

266 atoms.calc = calc 

267 

268 images.append(atoms) 

269 

270 if len(images) == 0: 

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

272 

273 return images[index]