Commit b4c0c0fb authored by Ask Hjorth Larsen's avatar Ask Hjorth Larsen
Browse files

get coordinates from parser log, or from coordinates file taken from parser...

get coordinates from parser log, or from coordinates file taken from parser log.  Get cell from static/info if it is listed
parent 9d18bf86
......@@ -392,33 +392,6 @@ def parse_input_file(fd):
return namespace.names
def kwargs2cell(kwargs):
# kwargs -> cell + remaining kwargs
# cell will be None if not ASE-compatible.
# Returns numbers verbatim; caller must convert units.
kwargs = normalize_keywords(kwargs)
if boxshape_is_ase_compatible(kwargs):
kwargs.pop('boxshape', None)
if 'lsize' in kwargs:
Lsize = kwargs.pop('lsize')
if not isinstance(Lsize, list):
Lsize = [[Lsize] * 3]
assert len(Lsize) == 1
cell = np.array([2 * float(l) for l in Lsize[0]])
elif 'latticeparameters' in kwargs:
# Eval latparam and latvec
latparam = np.array(kwargs.pop('latticeparameters'), float).T
cell = np.array(kwargs.pop('latticevectors', np.eye(3)), float)
for a, vec in zip(latparam, cell):
vec *= a
assert cell.shape == (3, 3)
cell = None
return cell, kwargs
def boxshape_is_ase_compatible(kwargs):
pdims = int(kwargs.get('periodicdimensions', 0))
default_boxshape = 'parallelepiped' if pdims > 0 else 'minimum'
......@@ -544,14 +517,14 @@ def kwargs2atoms(kwargs, directory=None):
if sum(atoms.pbc) != 0:
raise NotImplementedError('Periodic pdb not supported by ASE.')
if cell is None:
#if cell is None:
# cell could not be established from the file, so we set it on the
# Atoms now if possible:
cell, kwargs = kwargs2cell(kwargs)
if cell is not None:
cell *= length_unit
if cell is not None and atoms is not None:
atoms.cell = cell
# cell, kwargs = kwargs2cell(kwargs)
# if cell is not None:
# cell *= length_unit
# if cell is not None and atoms is not None:
# atoms.cell = cell
# In case of boxshape = sphere and similar, we still do not have
# a cell.
......@@ -568,10 +541,11 @@ def kwargs2atoms(kwargs, directory=None):
rcoords = kwargs.get('reducedcoordinates')
if rcoords is not None:
numbers, rpositions = get_positions_from_block('reducedcoordinates')
if cell is None:
raise ValueError('Cannot figure out what the cell is, '
'and thus cannot interpret reduced coordinates.')
atoms = Atoms(cell=cell, numbers=numbers, scaled_positions=rpositions)
#if cell is None:
# raise ValueError('Cannot figure out what the cell is, '
# 'and thus cannot interpret reduced coordinates.')
atoms = Atoms(cell=cell, numbers=numbers, positions=rpositions)
# Very invalid, but we get (scaled) positions from elsewhere anyway
if atoms is None:
raise OctopusParseError('Apparently there are no atoms.')
from __future__ import print_function
import re
import os
import sys
from glob import glob
from contextlib import contextmanager
import numpy as np
from ase.units import Bohr
import setup_paths
......@@ -13,7 +15,7 @@ from nomadcore.parser_backend import JsonParseEventsWriterBackend
from nomadcore.unit_conversion.unit_conversion import convert_unit, \
from aseoct import Octopus, parse_input_file, kwargs2cell
from aseoct import Octopus, parse_input_file
from octopus_info_parser import parse_infofile
from octopus_logfile_parser import parse_logfile
......@@ -40,6 +42,119 @@ is largely irrelevant.
def parse_gridinfo(metaInfoEnv, pew, fname):
results = {}
with open(fname) as fd:
for line in fd:
if '*** Grid ***' in line:
line = next(fd)
while '***********************' not in line:
if line.startswith('Simulation Box:'):
line = next(fd)
while line == '' or line[0].isspace():
m = re.match(r'\s*Type\s*=\s*(\S+)', line)
if m:
results['boxshape'] =
m = re.match(r'\s*Octopus will treat the system as periodic in (\S+) dim', line)
if m:
results['npbc'] = int(
m = re.match(r'\s*Lattice Vectors \[(\S+)\]', line)
if m:
results['unit'] =
cell = np.array([[float(x) for x in next(fd).split()] for _ in range(3)])
results['cell'] = cell
line = next(fd)
if line.startswith('Main mesh:'):
line = next(fd)
while line == '' or line[0].isspace():
m = re.match(r'\s*Spacing\s*\[(\S+)\]\s*=\s*\(\s*(\S+?),\s*(\S+?),\s*(\S+?)\)', line)
if m:
results['unit'] =
results['spacing'] = np.array(, 3, 4)).astype(float)
line = next(fd)
return results
def parse_coordinates_from_parserlog(fname):
results = {}
def buildblock(block):
imax = 1 + max(b[0] for b in block)
jmax = 1 + max(b[1] for b in block)
array = np.zeros((imax, jmax), object)
for i, j, val in block:
array[i, j] = val
labels = array[:, 0].astype(str)[1:-1]
numbers = array[:, 1:4].astype(float)
return labels, numbers
with open(fname) as fd:
blockname = None
block = None
for line in fd:
m = re.match(r"Opened block '(Coordinates|ReducedCoordinates)'", line)
if m:
blockname =
block = []
m = re.match(r'Closed block', line)
if m:
if blockname is None:
continue # This is some other block
labels, numbers = buildblock(block)
results['labels'] = labels
if blockname == 'Coordinates':
results['coords'] = numbers
assert blockname == 'ReducedCoordinates'
results['rcoords'] = numbers
blockname = None
block = None
if block is not None:
# Example: (0, 0) = "Ag"
paren_expression = r'\s*\((\d+),\s*(\d+)\)\s*=\s*(.+)'
m = re.match(paren_expression, line)
if not m:
# Example: Coordinates (0, 0) = "O"
m = re.match(r'\s*' + blockname + paren_expression, line)
if not m:
# Example: Coordinates[0][0] = "H"
m = re.match(r'\s*' + blockname + r'\s*\[(\d+)\]\[(\d+)\]\s*=\s*(.+)', line)
assert m is not None
i, j, val =, 2, 3)
block.append((int(i), int(j), val))
for name in ['PDBCoordinates',
m = re.match(name + r'\s*=\s*"(.+?)"', line)
if m:
assert block is None
results[name] =
m = re.match('XSFCoordinatesAnimStep\s*=\s*(\d+)', line)
if m:
results['XSFCoordinatesAnimStep'] = int(
return results
def normalize_names(names):
return [name.lower() for name in names]
......@@ -199,10 +314,13 @@ def register_octopus_keywords(pew, category, kwargs):
val = kwargs[keyword]
name, value = regularize_metadata_entry(normalized_name, val)
except (KeyError, ValueError): # unknown normalized_name or cannot convert
except Exception: # unknown normalized_name or cannot convert
# We can't crash on unknown keywords because we must support
# versions old and new alike.
# Some keywords (e.g. Spacing) specify float as type, but they
# can actually be blocks. (block is a type itself)
pew.addValue(name, value)
......@@ -275,23 +393,60 @@ def parse(fname, fd):
print('Add parsed values', file=fd)
with open_section('section_system') as system_gid:
# The Atoms object will always have a cell, even if it was not
# used in the Octopus calculation! Thus, to be more honest,
# we re-extract the cell at a level where we can distinguish:
cell, _unused = kwargs2cell(kwargs)
if cell is not None:
# ...and yet we add the ASE cell because that one is
# always measured in angstroms.
gridinfo = parse_gridinfo(metaInfoEnv, pew, fname)
cell_unit = gridinfo['unit']
if 'cell' in gridinfo:
cell = gridinfo['cell']
if cell_unit == 'A':
cell /= Bohr # Always Bohr now
assert lunit == 'b'
convert_unit(atoms.cell, 'angstrom'))
convert_unit(cell, 'bohr'))
# We will get the positions from the parser log.
coordinfo = parse_coordinates_from_parserlog(parser_log_path)
atoms = None
coords = None
if 'PDBCoordinates' in coordinfo:
atoms = read(coordinfo['PDBCoordinates'], format='proteindatabank')
elif 'XYZCoordinates' in coordinfo:
atoms = read(coordinfo['XYZCoordinates'], format='xyz')
elif 'XSFCoordinates' in coordinfo:
if 'XSFCoordinatesAnimStep' in coordinfo:
xxxx # read correct step. Take 1-indexation into account
atoms = read(coordinfo['XSFCoordinates'], format='xsf')
elif 'coords' in coordinfo:
coords = coordinfo['coords']
elif 'rcoords' in coordinfo:
# unit will be Bohr cf. handling of cell above
coords =['rcoords'], cell)
if atoms is not None:
coords = atoms.positions / Bohr
assert coords is not None, 'Cannot find coordinates'
labels = coordinfo.get('labels')
if labels is None:
labels = np.array(atoms.get_chemical_symbols())
pbc = np.zeros(3)
if 'npbc' in gridinfo:
pbc = np.zeros(3, bool)
pbc[:gridinfo['npbc']] = True
pew.addArrayValues('configuration_periodic_dimensions', pbc)
pew.addArrayValues('atom_labels', labels)
pew.addArrayValues('atom_positions', convert_unit(coords, 'bohr'))
# XXX FIXME atoms can be labeled in ways not compatible with ASE.
convert_unit(atoms.get_positions(), 'angstrom'))
# np.array(atoms.get_chemical_symbols()))
# convert_unit(atoms.get_positions(), 'angstrom'))
# np.array(atoms.pbc))
with open_section('section_single_configuration_calculation'):
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment