Coverage for /builds/kinetik161/ase/ase/io/octopus/input.py: 26.29%
369 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import os
2import re
4import numpy as np
6from ase import Atoms
7from ase.calculators.calculator import kpts2ndarray
8from ase.data import atomic_numbers
9from ase.io import read
10from ase.units import Bohr, Hartree
11from ase.utils import reader
13special_ase_keywords = {'kpts'}
16def process_special_kwargs(atoms, kwargs):
17 kwargs = kwargs.copy()
18 kpts = kwargs.pop('kpts', None)
19 if kpts is not None:
20 for kw in ['kpoints', 'reducedkpoints', 'kpointsgrid']:
21 if kw in kwargs:
22 raise ValueError('k-points specified multiple times')
24 kptsarray = kpts2ndarray(kpts, atoms)
25 nkpts = len(kptsarray)
26 fullarray = np.empty((nkpts, 4))
27 fullarray[:, 0] = 1.0 / nkpts # weights
28 fullarray[:, 1:4] = kptsarray
29 kwargs['kpointsreduced'] = fullarray.tolist()
31 # TODO xc=LDA/PBE etc.
33 # The idea is to get rid of the special keywords, since the rest
34 # will be passed to Octopus
35 # XXX do a better check of this
36 for kw in special_ase_keywords:
37 assert kw not in kwargs, kw
38 return kwargs
41class OctopusParseError(Exception):
42 pass # Cannot parse input file
45# Parse value as written in input file *or* something that one would be
46# passing to the ASE interface, i.e., this might already be a boolean
47def octbool2bool(value):
48 value = value.lower()
49 if isinstance(value, int):
50 return bool(value)
51 if value in ['true', 't', 'yes', '1']:
52 return True
53 elif value in ['no', 'f', 'false', '0']:
54 return False
55 else:
56 raise ValueError(f'Failed to interpret "{value}" as a boolean.')
59def list2block(name, rows):
60 """Construct 'block' of Octopus input.
62 convert a list of rows to a string with the format x | x | ....
63 for the octopus input file"""
64 lines = []
65 lines.append('%' + name)
66 for row in rows:
67 lines.append(' ' + ' | '.join(str(obj) for obj in row))
68 lines.append('%')
69 return lines
72def normalize_keywords(kwargs):
73 """Reduce keywords to unambiguous form (lowercase)."""
74 newkwargs = {}
75 for arg, value in kwargs.items():
76 lkey = arg.lower()
77 newkwargs[lkey] = value
78 return newkwargs
81def input_line_iter(lines):
82 """Convenient iterator for parsing input files 'cleanly'.
84 Discards comments etc."""
85 for line in lines:
86 line = line.split('#')[0].strip()
87 if not line or line.isspace():
88 continue
89 line = line.strip()
90 yield line
93def block2list(namespace, lines, header=None):
94 """Parse lines of block and return list of lists of strings."""
95 lines = iter(lines)
96 block = []
97 if header is None:
98 header = next(lines)
99 assert header.startswith('%'), header
100 name = header[1:].strip().lower()
101 for line in lines:
102 if line.startswith('%'): # Could also say line == '%' most likely.
103 break
104 tokens = [namespace.evaluate(token)
105 for token in line.strip().split('|')]
106 # XXX will fail for string literals containing '|'
107 block.append(tokens)
108 return name, block
111class OctNamespace:
112 def __init__(self):
113 self.names = {}
114 self.consts = {'pi': np.pi,
115 'angstrom': 1. / Bohr,
116 'ev': 1. / Hartree,
117 'yes': True,
118 'no': False,
119 't': True,
120 'f': False,
121 'i': 1j, # This will probably cause trouble
122 'true': True,
123 'false': False}
125 def evaluate(self, value):
126 value = value.strip()
128 for char in '"', "'": # String literal
129 if value.startswith(char):
130 assert value.endswith(char)
131 return value
133 value = value.lower()
135 if value in self.consts: # boolean or other constant
136 return self.consts[value]
138 if value in self.names: # existing variable
139 return self.names[value]
141 try: # literal integer
142 v = int(value)
143 except ValueError:
144 pass
145 else:
146 if v == float(v):
147 return v
149 try: # literal float
150 return float(value)
151 except ValueError:
152 pass
154 if ('*' in value or '/' in value
155 and not any(char in value for char in '()+')):
156 floatvalue = 1.0
157 op = '*'
158 for token in re.split(r'([\*/])', value):
159 if token in '*/':
160 op = token
161 continue
163 v = self.evaluate(token)
165 try:
166 v = float(v)
167 except TypeError:
168 try:
169 v = complex(v)
170 except ValueError:
171 break
172 except ValueError:
173 break # Cannot evaluate expression
174 else:
175 if op == '*':
176 floatvalue *= v
177 else:
178 assert op == '/', op
179 floatvalue /= v
180 else: # Loop completed successfully
181 return floatvalue
182 return value # unknown name, or complex arithmetic expression
184 def add(self, name, value):
185 value = self.evaluate(value)
186 self.names[name.lower().strip()] = value
189def parse_input_file(fd):
190 namespace = OctNamespace()
191 lines = input_line_iter(fd)
192 blocks = {}
193 while True:
194 try:
195 line = next(lines)
196 except StopIteration:
197 break
198 else:
199 if line.startswith('%'):
200 name, value = block2list(namespace, lines, header=line)
201 blocks[name] = value
202 else:
203 tokens = line.split('=', 1)
204 assert len(tokens) == 2, tokens
205 name, value = tokens
206 namespace.add(name, value)
208 namespace.names.update(blocks)
209 return namespace.names
212def kwargs2cell(kwargs):
213 # kwargs -> cell + remaining kwargs
214 # cell will be None if not ASE-compatible.
215 #
216 # Returns numbers verbatim; caller must convert units.
217 kwargs = normalize_keywords(kwargs)
219 if boxshape_is_ase_compatible(kwargs):
220 kwargs.pop('boxshape', None)
221 if 'lsize' in kwargs:
222 Lsize = kwargs.pop('lsize')
223 if not isinstance(Lsize, list):
224 Lsize = [[Lsize] * 3]
225 assert len(Lsize) == 1
226 cell = np.array([2 * float(x) for x in Lsize[0]])
227 elif 'latticeparameters' in kwargs:
228 # Eval latparam and latvec
229 latparam = np.array(kwargs.pop('latticeparameters'), float).T
230 cell = np.array(kwargs.pop('latticevectors', np.eye(3)), float)
231 for a, vec in zip(latparam, cell):
232 vec *= a
233 assert cell.shape == (3, 3)
234 else:
235 cell = None
236 return cell, kwargs
239def boxshape_is_ase_compatible(kwargs):
240 pdims = int(kwargs.get('periodicdimensions', 0))
241 default_boxshape = 'parallelepiped' if pdims > 0 else 'minimum'
242 boxshape = kwargs.get('boxshape', default_boxshape).lower()
243 # XXX add support for experimental keyword 'latticevectors'
244 return boxshape == 'parallelepiped'
247def kwargs2atoms(kwargs, directory=None):
248 """Extract atoms object from keywords and return remaining keywords.
250 Some keyword arguments may refer to files. The directory keyword
251 may be necessary to resolve the paths correctly, and is used for
252 example when running 'ase gui somedir/inp'."""
253 kwargs = normalize_keywords(kwargs)
255 # Only input units accepted nowadays are 'atomic'.
256 # But if we are loading an old file, and it specifies something else,
257 # we can be sure that the user wanted that back then.
259 coord_keywords = ['coordinates',
260 'xyzcoordinates',
261 'pdbcoordinates',
262 'reducedcoordinates',
263 'xsfcoordinates',
264 'xsfcoordinatesanimstep']
266 nkeywords = 0
267 for keyword in coord_keywords:
268 if keyword in kwargs:
269 nkeywords += 1
270 if nkeywords == 0:
271 raise OctopusParseError('No coordinates')
272 elif nkeywords > 1:
273 raise OctopusParseError('Multiple coordinate specifications present. '
274 'This may be okay in Octopus, but we do not '
275 'implement it.')
277 def get_positions_from_block(keyword):
278 # %Coordinates or %ReducedCoordinates -> atomic numbers, positions.
279 block = kwargs.pop(keyword)
280 positions = []
281 numbers = []
282 tags = []
283 types = {}
284 for row in block:
285 assert len(row) in [ndims + 1, ndims + 2]
286 row = row[:ndims + 1]
287 sym = row[0]
288 assert sym.startswith('"') or sym.startswith("'")
289 assert sym[0] == sym[-1]
290 sym = sym[1:-1]
291 pos0 = np.zeros(3)
292 ndim = int(kwargs.get('dimensions', 3))
293 pos0[:ndim] = [float(element) for element in row[1:]]
294 number = atomic_numbers.get(sym) # Use 0 ~ 'X' for unknown?
295 tag = 0
296 if number is None:
297 if sym not in types:
298 tag = len(types) + 1
299 types[sym] = tag
300 number = 0
301 tag = types[sym]
302 tags.append(tag)
303 numbers.append(number)
304 positions.append(pos0)
305 positions = np.array(positions)
306 tags = np.array(tags, int)
307 if types:
308 ase_types = {}
309 for sym, tag in types.items():
310 ase_types[('X', tag)] = sym
311 info = {'types': ase_types} # 'info' dict for Atoms object
312 else:
313 tags = None
314 info = None
315 return numbers, positions, tags, info
317 def read_atoms_from_file(fname, fmt):
318 assert fname.startswith('"') or fname.startswith("'")
319 assert fname[0] == fname[-1]
320 fname = fname[1:-1]
321 if directory is not None:
322 fname = os.path.join(directory, fname)
323 # XXX test xyz, pbd and xsf
324 if fmt == 'xsf' and 'xsfcoordinatesanimstep' in kwargs:
325 anim_step = kwargs.pop('xsfcoordinatesanimstep')
326 theslice = slice(anim_step, anim_step + 1, 1)
327 # XXX test animstep
328 else:
329 theslice = slice(None, None, 1)
330 images = read(fname, theslice, fmt)
331 if len(images) != 1:
332 raise OctopusParseError('Expected only one image. Don\'t know '
333 'what to do with %d images.' % len(images))
334 return images[0]
336 # We will attempt to extract cell and pbc from kwargs if 'lacking'.
337 # But they might have been left unspecified on purpose.
338 #
339 # We need to keep track of these two variables "externally"
340 # because the Atoms object assigns values when they are not given.
341 cell = None
342 pbc = None
343 adjust_positions_by_half_cell = False
345 atoms = None
346 xsfcoords = kwargs.pop('xsfcoordinates', None)
347 if xsfcoords is not None:
348 atoms = read_atoms_from_file(xsfcoords, 'xsf')
349 atoms.positions *= Bohr
350 atoms.cell *= Bohr
351 # As it turns out, non-periodic xsf is not supported by octopus.
352 # Also, it only supports fully periodic or fully non-periodic....
353 # So the only thing that we can test here is 3D fully periodic.
354 if sum(atoms.pbc) != 3:
355 raise NotImplementedError('XSF not fully periodic with Octopus')
356 cell = atoms.cell
357 pbc = atoms.pbc
358 # Position adjustment doesn't actually matter but this should work
359 # most 'nicely':
360 adjust_positions_by_half_cell = False
361 xyzcoords = kwargs.pop('xyzcoordinates', None)
362 if xyzcoords is not None:
363 atoms = read_atoms_from_file(xyzcoords, 'xyz')
364 atoms.positions *= Bohr
365 adjust_positions_by_half_cell = True
366 pdbcoords = kwargs.pop('pdbcoordinates', None)
367 if pdbcoords is not None:
368 atoms = read_atoms_from_file(pdbcoords, 'pdb')
369 pbc = atoms.pbc
370 adjust_positions_by_half_cell = True
371 # Due to an error in ASE pdb, we can only test the nonperiodic case.
372 # atoms.cell *= Bohr # XXX cell? Not in nonperiodic case...
373 atoms.positions *= Bohr
374 if sum(atoms.pbc) != 0:
375 raise NotImplementedError('Periodic pdb not supported by ASE.')
377 if cell is None:
378 # cell could not be established from the file, so we set it on the
379 # Atoms now if possible:
380 cell, kwargs = kwargs2cell(kwargs)
381 if cell is not None:
382 cell *= Bohr
383 if cell is not None and atoms is not None:
384 atoms.cell = cell
385 # In case of boxshape = sphere and similar, we still do not have
386 # a cell.
388 ndims = int(kwargs.get('dimensions', 3))
389 if ndims != 3:
390 raise NotImplementedError('Only 3D calculations supported.')
392 coords = kwargs.get('coordinates')
393 if coords is not None:
394 numbers, pos, tags, info = get_positions_from_block('coordinates')
395 pos *= Bohr
396 adjust_positions_by_half_cell = True
397 atoms = Atoms(cell=cell, numbers=numbers, positions=pos,
398 tags=tags, info=info)
399 rcoords = kwargs.get('reducedcoordinates')
400 if rcoords is not None:
401 numbers, spos, tags, info = get_positions_from_block(
402 'reducedcoordinates')
403 if cell is None:
404 raise ValueError('Cannot figure out what the cell is, '
405 'and thus cannot interpret reduced coordinates.')
406 atoms = Atoms(cell=cell, numbers=numbers, scaled_positions=spos,
407 tags=tags, info=info)
408 if atoms is None:
409 raise OctopusParseError('Apparently there are no atoms.')
411 # Either we have non-periodic BCs or the atoms object already
412 # got its BCs from reading the file. In the latter case
413 # we shall override only if PeriodicDimensions was given specifically:
415 if pbc is None:
416 pdims = int(kwargs.pop('periodicdimensions', 0))
417 pbc = np.zeros(3, dtype=bool)
418 pbc[:pdims] = True
419 atoms.pbc = pbc
421 if (cell is not None and cell.shape == (3,)
422 and adjust_positions_by_half_cell):
423 nonpbc = (atoms.pbc == 0)
424 atoms.positions[:, nonpbc] += np.array(cell)[None, nonpbc] / 2.0
426 return atoms, kwargs
429def generate_input(atoms, kwargs):
430 """Convert atoms and keyword arguments to Octopus input file."""
431 _lines = []
433 kwargs = {key.lower(): value for key, value in kwargs.items()}
435 def append(line):
436 _lines.append(line)
438 def extend(lines):
439 _lines.extend(lines)
440 append('')
442 def setvar(key, var):
443 append(f'{key} = {var}')
445 for kw in ['lsize', 'latticevectors', 'latticeparameters']:
446 assert kw not in kwargs
448 defaultboxshape = 'parallelepiped' if atoms.cell.rank > 0 else 'minimum'
449 boxshape = kwargs.pop('boxshape', defaultboxshape).lower()
450 use_ase_cell = (boxshape == 'parallelepiped')
451 atomskwargs = atoms2kwargs(atoms, use_ase_cell)
453 setvar('boxshape', boxshape)
455 if use_ase_cell:
456 if 'reducedcoordinates' in atomskwargs:
457 extend(list2block('LatticeParameters', [[1., 1., 1.]]))
458 block = list2block('LatticeVectors', atomskwargs['latticevectors'])
459 else:
460 assert 'lsize' in atomskwargs
461 block = list2block('LSize', atomskwargs['lsize'])
463 extend(block)
465 # Allow override or issue errors?
466 pdim = 'periodicdimensions'
467 if pdim in kwargs:
468 if int(kwargs[pdim]) != int(atomskwargs[pdim]):
469 raise ValueError('Cannot reconcile periodicity in input '
470 'with that of Atoms object')
471 setvar('periodicdimensions', atomskwargs[pdim])
473 # We should say that we want the forces if the user requests forces.
474 # However the Output keyword has changed between Octopus 10.1 and 11.4.
475 # We'll leave it to the user to manually set keywords correctly.
476 # The most complicated part is that the user probably needs to tell
477 # Octopus to save the forces in xcrysden format in order for us to
478 # parse them.
480 for key, val in kwargs.items():
481 # Most datatypes are straightforward but blocks require some attention.
482 if isinstance(val, list):
483 append('')
484 dict_data = list2block(key, val)
485 extend(dict_data)
486 else:
487 setvar(key, str(val))
488 append('')
490 if 'reducedcoordinates' in atomskwargs:
491 coord_block = list2block('ReducedCoordinates',
492 atomskwargs['reducedcoordinates'])
493 else:
494 coord_block = list2block('Coordinates', atomskwargs['coordinates'])
495 extend(coord_block)
496 return '\n'.join(_lines)
499def atoms2kwargs(atoms, use_ase_cell):
500 kwargs = {}
502 if any(atoms.pbc):
503 coordtype = 'reducedcoordinates'
504 coords = atoms.get_scaled_positions(wrap=False)
505 else:
506 coordtype = 'coordinates'
507 coords = atoms.positions / Bohr
509 if use_ase_cell:
510 cell = atoms.cell / Bohr
511 if coordtype == 'coordinates':
512 cell_offset = 0.5 * cell.sum(axis=0)
513 coords -= cell_offset
514 else:
515 assert coordtype == 'reducedcoordinates'
517 if atoms.cell.orthorhombic:
518 Lsize = 0.5 * np.diag(cell)
519 kwargs['lsize'] = [[repr(size) for size in Lsize]]
520 # ASE uses (0...cell) while Octopus uses -L/2...L/2.
521 # Lsize is really cell / 2, and we have to adjust our
522 # positions by subtracting Lsize (see construction of the coords
523 # block) in non-periodic directions.
524 else:
525 kwargs['latticevectors'] = cell.tolist()
527 types = atoms.info.get('types', {})
529 coord_block = []
530 for sym, pos, tag in zip(atoms.symbols, coords, atoms.get_tags()):
531 if sym == 'X':
532 sym = types.get((sym, tag))
533 if sym is None:
534 raise ValueError('Cannot represent atom X without tags and '
535 'species info in atoms.info')
536 coord_block.append([repr(sym)] + [repr(x) for x in pos])
538 kwargs[coordtype] = coord_block
539 npbc = sum(atoms.pbc)
540 for c in range(npbc):
541 if not atoms.pbc[c]:
542 msg = ('Boundary conditions of Atoms object inconsistent '
543 'with requirements of Octopus. pbc must be either '
544 '000, 100, 110, or 111.')
545 raise ValueError(msg)
546 kwargs['periodicdimensions'] = npbc
548 # TODO InitialSpins
549 #
550 # TODO can use maximumiterations + output/outputformat to extract
551 # things from restart file into output files without trouble.
552 #
553 # Velocities etc.?
554 return kwargs
557@reader
558def read_octopus_in(fd, get_kwargs=False):
559 kwargs = parse_input_file(fd)
561 # input files may contain internal references to other files such
562 # as xyz or xsf. We need to know the directory where the file
563 # resides in order to locate those. If fd is a real file
564 # object, it contains the path and we can use it. Else assume
565 # pwd.
566 #
567 # Maybe this is ugly; maybe it can lead to strange bugs if someone
568 # wants a non-standard file-like type. But it's probably better than
569 # failing 'ase gui somedir/inp'
570 try:
571 fname = fd.name
572 except AttributeError:
573 directory = None
574 else:
575 directory = os.path.split(fname)[0]
577 atoms, remaining_kwargs = kwargs2atoms(kwargs, directory=directory)
578 if get_kwargs:
579 return atoms, remaining_kwargs
580 else:
581 return atoms