parser_vasprun.py 56.6 KB
Newer Older
1
# Copyright 2016-2018 Fawzi Mohamed, Lauri Himanen, Danio Brambila, Ankit Kariryaa, Henning Glawe
Markus Scheidgen's avatar
Markus Scheidgen committed
2
#
3
4
5
#   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
Markus Scheidgen's avatar
Markus Scheidgen committed
6
#
7
#     http://www.apache.org/licenses/LICENSE-2.0
Markus Scheidgen's avatar
Markus Scheidgen committed
8
#
9
10
11
12
13
14
#   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.

15
16
17
from __future__ import division
from builtins import range
from builtins import object
18
import xml.etree.ElementTree
19
from xml.etree.ElementTree import ParseError
Markus Scheidgen's avatar
Markus Scheidgen committed
20
21
import sys
import bisect
22
from datetime import datetime
Markus Scheidgen's avatar
Markus Scheidgen committed
23
24
25
import os
import re
import traceback
26
import numpy as np
27
import ase.geometry
28
import ase.data
29
from math import pi
Markus Scheidgen's avatar
Markus Scheidgen committed
30
31
32
33
34
35
36
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
37

38
eV2J = convert_unit_function("eV", "J")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
39
40
eV2JV = np.vectorize(eV2J)

Markus Scheidgen's avatar
Markus Scheidgen committed
41

42

43
44
45
46
47
48
49
50
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
51
52
53
    # 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.
54
55
56
57
58
59
    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'
60
    # If a=b!=c, alpha=beta=gamma=90deg, tetragonal.
61
62
63
64
    elif abs(a - b) < eps and abs(angles - pi / 2).max() < eps:
        return 'tetragonal'
    elif abs(angles - pi / 2).max() < eps:
        return 'orthorhombic'
65
66
67
68
    # 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):
69
70
71
72
73
        return 'hexagonal'
    elif (c >= a and c >= b and alpha < pi / 2 and
          abs(angles[1:] - pi / 2).max() < eps):
        return 'monoclinic'
    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
74
        raise ValueError('Cannot find crystal structure')
75

76
77
78
79
80
vasp_to_metainfo_type_mapping = {
    'string': ['C'],
    'int': ['i'],
    'logical': ['b', 'C'],
    'float': ['f']}
81
82

