1import collections 

2from pathlib import Path 

3 

4import numpy as np 

5 

6from ase import Atoms 

7from ase.units import Bohr, Hartree 

8from ase.utils import reader, writer 

9 

10 

11elk_parameters = {'swidth': Hartree} 

12 

13 

14@reader 

15def read_elk(fd): 

16 """Import ELK atoms definition. 

17 

18 Reads unitcell, atom positions, magmoms from elk.in/GEOMETRY.OUT file. 

19 """ 

20 

21 lines = fd.readlines() 

22 

23 scale = np.ones(4) # unit cell scale 

24 positions = [] 

25 cell = [] 

26 symbols = [] 

27 magmoms = [] 

28 

29 # find cell scale 

30 for n, line in enumerate(lines): 

31 if line.split() == []: 

32 continue 

33 if line.strip() == 'scale': 

34 scale[0] = float(lines[n + 1]) 

35 elif line.startswith('scale'): 

36 scale[int(line.strip()[-1])] = float(lines[n + 1]) 

37 for n, line in enumerate(lines): 

38 if line.split() == []: 

39 continue 

40 if line.startswith('avec'): 

41 cell = np.array( 

42 [[float(v) * scale[1] for v in lines[n + 1].split()], 

43 [float(v) * scale[2] for v in lines[n + 2].split()], 

44 [float(v) * scale[3] for v in lines[n + 3].split()]]) 

45 if line.startswith('atoms'): 

46 lines1 = lines[n + 1:] # start subsearch 

47 spfname = [] 

48 natoms = [] 

49 atpos = [] 

50 bfcmt = [] 

51 for n1, line1 in enumerate(lines1): 

52 if line1.split() == []: 

53 continue 

54 if 'spfname' in line1: 

55 spfnamenow = lines1[n1].split()[0] 

56 spfname.append(spfnamenow) 

57 natomsnow = int(lines1[n1 + 1].split()[0]) 

58 natoms.append(natomsnow) 

59 atposnow = [] 

60 bfcmtnow = [] 

61 for l in lines1[n1 + 2:n1 + 2 + natomsnow]: 

62 atposnow.append([float(v) for v in l.split()[0:3]]) 

63 if len(l.split()) == 6: # bfcmt present 

64 bfcmtnow.append([float(v) for v in l.split()[3:]]) 

65 atpos.append(atposnow) 

66 bfcmt.append(bfcmtnow) 

67 # symbols, positions, magmoms based on ELK spfname, atpos, and bfcmt 

68 symbols = '' 

69 positions = [] 

70 magmoms = [] 

71 for n, s in enumerate(spfname): 

72 symbols += str(s[1:].split('.')[0]) * natoms[n] 

73 positions += atpos[n] # assumes fractional coordinates 

74 if len(bfcmt[n]) > 0: 

75 # how to handle cases of magmoms being one or three dim array? 

76 magmoms += [m[-1] for m in bfcmt[n]] 

77 atoms = Atoms(symbols, scaled_positions=positions, cell=[1, 1, 1]) 

78 if len(magmoms) > 0: 

79 atoms.set_initial_magnetic_moments(magmoms) 

80 # final cell scale 

81 cell = cell * scale[0] * Bohr 

82 

83 atoms.set_cell(cell, scale_atoms=True) 

84 atoms.pbc = True 

85 return atoms 

86 

87 

88@writer 

89def write_elk_in(fd, atoms, parameters=None): 

90 if parameters is None: 

91 parameters = {} 

92 

93 parameters = dict(parameters) 

94 species_path = parameters.pop('species_dir', None) 

95 

96 if parameters.get('spinpol') is None: 

97 if atoms.get_initial_magnetic_moments().any(): 

98 parameters['spinpol'] = True 

99 

100 if 'xctype' in parameters: 

101 if 'xc' in parameters: 

102 raise RuntimeError("You can't use both 'xctype' and 'xc'!") 

103 

104 if parameters.get('autokpt'): 

105 if 'kpts' in parameters: 

106 raise RuntimeError("You can't use both 'autokpt' and 'kpts'!") 

107 if 'ngridk' in parameters: 

108 raise RuntimeError( 

109 "You can't use both 'autokpt' and 'ngridk'!") 

110 if 'ngridk' in parameters: 

111 if 'kpts' in parameters: 

112 raise RuntimeError("You can't use both 'ngridk' and 'kpts'!") 

113 

114 if parameters.get('autoswidth'): 

115 if 'smearing' in parameters: 

116 raise RuntimeError( 

117 "You can't use both 'autoswidth' and 'smearing'!") 

118 if 'swidth' in parameters: 

119 raise RuntimeError( 

120 "You can't use both 'autoswidth' and 'swidth'!") 

