parser_vasprun.py

# Copyright 2016-2018 Fawzi Mohamed, Lauri Himanen, Danio Brambila, Ankit Kariryaa, Henning Glawe
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

from __future__ import division
from builtins import range
from builtins import object
import xml.etree.ElementTree
from xml.etree.ElementTree import ParseError
import sys
import bisect
from datetime import datetime
import os
import re
import traceback
import numpy as np
import ase.geometry
import ase.data
from math import pi
import xml.etree.ElementTree as ET
import logging

from nomadcore.parser_backend import JsonParseEventsWriterBackend
from nomadcore.local_meta_info import loadJsonFile, InfoKindEl
from nomadcore.unit_conversion.unit_conversion import convert_unit_function
from nomadcore.unit_conversion.unit_conversion import convert_unit

eV2J = convert_unit_function("eV", "J")
eV2JV = np.vectorize(eV2J)



def crystal_structure_from_cell(cell, eps=1e-4):
    """Return the crystal structure as a string calculated from the cell.
    """
    cellpar = ase.geometry.cell_to_cellpar(cell=cell)
    abc = cellpar[:3]
    angles = cellpar[3:] / 180 * pi
    a, b, c = abc
    alpha, beta, gamma = angles
    # According to:
    # https://www.physics-in-a-nutshell.com/article/6/symmetry-crystal-systems-and-bravais-lattices#the-seven-crystal-systems
    # If a=b=c and alpha=beta=gamma=90degrees we have cubic.
    if abc.ptp() < eps and abs(angles - pi / 2).max() < eps:
        return 'cubic'
    elif abc.ptp() < eps and abs(angles - pi / 3).max() < eps:
        return 'fcc'
    elif abc.ptp() < eps and abs(angles - np.arccos(-1 / 3)).max() < eps:
        return 'bcc'
    # If a=b!=c, alpha=beta=gamma=90deg, tetragonal.
    elif abs(a - b) < eps and abs(angles - pi / 2).max() < eps:
        return 'tetragonal'
    elif abs(angles - pi / 2).max() < eps:
        return 'orthorhombic'
    # if a = b != c , alpha = beta = 90deg, gamma = 120deg, hexagonal
    elif (abs(a - b) < eps
          and abs(gamma - pi / 3 * 2) < eps
          and abs(angles[:2] - pi / 2).max() < eps):
        return 'hexagonal'
    elif (c >= a and c >= b and alpha < pi / 2 and
          abs(angles[1:] - pi / 2).max() < eps):
        return 'monoclinic'
    else:
        raise ValueError('Cannot find crystal structure')

vasp_to_metainfo_type_mapping = {
    'string': ['C'],
    'int': ['i'],
    'logical': ['b', 'C'],
    'float': ['f']}

special_points = {
    'cubic': {'Γ': [0, 0, 0],
              'M': [1 / 2, 1 / 2, 0],
              'R': [1 / 2, 1 / 2, 1 / 2],
              'X': [0, 1 / 2, 0]},
    'fcc': {'Γ': [0, 0, 0],
            'K': [3 / 8, 3 / 8, 3 / 4],
            'L': [1 / 2, 1 / 2, 1 / 2],
            'U': [5 / 8, 1 / 4, 5 / 8],
            'W': [1 / 2, 1 / 4, 3 / 4],
            'X': [1 / 2, 0, 1 / 2]},
    'bcc': {'Γ': [0, 0, 0],
            'H': [1 / 2, -1 / 2, 1 / 2],
            'P': [1 / 4, 1 / 4, 1 / 4],
            'N': [0, 0, 1 / 2]},
    'tetragonal': {'Γ': [0, 0, 0],
                   'A': [1 / 2, 1 / 2, 1 / 2],
                   'M': [1 / 2, 1 / 2, 0],
                   'R': [0, 1 / 2, 1 / 2],
                   'X': [0, 1 / 2, 0],
                   'Z': [0, 0, 1 / 2]},
    'orthorhombic': {'Γ': [0, 0, 0],
                     'R': [1 / 2, 1 / 2, 1 / 2],
                     'S': [1 / 2, 1 / 2, 0],
                     'T': [0, 1 / 2, 1 / 2],
                     'U': [1 / 2, 0, 1 / 2],
                     'X': [1 / 2, 0, 0],
                     'Y': [0, 1 / 2, 0],
                     'Z': [0, 0, 1 / 2]},
    'hexagonal': {'Γ': [0, 0, 0],
                  'A': [0, 0, 1 / 2],
                  'H': [1 / 3, 1 / 3, 1 / 2],
                  'K': [1 / 3, 1 / 3, 0],
                  'L': [1 / 2, 0, 1 / 2],
                  'M': [1 / 2, 0, 0]}}


