parser_vasprun.py 50 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
import xml.etree.ElementTree as ET
import logging

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
36

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

40
41
42
43
44
vasp_to_metainfo_type_mapping = {
    'string': ['C'],
    'int': ['i'],
    'logical': ['b', 'C'],
    'float': ['f']}
45

46
47
48
49
50
51
52
53
54
55
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'}


56
def secondsFromEpoch(date):
Markus Scheidgen's avatar
Markus Scheidgen committed
57
58
    epoch = datetime(1970, 1, 1)
    ts = date-epoch
59
60
    return ts.seconds + ts.microseconds/1000.0

Markus Scheidgen's avatar
Markus Scheidgen committed
61
62
63
64
65
66
67

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


68
69
70
71
72
73
74
75
76
77
78
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(),
79
    'i': lambda x: int(float(x.strip())),
80
81
82
83
    'f': lambda x: float(x.strip()),
    'b': toBool,
}

84
85
86
87
88

class MyXMLParser(ET.XMLParser):

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

Markus Scheidgen's avatar
Markus Scheidgen committed
89
    def feed(self, data):
90
91
92
93
94
95
96
97
98
99
100
101
102
103
        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
104
105
106
        super(MyXMLParser, self).feed(mydata)


107
def transform2(y):
Markus Scheidgen's avatar
Markus Scheidgen committed
108
109
110
111
112
    if '**' in y:
        return float('nan')
    else:
        return y

113

Markus Scheidgen's avatar
Markus Scheidgen committed
114
def getVector(el, transform=float, field="v"):
115
116
    """ returns the vasp style vector contained in the element el (using field v).
    single elements are converted using the function convert"""
117
118
119
120
#
#    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)]
121

Markus Scheidgen's avatar
Markus Scheidgen committed
122

123
class VasprunContext(object):
124

125
    def __init__(self, logger=None):
126
        if logger is None:
127
            logger = logging.getLogger(__name__)
Markus Scheidgen's avatar
Markus Scheidgen committed
128
        self.logger = logger
129

130
131
132
133
        self.parser = None
        self.bands = None
        self.kpoints = None
        self.weights = None
134
135
        self.tetrahedrons = None
        self.tetrahedronVolume = None
136
        self.ispin = None
137
        self.ibrion = None
138
        self.lastSystemDescription = None
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
139
        self.labels = None
140
        self.singleConfCalcs = []
141
142
143
144
        self.vbTopE = None
        self.ebMinE = None
        self.eFermi = None
        self.cell = None
145
        self.angstrom_cell = None
146
147
        self.unknown_incars = {}

148
    sectionMap = {
149
        "modeling": ["section_run", "section_method"],
150
        "structure": ["section_system"],
151
        "calculation": ["section_single_configuration_calculation"]
152
153
    }

154
155
156
    def startedParsing(self, parser):
        self.parser = parser

157
    def onEnd_generator(self, parser, event, element, pathStr):
158
        backend = parser.backend
159
        program_name = g(element, "i/[@name='program']")
160
161
162
163
        if program_name.strip().upper() == "VASP":
            backend.addValue("program_name", "VASP")
        else:
            raise Exception("unexpected program name: %s" % program_name)
164
165
166
167
168
        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)
169
        backend.addValue("program_basis_set_type", "plane waves")
170
171
172
173
        date = g(element, "i/[@name='date']")
        pdate = None
        time = g(element, "i/[@name='time']")
        if date:
174
            pdate = datetime.strptime(date.strip(), "%Y %m %d")
175
        if pdate and time:
Markus Scheidgen's avatar
Markus Scheidgen committed
176
177
            pdate = datetime.combine(pdate.date(), datetime.strptime(
                time.strip(), "%H:%M:%S").timetz())
178
        if pdate:
Markus Scheidgen's avatar
Markus Scheidgen committed
179
180
            backend.addValue("program_compilation_datetime",
                             secondsFromEpoch(pdate))
181
        for i in element:
182
            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
183
184
                backend.pwarn("unexpected tag %s %s %r in generator" %
                              (i.tag, i.attrib, i.text))
