Coverage for /builds/ase/ase/ase/io/elk.py: 86.49%
259 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import collections
4from pathlib import Path
6import numpy as np
8from ase import Atoms
9from ase.units import Bohr, Hartree
10from ase.utils import reader, writer
12elk_parameters = {'swidth': Hartree}
15@reader
16def read_elk(fd):
17 """Import ELK atoms definition.
19 Reads unitcell, atom positions, magmoms from elk.in/GEOMETRY.OUT file.
20 """
22 lines = fd.readlines()
24 scale = np.ones(4) # unit cell scale
25 positions = []
26 cell = []
27 symbols = []
28 magmoms = []
30 # find cell scale
31 for n, line in enumerate(lines):
32 if line.split() == []:
33 continue
34 if line.strip() == 'scale':
35 scale[0] = float(lines[n + 1])
36 elif line.startswith('scale'):
37 scale[int(line.strip()[-1])] = float(lines[n + 1])
38 for n, line in enumerate(lines):
39 if line.split() == []:
40 continue
41 if line.startswith('avec'):
42 cell = np.array(
43 [[float(v) * scale[1] for v in lines[n + 1].split()],
44 [float(v) * scale[2] for v in lines[n + 2].split()],
45 [float(v) * scale[3] for v in lines[n + 3].split()]])
46 if line.startswith('atoms'):
47 lines1 = lines[n + 1:] # start subsearch
48 spfname = []
49 natoms = []
50 atpos = []
51 bfcmt = []
52 for n1, line1 in enumerate(lines1):
53 if line1.split() == []:
54 continue
55 if 'spfname' in line1:
56 spfnamenow = lines1[n1].split()[0]
57 spfname.append(spfnamenow)
58 natomsnow = int(lines1[n1 + 1].split()[0])
59 natoms.append(natomsnow)
60 atposnow = []
61 bfcmtnow = []
62 for line in lines1[n1 + 2:n1 + 2 + natomsnow]:
63 atposnow.append([float(v) for v in line.split()[0:3]])
64 if len(line.split()) == 6: # bfcmt present
65 bfcmtnow.append(
66 [float(v) for v in line.split()[3:]])
67 atpos.append(atposnow)
68 bfcmt.append(bfcmtnow)
69 # symbols, positions, magmoms based on ELK spfname, atpos, and bfcmt
70 symbols = ''
71 positions = []
72 magmoms = []
73 for n, s in enumerate(spfname):
74 symbols += str(s[1:].split('.')[0]) * natoms[n]
75 positions += atpos[n] # assumes fractional coordinates
76 if len(bfcmt[n]) > 0:
77 # how to handle cases of magmoms being one or three dim array?
78 magmoms += [m[-1] for m in bfcmt[n]]
79 atoms = Atoms(symbols, scaled_positions=positions, cell=[1, 1, 1])
80 if len(magmoms) > 0:
81 atoms.set_initial_magnetic_moments(magmoms)
82 # final cell scale
83 cell = cell * scale[0] * Bohr
85 atoms.set_cell(cell, scale_atoms=True)
86 atoms.pbc = True
87 return atoms
90@writer
91def write_elk_in(fd, atoms, parameters=None):
92 if parameters is None:
93 parameters = {}
95 parameters = dict(parameters)
97 species_path = parameters.pop('species_dir', None)
98 species_path = parameters.pop('sppath', species_path)
100 if parameters.get('spinpol') is None:
101 if atoms.get_initial_magnetic_moments().any():
102 parameters['spinpol'] = True
104 if 'xctype' in parameters:
105 if 'xc' in parameters:
106 raise RuntimeError("You can't use both 'xctype' and 'xc'!")
108 if parameters.get('autokpt'):
109 if 'kpts' in parameters:
110 raise RuntimeError("You can't use both 'autokpt' and 'kpts'!")
111 if 'ngridk' in parameters:
112 raise RuntimeError(
113 "You can't use both 'autokpt' and 'ngridk'!")
114 if 'ngridk' in parameters:
115 if 'kpts' in parameters:
116 raise RuntimeError("You can't use both 'ngridk' and 'kpts'!")
118 if parameters.get('autoswidth'):
119 if 'smearing' in parameters:
120 raise RuntimeError(
121 "You can't use both 'autoswidth' and 'smearing'!")
122 if 'swidth' in parameters:
123 raise RuntimeError(
124 "You can't use both 'autoswidth' and 'swidth'!")
126 inp = {}
127 inp.update(parameters)
129 if 'xc' in parameters:
130 xctype = {'LDA': 3, # PW92
131 'PBE': 20,
132 'REVPBE': 21,
133 'PBESOL': 22,
134 'WC06': 26,
135 'AM05': 30,
136 'mBJLDA': (100, 208, 12)}[parameters['xc']]
137 inp['xctype'] = xctype
138 del inp['xc']
140 if 'kpts' in parameters:
141 # XXX should generalize kpts handling.
142 from ase.calculators.calculator import kpts2mp
143 mp = kpts2mp(atoms, parameters['kpts'])
144 inp['ngridk'] = tuple(mp)
145 vkloff = [] # is this below correct?
146 for nk in mp:
147 if nk % 2 == 0: # shift kpoint away from gamma point
148 vkloff.append(0.5)
149 else:
150 vkloff.append(0)
151 inp['vkloff'] = vkloff
152 del inp['kpts']
154 if 'smearing' in parameters:
155 name = parameters.smearing[0].lower()
156 if name == 'methfessel-paxton':
157 stype = parameters.smearing[2]
158 else:
159 stype = {'gaussian': 0,
160 'fermi-dirac': 3,
161 }[name]
162 inp['stype'] = stype
163 inp['swidth'] = parameters.smearing[1]
164 del inp['smearing']
166 # convert keys to ELK units
167 for key, value in inp.items():
168 if key in elk_parameters:
169 inp[key] /= elk_parameters[key]
171 # write all keys
172 for key, value in inp.items():
173 fd.write(f'{key}\n')
174 if isinstance(value, bool):
175 fd.write(f'.{("false", "true")[value]}.\n\n')
176 elif isinstance(value, (int, float)):
177 fd.write(f'{value}\n\n')
178 else:
179 fd.write('%s\n\n' % ' '.join([str(x) for x in value]))
181 # cell
182 fd.write('avec\n')
183 for vec in atoms.cell:
184 fd.write('%.14f %.14f %.14f\n' % tuple(vec / Bohr))
185 fd.write('\n')
187 # atoms
188 species = {}
189 symbols = []
190 for a, (symbol, m) in enumerate(
191 zip(atoms.get_chemical_symbols(),
192 atoms.get_initial_magnetic_moments())):
193 if symbol in species:
194 species[symbol].append((a, m))
195 else:
196 species[symbol] = [(a, m)]
197 symbols.append(symbol)
198 fd.write('atoms\n%d\n' % len(species))
199 # scaled = atoms.get_scaled_positions(wrap=False)
200 scaled = np.linalg.solve(atoms.cell.T, atoms.positions.T).T
201 for symbol in symbols:
202 fd.write(f"'{symbol}.in' : spfname\n")
203 fd.write('%d\n' % len(species[symbol]))
204 for a, m in species[symbol]:
205 fd.write('%.14f %.14f %.14f 0.0 0.0 %.14f\n' %
206 (tuple(scaled[a]) + (m,)))
207 fd.write('\n')
209 # Elk seems to concatenate path and filename in such a way
210 # that we must put a / at the end:
211 if species_path is not None:
212 fd.write('sppath\n')
213 fd.write(f"'{species_path.rstrip('/')}/'\n\n")
216class ElkReader:
217 def __init__(self, path):
218 self.path = Path(path)
220 def _read_everything(self):
221 yield from self._read_energy()
223 with (self.path / 'INFO.OUT').open() as fd:
224 yield from parse_elk_info(fd)
226 with (self.path / 'EIGVAL.OUT').open() as fd:
227 yield from parse_elk_eigval(fd)
229 with (self.path / 'KPOINTS.OUT').open() as fd:
230 yield from parse_elk_kpoints(fd)
232 def read_everything(self):
233 dct = dict(self._read_everything())
235 # The eigenvalue/occupation tables do not say whether there are
236 # two spins, so we have to reshape them from 1 x K x SB to S x K x B:
237 spinpol = dct.pop('spinpol')
238 if spinpol:
239 for name in 'eigenvalues', 'occupations':
240 array = dct[name]
241 _, nkpts, nbands_double = array.shape
242 assert _ == 1
243 assert nbands_double % 2 == 0
244 nbands = nbands_double // 2
245 newarray = np.empty((2, nkpts, nbands))
246 newarray[0, :, :] = array[0, :, :nbands]
247 newarray[1, :, :] = array[0, :, nbands:]
248 if name == 'eigenvalues':
249 # Verify that eigenvalues are still sorted:
250 diffs = np.diff(newarray, axis=2)
251 assert all(diffs.flat[:] > 0)
252 dct[name] = newarray
253 return dct
255 def _read_energy(self):
256 txt = (self.path / 'TOTENERGY.OUT').read_text()
257 tokens = txt.split()
258 energy = float(tokens[-1]) * Hartree
259 yield 'free_energy', energy
260 yield 'energy', energy
263def parse_elk_kpoints(fd):
264 header = next(fd)
265 lhs, rhs = header.split(':', 1)
266 assert 'k-point, vkl, wkpt' in rhs, header
267 nkpts = int(lhs.strip())
269 kpts = np.empty((nkpts, 3))
270 weights = np.empty(nkpts)
272 for ikpt in range(nkpts):
273 line = next(fd)
274 tokens = line.split()
275 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
276 weights[ikpt] = float(tokens[4])
277 yield 'ibz_kpoints', kpts
278 yield 'kpoint_weights', weights
281def parse_elk_info(fd):
282 dct = collections.defaultdict(list)
283 fd = iter(fd)
285 spinpol = None
286 converged = False
287 actually_did_not_converge = False
288 # Legacy code kept track of both these things, which is strange.
289 # Why could a file both claim to converge *and* not converge?
291 # We loop over all lines and extract also data that occurs
292 # multiple times (e.g. in multiple SCF steps)
293 for line in fd:
294 # "name of quantity : 1 2 3"
295 tokens = line.split(':', 1)
296 if len(tokens) == 2:
297 lhs, rhs = tokens
298 dct[lhs.strip()].append(rhs.strip())
300 elif line.startswith('Convergence targets achieved'):
301 converged = True
302 elif 'reached self-consistent loops maximum' in line.lower():
303 actually_did_not_converge = True
305 if 'Spin treatment' in line:
306 # (Somewhat brittle doing multi-line stuff here.)
307 line = next(fd)
308 spinpol = line.strip() == 'spin-polarised'
310 yield 'converged', converged and not actually_did_not_converge
311 if spinpol is None:
312 raise RuntimeError('Could not determine spin treatment')
313 yield 'spinpol', spinpol
315 if 'Fermi' in dct:
316 yield 'fermi_level', float(dct['Fermi'][-1]) * Hartree
318 if 'total force' in dct:
319 forces = []
320 for line in dct['total force']:
321 forces.append(line.split())
322 yield 'forces', np.array(forces, float) * (Hartree / Bohr)
325def parse_elk_eigval(fd):
327 def match_int(line, word):
328 number, colon, word1 = line.split()
329 assert word1 == word
330 assert colon == ':'
331 return int(number)
333 def skip_spaces(line=''):
334 while not line.strip():
335 line = next(fd)
336 return line
338 line = skip_spaces()
339 nkpts = match_int(line, 'nkpt') # 10 : nkpts
340 line = next(fd)
341 nbands = match_int(line, 'nstsv') # 15 : nstsv
343 eigenvalues = np.empty((nkpts, nbands))
344 occupations = np.empty((nkpts, nbands))
345 kpts = np.empty((nkpts, 3))
347 for ikpt in range(nkpts):
348 line = skip_spaces()
349 tokens = line.split()
350 assert tokens[-1] == 'vkl', tokens
351 assert ikpt + 1 == int(tokens[0])
352 kpts[ikpt] = np.array(tokens[1:4]).astype(float)
354 line = next(fd) # "(state, eigenvalue and occupancy below)"
355 assert line.strip().startswith('(state,'), line
356 for iband in range(nbands):
357 line = next(fd)
358 tokens = line.split() # (band number, eigenval, occ)
359 assert iband + 1 == int(tokens[0])
360 eigenvalues[ikpt, iband] = float(tokens[1])
361 occupations[ikpt, iband] = float(tokens[2])
363 yield 'ibz_kpoints', kpts
364 yield 'eigenvalues', eigenvalues[None] * Hartree
365 yield 'occupations', occupations[None]