Coverage for ase / io / elk.py: 86.36%
264 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3import collections
4import numbers
5from pathlib import Path
7import numpy as np
9from ase import Atoms
10from ase.units import Bohr, Hartree
11from ase.utils import reader, writer
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 """
92 Writes an ASE Atoms object and optional parameters to an `elk.in` file.
94 The format of `elk.in` and the meaning of the parameters is documented under
95 https://elk.sourceforge.io/. All numeric parameters that are passed using
96 Elk-native keys must be in Elk units.
98 Parameters
99 ----------
100 fd : path or file object
101 A file path or an opened, writable file or file-like object
102 atoms : Atoms object
103 An ASE Atoms object with the atomic structure
104 parameters : dict
105 The keys of the dict are the names of the input blocks. The dict values
106 contain the contents of the input blocks. The contents can be several
107 different types or data structures: 1) bool, str, int, float for scalar
108 parameters; 2) iterables for input blocks containing several lines,
109 one element per line; 3) iterables for several parameters in a line,
110 where arrays are treated as iterables; 4) iterables at the innermost
111 level describing array parameters
113 Raises
114 ------
115 RuntimeError
116 This exception is raised in case of inconsistencies between parameters.
117 TypeError
118 This exception is raised in case of a wrong parameter type.
120 Returns
121 -------
122 None
123 The function writes directly to the provided file object.
125 See Also
126 --------
127 ase.calculators.elk.ELK : The ASE Elk calculator
129 Examples
130 --------
131 >>> import numpy as np
132 >>> from ase import Atoms
133 >>> from ase.io.elk import write_elk_in
134 >>> atoms = Atoms('FeAl', positions=[[0, 0, 0], [1.45]*3], cell=[2.9]*3)
135 >>> params = {'tasks': [[0], [10]], 'ngridk': [8, 8, 8], 'nempty': 8,
136 'bfieldc': np.array((0.0, 0.0, -0.01)), 'spinpol': True,
137 'dft+u': ((2, 1), (1, 2, 0.183, 0.034911967))}
138 >>> write_elk_in('/path/to/elk.in', atoms, parameters=params)
139 Note: ``np.array((0.0, 0.0, -0.01))``, ``[0.0, 0.0, -0.01]``,
140 ``[[0.0, 0.0, -0.01]]``, and any combinations of lists, tuples, and
141 ndarrays, are equivalent
142 """
144 if parameters is None:
145 parameters = {}
147 parameters = dict(parameters)
149 species_path = parameters.pop('species_dir', None)
150 species_path = parameters.pop('sppath', species_path)
152 if parameters.get('spinpol') is None:
153 if atoms.get_initial_magnetic_moments().any():
154 parameters['spinpol'] = True
156 if 'xctype' in parameters:
157 if 'xc' in parameters:
158 raise RuntimeError("You can't use both 'xctype' and 'xc'!")
160 if parameters.get('autokpt'):
161 if 'kpts' in parameters:
162 raise RuntimeError("You can't use both 'autokpt' and 'kpts'!")
163 if 'ngridk' in parameters:
164 raise RuntimeError(
165 "You can't use both 'autokpt' and 'ngridk'!")
166 if 'ngridk' in parameters:
167 if 'kpts' in parameters:
168 raise RuntimeError("You can't use both 'ngridk' and 'kpts'!")
170 if parameters.get('autoswidth'):
171 if 'smearing' in parameters:
172 raise RuntimeError(
173 "You can't use both 'autoswidth' and 'smearing'!")
174 if 'swidth' in parameters:
175 raise RuntimeError(
176 "You can't use both 'autoswidth' and 'swidth'!")
178 inp = {}
179 inp.update(parameters)
181 if 'xc' in parameters:
182 xctype = {'LDA': 3, # PW92
183 'PBE': 20,
184 'REVPBE': 21,
185 'PBESOL': 22,
186 'WC06': 26,
187 'AM05': 30,
188 'mBJLDA': (100, 208, 12)}[parameters['xc']]
189 inp['xctype'] = xctype
190 del inp['xc']
192 if 'kpts' in parameters:
193 # XXX should generalize kpts handling.
194 from ase.calculators.calculator import kpts2mp
195 mp = kpts2mp(atoms, parameters['kpts'])
196 inp['ngridk'] = tuple(mp)
197 vkloff = [] # is this below correct?
198 for nk in mp:
199 if nk % 2 == 0: # shift kpoint away from gamma point
200 vkloff.append(0.5)
201 else:
202 vkloff.append(0)
203 inp['vkloff'] = vkloff
204 del inp['kpts']
206 if 'smearing' in parameters:
207 name = parameters['smearing'][0].lower()
208 if name == 'methfessel-paxton':
209 stype = parameters['smearing'][2]
210 else:
211 stype = {'gaussian': 0,
212 'fermi-dirac': 3,
213 }[name]
214 inp['stype'] = stype
215 inp['swidth'] = parameters['smearing'][1] / Hartree
216 del inp['smearing']
218 # Elk seems to concatenate path and filename in such a way
219 # that we must put a / at the end:
220 if species_path is not None:
221 inp['sppath'] = f"'{species_path.rstrip('/')}/'"
223 def get_string_value(value_):
224 if isinstance(value_, bool):
225 return f'.{("false", "true")[value_]}.'
226 if isinstance(value_, (numbers.Real, str)):
227 return f'{value_}'
228 if isinstance(value_, collections.abc.Iterable):
229 if all(isinstance(v, (str, numbers.Real)) for v in value_):
230 return ' '.join(get_string_value(v) for v in value_)
231 lines = []
232 for v in value_:
233 if isinstance(v, (str, numbers.Real)):
234 lines.append(get_string_value(v))
235 else:
236 lines.append(' '.join(get_string_value(e) for e in v))
237 return '\n '.join(lines)
238 raise TypeError(f'Field type cannot be {type(value_)}')
240 # write all keys
241 for key, value in inp.items():
242 fd.write(f'{key}\n {get_string_value(value)}\n\n')
244 # cell
245 fd.write('avec\n')
246 for vec in atoms.cell:
247 fd.write(' %.14f %.14f %.14f\n' % tuple(vec / Bohr))
248 fd.write('\n')
250 # atoms
251 species = {}
252 symbols = []
253 for a, (symbol, m) in enumerate(
254 zip(atoms.get_chemical_symbols(),
255 atoms.get_initial_magnetic_moments())):
256 if symbol in species:
257 species[symbol].append((a, m))
258 else:
259 species[symbol] = [(a, m)]
260 symbols.append(symbol)
261 fd.write('atoms\n %d\n' % len(species))
262 # scaled = atoms.get_scaled_positions(wrap=False)
263 scaled = np.linalg.solve(atoms.cell.T, atoms.positions.T).T
264 for symbol in symbols:
265 fd.write(f" '{symbol}.in' : spfname\n")
266 fd.write(' %d\n' % len(species[symbol]))
267 for a, m in species[symbol]:
268 fd.write(' %.14f %.14f %.14f 0.0 0.0 %.14f\n' %
269 (tuple(scaled[a]) + (m,)))
272class ElkReader:
273 def __init__(self, path):
274 self.path = Path(path)
276 def _read_everything(self):
277 yield from self._read_energy()
279 with (self.path / 'INFO.OUT').open() as fd:
280 yield from parse_elk_info(fd)
282 with (self.path / 'EIGVAL.OUT').open() as fd:
283 yield from parse_elk_eigval(fd)
285 with (self.path / 'KPOINTS.OUT').open() as fd:
286 yield from parse_elk_kpoints(fd)
288 def read_everything(self):
289 dct = dict(self._read_everything())
291 # The eigenvalue/occupation tables do not say whether there are
292 # two spins, so we have to reshape them from 1 x K x SB to S x K x B:
293 spinpol = dct.pop('spinpol')
294 if spinpol:
295 for name in 'eigenvalues', 'occupations':
296 array = dct[name]
297 _, nkpts, nbands_double = array.shape
298 assert _ == 1
299 assert nbands_double % 2 == 0
300 nbands = nbands_double // 2
301 newarray = np.empty((2, nkpts, nbands))
302 newarray[0, :, :] = array[0, :, :nbands]
303 newarray[1, :, :] = array[0, :, nbands:]
304 if name == 'eigenvalues':
305 # Verify that eigenvalues are still sorted:
306 diffs = np.diff(newarray, axis=2)
307 assert all(diffs.flat[:] > 0)
308 dct[name] = newarray
309 return dct
311 def _read_energy(self):
312 txt = (self.path / 'TOTENERGY.OUT').read_text()
313 tokens = txt.split()
314 energy = float(tokens[-1]) * Hartree
315 yield 'free_energy', energy
316 yield 'energy', energy
319def parse_elk_kpoints(fd):
320 header = next(fd)
321 lhs, rhs = header.split(':', 1)
322 assert 'k-point, vkl, wkpt' in rhs, header
323 nkpts = int(lhs.strip())
325 kpts = np.empty((nkpts, 3))
326 weights = np.empty(nkpts)
328 for ikpt in range(nkpts):
329 line = next(fd)
330 tokens = line.split()
331 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
332 weights[ikpt] = float(tokens[4])
333 yield 'ibz_kpoints', kpts
334 yield 'kpoint_weights', weights
337def parse_elk_info(fd):
338 dct = collections.defaultdict(list)
339 fd = iter(fd)
341 spinpol = None
342 converged = False
343 actually_did_not_converge = False
344 # Legacy code kept track of both these things, which is strange.
345 # Why could a file both claim to converge *and* not converge?
347 # We loop over all lines and extract also data that occurs
348 # multiple times (e.g. in multiple SCF steps)
349 for line in fd:
350 # "name of quantity : 1 2 3"
351 tokens = line.split(':', 1)
352 if len(tokens) == 2:
353 lhs, rhs = tokens
354 dct[lhs.strip()].append(rhs.strip())
356 elif line.startswith('Convergence targets achieved'):
357 converged = True
358 elif 'reached self-consistent loops maximum' in line.lower():
359 actually_did_not_converge = True
361 if 'Spin treatment' in line:
362 # (Somewhat brittle doing multi-line stuff here.)
363 line = next(fd)
364 spinpol = line.strip() == 'spin-polarised'
366 yield 'converged', converged and not actually_did_not_converge
367 if spinpol is None:
368 raise RuntimeError('Could not determine spin treatment')
369 yield 'spinpol', spinpol
371 if 'Fermi' in dct:
372 yield 'fermi_level', float(dct['Fermi'][-1]) * Hartree
374 if 'total force' in dct:
375 forces = []
376 for line in dct['total force']:
377 forces.append(line.split())
378 yield 'forces', np.array(forces, float) * (Hartree / Bohr)
381def parse_elk_eigval(fd):
383 def match_int(line, word):
384 number, colon, word1 = line.split()
385 assert word1 == word
386 assert colon == ':'
387 return int(number)
389 def skip_spaces(line=''):
390 while not line.strip():
391 line = next(fd)
392 return line
394 line = skip_spaces()
395 nkpts = match_int(line, 'nkpt') # 10 : nkpts
396 line = next(fd)
397 nbands = match_int(line, 'nstsv') # 15 : nstsv
399 eigenvalues = np.empty((nkpts, nbands))
400 occupations = np.empty((nkpts, nbands))
401 kpts = np.empty((nkpts, 3))
403 for ikpt in range(nkpts):
404 line = skip_spaces()
405 tokens = line.split()
406 assert tokens[-1] == 'vkl', tokens
407 assert ikpt + 1 == int(tokens[0])
408 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
410 line = next(fd) # "(state, eigenvalue and occupancy below)"
411 assert line.strip().startswith('(state,'), line
412 for iband in range(nbands):
413 line = next(fd)
414 tokens = line.split() # (band number, eigenval, occ)
415 assert iband + 1 == int(tokens[0])
416 eigenvalues[ikpt, iband] = float(tokens[1])
417 occupations[ikpt, iband] = float(tokens[2])
419 yield 'ibz_kpoints', kpts
420 yield 'eigenvalues', eigenvalues[None] * Hartree
421 yield 'occupations', occupations[None]