special_points = {
83
    'cubic': {'Γ': [0, 0, 0],
84
85
86
              'M': [1 / 2, 1 / 2, 0],
              'R': [1 / 2, 1 / 2, 1 / 2],
              'X': [0, 1 / 2, 0]},
87
    'fcc': {'Γ': [0, 0, 0],
88
89
90
91
92
            '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]},
93
    'bcc': {'Γ': [0, 0, 0],
94
95
96
            'H': [1 / 2, -1 / 2, 1 / 2],
            'P': [1 / 4, 1 / 4, 1 / 4],
            'N': [0, 0, 1 / 2]},
97
    'tetragonal': {'Γ': [0, 0, 0],
98
99
100
101
102
                   '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]},
103
    'orthorhombic': {'Γ': [0, 0, 0],
104
105
106
107
108
109
110
                     '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]},
111
    'hexagonal': {'Γ': [0, 0, 0],
112
113
114
115
116
117
118
119
                  '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 = {
120
121
122
123
124
125
126
    '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'}
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176


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:
Daniel Speckhard's avatar
Daniel Speckhard committed
177
178
    eta = (1 - b * np.cos(alpha) / c) / (2 * np.sin(alpha)**2)
    nu = 1 / 2 - eta * c * np.cos(alpha) / b
179
    return {'Γ': [0, 0, 0],
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
            '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]}

Markus Scheidgen's avatar
Markus Scheidgen committed
196

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

Markus Scheidgen's avatar
Markus Scheidgen committed
203

204
def secondsFromEpoch(date):
Markus Scheidgen's avatar
Markus Scheidgen committed
205
206
    epoch = datetime(1970, 1, 1)
    ts = date-epoch
207
208
    return ts.seconds + ts.microseconds/1000.0

Markus Scheidgen's avatar
Markus Scheidgen committed
209
210
211
212
213
214
215

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*$")


216
217
218
219
220
221
222
223
224
225
226
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(),
227
    'i': lambda x: int(float(x.strip())),
228
229
230
231
    'f': lambda x: float(x.strip()),
    'b': toBool,
}

232
233
234
235
236

class MyXMLParser(ET.XMLParser):

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

Markus Scheidgen's avatar
Markus Scheidgen committed
237
    def feed(self, data):
238
239
240
241
242
243
244
245
246
247
248
249
250
251
        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
Markus Scheidgen's avatar
Markus Scheidgen committed
252
253
254
        super(MyXMLParser, self).feed(mydata)


255
def transform2(y):
Markus Scheidgen's avatar
Markus Scheidgen committed
256
257
258
259
260
    if '**' in y:
        return float('nan')
    else:
        return y

261

Markus Scheidgen's avatar
Markus Scheidgen committed
262
def getVector(el, transform=float, field="v"):
263
264
    """ returns the vasp style vector contained in the element el (using field v).
    single elements are converted using the function convert"""
265
266
267
268
#
#    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)]
269

Markus Scheidgen's avatar
Markus Scheidgen committed
270

271
class VasprunContext(object):
272

273
    def __init__(self, logger=None):
274
        if logger is None:
275
            logger = logging.getLogger(__name__)
Markus Scheidgen's avatar
Markus Scheidgen committed
276
        self.logger = logger
277

278
279
280
281
        self.parser = None
        self.bands = None
        self.kpoints = None
        self.weights = None
282
283
        self.tetrahedrons = None
        self.tetrahedronVolume = None
284
        self.ispin = None
285
        self.ibrion = None
286
        self.lastSystemDescription = None
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
287
        self.labels = None
288
        self.singleConfCalcs = []
289
290
291
292
        self.vbTopE = None
        self.ebMinE = None
        self.eFermi = None
        self.cell = None
293
        self.angstrom_cell = None
294
295
        self.unknown_incars = {}

296
    sectionMap = {
297
        "modeling": ["section_run", "section_method"],
298
        "structure": ["section_system"],
299
        "calculation": ["section_single_configuration_calculation"]
300
301
    }

302
303
304
    def startedParsing(self, parser):
        self.parser = parser

305
    def onEnd_generator(self, parser, event, element, pathStr):
306
        backend = parser.backend
307
        program_name = g(element, "i/[@name='program']")
308
309
310
311
        if program_name.strip().upper() == "VASP":
            backend.addValue("program_name", "VASP")
        else:
            raise Exception("unexpected program name: %s" % program_name)
312
313
314
315
316
        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)
317
        backend.addValue("program_basis_set_type", "plane waves")
318
319
320
321
        date = g(element, "i/[@name='date']")
        pdate = None
        time = g(element, "i/[@name='time']")
        if date:
322
            pdate = datetime.strptime(date.strip(), "%Y %m %d")
323
        if pdate and time:
Markus Scheidgen's avatar
Markus Scheidgen committed
324
325
            pdate = datetime.combine(pdate.date(), datetime.strptime(
                time.strip(), "%H:%M:%S").timetz())
326
        if pdate:
Markus Scheidgen's avatar
Markus Scheidgen committed
327
328
            backend.addValue("program_compilation_datetime",
                             secondsFromEpoch(pdate))
329
        for i in element:
330
            if i.tag != "i" or not i.attrib.get("name") in set(["program", "version", "subversion", "platform", "program_version", "date", "time"]):
Markus Scheidgen's avatar
Markus Scheidgen committed
331
332
                backend.pwarn("unexpected tag %s %s %r in generator" %
                              (i.tag, i.attrib, i.text))
333

334
    def onEnd_incar(self, parser, event, element, pathStr):
335
        backend = parser.backend
336
        metaEnv = parser.backend.metaInfoEnv()
337
        dft_plus_u = False
338
339
        ibrion = None
        nsw = 0
340
341
        for el in element:
            if el.tag == "v":
342
343
344
345
346
                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))
347
                    continue
348
                #- -
349
350
351
                vector_val = np.asarray(getVector(el))
                backend.addArrayValues(meta.get('name'), vector_val)
            elif el.tag == "i":
352
353
354
355
                name = el.attrib.get("name", None)
                meta = metaEnv['x_vasp_incar_' + name]
                valType = el.attrib.get("type")
                if not meta:
Markus Scheidgen's avatar
Markus Scheidgen committed
356
357
                    backend.pwarn("Unknown INCAR parameter (not registered in the meta data): %s %s %r" % (
                        el.tag, el.attrib, el.text))
358
359
360
361
                elif valType:
                    expectedMetaType = {
                        'string': ['C'],
                        'int': ['i'],
Markus Scheidgen's avatar
Markus Scheidgen committed
362
                        'logical': ['b', 'C']
363
364
                    }.get(valType)
                    if not expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
365
366
                        backend.pwarn("Unknown value type %s encountered in INCAR: %s %s %r" % (
                            valType, el.tag, el.attrib, el.text))
367
                    elif not meta.get('dtypeStr') in expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
368
369
                        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))