185

186
    def onEnd_incar(self, parser, event, element, pathStr):
187
        backend = parser.backend
188
        metaEnv = parser.backend.metaInfoEnv()
189
        dft_plus_u = False
190
191
        ibrion = None
        nsw = 0
192
193
        for el in element:
            if el.tag == "v":
194
195
196
197
198
                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))
199
                    continue
200
                #- -
201
202
203
                vector_val = np.asarray(getVector(el))
                backend.addArrayValues(meta.get('name'), vector_val)
            elif el.tag == "i":
204
205
206
207
                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
208
209
                    backend.pwarn("Unknown INCAR parameter (not registered in the meta data): %s %s %r" % (
                        el.tag, el.attrib, el.text))
210
211
212
213
                elif valType:
                    expectedMetaType = {
                        'string': ['C'],
                        'int': ['i'],
Markus Scheidgen's avatar
Markus Scheidgen committed
214
                        'logical': ['b', 'C']
215
216
                    }.get(valType)
                    if not expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
217
218
                        backend.pwarn("Unknown value type %s encountered in INCAR: %s %s %r" % (
                            valType, el.tag, el.attrib, el.text))
219
                    elif not meta.get('dtypeStr') in expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
220
221
                        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))
222
223
224
225
226
                    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
227
228
                            backend.pwarn(
                                "could not find converter for dtypeStr %s when handling meta info %s" % (dtypeStr, ))
229
230
                        elif shape:
                            vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
231
232
                            backend.addValue(
                                meta["name"], [converter(x) for x in vals])
233
234
235
                        else:
                            backend.addValue(meta["name"], converter(el.text))
                    if name == 'GGA':
236
237
                        # FIXME tmk: many options are not coded yet. See
                        # https://www.vasp.at/wiki/index.php/GGA
238
239
240
241
                        fMap = {
                            '91': ['GGA_X_PW91', 'GGA_C_PW91'],
                            'PE': ['GGA_X_PBE', 'GGA_C_PBE'],
                            'RP': ['GGA_X_RPBE', 'GGA_C_PBE'],
242
243
                            'PS': ['GGA_C_PBE_SOL', 'GGA_X_PBE_SOL'],
                            'MK': ['GGA_X_OPTB86_VDW']
244
245
246
                        }
                        functs = fMap.get(el.text.strip(), None)
                        if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
247
248
                            backend.pwarn("Unknown XC functional %s" %
                                          el.text.strip())
249
250
                        else:
                            for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
251
252
                                backend.openNonOverlappingSection(
                                    "section_XC_functionals")
253
                                backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
254
255
                                backend.closeNonOverlappingSection(
                                    "section_XC_functionals")
256
257
                    elif name == "ISPIN":
                        self.ispin = int(el.text.strip())
258
259
260
                    elif name == "LDAU":
                        if re.match(".?[Tt](?:[rR][uU][eE])?.?|[yY](?:[eE][sS])?|1", el.text.strip()):
                            dft_plus_u = True
261
262
263
264
                    elif name == "IBRION":
                        ibrion = int(el.text.strip())
                    elif name == "NSW":
                        nsw = int(el.text.strip())
265
266
267
            else:
                backend.pwarn("unexpected tag %s %s %r in incar" %
                              (el.tag, el.attrib, el.text))
268
269
        if ibrion is None:
            ibrion = -1 if nsw == 0 or nsw == 1 else 0
270
271
        if nsw == 0:
            ibrion = -1
272
        self.ibrion = ibrion
273
274
275
276
        if dft_plus_u:
            backend.addValue("electronic_structure_method", "DFT+U")
        else:
            backend.addValue("electronic_structure_method", "DFT")
277

278
    def onEnd_kpoints(self, parser, event, element, pathStr):
279
        backend = parser.backend
280
281
        self.bands = None
        self.kpoints = None
282
283
        self.weights = None
        for el in element:
284
            if el.tag == "generation":
285
                param = el.attrib.get("param", None) # eg. listgenerated, Monkhorst-Pack, Gamma
