Coverage for /builds/debichem-team/python-ase/ase/io/octopus/input.py : 25.07%
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(l) for l 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 def append(line):
435 _lines.append(line)
437 def extend(lines):
438 _lines.extend(lines)
439 append('')
441 def setvar(key, var):
442 append('%s = %s' % (key, var))
444 for kw in ['lsize', 'latticevectors', 'latticeparameters']:
445 assert kw not in kwargs
447 defaultboxshape = 'parallelepiped' if atoms.pbc.any() else 'minimum'
448 boxshape = kwargs.get('boxshape', defaultboxshape).lower()
449 use_ase_cell = (boxshape == 'parallelepiped')
450 atomskwargs = atoms2kwargs(atoms, use_ase_cell)
452 if use_ase_cell:
453 if 'lsize' in atomskwargs:
454 block = list2block('LSize', atomskwargs['lsize'])
455 elif 'latticevectors' in atomskwargs:
456 extend(list2block('LatticeParameters', [[1., 1., 1.]]))
457 block = list2block('LatticeVectors', atomskwargs['latticevectors'])
458 extend(block)
460 # Allow override or issue errors?
461 pdim = 'periodicdimensions'
462 if pdim in kwargs:
463 if int(kwargs[pdim]) != int(atomskwargs[pdim]):
464 raise ValueError('Cannot reconcile periodicity in input '
465 'with that of Atoms object')
466 setvar('periodicdimensions', atomskwargs[pdim])
468 # We like to output forces
469 if 'output' in kwargs:
470 output_string = kwargs.pop('output')
471 output_tokens = [token.strip()
472 for token in output_string.lower().split('+')]
473 else:
474 output_tokens = []
476 if 'forces' not in output_tokens:
477 output_tokens.append('forces')
478 setvar('output', ' + '.join(output_tokens))
479 # It is illegal to have output forces without any OutputFormat.
480 # Even though the forces are written in the same format no matter
481 # OutputFormat. Thus we have to make one up:
483 if 'outputformat' not in kwargs:
484 setvar('outputformat', 'xcrysden')
486 for key, val in kwargs.items():
487 # Most datatypes are straightforward but blocks require some attention.
488 if isinstance(val, list):
489 append('')
490 dict_data = list2block(key, val)
491 extend(dict_data)
492 else:
493 setvar(key, str(val))
494 append('')
496 coord_block = list2block('Coordinates', atomskwargs['coordinates'])
497 extend(coord_block)
498 return '\n'.join(_lines)
501def atoms2kwargs(atoms, use_ase_cell):
502 kwargs = {}
504 positions = atoms.positions / Bohr
506 if use_ase_cell:
507 cell = atoms.cell / Bohr
508 cell_offset = 0.5 * cell.sum(axis=0)
509 positions -= cell_offset
510 if atoms.cell.orthorhombic:
511 Lsize = 0.5 * np.diag(cell)
512 kwargs['lsize'] = [[repr(size) for size in Lsize]]
513 # ASE uses (0...cell) while Octopus uses -L/2...L/2.
514 # Lsize is really cell / 2, and we have to adjust our
515 # positions by subtracting Lsize (see construction of the coords
516 # block) in non-periodic directions.
517 else:
518 kwargs['latticevectors'] = cell.tolist()
520 types = atoms.info.get('types', {})
522 coord_block = []
523 for sym, pos, tag in zip(atoms.get_chemical_symbols(),
524 positions, atoms.get_tags()):
525 if sym == 'X':
526 sym = types.get((sym, tag))
527 if sym is None:
528 raise ValueError('Cannot represent atom X without tags and '
529 'species info in atoms.info')
530 coord_block.append([repr(sym)] + [repr(x) for x in pos])
532 kwargs['coordinates'] = coord_block
533 npbc = sum(atoms.pbc)
534 for c in range(npbc):
535 if not atoms.pbc[c]:
536 msg = ('Boundary conditions of Atoms object inconsistent '
537 'with requirements of Octopus. pbc must be either '
538 '000, 100, 110, or 111.')
539 raise ValueError(msg)
540 kwargs['periodicdimensions'] = npbc
542 # TODO InitialSpins
543 #
544 # TODO can use maximumiterations + output/outputformat to extract
545 # things from restart file into output files without trouble.
546 #
547 # Velocities etc.?
548 return kwargs
551@reader
552def read_octopus_in(fd, get_kwargs=False):
553 kwargs = parse_input_file(fd)
555 # input files may contain internal references to other files such
556 # as xyz or xsf. We need to know the directory where the file
557 # resides in order to locate those. If fd is a real file
558 # object, it contains the path and we can use it. Else assume
559 # pwd.
560 #
561 # Maybe this is ugly; maybe it can lead to strange bugs if someone
562 # wants a non-standard file-like type. But it's probably better than
563 # failing 'ase gui somedir/inp'
564 try:
565 fname = fd.name
566 except AttributeError:
567 directory = None
568 else:
569 directory = os.path.split(fname)[0]
571 atoms, remaining_kwargs = kwargs2atoms(kwargs, directory=directory)
572 if get_kwargs:
573 return atoms, remaining_kwargs
574 else:
575 return atoms