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_gpaw_out(fileobj, index): # -> Union[Atoms, List[Atoms]]: 

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

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

41 

42 blocks = [] 

43 i1 = 0 

44 for i2, line in enumerate(lines): 

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

46 if i1 > 0: 

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

48 i1 = i2 

49 blocks.append(lines[i1:]) 

50 

51 images: List[Atoms] = [] 

52 for lines in blocks: 

53 try: 

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

55 except ValueError: 

56 pass 

57 else: 

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

59 del lines[i + 2] # old format 

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

61 pbc = [] 

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

63 words = line.split() 

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

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

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

67 else: # new format with GUC 

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

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

70 

71 symbols = [] 

72 positions = [] 

73 magmoms = [] 

74 for line in lines[1:]: 

75 words = line.split() 

76 if len(words) < 5: 

77 break 

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

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

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

81 if len(words) > 5: 

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

83 magmoms.append(mom) 

84 if len(symbols): 

85 atoms = Atoms(symbols=symbols, positions=positions, 

86 cell=cell, pbc=pbc) 

87 else: 

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

89 

90 try: 

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

92 word = lines[ii].split() 

93 kx = int(word[2]) 

94 ky = int(word[4]) 

95 kz = int(word[6]) 

96 bz_kpts = (kx, ky, kz) 

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

98 except (ValueError, TypeError, IndexError): 

99 bz_kpts = None 

100 ibz_kpts = None 

101 

102 try: 

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

104 except ValueError: 

105 e = energy_contributions = None 

106 else: 

107 energy_contributions = {} 

108 for line in lines[i + 2:i + 8]: 

109 fields = line.split(':') 

110 energy_contributions[fields[0]] = float(fields[1]) 

111 line = lines[i + 10] 

112 assert (line.startswith('zero kelvin:') or 

113 line.startswith('extrapolated:')) 

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

115 

116 try: 

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

118 except ValueError: 

119 eFermi = None 

120 else: 

121 fields = lines[ii].split() 

122 try: 

123 def strip(string): 

124 for rubbish in '[],': 

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

126 return string 

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

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

129 except ValueError: 

130 eFermi = float(fields[-1]) 

131 

132 # read Eigenvalues and occupations 

133 ii1 = ii2 = 1e32 

134 try: 

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

136 except ValueError: 

137 pass 

138 try: 

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

140 except ValueError: 

141 pass 

142 ii = min(ii1, ii2) 

143 if ii == 1e32: 

144 kpts = None 

145 else: 

146 ii += 1 

147 words = lines[ii].split() 

148 vals = [] 

149 while(len(words) > 2): 

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

151 ii += 1 

152 words = lines[ii].split() 

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

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

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

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

157 if vals.shape[0] > 3: 

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

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

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

161 # read charge 

162 try: 

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

164 except ValueError: 

165 q = None 

166 else: 

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

168 # read dipole moment 

169 try: 

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

171 except ValueError: 

172 dipole = None 

173 else: 

174 line = lines[ii] 

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

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

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

178 

179 try: 

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

181 except ValueError: 

182 pass 

183 else: 

184 magmoms = [] 

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

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

187 magmoms.append(float(magmom)) 

188 

189 try: 

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

191 except ValueError: 

192 f = None 

193 else: 

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

195 

196 try: 

197 parameters = {} 

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

199 except ValueError: 

200 name = 'gpaw' 

201 else: 

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

203 # save uncorrected values 

204 parameters.update({ 

205 'calculator': 'gpaw', 

206 'uncorrected_energy': e, 

207 }) 

208 line = lines[ii + 1] 

209 assert line.startswith('energy:') 

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

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

212 

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

214 break 

215 

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

217 n = len(atoms) 

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

219 

220 if magmoms: 

221 atoms.set_initial_magnetic_moments(magmoms) 

222 else: 

223 magmoms = None 

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

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

226 dipole=dipole, magmoms=magmoms, 

227 efermi=eFermi, 

228 bzkpts=bz_kpts, ibzkpts=ibz_kpts) 

229 calc.name = name 

230 calc.parameters = parameters 

231 if energy_contributions is not None: 

232 calc.energy_contributions = energy_contributions 

233 if kpts is not None: 

234 calc.kpts = kpts 

235 atoms.calc = calc 

236 

237 images.append(atoms) 

238 

239 if len(images) == 0: 

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

241 

242 return images[index]