286
                if param:
Markus Scheidgen's avatar
Markus Scheidgen committed
287
288
                    backend.addValue(
                        "x_vasp_k_points_generation_method", param)
289
                if param == "listgenerated":
290
                    # This implies a path on k-space, potentially a bandstructure calculation
291
                    # Save k-path info into a dictionary
292
                    self.bands = {
Markus Scheidgen's avatar
Markus Scheidgen committed
293
                        "divisions": g(el, "i/[@name='divisions']", None),
294
295
                        "points": getVector(el)
                    }
296
297
298
299

                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
300
301
                    pass
                else:
302
                    backend.pwarn("Unknown k point generation method '%s'" %(param))
303
304
305
306
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "kpointlist":
                    self.kpoints = np.asarray(getVector(el))
307
                    backend.addArrayValues("k_mesh_points", self.kpoints)
308
309
                elif name == "weights":
                    self.weights = np.asarray(getVector(el))
Markus Scheidgen's avatar
Markus Scheidgen committed
310
311
                    backend.addArrayValues(
                        "k_mesh_weights", self.weights.flatten())
312
313
314
315
                elif name == "tetrahedronlist":
                    self.tetrahedrons = np.asarray(getVector(el), dtype=np.int)
                    backend.addArrayValues(
                        "x_vasp_tetrahedrons_list", self.tetrahedrons)
316
317
                else:
                    backend.pwarn("Unknown array %s in kpoints" % name)
318
319
320
321
322
323
324
325
326
327
328
            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)
329
            else:
330
                backend.pwarn("Unknown tag %s in kpoints" % el.tag)
331

332

333
    def onEnd_structure(self, parser, event, element, pathStr):
334
        backend = parser.backend
335
        gIndexes = parser.tagSections[pathStr]
336
        self.lastSystemDescription = gIndexes["section_system"]
337
        self.cell = None
338
339
340
341
342
343
        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":
344
                            conv = convert_unit_function("angstrom", "m")
Markus Scheidgen's avatar
Markus Scheidgen committed
345
346
                            self.cell = getVector(
                                cellEl, lambda x: conv(float(x)))
347
                            self.angstrom_cell = np.array(getVector(cellEl))
Markus Scheidgen's avatar
Markus Scheidgen committed
348
349
350
351
                            backend.addArrayValues(
                                "simulation_cell", np.asarray(self.cell))
                            backend.addArrayValues(
                                "configuration_periodic_dimensions", np.ones(3, dtype=bool))
352
                        elif name == "rec_basis":
353
354
                            pass
                        else:
Markus Scheidgen's avatar
Markus Scheidgen committed
355
356
                            backend.pwarn(
                                "Unexpected varray %s in crystal" % name)
357
358
                    elif cellEl.tag == "i":
                        if cellEl.attrib.get("name") != "volume":
Markus Scheidgen's avatar
Markus Scheidgen committed
359
360
                            backend.pwarn(
                                "Unexpected i value %s in crystal" % cellEl.attrib)
361
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
362
363
                        backend.pwarn("Unexpected tag %s %s %r in crystal" % (
                            cellEl.tag, cellEl.attrib, cellEl.text))
364
365
366
367
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "positions":
                    pos = getVector(el)
Markus Scheidgen's avatar
Markus Scheidgen committed
368
369
                    backend.addArrayValues(
                        "atom_positions", np.dot(np.asarray(pos), self.cell))
370
371
372
373
                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))
374
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
375
376
                    backend.pwarn(
                        "Unexpected varray in structure %s" % el.attrib)
377
378
379
            elif el.tag == "nose":
                nose = getVector(el)
                backend.addArrayValues("x_vasp_nose_thermostat", nose)
380
            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
381
                backend.pwarn("Unexpected tag in structure %s %s %r" %
382
                              (el.tag, el.attrib, el.text))
383
        if self.labels is not None:
384
            backend.addArrayValues("atom_labels", self.labels)
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
385

386
    def onEnd_eigenvalues(self, parser, event, element, pathStr):