370
371
372
373
374
                    else:
                        shape = meta.get("shape", None)
                        dtypeStr = meta.get("dtypeStr", None)
                        converter = metaTypeTransformers.get(dtypeStr)
                        if not converter:
Markus Scheidgen's avatar
Markus Scheidgen committed
375
376
                            backend.pwarn(
                                "could not find converter for dtypeStr %s when handling meta info %s" % (dtypeStr, ))
377
378
                        elif shape:
                            vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
379
380
                            backend.addValue(
                                meta["name"], [converter(x) for x in vals])
381
382
383
                        else:
                            backend.addValue(meta["name"], converter(el.text))
                    if name == 'GGA':
384
385
                        # FIXME tmk: many options are not coded yet. See
                        # https://www.vasp.at/wiki/index.php/GGA
386
387
388
389
                        fMap = {
                            '91': ['GGA_X_PW91', 'GGA_C_PW91'],
                            'PE': ['GGA_X_PBE', 'GGA_C_PBE'],
                            'RP': ['GGA_X_RPBE', 'GGA_C_PBE'],
390
391
                            'PS': ['GGA_C_PBE_SOL', 'GGA_X_PBE_SOL'],
                            'MK': ['GGA_X_OPTB86_VDW']
392
393
394
                        }
                        functs = fMap.get(el.text.strip(), None)
                        if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
395
396
                            backend.pwarn("Unknown XC functional %s" %
                                          el.text.strip())
397
398
                        else:
                            for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
399
400
                                backend.openNonOverlappingSection(
                                    "section_XC_functionals")
401
                                backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
402
403
                                backend.closeNonOverlappingSection(
                                    "section_XC_functionals")
404
405
                    elif name == "ISPIN":
                        self.ispin = int(el.text.strip())
406
407
408
                    elif name == "LDAU":
                        if re.match(".?[Tt](?:[rR][uU][eE])?.?|[yY](?:[eE][sS])?|1", el.text.strip()):
                            dft_plus_u = True
409
410
411
412
                    elif name == "IBRION":
                        ibrion = int(el.text.strip())
                    elif name == "NSW":
                        nsw = int(el.text.strip())
413
414
415
            else:
                backend.pwarn("unexpected tag %s %s %r in incar" %
                              (el.tag, el.attrib, el.text))
416
417
        if ibrion is None:
            ibrion = -1 if nsw == 0 or nsw == 1 else 0
418
419
        if nsw == 0:
            ibrion = -1
420
        self.ibrion = ibrion
421
422
423
424
        if dft_plus_u:
            backend.addValue("electronic_structure_method", "DFT+U")
        else:
            backend.addValue("electronic_structure_method", "DFT")
