parser_vasprun.py 52.5 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
Markus Scheidgen's avatar
Markus Scheidgen committed
19
20
import sys
import bisect
21
from datetime import datetime
Markus Scheidgen's avatar
Markus Scheidgen committed
22
23
24
import os
import re
import traceback
25
import numpy as np
26
import ase.geometry
27
import ase.data
28
from math import pi
Markus Scheidgen's avatar
Markus Scheidgen committed
29
30
31
32
33
34
35
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
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)

Markus Scheidgen's avatar
Markus Scheidgen committed
40

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


special_points = {
71
    'cubic': {'Γ': [0, 0, 0],
72
73
74
              'M': [1 / 2, 1 / 2, 0],
              'R': [1 / 2, 1 / 2, 1 / 2],
              'X': [0, 1 / 2, 0]},
75
    'fcc': {'Γ': [0, 0, 0],
76
77
78
79
80
            '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]},
81
    'bcc': {'Γ': [0, 0, 0],
82
83
84
            'H': [1 / 2, -1 / 2, 1 / 2],
            'P': [1 / 4, 1 / 4, 1 / 4],
            'N': [0, 0, 1 / 2]},
85
    'tetragonal': {'Γ': [0, 0, 0],
86
87
88
89
90
                   '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]},
91
    'orthorhombic': {'Γ': [0, 0, 0],
92
93
94
95
96
97
98
                     '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]},
99
    'hexagonal': {'Γ': [0, 0, 0],
100
101
102
103
104
105
106
107
                  '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 = {
108
109
110
111
112
113
114
    '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'}
115
116
117
118
119
120
121
122
123
124
125
126
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


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

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

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

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

    lattice = crystal_structure_from_cell(cell)

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

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

    # Here, we need the cell:
    eta = (1 - b * cos(alpha) / c) / (2 * sin(alpha)**2)
    nu = 1 / 2 - eta * c * cos(alpha) / b
167
    return {'Γ': [0, 0, 0],
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
            '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
184

185
186
187
188
189
190
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
191

192
def secondsFromEpoch(date):
Markus Scheidgen's avatar
Markus Scheidgen committed
193
194
    epoch = datetime(1970, 1, 1)
    ts = date-epoch
195
196
    return ts.seconds + ts.microseconds/1000.0

Markus Scheidgen's avatar
Markus Scheidgen committed
197
198
199
200
201
202
203

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


204
205
206
207
208
209
210
211
212
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

Markus Scheidgen's avatar
Markus Scheidgen committed
213

214
215
216
217
218
219
220
metaTypeTransformers = {
    'C': lambda x: x.strip(),
    'i': lambda x: int(x.strip()),
    'f': lambda x: float(x.strip()),
    'b': toBool,
}

221
222
223
224
225

class MyXMLParser(ET.XMLParser):

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

Markus Scheidgen's avatar
Markus Scheidgen committed
226
    def feed(self, data):
227
228
229
230
231
232
233
234
235
236
237
238
239
240
        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
241
242
243
        super(MyXMLParser, self).feed(mydata)


244
def transform2(y):
Markus Scheidgen's avatar
Markus Scheidgen committed
245
246
247
248
249
    if '**' in y:
        return float('nan')
    else:
        return y

250

Markus Scheidgen's avatar
Markus Scheidgen committed
251
def getVector(el, transform=float, field="v"):
252
253
    """ returns the vasp style vector contained in the element el (using field v).
    single elements are converted using the function convert"""
254
255
256
257
#
#    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)]
258

Markus Scheidgen's avatar
Markus Scheidgen committed
259

260
class VasprunContext(object):
Markus Scheidgen's avatar
Markus Scheidgen committed
261
262
    def __init__(self, logger=None):
        self.logger = logger
263
264
265
266
        self.parser = None
        self.bands = None
        self.kpoints = None
        self.weights = None
267
        self.ispin = None
268
        self.ibrion = None
269
        self.lastSystemDescription = None
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
270
        self.labels = None
271
        self.singleConfCalcs = []
272
273
274
275
        self.vbTopE = None
        self.ebMinE = None
        self.eFermi = None
        self.cell = None
276
        self.angstrom_cell = None
277

Markus Scheidgen's avatar
Markus Scheidgen committed
278
279
        self.unknown_incar_params = []

Markus Scheidgen's avatar
Markus Scheidgen committed
280
281
282
        if self.logger is None:
            logger = logging.getLogger(__name__)

283
    sectionMap = {
284
        "modeling": ["section_run", "section_method"],
285
        "structure": ["section_system"],
286
        "calculation": ["section_single_configuration_calculation"]
287
288
    }

289
290
291
    def startedParsing(self, parser):
        self.parser = parser

292
    def onEnd_generator(self, parser, event, element, pathStr):
293
        backend = parser.backend
294
        program_name = g(element, "i/[@name='program']")
295
296
297
298
        if program_name.strip().upper() == "VASP":
            backend.addValue("program_name", "VASP")
        else:
            raise Exception("unexpected program name: %s" % program_name)
299
300
301
302
303
        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)
304
        backend.addValue("program_basis_set_type", "plane waves")
305
306
307
308
        date = g(element, "i/[@name='date']")
        pdate = None
        time = g(element, "i/[@name='time']")
        if date:
309
            pdate = datetime.strptime(date.strip(), "%Y %m %d")
310
        if pdate and time:
Markus Scheidgen's avatar
Markus Scheidgen committed
311
312
            pdate = datetime.combine(pdate.date(), datetime.strptime(
                time.strip(), "%H:%M:%S").timetz())
313
        if pdate:
Markus Scheidgen's avatar
Markus Scheidgen committed
314
315
            backend.addValue("program_compilation_datetime",
                             secondsFromEpoch(pdate))
316
        for i in element:
317
            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
318
319
                backend.pwarn("unexpected tag %s %s %r in generator" %
                              (i.tag, i.attrib, i.text))
320

321
    def onEnd_incar(self, parser, event, element, pathStr):
322
        backend = parser.backend
323
        metaEnv = parser.backend.metaInfoEnv()
324
        dft_plus_u = False
325
326
        ibrion = None
        nsw = 0
327
328
        for el in element:
            if (el.tag != "i"):
Markus Scheidgen's avatar
Markus Scheidgen committed
329
330
                backend.pwarn("unexpected tag %s %s %r in incar" %
                              (el.tag, el.attrib, el.text))
331
332
333
334
335
            else:
                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
336
337
                    backend.pwarn("Unknown INCAR parameter (not registered in the meta data): %s %s %r" % (
                        el.tag, el.attrib, el.text))
338
339
340
341
                elif valType:
                    expectedMetaType = {
                        'string': ['C'],
                        'int': ['i'],
Markus Scheidgen's avatar
Markus Scheidgen committed
342
                        'logical': ['b', 'C']
343
344
                    }.get(valType)
                    if not expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
345
346
                        backend.pwarn("Unknown value type %s encountered in INCAR: %s %s %r" % (
                            valType, el.tag, el.attrib, el.text))
347
                    elif not meta.get('dtypeStr') in expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
348
349
                        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))
350
351
352
353
354
                    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
355
356
                            backend.pwarn(
                                "could not find converter for dtypeStr %s when handling meta info %s" % (dtypeStr, ))
357
358
                        elif shape:
                            vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
359
360
                            backend.addValue(
                                meta["name"], [converter(x) for x in vals])
361
362
363
364
365
366
367
368
369
370
371
                        else:
                            backend.addValue(meta["name"], converter(el.text))
                    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']
                        }
                        functs = fMap.get(el.text.strip(), None)
                        if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
372
373
                            backend.pwarn("Unknown XC functional %s" %
                                          el.text.strip())
374
375
                        else:
                            for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
376
377
                                backend.openNonOverlappingSection(
                                    "section_XC_functionals")
378
                                backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
379
380
                                backend.closeNonOverlappingSection(
                                    "section_XC_functionals")
381
382
                    elif name == "ISPIN":
                        self.ispin = int(el.text.strip())
383
384
385
                    elif name == "LDAU":
                        if re.match(".?[Tt](?:[rR][uU][eE])?.?|[yY](?:[eE][sS])?|1", el.text.strip()):
                            dft_plus_u = True
386
387
388
389
390
391
                    elif name == "IBRION":
                        ibrion = int(el.text.strip())
                    elif name == "NSW":
                        nsw = int(el.text.strip())
        if ibrion is None:
            ibrion = -1 if nsw == 0 or nsw == 1 else 0
392
393
        if nsw == 0:
            ibrion = -1
394
        self.ibrion = ibrion
395
396
397
398
        if dft_plus_u:
            backend.addValue("electronic_structure_method", "DFT+U")
        else:
            backend.addValue("electronic_structure_method", "DFT")
399

400
    def onEnd_kpoints(self, parser, event, element, pathStr):
401
        backend = parser.backend
402
403
404
405
406
407
408
        self.bands = None
        self.kpoints = None
        self.weights = None
        for el in element:
            if el.tag == "generation":
                param = el.attrib.get("param", None)
                if param:
Markus Scheidgen's avatar
Markus Scheidgen committed
409
410
                    backend.addValue(
                        "x_vasp_k_points_generation_method", param)
411
412
                if param == "listgenerated":
                    self.bands = {
Markus Scheidgen's avatar
Markus Scheidgen committed
413
                        "divisions": g(el, "i/[@name='divisions']", None),
414
415
416
417
418
419
420
421
422
423
                        "points": getVector(el)
                    }
                elif param == "Monkhorst-Pack":
                    pass
                else:
                    backend.pwarn("unknown k point generation method")
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "kpointlist":
                    self.kpoints = np.asarray(getVector(el))
424
                    backend.addArrayValues("k_mesh_points", self.kpoints)
425
426
                elif name == "weights":
                    self.weights = np.asarray(getVector(el))
Markus Scheidgen's avatar
Markus Scheidgen committed
427
428
                    backend.addArrayValues(
                        "k_mesh_weights", self.weights.flatten())
429
430
431
432
433
                else:
                    backend.pwarn("Unknown array %s in kpoints" % name)
            else:
                backend.pwarn("Unknown tag %s in kpoints" % el.tag)

434
    def onEnd_structure(self, parser, event, element, pathStr):
435
        backend = parser.backend
436
        gIndexes = parser.tagSections[pathStr]
437
        self.lastSystemDescription = gIndexes["section_system"]
438
        self.cell = None
439
440
441
442
443
444
        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":
445
                            conv = convert_unit_function("angstrom", "m")
Markus Scheidgen's avatar
Markus Scheidgen committed
446
447
                            self.cell = getVector(
                                cellEl, lambda x: conv(float(x)))
448
                            self.angstrom_cell = np.array(getVector(cellEl))
Markus Scheidgen's avatar
Markus Scheidgen committed
449
450
451
452
                            backend.addArrayValues(
                                "simulation_cell", np.asarray(self.cell))
                            backend.addArrayValues(
                                "configuration_periodic_dimensions", np.ones(3, dtype=bool))
453
                        elif name == "rec_basis":
454
455
                            pass
                        else:
Markus Scheidgen's avatar
Markus Scheidgen committed
456
457
                            backend.pwarn(
                                "Unexpected varray %s in crystal" % name)
458
459
                    elif cellEl.tag == "i":
                        if cellEl.attrib.get("name") != "volume":
Markus Scheidgen's avatar
Markus Scheidgen committed
460
461
                            backend.pwarn(
                                "Unexpected i value %s in crystal" % cellEl.attrib)
462
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
463
464
                        backend.pwarn("Unexpected tag %s %s %r in crystal" % (
                            cellEl.tag, cellEl.attrib, cellEl.text))
465
466
467
468
            elif el.tag == "varray":
                name = el.attrib.get("name", None)
                if name == "positions":
                    pos = getVector(el)
Markus Scheidgen's avatar
Markus Scheidgen committed
469
470
                    backend.addArrayValues(
                        "atom_positions", np.dot(np.asarray(pos), self.cell))
471
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
472
473
                    backend.pwarn(
                        "Unexpected varray in structure %s" % el.attrib)
474
            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
475
476
                backend.pwarn("Unexpected tag in structure %s %s %r" %
                              el.tag, el.attrib, el.text)
477
        if self.labels is not None:
478
            backend.addArrayValues("atom_labels", self.labels)
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
479

480
    def onEnd_eigenvalues(self, parser, event, element, pathStr):
481
482
        if pathStr != "modeling/calculation/eigenvalues":
            return True
483
484
485
        backend = parser.backend
        eigenvalues = None
        occupation = None
486
487
488
489
490
491
492
493
        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":
494
                        isp = -1
495
496
                        for spinEl in arrEl:
                            if spinEl.tag == "set":
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
497
                                ik = -1
498
                                isp += 1
499
500
                                for kEl in spinEl:
                                    if kEl.tag == "set":
501
                                        ik += 1
Markus Scheidgen's avatar
Markus Scheidgen committed
502
503
                                        bands = np.asarray(
                                            getVector(kEl, field="r"))
504
                                        if eigenvalues is None:
Markus Scheidgen's avatar
Markus Scheidgen committed
505
506
507
508
509
510
                                            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]