special_paths = {
    'cubic': 'ΓXMΓRX,MR',
    'fcc': 'ΓXWKΓLUWLK,UX',
    'bcc': 'ΓHNΓPH,PN',
    'tetragonal': 'ΓXMΓZRAZXR,MA',
    'orthorhombic': 'ΓXSYΓZURTZ,YT,UX,SR',
    'hexagonal': 'ΓMKΓALHA,LM,KH',
    'monoclinic': 'ΓYHCEM1AXH1,MDZ,YD'}


def get_special_points(cell, eps=1e-4):
    """Return dict of special points.

    The definitions are from a paper by Wahyu Setyawana and Stefano
    Curtarolo::

        http://dx.doi.org/10.1016/j.commatsci.2010.05.010

    lattice: str
        One of the following: cubic, fcc, bcc, orthorhombic, tetragonal,
        hexagonal or monoclinic.
    cell: 3x3 ndarray
        Unit cell.
    eps: float
        Tolerance for cell-check.
    """

    lattice = crystal_structure_from_cell(cell)

    cellpar = ase.geometry.cell_to_cellpar(cell=cell)
    abc = cellpar[:3]
    angles = cellpar[3:] / 180 * pi
    a, b, c = abc
    alpha, beta, gamma = angles

    # Check that the unit-cells are as in the Setyawana-Curtarolo paper:
    if lattice == 'cubic':
        assert abc.ptp() < eps and abs(angles - pi / 2).max() < eps
    elif lattice == 'fcc':
        assert abc.ptp() < eps and abs(angles - pi / 3).max() < eps
    elif lattice == 'bcc':
        angle = np.arccos(-1 / 3)
        assert abc.ptp() < eps and abs(angles - angle).max() < eps
    elif lattice == 'tetragonal':
        assert abs(a - b) < eps and abs(angles - pi / 2).max() < eps
    elif lattice == 'orthorhombic':
        assert abs(angles - pi / 2).max() < eps
    elif lattice == 'hexagonal':
        assert abs(a - b) < eps
        assert abs(gamma - pi / 3 * 2) < eps
        assert abs(angles[:2] - pi / 2).max() < eps
    elif lattice == 'monoclinic':
        assert c >= a and c >= b
        assert alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps
    if lattice != 'monoclinic':
        return special_points[lattice]

    # Here, we need the cell:
    eta = (1 - b * np.cos(alpha) / c) / (2 * np.sin(alpha)**2)
    nu = 1 / 2 - eta * c * np.cos(alpha) / b
    return {'Γ': [0, 0, 0],
            'A': [1 / 2, 1 / 2, 0],
            'C': [0, 1 / 2, 1 / 2],
            'D': [1 / 2, 0, 1 / 2],
            'D1': [1 / 2, 0, -1 / 2],
            'E': [1 / 2, 1 / 2, 1 / 2],
            'H': [0, eta, 1 - nu],
            'H1': [0, 1 - eta, nu],
            'H2': [0, eta, -nu],
            'M': [1 / 2, eta, 1 - nu],
            'M1': [1 / 2, 1 - eta, nu],
            'M2': [1 / 2, eta, -nu],
            'X': [0, 1 / 2, 0],
            'Y': [0, 0, 1 / 2],
            'Y1': [0, 0, -1 / 2],
            'Z': [1 / 2, 0, 0]}


def findLabel(labels, value):
    for k, v in labels.items():
        if np.all(np.abs(v-value) < 1.e-5):
            return k
    return "?"


def secondsFromEpoch(date):
    epoch = datetime(1970, 1, 1)
    ts = date-epoch
    return ts.seconds + ts.microseconds/1000.0


trueRe = re.compile(
    r"\s*(?:\.?[Tt](?:[Rr][Uu][Ee])?\.?|1|[Yy](?:[Ee][Ss])?|[Jj][Aa]?)\s*$")
falseRe = re.compile(
    r"\s*(?:\.?[fF](?:[Aa][Ll][Ss][Ee])?\.?|0|[Nn](?:[Oo]|[Ee][Ii][Nn])?)\s*$")


def toBool(value):
    if falseRe.match(value):
        return False
    elif trueRe.match(value):
        return True
    else:
        backend.pwarn("Unexpected value for boolean field: %s" % (value))
        return None

metaTypeTransformers = {
    'C': lambda x: x.strip(),
    'i': lambda x: int(float(x.strip())),
    'f': lambda x: float(x.strip()),
    'b': toBool,
}