387
388
        if pathStr != "modeling/calculation/eigenvalues":
            return True
389
390
391
        backend = parser.backend
        eigenvalues = None
        occupation = None
392
393
394
395
396
397
398
399
        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":
400
                        isp = -1
401
402
                        for spinEl in arrEl:
                            if spinEl.tag == "set":
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
403
                                ik = -1
404
                                isp += 1
405
406
                                for kEl in spinEl:
                                    if kEl.tag == "set":
407
                                        ik += 1
Markus Scheidgen's avatar
Markus Scheidgen committed
408
409
                                        bands = np.asarray(
                                            getVector(kEl, field="r"))
410
                                        if eigenvalues is None:
Markus Scheidgen's avatar
Markus Scheidgen committed
411
412
413
414
415
416
                                            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]
417
                                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
418
419
                                        backend.pwarn(
                                            "unexpected tag %s in k array of the eigenvalues" % kEl.tag)
420
                            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
421
422
                                backend.pwarn(
                                    "unexpected tag %s in spin array of the eigenvalues" % spinEl.tag)
423
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
424
425
                        backend.pwarn(
                            "unexpected tag %s in array of the eigenvalues" % arrEl.tag)
426
                if eigenvalues is not None:
427

Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
428
                    ev = eV2JV(eigenvalues)
429
430
                    vbTopE = []
                    ebMinE = []
431
                    for ispin in range(occupation.shape[0]):
432
433
                        vbTopE.append(float('-inf'))
                        ebMinE.append(float('inf'))
434
                        for ik in range(occupation.shape[1]):
Markus Scheidgen's avatar
Markus Scheidgen committed
435
436
437
                            ebIndex = bisect.bisect_right(
                                -occupation[ispin, ik, :], -0.5) - 1
                            vbTopIndex = ebIndex - 1
438
439
                            if vbTopIndex >= 0:
                                vbTopK = ev[ispin, ik, vbTopIndex]
440
441
                                if vbTopK > vbTopE[ispin]:
                                    vbTopE[ispin] = vbTopK
442
443
                            if ebIndex < ev.shape[2]:
                                ebMinK = ev[ispin, ik, ebIndex]
444
445
                                if ebMinK < ebMinE[ispin]:
                                    ebMinE[ispin] = ebMinK
446
447
                    self.vbTopE = vbTopE
                    self.ebMinE = ebMinE
Markus Scheidgen's avatar
Markus Scheidgen committed
448
449
450
451
                    backend.addArrayValues(
                        "energy_reference_highest_occupied", np.array(vbTopE))
                    backend.addArrayValues(
                        "energy_reference_lowest_unoccupied", np.array(ebMinE))
452
453
454
                    if self.bands:
                        divisions = int(self.bands['divisions'])
                        backend.openNonOverlappingSection("section_k_band")
455
                        nsegments = self.kpoints.shape[0] // divisions
Markus Scheidgen's avatar
Markus Scheidgen committed
456
457
458
                        kpt = np.reshape(
                            self.kpoints, (nsegments, divisions, 3))
                        energies = np.reshape(
459
                            ev, (self.ispin, nsegments, divisions, bands.shape[0]))
Markus Scheidgen's avatar
Markus Scheidgen committed
460
461
                        occ = np.reshape(
                            occupation, (self.ispin, nsegments, divisions, bands.shape[0]))
462
                        for isegment in range(nsegments):
Markus Scheidgen's avatar
Markus Scheidgen committed
463
464
465
466
467
468
469
470
                            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])
471
                            # "band_segm_labels"