511
                                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
512
513
                                        backend.pwarn(
                                            "unexpected tag %s in k array of the eigenvalues" % kEl.tag)
514
                            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
515
516
                                backend.pwarn(
                                    "unexpected tag %s in spin array of the eigenvalues" % spinEl.tag)
517
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
518
519
                        backend.pwarn(
                            "unexpected tag %s in array of the eigenvalues" % arrEl.tag)
520
                if eigenvalues is not None:
521

Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
522
                    ev = eV2JV(eigenvalues)
523
524
                    vbTopE = []
                    ebMinE = []
525
                    for ispin in range(occupation.shape[0]):
526
527
                        vbTopE.append(float('-inf'))
                        ebMinE.append(float('inf'))
528
                        for ik in range(occupation.shape[1]):
Markus Scheidgen's avatar
Markus Scheidgen committed
529
530
531
                            ebIndex = bisect.bisect_right(
                                -occupation[ispin, ik, :], -0.5) - 1
                            vbTopIndex = ebIndex - 1
532
533
                            if vbTopIndex >= 0:
                                vbTopK = ev[ispin, ik, vbTopIndex]
534
535
                                if vbTopK > vbTopE[ispin]:
                                    vbTopE[ispin] = vbTopK
536
537
                            if ebIndex < ev.shape[2]:
                                ebMinK = ev[ispin, ik, ebIndex]
