Coverage for /builds/kinetik161/ase/ase/io/elk.py: 86.33%
256 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
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
10elk_parameters = {'swidth': Hartree}
13@reader
14def read_elk(fd):
15 """Import ELK atoms definition.
17 Reads unitcell, atom positions, magmoms from elk.in/GEOMETRY.OUT file.
18 """
20 lines = fd.readlines()
22 scale = np.ones(4) # unit cell scale
23 positions = []
24 cell = []
25 symbols = []
26 magmoms = []
28 # find cell scale
29 for n, line in enumerate(lines):
30 if line.split() == []:
31 continue
32 if line.strip() == 'scale':
33 scale[0] = float(lines[n + 1])
34 elif line.startswith('scale'):
35 scale[int(line.strip()[-1])] = float(lines[n + 1])
36 for n, line in enumerate(lines):
37 if line.split() == []:
38 continue
39 if line.startswith('avec'):
40 cell = np.array(
41 [[float(v) * scale[1] for v in lines[n + 1].split()],
42 [float(v) * scale[2] for v in lines[n + 2].split()],
43 [float(v) * scale[3] for v in lines[n + 3].split()]])
44 if line.startswith('atoms'):
45 lines1 = lines[n + 1:] # start subsearch
46 spfname = []
47 natoms = []
48 atpos = []
49 bfcmt = []
50 for n1, line1 in enumerate(lines1):
51 if line1.split() == []:
52 continue
53 if 'spfname' in line1:
54 spfnamenow = lines1[n1].split()[0]
55 spfname.append(spfnamenow)
56 natomsnow = int(lines1[n1 + 1].split()[0])
57 natoms.append(natomsnow)
58 atposnow = []
59 bfcmtnow = []
60 for line in lines1[n1 + 2:n1 + 2 + natomsnow]:
61 atposnow.append([float(v) for v in line.split()[0:3]])
62 if len(line.split()) == 6: # bfcmt present
63 bfcmtnow.append(
64 [float(v) for v in line.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
83 atoms.set_cell(cell, scale_atoms=True)
84 atoms.pbc = True
85 return atoms
88@writer
89def write_elk_in(fd, atoms, parameters=None):
90 if parameters is None:
91 parameters = {}
93 parameters = dict(parameters)
94 species_path = parameters.pop('species_dir', None)
96 if parameters.get('spinpol') is None:
97 if atoms.get_initial_magnetic_moments().any():
98 parameters['spinpol'] = True
100 if 'xctype' in parameters:
101 if 'xc' in parameters:
102 raise RuntimeError("You can't use both 'xctype' and 'xc'!")
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'!")
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'!")
122 inp = {}
123 inp.update(parameters)
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']
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']
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']
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]
167 # write all keys
168 for key, value in inp.items():
169 fd.write(f'{key}\n')
170 if isinstance(value, bool):
171 fd.write(f'.{("false", "true")[value]}.\n\n')
172 elif isinstance(value, (int, float)):
173 fd.write(f'{value}\n\n')
174 else:
175 fd.write('%s\n\n' % ' '.join([str(x) for x in value]))
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')
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(f"'{symbol}.in' : spfname\n")
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,)))
204 # if sppath is present in elk.in it overwrites species blocks!
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")
212class ElkReader:
213 def __init__(self, path):
214 self.path = Path(path)
216 def _read_everything(self):
217 yield from self._read_energy()
219 with (self.path / 'INFO.OUT').open() as fd:
220 yield from parse_elk_info(fd)
222 with (self.path / 'EIGVAL.OUT').open() as fd:
223 yield from parse_elk_eigval(fd)
225 with (self.path / 'KPOINTS.OUT').open() as fd:
226 yield from parse_elk_kpoints(fd)
228 def read_everything(self):
229 dct = dict(self._read_everything())
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
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
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())
265 kpts = np.empty((nkpts, 3))
266 weights = np.empty(nkpts)
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
277def parse_elk_info(fd):
278 dct = collections.defaultdict(list)
279 fd = iter(fd)
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?
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())
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
301 if 'Spin treatment' in line:
302 # (Somewhat brittle doing multi-line stuff here.)
303 line = next(fd)
304 spinpol = line.strip() == 'spin-polarised'
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
311 if 'Fermi' in dct:
312 yield 'fermi_level', float(dct['Fermi'][-1]) * Hartree
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)
321def parse_elk_eigval(fd):
323 def match_int(line, word):
324 number, colon, word1 = line.split()
325 assert word1 == word
326 assert colon == ':'
327 return int(number)
329 def skip_spaces(line=''):
330 while not line.strip():
331 line = next(fd)
332 return line
334 line = skip_spaces()
335 nkpts = match_int(line, 'nkpt') # 10 : nkpts
336 line = next(fd)
337 nbands = match_int(line, 'nstsv') # 15 : nstsv
339 eigenvalues = np.empty((nkpts, nbands))
340 occupations = np.empty((nkpts, nbands))
341 kpts = np.empty((nkpts, 3))
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)
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])
359 yield 'ibz_kpoints', kpts
360 yield 'eigenvalues', eigenvalues[None] * Hartree
361 yield 'occupations', occupations[None]