425

426
    def onEnd_kpoints(self, parser, event, element, pathStr):
427
        backend = parser.backend
428
429
        self.bands = None
        self.kpoints = None
430
431
        self.weights = None
        for el in element:
432
            if el.tag == "generation":
433
                param = el.attrib.get("param", None) # eg. listgenerated, Monkhorst-Pack, Gamma
434
                if param:
Markus Scheidgen's avatar
Markus Scheidgen committed
435
436
                    backend.addValue(
                        "x_vasp_k_points_generation_method", param)
437
                if param == "listgenerated":
438
                    # This implies a path on k-space, potentially a bandstructure calculation
439
                    # Save k-path info into a dictionary
440
                    self.bands = {
Markus Scheidgen's avatar
Markus Scheidgen committed
441
                        "divisions": g(el, "i/[@name='divisions']", None),
442
443
                        "points": getVector(el)
                    }
444
445
446
447

                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
448
449
                    pass
                else:
450
                    backend.pwarn("Unknown k point generation method '%s'" %(param))
451
452
453
454
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "kpointlist":
                    self.kpoints = np.asarray(getVector(el))
455
                    backend.addArrayValues("k_mesh_points", self.kpoints)
456
457
                elif name == "weights":
                    self.weights = np.asarray(getVector(el))
Markus Scheidgen's avatar
Markus Scheidgen committed
458
459
                    backend.addArrayValues(
                        "k_mesh_weights", self.weights.flatten())
460
461
462
463
                elif name == "tetrahedronlist":
                    self.tetrahedrons = np.asarray(getVector(el), dtype=np.int)
                    backend.addArrayValues(
                        "x_vasp_tetrahedrons_list", self.tetrahedrons)
464
465
                else:
                    backend.pwarn("Unknown array %s in kpoints" % name)
466
467
468
469
470
471
472
473
474
475
476
            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)
477
            else:
478
                backend.pwarn("Unknown tag %s in kpoints" % el.tag)
479

480

481
    def onEnd_structure(self, parser, event, element, pathStr):
482
        backend = parser.backend
483
        gIndexes = parser.tagSections[pathStr]
484
        self.lastSystemDescription = gIndexes["section_system"]
485
        self.cell = None
486
487
488
489
490
491
        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":
492
                            conv = convert_unit_function("angstrom", "m")
Markus Scheidgen's avatar
Markus Scheidgen committed
493
494
                            self.cell = getVector(
                                cellEl, lambda x: conv(float(x)))
495
                            self.angstrom_cell = np.array(getVector(cellEl))
Markus Scheidgen's avatar
Markus Scheidgen committed
496
497
498
499
                            backend.addArrayValues(
                                "simulation_cell", np.asarray(self.cell))
                            backend.addArrayValues(
                                "configuration_periodic_dimensions", np.ones(3, dtype=bool))
500
                        elif name == "rec_basis":
501
502
                            pass
                        else:
Markus Scheidgen's avatar
Markus Scheidgen committed
503
504
                            backend.pwarn(
                                "Unexpected varray %s in crystal" % name)
505
506
                    elif cellEl.tag == "i":
                        if cellEl.attrib.get("name") != "volume":
Markus Scheidgen's avatar
Markus Scheidgen committed
507
508
                            backend.pwarn(
                                "Unexpected i value %s in crystal" % cellEl.attrib)
509
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
510
511
                        backend.pwarn("Unexpected tag %s %s %r in crystal" % (
                            cellEl.tag, cellEl.attrib, cellEl.text))
512
513
514
515
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "positions":
                    pos = getVector(el)
Markus Scheidgen's avatar
Markus Scheidgen committed
516
517
                    backend.addArrayValues(
                        "atom_positions", np.dot(np.asarray(pos), self.cell))
518
519
520
521
                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))
522
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
523
524
                    backend.pwarn(
                        "Unexpected varray in structure %s" % el.attrib)