538
539
                                if ebMinK < ebMinE[ispin]:
                                    ebMinE[ispin] = ebMinK
540
541
                    self.vbTopE = vbTopE
                    self.ebMinE = ebMinE
Markus Scheidgen's avatar
Markus Scheidgen committed
542
543
544
545
                    backend.addArrayValues(
                        "energy_reference_highest_occupied", np.array(vbTopE))
                    backend.addArrayValues(
                        "energy_reference_lowest_unoccupied", np.array(ebMinE))
546
547
548
                    if self.bands:
                        divisions = int(self.bands['divisions'])
                        backend.openNonOverlappingSection("section_k_band")
549
                        nsegments = self.kpoints.shape[0] // divisions
Markus Scheidgen's avatar
Markus Scheidgen committed
550
551
552
553
554
555
                        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]))
556
                        for isegment in range(nsegments):
Markus Scheidgen's avatar
Markus Scheidgen committed
557
558
559
560
561
562
563
564
                            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])
565
                            # "band_segm_labels"
Markus Scheidgen's avatar
Markus Scheidgen committed
566
567
568
569
                            backend.addArrayValues("band_segm_start_end", np.asarray(
                                [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment")
570
                        backend.closeNonOverlappingSection("section_k_band")
Markus Scheidgen's avatar
Markus Scheidgen committed
571
572
                        backend.openNonOverlappingSection(
                            "section_k_band_normalized")
573
574
                        specialPoints = {}
                        try:
Markus Scheidgen's avatar
Markus Scheidgen committed
575
576
577
578
                            specialPoints = get_special_points(
                                convert_unit_function("m", "angstrom")(self.cell))
                        except Exception as e:
                            self.logger.error("failed to get special points", exc_info=e)
579
                        for isegment in range(nsegments):
Markus Scheidgen's avatar
Markus Scheidgen committed
580
581
582
583
584
585
586
587
588
589
                            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]]))