Markus Scheidgen's avatar
Markus Scheidgen committed
472
473
474
475
                            backend.addArrayValues("band_segm_start_end", np.asarray(
                                [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment")
476
                        backend.closeNonOverlappingSection("section_k_band")
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
                        backend.openNonOverlappingSection(
                            "section_k_band_normalized")
                        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.closeNonOverlappingSection(
                                "section_k_band_segment_normalized")

                        backend.closeNonOverlappingSection(
                            "section_k_band_normalized")
495
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
496
497
                        backend.openNonOverlappingSection(
                            "section_eigenvalues")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
498
                        backend.addArrayValues("eigenvalues_values", ev)
Markus Scheidgen's avatar
Markus Scheidgen committed
499
500
501
502
                        backend.addArrayValues(
                            "eigenvalues_occupation", occupation)
                        backend.closeNonOverlappingSection(
                            "section_eigenvalues")
503
504
505
            else:
                backend.pwarn("unexpected tag %s in the eigenvalues" % el.tag)

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

509
    def onStart_calculation(self, parser, event, element, pathStr):
510
        backend = parser.backend
511
        gIndexes = parser.tagSections[pathStr]
Markus Scheidgen's avatar
Markus Scheidgen committed
512
513
        self.singleConfCalcs.append(
            gIndexes["section_single_configuration_calculation"])
514
515
        if self.waveCut:
            backend.openNonOverlappingSection("section_basis_set")
Markus Scheidgen's avatar
Markus Scheidgen committed
516
517
            backend.addValue(
                "mapping_section_basis_set_cell_dependent", self.waveCut)
518
519
            backend.closeNonOverlappingSection("section_basis_set")

520
521
    def onEnd_modeling(self, parser, event, element, pathStr):
        backend = parser.backend
522
        backend.addValue("x_vasp_unknown_incars", self.unknown_incars)
523
524
525
526
527
528
529
530
531
532
533
        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
534
535
        backend.addArrayValues(
            "frame_sequence_local_frames_ref", np.asarray(self.singleConfCalcs))
536
        backend.closeSection("section_frame_sequence", frameSequenceGIndex)
537

538
539

    def onEnd_calculation(self, parser, event, element, pathStr):
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
540
        eConv = eV2J
541
542
543
        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
544
545
        backend.addValue(
            "single_configuration_calculation_to_system_ref", self.lastSystemDescription)
546
        gIndexes = parser.tagSections["/modeling"]
Markus Scheidgen's avatar
Markus Scheidgen committed
547
548
        backend.addValue(
            "single_configuration_to_calculation_method_ref", gIndexes["section_method"])
549
550
551
552
553
554
555
556
        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)
557
                        elif name == "e_wo_entrp":
558
559
560
561
562
563
                            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
564
565
                            backend.pwarn(
                                "Unexpected i tag with name %s in energy section" % name)
566
567
568
569
570
571
572
573
574
                    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)

575
    def onEnd_atominfo(self, parser, event, element, pathStr):
576
577
578
579
        nAtoms = None
        nAtomTypes = None
        atomTypes = []
        labels = []
580
        labels2 = None
581
        atomTypesDesc = []
582
        backend = parser.backend
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
        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
599
600
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
601
602
603
604
605
                                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
606
607
                            backend.pwarn(
                                "unexpected tag %s in atoms array in atominfo" % atomsEl.tag)
608
                elif name == "atomtypes":
609
610
611
612
613
614
615
                    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
616
617
                            fieldTypes.append(
                                atomsEl.attrib.get("type", "float"))
618
                        elif atomsEl.tag == "set":
Markus Scheidgen's avatar
Markus Scheidgen committed
619
620
                            expectedKeys = ["atomspertype", "element",
                                            "mass", "valence", "pseudopotential"]
621
                            if keys != expectedKeys:
Markus Scheidgen's avatar
Markus Scheidgen committed
622
623
                                backend.pwarn(
                                    "unexpected fields in atomtype: %s vs %s" % (keys, expectedKeys))
624
625
                            for atomsLine in atomsEl:
                                if atomsLine.tag != "rc":
Markus Scheidgen's avatar
Markus Scheidgen committed
626
627
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
                                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
643
644
                            backend.pwarn(
                                "unexpected tag %s in atomtypes array in atominfo" % atomsEl.tag)
645
                    kindIds = []
Markus Scheidgen's avatar
Markus Scheidgen committed
646
                    nEl = {}
647
                    kindLabels = []
648
                    for atomDesc in atomTypesDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
649
650
                        kindId = backend.openSection(
                            "section_method_atom_kind")
651
                        if 'element' in atomDesc:
652
                            elName = atomDesc['element'].strip()
653
                            try:
654
                                elNr = ase.data.chemical_symbols.index(elName)
Markus Scheidgen's avatar
Markus Scheidgen committed
655
656
657
658
659
660
661
                                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)