525
526
527
            elif el.tag == "nose":
                nose = getVector(el)
                backend.addArrayValues("x_vasp_nose_thermostat", nose)
528
            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
529
                backend.pwarn("Unexpected tag in structure %s %s %r" %
530
                              (el.tag, el.attrib, el.text))
531
        if self.labels is not None:
532
            backend.addArrayValues("atom_labels", self.labels)
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
533

534
    def onEnd_eigenvalues(self, parser, event, element, pathStr):
535
536
        if pathStr != "modeling/calculation/eigenvalues":
            return True
537
538
539
        backend = parser.backend
        eigenvalues = None
        occupation = None
540
541
542
543
544
545
546
547
        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":
548
                        isp = -1
549
550
                        for spinEl in arrEl:
                            if spinEl.tag == "set":
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
551
                                ik = -1
552
                                isp += 1
553
554
                                for kEl in spinEl:
                                    if kEl.tag == "set":
555
                                        ik += 1
Markus Scheidgen's avatar
Markus Scheidgen committed
556
557
                                        bands = np.asarray(
                                            getVector(kEl, field="r"))
558
                                        if eigenvalues is None:
Markus Scheidgen's avatar
Markus Scheidgen committed
559
560
561
562
563
564
                                            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]
565
                                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
566
567
                                        backend.pwarn(
                                            "unexpected tag %s in k array of the eigenvalues" % kEl.tag)
568
                            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
569
570
                                backend.pwarn(
                                    "unexpected tag %s in spin array of the eigenvalues" % spinEl.tag)
571
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
572
573
                        backend.pwarn(
                            "unexpected tag %s in array of the eigenvalues" % arrEl.tag)
574
                if eigenvalues is not None:
575

Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
576
                    ev = eV2JV(eigenvalues)
577
578
                    vbTopE = []
                    ebMinE = []
579
                    for ispin in range(occupation.shape[0]):
580
581
                        vbTopE.append(float('-inf'))
                        ebMinE.append(float('inf'))
582
                        for ik in range(occupation.shape[1]):
Markus Scheidgen's avatar
Markus Scheidgen committed
583
584
585
                            ebIndex = bisect.bisect_right(
                                -occupation[ispin, ik, :], -0.5) - 1
                            vbTopIndex = ebIndex - 1
586
587
                            if vbTopIndex >= 0:
                                vbTopK = ev[ispin, ik, vbTopIndex]
588
589
                                if vbTopK > vbTopE[ispin]:
                                    vbTopE[ispin] = vbTopK
590
591
                            if ebIndex < ev.shape[2]:
                                ebMinK = ev[ispin, ik, ebIndex]
592
593
                                if ebMinK < ebMinE[ispin]:
                                    ebMinE[ispin] = ebMinK
594
595
                    self.vbTopE = vbTopE
                    self.ebMinE = ebMinE
Markus Scheidgen's avatar
Markus Scheidgen committed
596
597
598
599
                    backend.addArrayValues(
                        "energy_reference_highest_occupied", np.array(vbTopE))
                    backend.addArrayValues(
                        "energy_reference_lowest_unoccupied", np.array(ebMinE))
600
601
602
                    if self.bands:
                        divisions = int(self.bands['divisions'])
                        backend.openNonOverlappingSection("section_k_band")
603
                        nsegments = self.kpoints.shape[0] // divisions
Markus Scheidgen's avatar
Markus Scheidgen committed
604
605
606
607
608
609
                        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]))
610
                        for isegment in range(nsegments):
Markus Scheidgen's avatar
Markus Scheidgen committed
611
612
613
614
615
616
617
618
                            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])
619
                            # "band_segm_labels"