590
591
592
                            backend.addValue("band_segm_labels_normalized",
                                             [findLabel(specialPoints, kpt[isegment, 0]),
                                              findLabel(specialPoints, kpt[isegment, divisions - 1])])
Markus Scheidgen's avatar
Markus Scheidgen committed
593
594
                            backend.closeNonOverlappingSection(
                                "section_k_band_segment_normalized")
595

Markus Scheidgen's avatar
Markus Scheidgen committed
596
597
                        backend.closeNonOverlappingSection(
                            "section_k_band_normalized")
598
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
599
600
                        backend.openNonOverlappingSection(
                            "section_eigenvalues")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
601
                        backend.addArrayValues("eigenvalues_values", ev)
Markus Scheidgen's avatar
Markus Scheidgen committed
602
603
604
605
                        backend.addArrayValues(
                            "eigenvalues_occupation", occupation)
                        backend.closeNonOverlappingSection(
                            "section_eigenvalues")
606
607
608
            else:
                backend.pwarn("unexpected tag %s in the eigenvalues" % el.tag)

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

612
613
    def onStart_calculation(self, parser, event, element, pathStr):
        gIndexes = parser.tagSections[pathStr]
Markus Scheidgen's avatar
Markus Scheidgen committed
614
615
        self.singleConfCalcs.append(
            gIndexes["section_single_configuration_calculation"])