121 

122 inp = {} 

123 inp.update(parameters) 

124 

125 if 'xc' in parameters: 

126 xctype = {'LDA': 3, # PW92 

127 'PBE': 20, 

128 'REVPBE': 21, 

129 'PBESOL': 22, 

130 'WC06': 26, 

131 'AM05': 30, 

132 'mBJLDA': (100, 208, 12)}[parameters['xc']] 

133 inp['xctype'] = xctype 

134 del inp['xc'] 

135 

136 if 'kpts' in parameters: 

137 # XXX should generalize kpts handling. 

138 from ase.calculators.calculator import kpts2mp 

139 mp = kpts2mp(atoms, parameters['kpts']) 

140 inp['ngridk'] = tuple(mp) 

141 vkloff = [] # is this below correct? 

142 for nk in mp: 

143 if nk % 2 == 0: # shift kpoint away from gamma point 

144 vkloff.append(0.5) 

145 else: 

146 vkloff.append(0) 

147 inp['vkloff'] = vkloff 

148 del inp['kpts'] 

149 

150 if 'smearing' in parameters: 

151 name = parameters.smearing[0].lower() 

152 if name == 'methfessel-paxton': 

153 stype = parameters.smearing[2] 

154 else: 

155 stype = {'gaussian': 0, 

156 'fermi-dirac': 3, 

157 }[name] 

158 inp['stype'] = stype 

159 inp['swidth'] = parameters.smearing[1] 

160 del inp['smearing'] 

161 

162 # convert keys to ELK units 

163 for key, value in inp.items(): 

164 if key in elk_parameters: 

165 inp[key] /= elk_parameters[key] 

166 

167 # write all keys 

168 for key, value in inp.items(): 

169 fd.write('%s\n' % key) 

170 if isinstance(value, bool): 

171 fd.write('.%s.\n\n' % ('false', 'true')[value]) 

172 elif isinstance(value, (int, float)): 

173 fd.write('%s\n\n' % value) 

174 else: 

175 fd.write('%s\n\n' % ' '.join([str(x) for x in value])) 

176 

177 # cell 

178 fd.write('avec\n') 

179 for vec in atoms.cell: 

180 fd.write('%.14f %.14f %.14f\n' % tuple(vec / Bohr)) 

181 fd.write('\n') 

182 

183 # atoms 

184 species = {} 

185 symbols = [] 

186 for a, (symbol, m) in enumerate( 

187 zip(atoms.get_chemical_symbols(), 

188 atoms.get_initial_magnetic_moments())): 

189 if symbol in species: 

190 species[symbol].append((a, m)) 

191 else: 

192 species[symbol] = [(a, m)] 

193 symbols.append(symbol) 

194 fd.write('atoms\n%d\n' % len(species)) 

195 # scaled = atoms.get_scaled_positions(wrap=False) 

196 scaled = np.linalg.solve(atoms.cell.T, atoms.positions.T).T 

197 for symbol in symbols: 

198 fd.write("'%s.in' : spfname\n" % symbol) 

199 fd.write('%d\n' % len(species[symbol])) 

200 for a, m in species[symbol]: 

201 fd.write('%.14f %.14f %.14f 0.0 0.0 %.14f\n' % 

202 (tuple(scaled[a]) + (m,))) 

203 

204 # if sppath is present in elk.in it overwrites species blocks! 

205 

206 # Elk seems to concatenate path and filename in such a way 

207 # that we must put a / at the end: 

208 if species_path is not None: 

209 fd.write(f"sppath\n'{species_path}/'\n\n") 

210 

211 

212class ElkReader: 

213 def __init__(self, path): 

214 self.path = Path(path) 

215 

216 def _read_everything(self): 

217 yield from self._read_energy() 

218 

219 with (self.path / 'INFO.OUT').open() as fd: 

220 yield from parse_elk_info(fd) 

221 

222 with (self.path / 'EIGVAL.OUT').open() as fd: 

223 yield from parse_elk_eigval(fd) 

224 

225 with (self.path / 'KPOINTS.OUT').open() as fd: 

226 yield from parse_elk_kpoints(fd) 

227 

228 def read_everything(self): 

229 dct = dict(self._read_everything()) 

230 

231 # The eigenvalue/occupation tables do not say whether there are 

232 # two spins, so we have to reshape them from 1 x K x SB to S x K x B: 

233 spinpol = dct.pop('spinpol') 

234 if spinpol: 

235 for name in 'eigenvalues', 'occupations': 

236 array = dct[name] 

237 _, nkpts, nbands_double = array.shape 

238 assert _ == 1 

239 assert nbands_double % 2 == 0 

240 nbands = nbands_double // 2 

241 newarray = np.empty((2, nkpts, nbands)) 