Markus Scheidgen's avatar
Markus Scheidgen committed
620
621
622
623
                            backend.addArrayValues("band_segm_start_end", np.asarray(
                                [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment")
624
                        backend.closeNonOverlappingSection("section_k_band")
Markus Scheidgen's avatar
Markus Scheidgen committed
625
626
                        backend.openNonOverlappingSection(
                            "section_k_band_normalized")
627
628
                        specialPoints = {}
                        try:
Markus Scheidgen's avatar
Markus Scheidgen committed
629
630
631
                            specialPoints = get_special_points(
                                convert_unit_function("m", "angstrom")(self.cell))
                        except Exception as e:
Daniel Speckhard's avatar
Daniel Speckhard committed
632
                            self.logger.warn("failed to get special points/xtal structure", exc_info=e)
633
                        for isegment in range(nsegments):
Markus Scheidgen's avatar
Markus Scheidgen committed
634
635
636
637
638
639
640
641
642
643
                            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]]))
644
645
646
                            backend.addValue("band_segm_labels_normalized",
                                             [findLabel(specialPoints, kpt[isegment, 0]),
                                              findLabel(specialPoints, kpt[isegment, divisions - 1])])
Markus Scheidgen's avatar
Markus Scheidgen committed
647
648
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment_normalized")
649

Markus Scheidgen's avatar
Markus Scheidgen committed
650
651
                        backend.closeNonOverlappingSection(
                            "section_k_band_normalized")
652
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
653
654
                        backend.openNonOverlappingSection(
                            "section_eigenvalues")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
655
                        backend.addArrayValues("eigenvalues_values", ev)
Markus Scheidgen's avatar
Markus Scheidgen committed
656
657
658
659
                        backend.addArrayValues(
                            "eigenvalues_occupation", occupation)
                        backend.closeNonOverlappingSection(
                            "section_eigenvalues")
660
661
662
            else:
                backend.pwarn("unexpected tag %s in the eigenvalues" % el.tag)

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

666
    def onStart_calculation(self, parser, event, element, pathStr):
667
        backend = parser.backend
668
        gIndexes = parser.tagSections[pathStr]
Markus Scheidgen's avatar
Markus Scheidgen committed
669
670
        self.singleConfCalcs.append(
            gIndexes["section_single_configuration_calculation"])
671
672
        if self.waveCut:
            backend.openNonOverlappingSection("section_basis_set")
Markus Scheidgen's avatar
Markus Scheidgen committed
673
674
            backend.addValue(
                "mapping_section_basis_set_cell_dependent", self.waveCut)
675
676
            backend.closeNonOverlappingSection("section_basis_set")

677
678
    def onEnd_modeling(self, parser, event, element, pathStr):
        backend = parser.backend
679
        backend.addValue("x_vasp_unknown_incars", self.unknown_incars)
680
681
682
683
684
685
686
687
688
689
690
        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)
Markus Scheidgen's avatar
Markus Scheidgen committed
691
692
        backend.addArrayValues(
            "frame_sequence_local_frames_ref", np.asarray(self.singleConfCalcs))
693
        backend.closeSection("section_frame_sequence", frameSequenceGIndex)
694

695
696

    def onEnd_calculation(self, parser, event, element, pathStr):
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
697
        eConv = eV2J
698
699
700
        fConv = convert_unit_function("eV/angstrom", "N")
        pConv = convert_unit_function("eV/angstrom^3", "Pa")
        backend = parser.backend
Markus Scheidgen's avatar
Markus Scheidgen committed
701
702
        backend.addValue(
            "single_configuration_calculation_to_system_ref", self.lastSystemDescription)
703
        gIndexes = parser.tagSections["/modeling"]
Markus Scheidgen's avatar
Markus Scheidgen committed
704
705
        backend.addValue(
            "single_configuration_to_calculation_method_ref", gIndexes["section_method"])
706
707
708
709
710
711
712
713
        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)
714
                        elif name == "e_wo_entrp":
715
716
717
718
719
720
                            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:
Markus Scheidgen's avatar
Markus Scheidgen committed
721
722
                            backend.pwarn(
                                "Unexpected i tag with name %s in energy section" % name)
723
724
725
726
727
728
729
730
731
                    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)