616
617
        if self.waveCut:
            backend.openNonOverlappingSection("section_basis_set")
Markus Scheidgen's avatar
Markus Scheidgen committed
618
619
            backend.addValue(
                "mapping_section_basis_set_cell_dependent", self.waveCut)
620
621
            backend.closeNonOverlappingSection("section_basis_set")

622
623
624
625
626
627
628
629
630
631
632
633
634
    def onEnd_modeling(self, parser, event, element, pathStr):
        backend = parser.backend
        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
635
636
        backend.addArrayValues(
            "frame_sequence_local_frames_ref", np.asarray(self.singleConfCalcs))
637
638
639
        backend.closeSection("section_frame_sequence", frameSequenceGIndex)

    def onEnd_calculation(self, parser, event, element, pathStr):
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
640
        eConv = eV2J
641
642
643
        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
644
645
        backend.addValue(
            "single_configuration_calculation_to_system_ref", self.lastSystemDescription)
646
        gIndexes = parser.tagSections["/modeling"]
Markus Scheidgen's avatar
Markus Scheidgen committed
647
648
        backend.addValue(
            "single_configuration_to_calculation_method_ref", gIndexes["section_method"])
649
650
651
652
653
654
655
656
        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)
657
                        elif name == "e_wo_entrp":
658
659
660
661
662
663
                            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
664
665
                            backend.pwarn(
                                "Unexpected i tag with name %s in energy section" % name)
666
667
668
669
670
671
672
673
674
                    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)

675
    def onEnd_atominfo(self, parser, event, element, pathStr):