class MyXMLParser(ET.XMLParser):

    rx = re.compile("&#([0-9]+);|&#x([0-9a-fA-F]+);")

    def feed(self, data):
        m = self.rx.search(data)
        if m is not None:
            target = m.group(1)
            if target:
                num = int(target)
            else:
                num = int(m.group(2), 16)
            if not(num in (0x9, 0xA, 0xD) or 0x20 <= num <= 0xD7FF
                   or 0xE000 <= num <= 0xFFFD or 0x10000 <= num <= 0x10FFFF):
                # is invalid xml character, cut it out of the stream
                mstart, mend = m.span()
                mydata = data[:mstart] + data[mend:]
        else:
            mydata = data
        super(MyXMLParser, self).feed(mydata)


def transform2(y):
    if '**' in y:
        return float('nan')
    else:
        return y


def getVector(el, transform=float, field="v"):
    """ returns the vasp style vector contained in the element el (using field v).
    single elements are converted using the function convert"""
#
#    for x in el.findall(field):
#        for y in re.split(r"\s+", x.text.strip()):
    return [[transform(transform2(y)) for y in re.split(r"\s+", x.text.strip())] for x in el.findall(field)]


class VasprunContext(object):

    def __init__(self, logger=None):
        if logger is None:
            logger = logging.getLogger(__name__)
        self.logger = logger

        self.parser = None
        self.bands = None
        self.kpoints = None
        self.weights = None
        self.tetrahedrons = None
        self.tetrahedronVolume = None
        self.ispin = None
        self.ibrion = None
        self.lastSystemDescription = None
        self.labels = None
        self.singleConfCalcs = []
        self.vbTopE = None
        self.ebMinE = None
        self.eFermi = None
        self.cell = None
        self.angstrom_cell = None
        self.unknown_incars = {}

    sectionMap = {
        "modeling": ["section_run", "section_method"],
        "structure": ["section_system"],
        "calculation": ["section_single_configuration_calculation"]
    }

    def startedParsing(self, parser):
        self.parser = parser

    def onEnd_generator(self, parser, event, element, pathStr):
        backend = parser.backend
        program_name = g(element, "i/[@name='program']")
        if program_name.strip().upper() == "VASP":
            backend.addValue("program_name", "VASP")
        else:
            raise Exception("unexpected program name: %s" % program_name)
        version = (g(element, "i/[@name='version']", "") + " " +
                   g(element, "i/[@name='subversion']", "") + " " +
                   g(element, "i/[@name='platform']", ""))
        if not version.isspace():
            backend.addValue("program_version", version)
        backend.addValue("program_basis_set_type", "plane waves")
        date = g(element, "i/[@name='date']")
        pdate = None
        time = g(element, "i/[@name='time']")
        if date:
            pdate = datetime.strptime(date.strip(), "%Y %m %d")
        if pdate and time:
            pdate = datetime.combine(pdate.date(), datetime.strptime(
                time.strip(), "%H:%M:%S").timetz())
        if pdate:
            backend.addValue("program_compilation_datetime",
                             secondsFromEpoch(pdate))
        for i in element:
            if i.tag != "i" or not i.attrib.get("name") in set(["program", "version", "subversion", "platform", "program_version", "date", "time"]):
                backend.pwarn("unexpected tag %s %s %r in generator" %
                              (i.tag, i.attrib, i.text))

    def onEnd_incar(self, parser, event, element, pathStr):
        backend = parser.backend
        metaEnv = parser.backend.metaInfoEnv()
        dft_plus_u = False
        ibrion = None
        nsw = 0
        for el in element:
            if el.tag == "v":
                name = el.attrib.get("name", None)
                meta = metaEnv['x_vasp_incar_' + name]
                if not meta:
                    backend.pwarn("Unknown INCAR parameter (not registered in the meta data): %s %s %r" % (
                        el.tag, el.attrib, el.text))
                    continue
                #- -
                vector_val = np.asarray(getVector(el))
                backend.addArrayValues(meta.get('name'), vector_val)
            elif el.tag == "i":
                name = el.attrib.get("name", None)
                meta = metaEnv['x_vasp_incar_' + name]
                valType = el.attrib.get("type")
                if not meta:
                    backend.pwarn("Unknown INCAR parameter (not registered in the meta data): %s %s %r" % (
                        el.tag, el.attrib, el.text))
                elif valType:
                    expectedMetaType = {
                        'string': ['C'],
                        'int': ['i'],
                        'logical': ['b', 'C']
                    }.get(valType)
                    if not expectedMetaType:
                        backend.pwarn("Unknown value type %s encountered in INCAR: %s %s %r" % (
                            valType, el.tag, el.attrib, el.text))
                    elif not meta.get('dtypeStr') in expectedMetaType:
                        backend.pwarn("type mismatch between meta data %s and INCAR type %s for %s %s %r" % (
                            meta.get('dtypeStr'), valType, el.tag, el.attrib, el.text))
                    else:
                        shape = meta.get("shape", None)
                        dtypeStr = meta.get("dtypeStr", None)
                        converter = metaTypeTransformers.get(dtypeStr)
                        if not converter:
                            backend.pwarn(
                                "could not find converter for dtypeStr %s when handling meta info %s" % (dtypeStr, ))
                        elif shape:
                            vals = re.split("\s+", el.text.strip())
                            backend.addValue(
                                meta["name"], [converter(x) for x in vals])
                        else:
                            backend.addValue(meta["name"], converter(el.text))
                    if name == 'GGA':
                        # FIXME tmk: many options are not coded yet. See
                        # https://www.vasp.at/wiki/index.php/GGA
                        fMap = {
                            '91': ['GGA_X_PW91', 'GGA_C_PW91'],
                            'PE': ['GGA_X_PBE', 'GGA_C_PBE'],
                            'RP': ['GGA_X_RPBE', 'GGA_C_PBE'],
                            'PS': ['GGA_C_PBE_SOL', 'GGA_X_PBE_SOL'],
                            'MK': ['GGA_X_OPTB86_VDW']
                        }
                        functs = fMap.get(el.text.strip(), None)
                        if not functs:
                            backend.pwarn("Unknown XC functional %s" %
                                          el.text.strip())
                        else:
                            for f in functs:
                                backend.openNonOverlappingSection(
                                    "section_XC_functionals")
                                backend.addValue("XC_functional_name", f)
                                backend.closeNonOverlappingSection(
                                    "section_XC_functionals")
                    elif name == "ISPIN":
                        self.ispin = int(el.text.strip())
                    elif name == "LDAU":
                        if re.match(".?[Tt](?:[rR][uU][eE])?.?|[yY](?:[eE][sS])?|1", el.text.strip()):
                            dft_plus_u = True
                    elif name == "IBRION":
                        ibrion = int(el.text.strip())
                    elif name == "NSW":
                        nsw = int(el.text.strip())
            else:
                backend.pwarn("unexpected tag %s %s %r in incar" %
                              (el.tag, el.attrib, el.text))
        if ibrion is None:
            ibrion = -1 if nsw == 0 or nsw == 1 else 0
        if nsw == 0:
            ibrion = -1
        self.ibrion = ibrion
        if dft_plus_u:
            backend.addValue("electronic_structure_method", "DFT+U")
        else:
            backend.addValue("electronic_structure_method", "DFT")

    def onEnd_kpoints(self, parser, event, element, pathStr):
        backend = parser.backend
        self.bands = None
        self.kpoints = None
        self.weights = None
        for el in element:
            if el.tag == "generation":
                param = el.attrib.get("param", None) # eg. listgenerated, Monkhorst-Pack, Gamma
                if param:
                    backend.addValue(
                        "x_vasp_k_points_generation_method", param)
                if param == "listgenerated":
                    # This implies a path on k-space, potentially a bandstructure calculation
                    # Save k-path info into a dictionary
                    self.bands = {
                        "divisions": g(el, "i/[@name='divisions']", None),
                        "points": getVector(el)
                    }

                elif param in ["Monkhorst-Pack", "Gamma"]:
                    # This implies a (2D|3D) mesh on k-space, i.e., not a badstructure calculation
                    # Hence, do nothing: k-points will be stored in the `varray` if-block
                    pass
                else:
                    backend.pwarn("Unknown k point generation method '%s'" %(param))
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "kpointlist":
                    self.kpoints = np.asarray(getVector(el))
                    backend.addArrayValues("k_mesh_points", self.kpoints)
                elif name == "weights":
                    self.weights = np.asarray(getVector(el))
                    backend.addArrayValues(
                        "k_mesh_weights", self.weights.flatten())
                elif name == "tetrahedronlist":
                    self.tetrahedrons = np.asarray(getVector(el), dtype=np.int)
                    backend.addArrayValues(
                        "x_vasp_tetrahedrons_list", self.tetrahedrons)
                else:
                    backend.pwarn("Unknown array %s in kpoints" % name)
            elif el.tag == "i":
                name = el.attrib.get("name", None)
                if name == "volumeweight":
                    ang2m = convert_unit_function("angstrom", "m")

                    # get volume and transform to meters^3
                    vol_cubic_angs = float(el.text.strip())
                    vol_cubic_meters = ang2m(ang2m(ang2m(vol_cubic_angs)))

                    backend.addArrayValues("x_vasp_tetrahedron_volume",
                        vol_cubic_meters)
            else:
                backend.pwarn("Unknown tag %s in kpoints" % el.tag)


    def onEnd_structure(self, parser, event, element, pathStr):
        backend = parser.backend
        gIndexes = parser.tagSections[pathStr]
        self.lastSystemDescription = gIndexes["section_system"]
        self.cell = None
        for el in element:
            if (el.tag == "crystal"):
                for cellEl in el:
                    if cellEl.tag == "varray":
                        name = cellEl.attrib.get("name", None)
                        if name == "basis":
                            conv = convert_unit_function("angstrom", "m")
                            self.cell = getVector(
                                cellEl, lambda x: conv(float(x)))
                            self.angstrom_cell = np.array(getVector(cellEl))
                            backend.addArrayValues(
                                "simulation_cell", np.asarray(self.cell))
                            backend.addArrayValues(
                                "configuration_periodic_dimensions", np.ones(3, dtype=bool))
                        elif name == "rec_basis":
                            pass
                        else:
                            backend.pwarn(
                                "Unexpected varray %s in crystal" % name)
                    elif cellEl.tag == "i":
                        if cellEl.attrib.get("name") != "volume":
                            backend.pwarn(
                                "Unexpected i value %s in crystal" % cellEl.attrib)
                    else:
                        backend.pwarn("Unexpected tag %s %s %r in crystal" % (
                            cellEl.tag, cellEl.attrib, cellEl.text))
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "positions":
                    pos = getVector(el)
                    backend.addArrayValues(
                        "atom_positions", np.dot(np.asarray(pos), self.cell))
                elif name == "selective":
                    atom_sel = getVector(el, transform=lambda item: item == 'T')
                    backend.addArrayValues(
                        "x_vasp_selective_dynamics", np.asarray(atom_sel, dtype=np.bool))
                else:
                    backend.pwarn(
                        "Unexpected varray in structure %s" % el.attrib)
            elif el.tag == "nose":
                nose = getVector(el)
                backend.addArrayValues("x_vasp_nose_thermostat", nose)
            else:
                backend.pwarn("Unexpected tag in structure %s %s %r" %
                              (el.tag, el.attrib, el.text))
        if self.labels is not None:
            backend.addArrayValues("atom_labels", self.labels)

    def onEnd_eigenvalues(self, parser, event, element, pathStr):
        if pathStr != "modeling/calculation/eigenvalues":
            return True
        backend = parser.backend
        eigenvalues = None
        occupation = None
        for el in element:
            if el.tag == "array":
                for arrEl in el:
                    if arrEl.tag == "dimension":
                        pass
                    elif arrEl.tag == "field":
                        pass
                    elif arrEl.tag == "set":
                        isp = -1
                        for spinEl in arrEl:
                            if spinEl.tag == "set":
                                ik = -1
                                isp += 1
                                for kEl in spinEl:
                                    if kEl.tag == "set":
                                        ik += 1
                                        bands = np.asarray(
                                            getVector(kEl, field="r"))
                                        if eigenvalues is None:
                                            eigenvalues = np.zeros(
                                                (self.ispin, self.kpoints.shape[0],  bands.shape[0]), dtype=float)
                                            occupation = np.zeros(
                                                (self.ispin, self.kpoints.shape[0],  bands.shape[0]), dtype=float)
                                        eigenvalues[isp, ik] = bands[:, 0]
                                        occupation[isp, ik] = bands[:, 1]
                                    else:
                                        backend.pwarn(
                                            "unexpected tag %s in k array of the eigenvalues" % kEl.tag)
                            else:
                                backend.pwarn(
                                    "unexpected tag %s in spin array of the eigenvalues" % spinEl.tag)
                    else:
                        backend.pwarn(
                            "unexpected tag %s in array of the eigenvalues" % arrEl.tag)
                if eigenvalues is not None:

                    ev = eV2JV(eigenvalues)
                    vbTopE = []
                    ebMinE = []
                    for ispin in range(occupation.shape[0]):
                        vbTopE.append(float('-inf'))
                        ebMinE.append(float('inf'))
                        for ik in range(occupation.shape[1]):
                            ebIndex = bisect.bisect_right(
                                -occupation[ispin, ik, :], -0.5) - 1
                            vbTopIndex = ebIndex - 1
                            if vbTopIndex >= 0:
                                vbTopK = ev[ispin, ik, vbTopIndex]
                                if vbTopK > vbTopE[ispin]:
                                    vbTopE[ispin] = vbTopK
                            if ebIndex < ev.shape[2]:
                                ebMinK = ev[ispin, ik, ebIndex]
                                if ebMinK < ebMinE[ispin]:
                                    ebMinE[ispin] = ebMinK
                    self.vbTopE = vbTopE
                    self.ebMinE = ebMinE
                    backend.addArrayValues(
                        "energy_reference_highest_occupied", np.array(vbTopE))
                    backend.addArrayValues(
                        "energy_reference_lowest_unoccupied", np.array(ebMinE))
                    if self.bands:
                        divisions = int(self.bands['divisions'])
                        backend.openNonOverlappingSection("section_k_band")
                        nsegments = self.kpoints.shape[0] // divisions
                        kpt = np.reshape(
                            self.kpoints, (nsegments, divisions, 3))
                        energies = np.reshape(
                            ev, (self.ispin, nsegments, divisions,  bands.shape[0]))
                        occ = np.reshape(
                            occupation, (self.ispin, nsegments, divisions, bands.shape[0]))
                        for isegment in range(nsegments):
                            backend.openNonOverlappingSection(
                                "section_k_band_segment")
                            backend.addArrayValues(
                                "band_energies", energies[:, isegment, :, :])
                            backend.addArrayValues(
                                "band_occupations", occ[:, isegment, :, :])
                            backend.addArrayValues(
                                "band_k_points", kpt[isegment])
                            # "band_segm_labels"
                            backend.addArrayValues("band_segm_start_end", np.asarray(
                                [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment")
                        backend.closeNonOverlappingSection("section_k_band")
                        backend.openNonOverlappingSection(
                            "section_k_band_normalized")
                        specialPoints = {}
                        try:
                            specialPoints = get_special_points(
                                convert_unit_function("m", "angstrom")(self.cell))
                        except Exception as e:
                            self.logger.warn("failed to get special points/xtal structure", exc_info=e)
                        for isegment in range(nsegments):
                            backend.openNonOverlappingSection(
                                "section_k_band_segment_normalized")
                            backend.addArrayValues(
                                "band_energies_normalized", energies[:, isegment, :, :] - max(self.vbTopE))
                            backend.addArrayValues(
                                "band_occupations_normalized", occ[:, isegment, :, :])
                            backend.addArrayValues(
                                "band_k_points_normalized", kpt[isegment])
                            backend.addArrayValues("band_segm_start_end_normalized", np.asarray(
                                [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
                            backend.addValue("band_segm_labels_normalized",
                                             [findLabel(specialPoints, kpt[isegment, 0]),
                                              findLabel(specialPoints, kpt[isegment, divisions - 1])])
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment_normalized")

                        backend.closeNonOverlappingSection(
                            "section_k_band_normalized")
                    else:
                        backend.openNonOverlappingSection(
                            "section_eigenvalues")
                        backend.addArrayValues("eigenvalues_values", ev)
                        backend.addArrayValues(
                            "eigenvalues_occupation", occupation)
                        backend.closeNonOverlappingSection(
                            "section_eigenvalues")
            else:
                backend.pwarn("unexpected tag %s in the eigenvalues" % el.tag)

    def onEnd_scstep(self, parser, event, element, pathStr):
        pass

    def onStart_calculation(self, parser, event, element, pathStr):
        backend = parser.backend
        gIndexes = parser.tagSections[pathStr]
        self.singleConfCalcs.append(
            gIndexes["section_single_configuration_calculation"])
        if self.waveCut:
            backend.openNonOverlappingSection("section_basis_set")
            backend.addValue(
                "mapping_section_basis_set_cell_dependent", self.waveCut)
            backend.closeNonOverlappingSection("section_basis_set")

    def onEnd_modeling(self, parser, event, element, pathStr):
        backend = parser.backend
        backend.addValue("x_vasp_unknown_incars", self.unknown_incars)
        if self.ibrion is None or self.ibrion == -1:
            return
        samplingGIndex = backend.openSection("section_sampling_method")
        if self.ibrion == 0:
            sampling_method = "molecular_dynamics"
        else:
            sampling_method = "geometry_optimization"
        backend.addValue("sampling_method", sampling_method)
        backend.closeSection("section_sampling_method", samplingGIndex)
        frameSequenceGIndex = backend.openSection("section_frame_sequence")
        backend.addValue("frame_sequence_to_sampling_ref", samplingGIndex)
        backend.addArrayValues(
            "frame_sequence_local_frames_ref", np.asarray(self.singleConfCalcs))
        backend.closeSection("section_frame_sequence", frameSequenceGIndex)


    def onEnd_calculation(self, parser, event, element, pathStr):
        eConv = eV2J
        fConv = convert_unit_function("eV/angstrom", "N")
        pConv = convert_unit_function("eV/angstrom^3", "Pa")
        backend = parser.backend
        backend.addValue(
            "single_configuration_calculation_to_system_ref", self.lastSystemDescription)
        gIndexes = parser.tagSections["/modeling"]
        backend.addValue(
            "single_configuration_to_calculation_method_ref", gIndexes["section_method"])
        for el in element:
            if el.tag == "energy":
                for enEl in el:
                    if enEl.tag == "i":
                        name = enEl.attrib.get("name", None)
                        if name == "e_fr_energy":
                            value = eConv(float(enEl.text.strip()))
                            backend.addValue("energy_free", value)
                        elif name == "e_wo_entrp":
                            value = eConv(float(enEl.text.strip()))
                            backend.addValue("energy_total", value)
                        elif name == "e_0_energy":
                            value = eConv(float(enEl.text.strip()))
                            backend.addValue("energy_total_T0", value)
                        else:
                            backend.pwarn(
                                "Unexpected i tag with name %s in energy section" % name)
                    elif enEl.tag == "varray":
                        name = enEl.attrib.get("name", None)
                        if name == "forces":
                            f = getVector(enEl, lambda x: fConv(float(x)))
                            backend.addValue("atom_forces", f)
                        elif name == 'stress':
                            f = getVector(enEl, lambda x: pConv(float(x)))
                            backend.addValue("stress_tensor", f)

    def onEnd_atominfo(self, parser, event, element, pathStr):
        nAtoms = None
        nAtomTypes = None
        atomTypes = []
        labels = []
        labels2 = None
        atomTypesDesc = []
        backend = parser.backend
        for el in element:
            if el.tag == "atoms":
                nAtoms = int(el.text.strip())
            elif el.tag == "types":
                nAtomTypes = int(el.text.strip())
            elif el.tag == "array":
                name = el.attrib.get("name", None)
                if name == "atoms":
                    for atomsEl in el:
                        if atomsEl.tag == "dimension":
                            pass
                        elif atomsEl.tag == "field":
                            pass
                        elif atomsEl.tag == "set":
                            for atomsLine in atomsEl:
                                if atomsLine.tag != "rc":
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
                                else:
                                    line = atomsLine.findall("c")
                                    labels.append(line[0].text.strip())
                                    atomTypes.append(int(line[1].text.strip()))
                        else:
                            backend.pwarn(
                                "unexpected tag %s in atoms array in atominfo" % atomsEl.tag)
                elif name == "atomtypes":
                    keys = []
                    fieldTypes = []
                    for atomsEl in el:
                        if atomsEl.tag == "dimension":
                            pass
                        elif atomsEl.tag == "field":
                            keys.append(atomsEl.text.strip())
                            fieldTypes.append(
                                atomsEl.attrib.get("type", "float"))
                        elif atomsEl.tag == "set":
                            expectedKeys = ["atomspertype", "element",
                                            "mass", "valence", "pseudopotential"]
                            if keys != expectedKeys:
                                backend.pwarn(
                                    "unexpected fields in atomtype: %s vs %s" % (keys, expectedKeys))
                            for atomsLine in atomsEl:
                                if atomsLine.tag != "rc":
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
                                else:
                                    line = atomsLine.findall("c")
                                    typeDesc = {}
                                    for i, k in enumerate(keys):
                                        fieldType = fieldTypes[i]
                                        value = line[i].text
                                        if fieldType == "float":
                                            value = float(value)
                                        elif fieldType == "int":
                                            value = int(value)
                                        else:
                                            pass
                                        typeDesc[k] = value
                                    atomTypesDesc.append(typeDesc)
                        else:
                            backend.pwarn(
                                "unexpected tag %s in atomtypes array in atominfo" % atomsEl.tag)
                    kindIds = []
                    nEl = {}
                    kindLabels = []
                    for atomDesc in atomTypesDesc:
                        kindId = backend.openSection(
                            "section_method_atom_kind")
                        if 'element' in atomDesc:
                            elName = atomDesc['element'].strip()
                            try:
                                elNr = ase.data.chemical_symbols.index(elName)
                                backend.addValue(
                                    "method_atom_kind_atom_number", elNr)
                            except Exception as e:
                                self.logger.error(
                                    "error finding element number for %r" % atomDesc['element'].strip(),
                                    exc_info=e)
                            nElNow = 1 + nEl.get(elName, 0)
                            nEl[elName] = nElNow
                            elLabel = elName + \
                                (str(nElNow) if nElNow > 1 else "")
                            kindLabels.append(elLabel)
                            backend.addValue("method_atom_kind_label", elLabel)
                            if "mass" in atomDesc:
                                backend.addValue(
                                    "method_atom_kind_mass", atomDesc["mass"])
                            if "valence" in atomDesc:
                                backend.addValue(
                                    "method_atom_kind_explicit_electrons", atomDesc["valence"])
                            if "pseudopotential" in atomDesc:
                                backend.addValue(
                                    "method_atom_kind_pseudopotential_name", atomDesc["pseudopotential"])
                        kindIds.append(kindId)
                        backend.closeSection(
                            "section_method_atom_kind", kindId)
                    backend.addArrayValues("x_vasp_atom_kind_refs", np.asarray(
                        [kindIds[i-1] for i in atomTypes]))
                    labels2 = [kindLabels[i-1] for i in atomTypes]
                else:
                    backend.pwarn(
                        "unexpected array named %s in atominfo" % name)
            else:
                backend.pwarn("unexpected tag %s in atominfo" % el.tag)
        self.labels = np.asarray(labels2) if labels2 else np.asarray(labels)

    def incarOutTag(self, el):
        backend = self.parser.backend
        metaEnv = self.parser.backend.metaInfoEnv()
        if (el.tag != "i"):
            backend.pwarn("unexpected tag %s %s %r in incar" %
                          (el.tag, el.attrib, el.text))
        else:
            name    = el.attrib.get("name", None)
            valType = el.attrib.get("type")
            meta = metaEnv['x_vasp_incarOut_' + name]


            if not meta:
                # Unknown_Incars_Begin: storage into a dictionary
                if not valType:
                    # On vasp's xml files, valType *could* be absent if incar value is float
                    valType = 'float'

                # map vasp's datatype to nomad's datatype [b, f, i, C, D, R]
                nomad_dtypeStr = vasp_to_metainfo_type_mapping[valType][0]

                converter = metaTypeTransformers.get(nomad_dtypeStr)
                text_value = el.text.strip() # text representation of incar value
                try:
                    pyvalue = converter(text_value) # python data type
                except Exception:
                    pyvalue = text_value

                # save (name, pyvalue) into a dict
                self.unknown_incars[name] = pyvalue
                # Unknown_Incars_end
            else:
                if not valType:
                    valType = 'float'

                vasp_metainfo_type = vasp_to_metainfo_type_mapping.get(valType)[0]
                metainfo_type = meta.get('dtypeStr')
                if not vasp_metainfo_type:
                    backend.pwarn("Unknown value type %s encountered in INCAR out: %s %s %r" % (
                        valType, el.tag, el.attrib, el.text))

                elif metainfo_type != vasp_metainfo_type:
                    if  (metainfo_type == 'C' and vasp_metainfo_type == 'b'):
                        pass
                    elif  (metainfo_type == 'i' and vasp_metainfo_type == 'f'):
                        pass
                    else:
                        backend.pwarn("Data type mismatch: %s. Vasp_type: %s, metainfo_type: %s " %
                        (name, vasp_metainfo_type, metainfo_type))
                try:
                    shape = meta.get("shape", None)
                    converter = metaTypeTransformers.get(metainfo_type)
                    if not converter:
                        backend.pwarn(
                            "could not find converter for dtypeStr %s when handling meta info %s" %
                            (metainfo_type, meta ))
                    elif shape:
                        vals = re.split("\s+", el.text.strip())
                        backend.addValue(
                            meta["name"], [converter(x) for x in vals])
                    else:
                        # If-block to handle incars without value
                        if el.text == None:
                            el.text = ''
                        backend.addValue(meta["name"], converter(el.text))

                except:
                    backend.pwarn("Exception trying to handle incarOut %s: %s" % (
                        name, traceback.format_exc()))

                if name == 'ENMAX' or name == 'PREC':
                    if name == 'ENMAX':
                        self.enmax = converter(el.text)
                    if name == 'PREC':
                        if 'acc' in converter(el.text):
                            self.prec = 1.3
                        else:
                            self.prec = 1.0
                if name == 'GGA':
                    fMap = {
                        '91': ['GGA_X_PW91', 'GGA_C_PW91'],
                        'PE': ['GGA_X_PBE', 'GGA_C_PBE'],
                        'RP': ['GGA_X_RPBE', 'GGA_C_PBE'],
                        'PS': ['GGA_C_PBE_SOL', 'GGA_X_PBE_SOL'],
                        'MK': ['GGA_X_OPTB86_VDW'],
                        '--': ['GGA_X_PBE', 'GGA_C_PBE']  # should check potcar
                    }
                    functs = fMap.get(el.text.strip(), None)
                    if not functs:
                        backend.pwarn("Unknown XC functional %s" %
                                      el.text.strip())
                    else:
                        for f in functs:
                            backend.openNonOverlappingSection(
                                "section_XC_functionals")
                            backend.addValue("XC_functional_name", f)
                            backend.closeNonOverlappingSection(
                                "section_XC_functionals")
                elif name == "ISPIN":
                    self.ispin = int(el.text.strip())


    def separatorScan(self, element, backend, depth=0):
        for separators in element:
            if separators.tag == "separator":
                separatorName = separators.attrib.get("name")
                for el in separators:
                    if el.tag == "i":
                        self.incarOutTag(el)
                    elif el.tag == "separator":</