662
                            nEl[elName] = nElNow
Markus Scheidgen's avatar
Markus Scheidgen committed
663
664
                            elLabel = elName + \
                                (str(nElNow) if nElNow > 1 else "")
665
666
                            kindLabels.append(elLabel)
                            backend.addValue("method_atom_kind_label", elLabel)
667
                            if "mass" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
668
669
                                backend.addValue(
                                    "method_atom_kind_mass", atomDesc["mass"])
670
                            if "valence" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
671
672
                                backend.addValue(
                                    "method_atom_kind_explicit_electrons", atomDesc["valence"])
673
                            if "pseudopotential" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
674
675
                                backend.addValue(
                                    "method_atom_kind_pseudopotential_name", atomDesc["pseudopotential"])
676
                        kindIds.append(kindId)
Markus Scheidgen's avatar
Markus Scheidgen committed
677
678
679
680
                        backend.closeSection(
                            "section_method_atom_kind", kindId)
                    backend.addArrayValues("x_vasp_atom_kind_refs", np.asarray(
                        [kindIds[i-1] for i in atomTypes]))
681
                    labels2 = [kindLabels[i-1] for i in atomTypes]
682
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
683
684
                    backend.pwarn(
                        "unexpected array named %s in atominfo" % name)
685
686
            else:
                backend.pwarn("unexpected tag %s in atominfo" % el.tag)
687
        self.labels = np.asarray(labels2) if labels2 else np.asarray(labels)
688

689
    def incarOutTag(self, el):
690
691
        backend = self.parser.backend
        metaEnv = self.parser.backend.metaInfoEnv()
692
        if (el.tag != "i"):
Markus Scheidgen's avatar
Markus Scheidgen committed
693
694
            backend.pwarn("unexpected tag %s %s %r in incar" %
                          (el.tag, el.attrib, el.text))
695
        else:
696
            name    = el.attrib.get("name", None)
697
            valType = el.attrib.get("type")
698
            meta = metaEnv['x_vasp_incarOut_' + name]
699

700
701
702

            if not meta:
                # Unknown_Incars_Begin: storage into a dictionary
703
                if not valType:
704
                    # On vasp's xml files, valType *could* be absent if incar value is float
705
706
707
708
709
710
                    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)
711
712
                text_value = el.text.strip() # text representation of incar value
                try:
713
                    pyvalue = converter(text_value) # python data type
714
                except Exception:
715
                    pyvalue = text_value
716

717
718
719
                # save (name, pyvalue) into a dict
                self.unknown_incars[name] = pyvalue
                # Unknown_Incars_end
720
            else:
721
722
723
724
                if not valType:
                    valType = 'float'

                vasp_metainfo_type = vasp_to_metainfo_type_mapping.get(valType)[0]
725
                metainfo_type = meta.get('dtypeStr')
726
727
728
729
                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))

730
731
                elif metainfo_type != vasp_metainfo_type:
                    if  (metainfo_type == 'C' and vasp_metainfo_type == 'b'):
732
                        pass
733
                    elif  (metainfo_type == 'i' and vasp_metainfo_type == 'f'):
734
735
736
737
                        pass
                    else:
                        backend.pwarn("Data type mismatch: %s. Vasp_type: %s, metainfo_type: %s " %
                        (name, vasp_metainfo_type, metainfo_type))
738
739
740
                try:
                    shape = meta.get("shape", None)
                    converter = metaTypeTransformers.get(metainfo_type)
741
                    if not converter:
Markus Scheidgen's avatar
Markus Scheidgen committed
742
                        backend.pwarn(
743
                            "could not find converter for dtypeStr %s when handling meta info %s" %
744
                            (metainfo_type, meta ))
745
746
                    elif shape:
                        vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
