diff --git a/vaspparser/metainfo/vasp_incars.py b/vaspparser/metainfo/vasp_incars.py
index f4d958d6a530aa48d944c6eed434b0494662ad0c..21ee88c0e62de2a187f3717c954833aaa6b65c0a 100644
--- a/vaspparser/metainfo/vasp_incars.py
+++ b/vaspparser/metainfo/vasp_incars.py
@@ -856,7 +856,7 @@ class section_method(public.section_method):
x_vasp_incar_KPOINT_BSE = Quantity(
type=np.dtype(np.int32),
- shape=['0..3'],
+ shape=[],
description='''
The flag KPOINT_BSE allows to calculate the dielectric matrix at one of the
kpoints used to sample the Brillouin zone. NOTE: Either specify one or three
diff --git a/vaspparser/parser_vasprun.py b/vaspparser/parser_vasprun.py
index b4c96b055406668492eb6bd6966266de7157d25d..50ce84c4b52855f5d94bc49d342e8a78e675f419 100644
--- a/vaspparser/parser_vasprun.py
+++ b/vaspparser/parser_vasprun.py
@@ -37,12 +37,161 @@ from nomadcore.unit_conversion.unit_conversion import convert_unit
eV2J = convert_unit_function("eV", "J")
eV2JV = np.vectorize(eV2J)
+
+
+def crystal_structure_from_cell(cell, eps=1e-4):
+ """Return the crystal structure as a string calculated from the cell.
+ """
+ cellpar = ase.geometry.cell_to_cellpar(cell=cell)
+ abc = cellpar[:3]
+ angles = cellpar[3:] / 180 * pi
+ a, b, c = abc
+ alpha, beta, gamma = angles
+ # According to:
+ # https://www.physics-in-a-nutshell.com/article/6/symmetry-crystal-systems-and-bravais-lattices#the-seven-crystal-systems
+ # If a=b=c and alpha=beta=gamma=90degrees we have cubic.
+ if abc.ptp() < eps and abs(angles - pi / 2).max() < eps:
+ return 'cubic'
+ elif abc.ptp() < eps and abs(angles - pi / 3).max() < eps:
+ return 'fcc'
+ elif abc.ptp() < eps and abs(angles - np.arccos(-1 / 3)).max() < eps:
+ return 'bcc'
+ # If a=b!=c, alpha=beta=gamma=90deg, tetragonal.
+ elif abs(a - b) < eps and abs(angles - pi / 2).max() < eps:
+ return 'tetragonal'
+ elif abs(angles - pi / 2).max() < eps:
+ return 'orthorhombic'
+ # if a = b != c , alpha = beta = 90deg, gamma = 120deg, hexagonal
+ elif (abs(a - b) < eps
+ and abs(gamma - pi / 3 * 2) < eps
+ and abs(angles[:2] - pi / 2).max() < eps):
+ return 'hexagonal'
+ elif (c >= a and c >= b and alpha < pi / 2 and
+ abs(angles[1:] - pi / 2).max() < eps):
+ return 'monoclinic'
+ else:
+ raise ValueError('Cannot find crystal structure')
+
vasp_to_metainfo_type_mapping = {
'string': ['C'],
'int': ['i'],
'logical': ['b', 'C'],
'float': ['f']}
+special_points = {
+ 'cubic': {'Γ': [0, 0, 0],
+ 'M': [1 / 2, 1 / 2, 0],
+ 'R': [1 / 2, 1 / 2, 1 / 2],
+ 'X': [0, 1 / 2, 0]},
+ 'fcc': {'Γ': [0, 0, 0],
+ 'K': [3 / 8, 3 / 8, 3 / 4],
+ 'L': [1 / 2, 1 / 2, 1 / 2],
+ 'U': [5 / 8, 1 / 4, 5 / 8],
+ 'W': [1 / 2, 1 / 4, 3 / 4],
+ 'X': [1 / 2, 0, 1 / 2]},
+ 'bcc': {'Γ': [0, 0, 0],
+ 'H': [1 / 2, -1 / 2, 1 / 2],
+ 'P': [1 / 4, 1 / 4, 1 / 4],
+ 'N': [0, 0, 1 / 2]},
+ 'tetragonal': {'Γ': [0, 0, 0],
+ 'A': [1 / 2, 1 / 2, 1 / 2],
+ 'M': [1 / 2, 1 / 2, 0],
+ 'R': [0, 1 / 2, 1 / 2],
+ 'X': [0, 1 / 2, 0],
+ 'Z': [0, 0, 1 / 2]},
+ 'orthorhombic': {'Γ': [0, 0, 0],
+ 'R': [1 / 2, 1 / 2, 1 / 2],
+ 'S': [1 / 2, 1 / 2, 0],
+ 'T': [0, 1 / 2, 1 / 2],
+ 'U': [1 / 2, 0, 1 / 2],
+ 'X': [1 / 2, 0, 0],
+ 'Y': [0, 1 / 2, 0],
+ 'Z': [0, 0, 1 / 2]},
+ 'hexagonal': {'Γ': [0, 0, 0],
+ 'A': [0, 0, 1 / 2],
+ 'H': [1 / 3, 1 / 3, 1 / 2],
+ 'K': [1 / 3, 1 / 3, 0],
+ 'L': [1 / 2, 0, 1 / 2],
+ 'M': [1 / 2, 0, 0]}}
+
+
+special_paths = {
+ 'cubic': 'ΓXMΓRX,MR',
+ 'fcc': 'ΓXWKΓLUWLK,UX',
+ 'bcc': 'ΓHNΓPH,PN',
+ 'tetragonal': 'ΓXMΓZRAZXR,MA',
+ 'orthorhombic': 'ΓXSYΓZURTZ,YT,UX,SR',
+ 'hexagonal': 'ΓMKΓALHA,LM,KH',
+ 'monoclinic': 'ΓYHCEM1AXH1,MDZ,YD'}
+
+
+def get_special_points(cell, eps=1e-4):
+ """Return dict of special points.
+
+ The definitions are from a paper by Wahyu Setyawana and Stefano
+ Curtarolo::
+
+ http://dx.doi.org/10.1016/j.commatsci.2010.05.010
+
+ lattice: str
+ One of the following: cubic, fcc, bcc, orthorhombic, tetragonal,
+ hexagonal or monoclinic.
+ cell: 3x3 ndarray
+ Unit cell.
+ eps: float
+ Tolerance for cell-check.
+ """
+
+ lattice = crystal_structure_from_cell(cell)
+
+ cellpar = ase.geometry.cell_to_cellpar(cell=cell)
+ abc = cellpar[:3]
+ angles = cellpar[3:] / 180 * pi
+ a, b, c = abc
+ alpha, beta, gamma = angles
+
+ # Check that the unit-cells are as in the Setyawana-Curtarolo paper:
+ if lattice == 'cubic':
+ assert abc.ptp() < eps and abs(angles - pi / 2).max() < eps
+ elif lattice == 'fcc':
+ assert abc.ptp() < eps and abs(angles - pi / 3).max() < eps
+ elif lattice == 'bcc':
+ angle = np.arccos(-1 / 3)
+ assert abc.ptp() < eps and abs(angles - angle).max() < eps
+ elif lattice == 'tetragonal':
+ assert abs(a - b) < eps and abs(angles - pi / 2).max() < eps
+ elif lattice == 'orthorhombic':
+ assert abs(angles - pi / 2).max() < eps
+ elif lattice == 'hexagonal':
+ assert abs(a - b) < eps
+ assert abs(gamma - pi / 3 * 2) < eps
+ assert abs(angles[:2] - pi / 2).max() < eps
+ elif lattice == 'monoclinic':
+ assert c >= a and c >= b
+ assert alpha < pi / 2 and abs(angles[1:] - pi / 2).max() < eps
+ if lattice != 'monoclinic':
+ return special_points[lattice]
+
+ # Here, we need the cell:
+ eta = (1 - b * np.cos(alpha) / c) / (2 * np.sin(alpha)**2)
+ nu = 1 / 2 - eta * c * np.cos(alpha) / b
+ return {'Γ': [0, 0, 0],
+ 'A': [1 / 2, 1 / 2, 0],
+ 'C': [0, 1 / 2, 1 / 2],
+ 'D': [1 / 2, 0, 1 / 2],
+ 'D1': [1 / 2, 0, -1 / 2],
+ 'E': [1 / 2, 1 / 2, 1 / 2],
+ 'H': [0, eta, 1 - nu],
+ 'H1': [0, 1 - eta, nu],
+ 'H2': [0, eta, -nu],
+ 'M': [1 / 2, eta, 1 - nu],
+ 'M1': [1 / 2, 1 - eta, nu],
+ 'M2': [1 / 2, eta, -nu],
+ 'X': [0, 1 / 2, 0],
+ 'Y': [0, 0, 1 / 2],
+ 'Y1': [0, 0, -1 / 2],
+ 'Z': [1 / 2, 0, 0]}
+
def findLabel(labels, value):
for k, v in labels.items():
@@ -466,11 +615,39 @@ class VasprunContext(object):
"band_occupations", occ[:, isegment, :, :])
backend.addArrayValues(
"band_k_points", kpt[isegment])
+ # "band_segm_labels"
backend.addArrayValues("band_segm_start_end", np.asarray(
[kpt[isegment, 0], kpt[isegment, divisions - 1]]))
backend.closeNonOverlappingSection(
"section_k_band_segment")
backend.closeNonOverlappingSection("section_k_band")
+ backend.openNonOverlappingSection(
+ "section_k_band_normalized")
+ specialPoints = {}
+ try:
+ specialPoints = get_special_points(
+ convert_unit_function("m", "angstrom")(self.cell))
+ except Exception as e:
+ self.logger.warn("failed to get special points/xtal structure", exc_info=e)
+ for isegment in range(nsegments):
+ backend.openNonOverlappingSection(
+ "section_k_band_segment_normalized")
+ backend.addArrayValues(
+ "band_energies_normalized", energies[:, isegment, :, :] - max(self.vbTopE))
+ backend.addArrayValues(
+ "band_occupations_normalized", occ[:, isegment, :, :])
+ backend.addArrayValues(
+ "band_k_points_normalized", kpt[isegment])
+ backend.addArrayValues("band_segm_start_end_normalized", np.asarray(
+ [kpt[isegment, 0], kpt[isegment, divisions - 1]]))
+ backend.addValue("band_segm_labels_normalized",
+ [findLabel(specialPoints, kpt[isegment, 0]),
+ findLabel(specialPoints, kpt[isegment, divisions - 1])])
+ backend.closeNonOverlappingSection(
+ "section_k_band_segment_normalized")
+
+ backend.closeNonOverlappingSection(
+ "section_k_band_normalized")
else:
backend.openNonOverlappingSection(
"section_eigenvalues")