676
677
678
679
        nAtoms = None
        nAtomTypes = None
        atomTypes = []
        labels = []
680
        labels2 = None
681
        atomTypesDesc = []
682
        backend = parser.backend
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
        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
699
700
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
701
702
703
704
705
                                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
706
707
                            backend.pwarn(
                                "unexpected tag %s in atoms array in atominfo" % atomsEl.tag)
708
                elif name == "atomtypes":
709
710
711
712
713
714
715
                    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
716
717
                            fieldTypes.append(
                                atomsEl.attrib.get("type", "float"))
718
                        elif atomsEl.tag == "set":
Markus Scheidgen's avatar
Markus Scheidgen committed
719
720
                            expectedKeys = ["atomspertype", "element",
                                            "mass", "valence", "pseudopotential"]
721
                            if keys != expectedKeys:
Markus Scheidgen's avatar
Markus Scheidgen committed
722
723
                                backend.pwarn(
                                    "unexpected fields in atomtype: %s vs %s" % (keys, expectedKeys))
724
725
                            for atomsLine in atomsEl:
                                if atomsLine.tag != "rc":
Markus Scheidgen's avatar
Markus Scheidgen committed
726
727
                                    backend.pwarn(
                                        "unexpected tag %s in atoms array in atominfo" % atomsLine.tag)
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
                                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
743
744
                            backend.pwarn(
                                "unexpected tag %s in atomtypes array in atominfo" % atomsEl.tag)
745
                    kindIds = []
Markus Scheidgen's avatar
Markus Scheidgen committed
746
                    nEl = {}
747
                    kindLabels = []
748
                    for atomDesc in atomTypesDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
749
750
                        kindId = backend.openSection(
                            "section_method_atom_kind")
751
                        if 'element' in atomDesc:
752
                            elName = atomDesc['element'].strip()
753
                            try:
754
                                elNr = ase.data.chemical_symbols.index(elName)
Markus Scheidgen's avatar
Markus Scheidgen committed
755
756
757
758
759
760
761
                                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)
762
                            nEl[elName] = nElNow
Markus Scheidgen's avatar
Markus Scheidgen committed
763
764
                            elLabel = elName + \
                                (str(nElNow) if nElNow > 1 else "")
765
766
                            kindLabels.append(elLabel)
                            backend.addValue("method_atom_kind_label", elLabel)
767
                            if "mass" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
768
769
                                backend.addValue(
                                    "method_atom_kind_mass", atomDesc["mass"])
770
                            if "valence" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
771
772
                                backend.addValue(
                                    "method_atom_kind_explicit_electrons", atomDesc["valence"])
773
                            if "pseudopotential" in atomDesc:
Markus Scheidgen's avatar
Markus Scheidgen committed
774
775
                                backend.addValue(
                                    "method_atom_kind_pseudopotential_name", atomDesc["pseudopotential"])
776
                        kindIds.append(kindId)
Markus Scheidgen's avatar
Markus Scheidgen committed
777
778
779
780
                        backend.closeSection(
                            "section_method_atom_kind", kindId)
                    backend.addArrayValues("x_vasp_atom_kind_refs", np.asarray(
                        [kindIds[i-1] for i in atomTypes]))
781
                    labels2 = [kindLabels[i-1] for i in atomTypes]
782
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
783
784
                    backend.pwarn(
                        "unexpected array named %s in atominfo" % name)
785
786
            else:
                backend.pwarn("unexpected tag %s in atominfo" % el.tag)
787
        self.labels = np.asarray(labels2) if labels2 else np.asarray(labels)
788

789
    def incarOutTag(self, el):
790
791
        backend = self.parser.backend
        metaEnv = self.parser.backend.metaInfoEnv()
792
        if (el.tag != "i"):
Markus Scheidgen's avatar
Markus Scheidgen committed
793
794
            backend.pwarn("unexpected tag %s %s %r in incar" %
                          (el.tag, el.attrib, el.text))