242 newarray[0, :, :] = array[0, :, :nbands] 

243 newarray[1, :, :] = array[0, :, nbands:] 

244 if name == 'eigenvalues': 

245 # Verify that eigenvalues are still sorted: 

246 diffs = np.diff(newarray, axis=2) 

247 assert all(diffs.flat[:] > 0) 

248 dct[name] = newarray 

249 return dct 

250 

251 def _read_energy(self): 

252 txt = (self.path / 'TOTENERGY.OUT').read_text() 

253 tokens = txt.split() 

254 energy = float(tokens[-1]) * Hartree 

255 yield 'free_energy', energy 

256 yield 'energy', energy 

257 

258 

259def parse_elk_kpoints(fd): 

260 header = next(fd) 

261 lhs, rhs = header.split(':', 1) 

262 assert 'k-point, vkl, wkpt' in rhs, header 

263 nkpts = int(lhs.strip()) 

264 

265 kpts = np.empty((nkpts, 3)) 

266 weights = np.empty(nkpts) 

267 

268 for ikpt in range(nkpts): 

269 line = next(fd) 

270 tokens = line.split() 

271 kpts[ikpt] = np.array(tokens[1:4]).astype(float) 

272 weights[ikpt] = float(tokens[4]) 

273 yield 'ibz_kpoints', kpts 

274 yield 'kpoint_weights', weights 

275 

276 

277def parse_elk_info(fd): 

278 dct = collections.defaultdict(list) 

279 fd = iter(fd) 

280 

281 spinpol = None 

282 converged = False 

283 actually_did_not_converge = False 

284 # Legacy code kept track of both these things, which is strange. 

285 # Why could a file both claim to converge *and* not converge? 

286 

287 # We loop over all lines and extract also data that occurs 

288 # multiple times (e.g. in multiple SCF steps) 

289 for line in fd: 

290 # "name of quantity : 1 2 3" 

291 tokens = line.split(':', 1) 

292 if len(tokens) == 2: 

293 lhs, rhs = tokens 

294 dct[lhs.strip()].append(rhs.strip()) 

295 

296 elif line.startswith('Convergence targets achieved'): 

297 converged = True 

298 elif 'reached self-consistent loops maximum' in line.lower(): 

299 actually_did_not_converge = True 

300 

301 if 'Spin treatment' in line: 

302 # (Somewhat brittle doing multi-line stuff here.) 

303 line = next(fd) 

304 spinpol = line.strip() == 'spin-polarised' 

305 

306 yield 'converged', converged and not actually_did_not_converge 

307 if spinpol is None: 

308 raise RuntimeError('Could not determine spin treatment') 

309 yield 'spinpol', spinpol 

310 

311 if 'Fermi' in dct: 

312 yield 'fermi_level', float(dct['Fermi'][-1]) * Hartree 

313 

314 if 'total force' in dct: 

315 forces = [] 

316 for line in dct['total force']: 

317 forces.append(line.split()) 

318 yield 'forces', np.array(forces, float) * (Hartree / Bohr) 

319 

320 

321def parse_elk_eigval(fd): 

322 

323 def match_int(line, word): 

324 number, colon, word1 = line.split() 

325 assert word1 == word 

326 assert colon == ':' 

327 return int(number) 

328 

329 def skip_spaces(line=''): 

330 while not line.strip(): 

331 line = next(fd) 

332 return line 

333 

334 line = skip_spaces() 

335 nkpts = match_int(line, 'nkpt') # 10 : nkpts 

336 line = next(fd) 

337 nbands = match_int(line, 'nstsv') # 15 : nstsv 

338 

339 eigenvalues = np.empty((nkpts, nbands)) 

340 occupations = np.empty((nkpts, nbands)) 

341 kpts = np.empty((nkpts, 3)) 

342 

343 for ikpt in range(nkpts): 

344 line = skip_spaces() 

345 tokens = line.split() 

346 assert tokens[-1] == 'vkl', tokens 

347 assert ikpt + 1 == int(tokens[0]) 

348 kpts[ikpt] = np.array(tokens[1:4]).astype(float) 

349 

350 line = next(fd) # "(state, eigenvalue and occupancy below)" 

351 assert line.strip().startswith('(state,'), line 

352 for iband in range(nbands): 

353 line = next(fd) 

354 tokens = line.split() # (band number, eigenval, occ) 

355 assert iband + 1 == int(tokens[0]) 

356 eigenvalues[ikpt, iband] = float(tokens[1]) 

357 occupations[ikpt, iband] = float(tokens[2]) 

358 

359 yield 'ibz_kpoints', kpts 

360 yield 'eigenvalues', eigenvalues[None] * Hartree 

361 yield 'occupations', occupations[None]