732
    def onEnd_atominfo(self, parser, event, element, pathStr):
733
734
735
736
        nAtoms = None
        nAtomTypes = None
        atomTypes = []
        labels = []
737
        labels2 = None
738
        atomTypesDesc = []
739
        backend = parser.backend
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
        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":
Markus Scheidgen's avatar
Markus Scheidgen committed
756
757
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
758
759
760
761
762
                                else:
                                    line = atomsLine.findall("c")
                                    labels.append(line[0].text.strip())
                                    atomTypes.append(int(line[1].text.strip()))
                        else:
Markus Scheidgen's avatar
Markus Scheidgen committed
763
764
                            backend.pwarn(
                                "unexpected tag %s in atoms array in atominfo" % atomsEl.tag)
765
                elif name == "atomtypes":
766
767
768
769
770
771
772
                    keys = []
                    fieldTypes = []
                    for atomsEl in el:
                        if atomsEl.tag == "dimension":
                            pass
                        elif atomsEl.tag == "field":
                            keys.append(atomsEl.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
773
774
                            fieldTypes.append(
                                atomsEl.attrib.get("type", "float"))
775
                        elif atomsEl.tag == "set":
Markus Scheidgen's avatar
Markus Scheidgen committed
776
777
                            expectedKeys = ["atomspertype", "element",
                                            "mass", "valence", "pseudopotential"]
778
                            if keys != expectedKeys:
Markus Scheidgen's avatar
Markus Scheidgen committed
779
780
                                backend.pwarn(
                                    "unexpected fields in atomtype: %s vs %s" % (keys, expectedKeys))
781
782
                            for atomsLine in atomsEl:
                                if atomsLine.tag != "rc":
Markus Scheidgen's avatar
Markus Scheidgen committed
783
784
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
                                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:
Markus Scheidgen's avatar
Markus Scheidgen committed
800
801
                            backend.pwarn(
                                "unexpected tag %s in atomtypes array in atominfo" % atomsEl.tag)
802
                    kindIds = []
Markus Scheidgen's avatar
Markus Scheidgen committed
803
                    nEl = {}
804
                    kindLabels = []
805
                    for atomDesc in atomTypesDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
806
807
                        kindId = backend.openSection(
                            "section_method_atom_kind")
808
                        if 'element' in atomDesc:
809
                            elName = atomDesc['element'].strip()
810
                            try:
811
                                elNr = ase.data.chemical_symbols.index(elName)
Markus Scheidgen's avatar
Markus Scheidgen committed
812
813
814
815
816
817
818
                                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)
819
                            nEl[elName] = nElNow
Markus Scheidgen's avatar
Markus Scheidgen committed
820
821
                            elLabel = elName + \
                                (str(nElNow) if nElNow > 1 else "")
822
823
                            kindLabels.append(elLabel)
                            backend.addValue("method_atom_kind_label", elLabel)
824
                            if "mass" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
825
826
                                backend.addValue(
                                    "method_atom_kind_mass", atomDesc["mass"])
827
                            if "valence" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
828
829
                                backend.addValue(
                                    "method_atom_kind_explicit_electrons", atomDesc["valence"])
830
                            if "pseudopotential" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
831
832
                                backend.addValue(
                                    "method_atom_kind_pseudopotential_name", atomDesc["pseudopotential"])
833
                        kindIds.append(kindId)
Markus Scheidgen's avatar
Markus Scheidgen committed
834
835
836
837
                        backend.closeSection(
                            "section_method_atom_kind", kindId)
                    backend.addArrayValues("x_vasp_atom_kind_refs", np.asarray(
                        [kindIds[i-1] for i in atomTypes]))
838
                    labels2 = [kindLabels[i-1] for i in atomTypes]
839
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
840
841
                    backend.pwarn(
                        "unexpected array named %s in atominfo" % name)
842
843
            else:
                backend.pwarn("unexpected tag %s in atominfo" % el.tag)