795
796
797
798
799
        else:
            name = el.attrib.get("name", None)
            meta = metaEnv['x_vasp_incarOut_' + name]
            valType = el.attrib.get("type")
            if not meta:
Markus Scheidgen's avatar
Markus Scheidgen committed
800
                self.unknown_incar_params.append((el.tag, el.attrib, el.text))
801
802
803
804
805
            else:
                if valType:
                    expectedMetaType = {
                        'string': ['C'],
                        'int': ['i'],
Markus Scheidgen's avatar
Markus Scheidgen committed
806
                        'logical': ['b', 'C']
807
808
                    }.get(valType)
                    if not expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
809
810
                        backend.pwarn("Unknown value type %s encountered in INCAR out: %s %s %r" % (
                            valType, el.tag, el.attrib, el.text))
811
                    elif not meta.get('dtypeStr') in expectedMetaType:
Markus Scheidgen's avatar
Markus Scheidgen committed
812
813
                        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))
814
                try:
815
816
817
818
                    shape = meta.get("shape", None)
                    dtypeStr = meta.get("dtypeStr", None)
                    converter = metaTypeTransformers.get(dtypeStr)
                    if not converter:
Markus Scheidgen's avatar
Markus Scheidgen committed
819
820
                        backend.pwarn(
                            "could not find converter for dtypeStr %s when handling meta info %s" % (dtypeStr, ))
821
822
                    elif shape:
                        vals = re.split("\s+", el.text.strip())
Markus Scheidgen's avatar
Markus Scheidgen committed
823
824
                        backend.addValue(
                            meta["name"], [converter(x) for x in vals])
825
826
                    else:
                        backend.addValue(meta["name"], converter(el.text))
827
                except:
Markus Scheidgen's avatar
Markus Scheidgen committed
828
829
                    backend.pwarn("Exception trying to handle incarOut %s: %s" % (
                        name, traceback.format_exc()))
830
                if name == 'ENMAX' or name == 'PREC':
Markus Scheidgen's avatar
Markus Scheidgen committed
831
832
833
834
835
836
837
                    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
838
839
840
841
842
843
                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'],
Markus Scheidgen's avatar
Markus Scheidgen committed
844
                        '--': ['GGA_X_PBE', 'GGA_C_PBE']  # should check potcar
845
846
847
                    }
                    functs = fMap.get(el.text.strip(), None)
                    if not functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
848
849
                        backend.pwarn("Unknown XC functional %s" %
                                      el.text.strip())
850
851
                    else:
                        for f in functs:
Markus Scheidgen's avatar
Markus Scheidgen committed
852
853
                            backend.openNonOverlappingSection(
                                "section_XC_functionals")
854
                            backend.addValue("XC_functional_name", f)
Markus Scheidgen's avatar
Markus Scheidgen committed
855
856
                            backend.closeNonOverlappingSection(
                                "section_XC_functionals")
857
858
859
                elif name == "ISPIN":
                    self.ispin = int(el.text.strip())

Markus Scheidgen's avatar
Markus Scheidgen committed
860
    def separatorScan(self, element, backend, depth=0):
861
862
863
864
865
866
867
        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
868
                        self.separatorScan(el, backend, depth + 1)
869
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
870
871
872
                        # backend.pwarn("unexpected tag %s %s in parameters separator %s at depth %d" % (
                        #     el.tag, el.attrib, separatorName, depth))
                        pass
873
874
875
            elif separators.tag == "i":
                self.incarOutTag(separators)
            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
876
877
878
                # backend.pwarn("unexpected tag %s %s in parameters at depth %d" % (
                #     separators.tag, separators.attrib, depth))
                pass
879

880
    def onEnd_parameters(self, parser, event, element, pathStr):
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
881
        self.separatorScan(element, parser.backend)
882
        backend = parser.backend
883
        try:
Markus Scheidgen's avatar
Markus Scheidgen committed
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
            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:
                backend.pwarn(
                    "Missing ENMAX for calculating plane wave basis cut off ")
900
        except AttributeError:
Markus Scheidgen's avatar
Markus Scheidgen committed
901
902
            backend.pwarn(
                "Missing PREC for calculating plane wave basis cut off ")
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
903
904
905

    def onEnd_dos(self, parser, event, element, pathStr):
        "density of states"
906
        backend = parser.backend
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
907
908
909
910
        backend.openNonOverlappingSection("section_dos")
        for el in element:
            if el.tag == "i":
                if el.attrib.get("name") == "efermi":
911
912
                    self.eFermi = eV2J(float(el.text.strip()))
                    backend.addValue("dos_fermi_energy", self.eFermi)
Markus Scheidgen's avatar
Markus Scheidgen committed
913
914
                    backend.addArrayValues(
                        "energy_reference_fermi", np.array([self.eFermi]*self.ispin))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
915
                else:
Markus Scheidgen's avatar
Markus Scheidgen committed
916
917
                    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
918
919
920
921
922
923
924
925
926
927
            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
928
929
                                        dosL.append(
                                            getVector(spinComponent, field="r"))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
930
                                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
931
932
                                        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
933
934
                                dosA = np.asarray(dosL)
                                if len(dosA.shape) != 3:
Markus Scheidgen's avatar
Markus Scheidgen committed
935
936
                                    raise Exception("unexpected shape %s (%s) for total dos (ragged arrays?)" % (
                                        dosA.shape), dosA.dtype)
937
938
939
940
941
                                dosE = eV2JV(dosA[0, :, 0])
                                dosI = dosA[:, :, 2]
                                dosV = dosA[:, :, 1]

                                # Convert the DOS values to SI. VASP uses the
942
943
944
945
                                # following units in the output:
                                # states/eV/cell. This means that the volume
                                # dependence has been introduced by multiplying
                                # by the cell volume
946
947
                                # the integrated dos value is the number of electrons until that energy level
                                # and thus not directly energy dependent anymore
948
949
                                joule_in_ev = convert_unit(1, "eV", "J")
                                dosV = dosV / joule_in_ev
950

951
952
                                if self.vbTopE:
                                    eRef = max(self.vbTopE)
953
954
                                else:
                                    eRef = self.eFermi
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
955
                                backend.addArrayValues("dos_energies", dosE)
Markus Scheidgen's avatar
Markus Scheidgen committed
956
957
                                backend.addArrayValues(
                                    "dos_energies_normalized", dosE - eRef)
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
958
                                backend.addArrayValues("dos_values", dosV)
Markus Scheidgen's avatar
Markus Scheidgen committed
959
960
                                backend.addArrayValues(
                                    "dos_integrated_values", dosI)
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
961
                            else:
Markus Scheidgen's avatar
Markus Scheidgen committed
962
963
                                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
964
                    else:
Markus Scheidgen's avatar
Markus Scheidgen committed
965
966
                        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
967
968
969
            elif el.tag == "partial":
                for el1 in el:
                    if el1.tag == "array":
Markus Scheidgen's avatar
Markus Scheidgen committed
970
                        lm = []
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
971
972
973
974
975
976
977
978
                        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
979
980
981
982
983
984
985
986
987
988
989
990
                                        "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
991
                                    }
Markus Scheidgen's avatar
Markus Scheidgen committed
992
993
                                    lm.append(
                                        strLm.get(el2.text.strip(), [-1, -1]))
Mohamed, Fawzi Roberto (fawzi)'s avatar
add dos    
Mohamed, Fawzi Roberto (fawzi) committed
994
995
996
997
998
999
1000
                            elif el2.tag == "set":
                                dosL = []
                                for atom in el2:
                                    if atom.tag == "set":
                                        atomL = []
                                        dosL.append(atomL)
                                        for spinComponent in atom:
For faster browsing, not all history is shown. View entire blame