747
748
                        backend.addValue(
                            meta["name"], [converter(x) for x in vals])
749
                    else:
750
751
                        # If-block to handle incars without value
                        if el.text == None:
752
                            el.text = ''
753
                        backend.addValue(meta["name"], converter(el.text))
754

755
                except:
Markus Scheidgen's avatar
Markus Scheidgen committed
756
757
                    backend.pwarn("Exception trying to handle incarOut %s: %s" % (
                        name, traceback.format_exc()))
758

759
                if name == 'ENMAX' or name == 'PREC':
Markus Scheidgen's avatar
Markus Scheidgen committed
760
761
762
763
764
765
766
                    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
767
768
769
770
771
772
                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'],
773
                        'MK': ['GGA_X_OPTB86_VDW'],
Markus Scheidgen's avatar
Markus Scheidgen committed
774
                        '--': ['GGA_X_PBE', 'GGA_C_PBE']  # should check potcar
775
776
777
                    }
                    functs = fMap.get(el.text.strip(), None)
                    if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
778
779
                        backend.pwarn("Unknown XC functional %s" %
                                      el.text.strip())
780
781
                    else:
                        for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
782
783
                            backend.openNonOverlappingSection(
                                "section_XC_functionals")
784
                            backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
785
786
                            backend.closeNonOverlappingSection(
                                "section_XC_functionals")
787
788
                elif name == "ISPIN":
                    self.ispin = int(el.text.strip())
789

790

Markus Scheidgen's avatar
Markus Scheidgen committed
791
    def separatorScan(self, element, backend, depth=0):
792
793
794
795
796
797
798
        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":
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
799
                        self.separatorScan(el, backend, depth + 1)
800
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
801
802
803
                        # backend.pwarn("unexpected tag %s %s in parameters separator %s at depth %d" % (
                        #     el.tag, el.attrib, separatorName, depth))
                        pass
804
805
806
            elif separators.tag == "i":
                self.incarOutTag(separators)
            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
807
808
809
                # backend.pwarn("unexpected tag %s %s in parameters at depth %d" % (
                #     separators.tag, separators.attrib, depth))
                pass
810

811
    def onEnd_parameters(self, parser, event, element, pathStr):
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
812
        self.separatorScan(element, parser.backend)
813
        backend = parser.backend
814
        try:
Markus Scheidgen's avatar
Markus Scheidgen committed
815
816
817
818
819
820
821
822
823
824
825
826
827
828
            self.prec
            try:
                self.enmax
                self.waveCut = backend.openNonOverlappingSection(
                    "section_basis_set_cell_dependent")
                backend.addValue("basis_set_planewave_cutoff",
                                 eV2J(self.enmax*self.prec))
                backend.closeNonOverlappingSection(
                    "section_basis_set_cell_dependent")
                backend.openNonOverlappingSection("section_method_basis_set")
                backend.addValue(
                    "mapping_section_method_basis_set_cell_associated", self.waveCut)
                backend.closeNonOverlappingSection("section_method_basis_set")
            except AttributeError:
829
                import traceback
830
                traceback.print_exc()
Markus Scheidgen's avatar
Markus Scheidgen committed
831
832
                backend.pwarn(
                    "Missing ENMAX for calculating plane wave basis cut off ")
833
        except AttributeError:
Markus Scheidgen's avatar
Markus Scheidgen committed
834
835
            backend.pwarn(
                "Missing PREC for calculating plane wave basis cut off ")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
836
837
838

    def onEnd_dos(self, parser, event, element, pathStr):
        "density of states"
839
        backend = parser.backend
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
840
841
842
843
        backend.openNonOverlappingSection("section_dos")
        for el in element:
            if el.tag == "i":
                if el.attrib.get("name") == "efermi":
844
                    self.eFermi = eV2J(float(el.text.strip()))
Markus Scheidgen's avatar
Markus Scheidgen committed
845
                    backend.addArrayValues(
846
                        "energy_reference_fermi", np.array([self.eFermi] * self.ispin))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
