Coverage for /builds/ase/ase/ase/io/elk.py : 86.33%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import collections
2from pathlib import Path
4import numpy as np
6from ase import Atoms
7from ase.units import Bohr, Hartree
8from ase.utils import reader, writer
11elk_parameters = {'swidth': Hartree}
14@reader
15def read_elk(fd):
16 """Import ELK atoms definition.
18 Reads unitcell, atom positions, magmoms from elk.in/GEOMETRY.OUT file.
19 """
21 lines = fd.readlines()
23 scale = np.ones(4) # unit cell scale
24 positions = []
25 cell = []
26 symbols = []
27 magmoms = []
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 line in lines1[n1 + 2:n1 + 2 + natomsnow]:
62 atposnow.append([float(v) for v in line.split()[0:3]])
63 if len(line.split()) == 6: # bfcmt present
64 bfcmtnow.append(
65 [float(v) for v in line.split()[3:]])
66 atpos.append(atposnow)
67 bfcmt.append(bfcmtnow)
68 # symbols, positions, magmoms based on ELK spfname, atpos, and bfcmt
69 symbols = ''
70 positions = []
71 magmoms = []
72 for n, s in enumerate(spfname):
73 symbols += str(s[1:].split('.')[0]) * natoms[n]
74 positions += atpos[n] # assumes fractional coordinates
75 if len(bfcmt[n]) > 0:
76 # how to handle cases of magmoms being one or three dim array?
77 magmoms += [m[-1] for m in bfcmt[n]]
78 atoms = Atoms(symbols, scaled_positions=positions, cell=[1, 1, 1])
79 if len(magmoms) > 0:
80 atoms.set_initial_magnetic_moments(magmoms)
81 # final cell scale
82 cell = cell * scale[0] * Bohr
84 atoms.set_cell(cell, scale_atoms=True)
85 atoms.pbc = True
86 return atoms
89@writer
90def write_elk_in(fd, atoms, parameters=None):
91 if parameters is None:
92 parameters = {}
94 parameters = dict(parameters)
95 species_path = parameters.pop('species_dir', None)
97 if parameters.get('spinpol') is None:
98 if atoms.get_initial_magnetic_moments().any():
99 parameters['spinpol'] = True
101 if 'xctype' in parameters:
102 if 'xc' in parameters:
103 raise RuntimeError("You can't use both 'xctype' and 'xc'!")
105 if parameters.get('autokpt'):
106 if 'kpts' in parameters:
107 raise RuntimeError("You can't use both 'autokpt' and 'kpts'!")
108 if 'ngridk' in parameters:
109 raise RuntimeError(
110 "You can't use both 'autokpt' and 'ngridk'!")
111 if 'ngridk' in parameters:
112 if 'kpts' in parameters:
113 raise RuntimeError("You can't use both 'ngridk' and 'kpts'!")
115 if parameters.get('autoswidth'):
116 if 'smearing' in parameters:
117 raise RuntimeError(
118 "You can't use both 'autoswidth' and 'smearing'!")
119 if 'swidth' in parameters:
120 raise RuntimeError(
121 "You can't use both 'autoswidth' and 'swidth'!")
123 inp = {}
124 inp.update(parameters)
126 if 'xc' in parameters:
127 xctype = {'LDA': 3, # PW92
128 'PBE': 20,
129 'REVPBE': 21,
130 'PBESOL': 22,
131 'WC06': 26,
132 'AM05': 30,
133 'mBJLDA': (100, 208, 12)}[parameters['xc']]
134 inp['xctype'] = xctype
135 del inp['xc']
137 if 'kpts' in parameters:
138 # XXX should generalize kpts handling.
139 from ase.calculators.calculator import kpts2mp
140 mp = kpts2mp(atoms, parameters['kpts'])
141 inp['ngridk'] = tuple(mp)
142 vkloff = [] # is this below correct?
143 for nk in mp:
144 if nk % 2 == 0: # shift kpoint away from gamma point
145 vkloff.append(0.5)
146 else:
147 vkloff.append(0)
148 inp['vkloff'] = vkloff
149 del inp['kpts']
151 if 'smearing' in parameters:
152 name = parameters.smearing[0].lower()
153 if name == 'methfessel-paxton':
154 stype = parameters.smearing[2]
155 else:
156 stype = {'gaussian': 0,
157 'fermi-dirac': 3,
158 }[name]
159 inp['stype'] = stype
160 inp['swidth'] = parameters.smearing[1]
161 del inp['smearing']
163 # convert keys to ELK units
164 for key, value in inp.items():
165 if key in elk_parameters:
166 inp[key] /= elk_parameters[key]
168 # write all keys
169 for key, value in inp.items():
170 fd.write('%s\n' % key)
171 if isinstance(value, bool):
172 fd.write('.%s.\n\n' % ('false', 'true')[value])
173 elif isinstance(value, (int, float)):
174 fd.write('%s\n\n' % value)
175 else:
176 fd.write('%s\n\n' % ' '.join([str(x) for x in value]))
178 # cell
179 fd.write('avec\n')
180 for vec in atoms.cell:
181 fd.write('%.14f %.14f %.14f\n' % tuple(vec / Bohr))
182 fd.write('\n')
184 # atoms
185 species = {}
186 symbols = []
187 for a, (symbol, m) in enumerate(
188 zip(atoms.get_chemical_symbols(),
189 atoms.get_initial_magnetic_moments())):
190 if symbol in species:
191 species[symbol].append((a, m))
192 else:
193 species[symbol] = [(a, m)]
194 symbols.append(symbol)
195 fd.write('atoms\n%d\n' % len(species))
196 # scaled = atoms.get_scaled_positions(wrap=False)
197 scaled = np.linalg.solve(atoms.cell.T, atoms.positions.T).T
198 for symbol in symbols:
199 fd.write("'%s.in' : spfname\n" % symbol)
200 fd.write('%d\n' % len(species[symbol]))
201 for a, m in species[symbol]:
202 fd.write('%.14f %.14f %.14f 0.0 0.0 %.14f\n' %
203 (tuple(scaled[a]) + (m,)))
205 # if sppath is present in elk.in it overwrites species blocks!
207 # Elk seems to concatenate path and filename in such a way
208 # that we must put a / at the end:
209 if species_path is not None:
210 fd.write(f"sppath\n'{species_path}/'\n\n")
213class ElkReader:
214 def __init__(self, path):
215 self.path = Path(path)
217 def _read_everything(self):
218 yield from self._read_energy()
220 with (self.path / 'INFO.OUT').open() as fd:
221 yield from parse_elk_info(fd)
223 with (self.path / 'EIGVAL.OUT').open() as fd:
224 yield from parse_elk_eigval(fd)
226 with (self.path / 'KPOINTS.OUT').open() as fd:
227 yield from parse_elk_kpoints(fd)
229 def read_everything(self):
230 dct = dict(self._read_everything())
232 # The eigenvalue/occupation tables do not say whether there are
233 # two spins, so we have to reshape them from 1 x K x SB to S x K x B:
234 spinpol = dct.pop('spinpol')
235 if spinpol:
236 for name in 'eigenvalues', 'occupations':
237 array = dct[name]
238 _, nkpts, nbands_double = array.shape
239 assert _ == 1
240 assert nbands_double % 2 == 0
241 nbands = nbands_double // 2
242 newarray = np.empty((2, nkpts, nbands))
243 newarray[0, :, :] = array[0, :, :nbands]
244 newarray[1, :, :] = array[0, :, nbands:]
245 if name == 'eigenvalues':
246 # Verify that eigenvalues are still sorted:
247 diffs = np.diff(newarray, axis=2)
248 assert all(diffs.flat[:] > 0)
249 dct[name] = newarray
250 return dct
252 def _read_energy(self):
253 txt = (self.path / 'TOTENERGY.OUT').read_text()
254 tokens = txt.split()
255 energy = float(tokens[-1]) * Hartree
256 yield 'free_energy', energy
257 yield 'energy', energy
260def parse_elk_kpoints(fd):
261 header = next(fd)
262 lhs, rhs = header.split(':', 1)
263 assert 'k-point, vkl, wkpt' in rhs, header
264 nkpts = int(lhs.strip())
266 kpts = np.empty((nkpts, 3))
267 weights = np.empty(nkpts)
269 for ikpt in range(nkpts):
270 line = next(fd)
271 tokens = line.split()
272 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
273 weights[ikpt] = float(tokens[4])
274 yield 'ibz_kpoints', kpts
275 yield 'kpoint_weights', weights
278def parse_elk_info(fd):
279 dct = collections.defaultdict(list)
280 fd = iter(fd)
282 spinpol = None
283 converged = False
284 actually_did_not_converge = False
285 # Legacy code kept track of both these things, which is strange.
286 # Why could a file both claim to converge *and* not converge?
288 # We loop over all lines and extract also data that occurs
289 # multiple times (e.g. in multiple SCF steps)
290 for line in fd:
291 # "name of quantity : 1 2 3"
292 tokens = line.split(':', 1)
293 if len(tokens) == 2:
294 lhs, rhs = tokens
295 dct[lhs.strip()].append(rhs.strip())
297 elif line.startswith('Convergence targets achieved'):
298 converged = True
299 elif 'reached self-consistent loops maximum' in line.lower():
300 actually_did_not_converge = True
302 if 'Spin treatment' in line:
303 # (Somewhat brittle doing multi-line stuff here.)
304 line = next(fd)
305 spinpol = line.strip() == 'spin-polarised'
307 yield 'converged', converged and not actually_did_not_converge
308 if spinpol is None:
309 raise RuntimeError('Could not determine spin treatment')
310 yield 'spinpol', spinpol
312 if 'Fermi' in dct:
313 yield 'fermi_level', float(dct['Fermi'][-1]) * Hartree
315 if 'total force' in dct:
316 forces = []
317 for line in dct['total force']:
318 forces.append(line.split())
319 yield 'forces', np.array(forces, float) * (Hartree / Bohr)
322def parse_elk_eigval(fd):
324 def match_int(line, word):
325 number, colon, word1 = line.split()
326 assert word1 == word
327 assert colon == ':'
328 return int(number)
330 def skip_spaces(line=''):
331 while not line.strip():
332 line = next(fd)
333 return line
335 line = skip_spaces()
336 nkpts = match_int(line, 'nkpt') # 10 : nkpts
337 line = next(fd)
338 nbands = match_int(line, 'nstsv') # 15 : nstsv
340 eigenvalues = np.empty((nkpts, nbands))
341 occupations = np.empty((nkpts, nbands))
342 kpts = np.empty((nkpts, 3))
344 for ikpt in range(nkpts):
345 line = skip_spaces()
346 tokens = line.split()
347 assert tokens[-1] == 'vkl', tokens
348 assert ikpt + 1 == int(tokens[0])
349 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
351 line = next(fd) # "(state, eigenvalue and occupancy below)"
352 assert line.strip().startswith('(state,'), line
353 for iband in range(nbands):
354 line = next(fd)
355 tokens = line.split() # (band number, eigenval, occ)
356 assert iband + 1 == int(tokens[0])
357 eigenvalues[ikpt, iband] = float(tokens[1])
358 occupations[ikpt, iband] = float(tokens[2])
360 yield 'ibz_kpoints', kpts
361 yield 'eigenvalues', eigenvalues[None] * Hartree
362 yield 'occupations', occupations[None]