844
        self.labels = np.asarray(labels2) if labels2 else np.asarray(labels)
845

846
    def incarOutTag(self, el):
847
848
        backend = self.parser.backend
        metaEnv = self.parser.backend.metaInfoEnv()
849
        if (el.tag != "i"):
Markus Scheidgen's avatar
Markus Scheidgen committed
850
851
            backend.pwarn("unexpected tag %s %s %r in incar" %
                          (el.tag, el.attrib, el.text))
852
        else:
853
            name    = el.attrib.get("name", None)
854
            valType = el.attrib.get("type")
855
            meta = metaEnv['x_vasp_incarOut_' + name]
856

857
858
859

            if not meta:
                # Unknown_Incars_Begin: storage into a dictionary
860
                if not valType:
861
                    # On vasp's xml files, valType *could* be absent if incar value is float
862
863
864
865
866
867
                    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)
868
869
                text_value = el.text.strip() # text representation of incar value
                try:
870
                    pyvalue = converter(text_value) # python data type
871
                except Exception:
872
                    pyvalue = text_value
873

874
875
876
                # save (name, pyvalue) into a dict
                self.unknown_incars[name] = pyvalue
                # Unknown_Incars_end
877
            else:
878
879
880
881
                if not valType:
                    valType = 'float'

                vasp_metainfo_type = vasp_to_metainfo_type_mapping.get(valType)[0]
882
                metainfo_type = meta.get('dtypeStr')
883
884
885
886
                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))

887
888
                elif metainfo_type != vasp_metainfo_type:
                    if  (metainfo_type == 'C' and vasp_metainfo_type == 'b'):
889
                        pass
890
                    elif  (metainfo_type == 'i' and vasp_metainfo_type == 'f'):
891
892
893
894
                        pass
                    else:
                        backend.pwarn("Data type mismatch: %s. Vasp_type: %s, metainfo_type: %s " %
                        (name, vasp_metainfo_type, metainfo_type))
895
896
897
                try:
                    shape = meta.get("shape", None)
                    converter = metaTypeTransformers.get(metainfo_type)
898
                    if not converter:
Markus Scheidgen's avatar
Markus Scheidgen committed
899
                        backend.pwarn(
900
                            "could not find converter for dtypeStr %s when handling meta info %s" %
901
                            (metainfo_type, meta ))
902
903
                    elif shape:
                        vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
904
905
                        backend.addValue(
                            meta["name"], [converter(x) for x in vals])
906
                    else:
907
908
                        # If-block to handle incars without value
                        if el.text == None:
909
                            el.text = ''
910
                        backend.addValue(meta["name"], converter(el.text))
911

912
                except:
Markus Scheidgen's avatar
Markus Scheidgen committed
913
914
                    backend.pwarn("Exception trying to handle incarOut %s: %s" % (
                        name, traceback.format_exc()))
915

916
                if name == 'ENMAX' or name == 'PREC':
Markus Scheidgen's avatar
Markus Scheidgen committed
917
918
919
920
921
922
923
                    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
924
925
926
927
928
929
                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'],
930
                        'MK': ['GGA_X_OPTB86_VDW'],
Markus Scheidgen's avatar
Markus Scheidgen committed
931
                        '--': ['GGA_X_PBE', 'GGA_C_PBE']  # should check potcar
932
933
934
                    }
                    functs = fMap.get(el.text.strip(), None)
                    if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
935
936
                        backend.pwarn("Unknown XC functional %s" %
                                      el.text.strip())
937
938
                    else:
                        for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
939
940
                            backend.openNonOverlappingSection(
                                "section_XC_functionals")
941
                            backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
942
943
                            backend.closeNonOverlappingSection(
                                "section_XC_functionals")
944
945
                elif name == "ISPIN":
                    self.ispin = int(el.text.strip())
946

947

Markus Scheidgen's avatar
Markus Scheidgen committed
948
    def separatorScan(self, element, backend, depth=0):
949
950
951
952
953
954
955
        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":</