Coverage for /builds/kinetik161/ase/ase/io/gpaw_out.py: 84.47%

1import re 

2from typing import List, Tuple, Union 

3 

4import numpy as np 

5 

6from ase.atoms import Atoms 

7from ase.calculators.singlepoint import (SinglePointDFTCalculator, 

8 SinglePointKPoint) 

9 

10 

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

12 for i, line in enumerate(lines): 

13 if line.startswith(string): 

14 return i 

15 raise ValueError 

16 

17 

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

19 repat = re.compile(pattern) 

20 for i, line in enumerate(lines): 

21 if repat.match(line): 

22 return i 

23 raise ValueError 

24 

25 

26def read_forces(lines: List[str], 

27 ii: int, 

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

29 f = [] 

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

31 try: 

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

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

34 except (ValueError, IndexError) as m: 

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

36 return f, i 

37 

38 

39def read_stresses(lines: List[str], 

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

41 s = [] 

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

43 try: 

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

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

46 except (ValueError, IndexError) as m: 

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

48 return s, i 

49 

50 

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

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

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

54 

55 blocks = [] 

56 i1 = 0 

57 for i2, line in enumerate(lines): 

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

59 if i1 > 0: 

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

61 i1 = i2 

62 blocks.append(lines[i1:]) 

63 

64 images: List[Atoms] = [] 

65 for lines in blocks: 

66 try: 

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

68 except ValueError: 

69 pass 

70 else: 

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

72 del lines[i + 2] # old format 

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

74 pbc = [] 

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

76 words = line.split() 

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

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

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

80 else: # new format with GUC 

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

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

83 

84 symbols = [] 

85 positions = [] 

86 magmoms = [] 

87 for line in lines[1:]: 

88 words = line.split() 

89 if len(words) < 5: 

90 break 

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

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

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

94 if len(words) > 5: 

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

96 magmoms.append(mom) 

97 if len(symbols): 

98 atoms = Atoms(symbols=symbols, positions=positions, 

99 cell=cell, pbc=pbc) 

100 else: 

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

102 

103 try: 

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

105 word = lines[ii].split() 

106 kx = int(word[2]) 

107 ky = int(word[4]) 

108 kz = int(word[6]) 

109 bz_kpts = (kx, ky, kz) 

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

111 except (ValueError, TypeError, IndexError): 

112 bz_kpts = None 

113 ibz_kpts = None 

114 

115 try: 

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

117 except ValueError: 

118 e = energy_contributions = None 

119 else: 

120 energy_contributions = {} 

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

122 fields = line.split(':') 

123 if len(fields) == 2: 

124 name = fields[0] 

125 energy = float(fields[1]) 

126 energy_contributions[name] = energy 

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

128 e = energy 

129 break 

130 else: # no break 

131 raise ValueError 

132 

133 try: 

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

135 except ValueError: 

136 eFermi = None 

137 else: 

138 fields = lines[ii].split() 

139 try: 

140 def strip(string): 

141 for rubbish in '[],': 

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

143 return string 

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

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

146 except ValueError: 

147 eFermi = float(fields[-1]) 

148 

149 # read Eigenvalues and occupations 

150 ii1 = ii2 = 1e32 

151 try: 

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

153 except ValueError: 

154 pass 

155 try: 

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

157 except ValueError: 

158 pass 

159 ii = min(ii1, ii2) 

160 if ii == 1e32: 

161 kpts = None 

162 else: 

163 ii += 1 

164 words = lines[ii].split() 

165 vals = [] 

166 while len(words) > 2: 

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

168 ii += 1 

169 words = lines[ii].split() 

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

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

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

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

174 if vals.shape[0] > 3: 

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

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

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

178 # read charge 

179 try: 

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

181 except ValueError: 

182 q = None 

183 else: 

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

185 # read dipole moment 

186 try: 

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

188 except ValueError: 

189 dipole = None 

190 else: 

191 line = lines[ii] 

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

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

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

195 

196 try: 

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

198 except ValueError: 

199 pass 

200 else: 

201 magmoms = [] 

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

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

204 magmoms.append(float(magmom)) 

205 

206 try: 

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

208 except ValueError: 

209 f = None 

210 else: 

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

212 

213 try: 

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

215 except ValueError: 

216 stress_tensor = None 

217 else: 

218 stress_tensor, i = read_stresses(lines, ii) 

219 

220 try: 

221 parameters = {} 

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

223 except ValueError: 

224 name = 'gpaw' 

225 else: 

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

227 # save uncorrected values 

228 parameters.update({ 

229 'calculator': 'gpaw', 

230 'uncorrected_energy': e, 

231 }) 

232 line = lines[ii + 1] 

233 assert line.startswith('energy:') 

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

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

236 

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

238 break 

239 

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

241 n = len(atoms) 

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

243 

244 if magmoms: 

245 atoms.set_initial_magnetic_moments(magmoms) 

246 else: 

247 magmoms = None 

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

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

250 dipole=dipole, magmoms=magmoms, 

251 efermi=eFermi, 

252 bzkpts=bz_kpts, ibzkpts=ibz_kpts, 

253 stress=stress_tensor) 

254 calc.name = name 

255 calc.parameters = parameters 

256 if energy_contributions is not None: 

257 calc.energy_contributions = energy_contributions 

258 if kpts is not None: 

259 calc.kpts = kpts 

260 atoms.calc = calc 

261 

262 images.append(atoms) 

263 

264 if len(images) == 0: 

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

266 

267 return images[index]