847
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
848
849
                    backend.pwarn("unexpected tag %s %s in dos" %
                                  (el.tag, el.attrib))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
850
851
852
853
854
855
856
857
858
859
            elif el.tag == "total":
                for el1 in el:
                    if el1.tag == "array":
                        for el2 in el1:
                            if el2.tag == "dimension" or el2.tag == "field":
                                pass
                            elif el2.tag == "set":
                                dosL = []
                                for spinComponent in el2:
                                    if spinComponent.tag == "set":
Markus Scheidgen's avatar
Markus Scheidgen committed
860
861
                                        dosL.append(
                                            getVector(spinComponent, field="r"))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
862
                                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
863
864
                                        backend.pwarn("unexpected tag %s %s in dos total array set" % (
                                            spinComponent.tag, spinComponent.attrib))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
865
866
                                dosA = np.asarray(dosL)
                                if len(dosA.shape) != 3:
Markus Scheidgen's avatar
Markus Scheidgen committed
867
868
                                    raise Exception("unexpected shape %s (%s) for total dos (ragged arrays?)" % (
                                        dosA.shape), dosA.dtype)
869
870
871
872
873
                                dosE = eV2JV(dosA[0, :, 0])
                                dosI = dosA[:, :, 2]
                                dosV = dosA[:, :, 1]

                                # Convert the DOS values to SI. VASP uses the
874
875
876
877
                                # following units in the output:
                                # states/eV/cell. This means that the volume
                                # dependence has been introduced by multiplying
                                # by the cell volume
878
879
                                # the integrated dos value is the number of electrons until that energy level
                                # and thus not directly energy dependent anymore
880
881
                                joule_in_ev = convert_unit(1, "eV", "J")
                                dosV = dosV / joule_in_ev
882

Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
883
                                backend.addArrayValues("dos_energies", dosE)
884
885
                                cell_volume = np.abs(np.linalg.det(self.cell))
                                backend.addArrayValues("dos_values", dosV * cell_volume)
Markus Scheidgen's avatar
Markus Scheidgen committed
886
887
                                backend.addArrayValues(
                                    "dos_integrated_values", dosI)
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
888
                            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
889
890
                                backend.pwarn("unexpected tag %s %s in dos total array" % (
                                    el2.tag, el2.attrib))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
891
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
892
893
                        backend.pwarn("unexpected tag %s %s in dos total" % (
                            el2.tag, el2.attrib))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
894
895
896
            elif el.tag == "partial":
                for el1 in el:
                    if el1.tag == "array":
Markus Scheidgen's avatar
Markus Scheidgen committed
897
                        lm = []
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
898
899
900
901
902
903
904
905
                        for el2 in el1:
                            if el2.tag == "dimension":
                                pass
                            elif el2.tag == "field":
                                if el2.text.strip() == "energy":
                                    pass
                                else:
                                    strLm = {
Markus Scheidgen's avatar
Markus Scheidgen committed
906
907
908
909
910
911
912
913
914
915
916
917
                                        "s": [0, 0],
                                        "p": [1, -1],
                                        "px": [1, 0],
                                        "py": [1, 1],
                                        "pz": [1, 2],
                                        "d": [2, -1],
                                        "dx2": [2, 0],
                                        "dxy": [2, 1],
                                        "dxz": [2, 2],
                                        "dy2": [2, 3],
                                        "dyz": [2, 4],
                                        "dz2": [2, 5]
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
918
                                    }
Markus Scheidgen's avatar
Markus Scheidgen committed
919
920
                                    lm.append(
                                        strLm.get(el2.text.strip(), [-1, -1]))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
921
922
923
924
925
926
927
928
                            elif el2.tag == "set":
                                dosL = []
                                for atom in el2:
                                    if atom.tag == "set":
                                        atomL = []
                                        dosL.append(atomL)
                                        for spinComponent in atom:
                                            if spinComponent.tag == "set":
Markus Scheidgen's avatar
Markus Scheidgen committed
929
930
                                                atomL.append(
                                                    getVector(spinComponent, field="r"))
Mohamed, Fawzi Roberto (fawzi)'s avatar