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