Commit 7766a086 authored by Lauri Himanen's avatar Lauri Himanen
Browse files

Moved band segment labeling from the VASP parser into BandStructureNormalizer,...

Moved band segment labeling from the VASP parser into BandStructureNormalizer, added proper tests for label detection, fixed issues with detecting band structure calculations.
parent 3d76aad4
Pipeline #73910 passed with stages
in 37 minutes and 48 seconds
Subproject commit 1e26aeaf2bbb27826075d02c0a80c2b04f482d20
Subproject commit d42a5739be64bfc74b383122b10bf0d1b0a3073d
......@@ -14,17 +14,19 @@
import json
import numpy as np
import ase
from nomad.datamodel.metainfo.public import section_k_band, section_band_gap, section_system
from nomad.normalizing.normalizer import Normalizer
from nomad import config, atomutils
from nomad.constants import pi
class BandStructureNormalizer(Normalizer):
"""Normalizer with the following responsibilities:
- Calculates band gap(s) if present (section_band_gap, section_band_gap_spin_up, section_band_gap_spin_down).
- TODO: Creates labels for special points within the band path (band_path_labels)
- TODO: Creates labels for special points within the band path (band_segm_labels)
- TODO: Determines if the path is a standard one or not (is_standard)
"""
def normalize(self, logger=None) -> None:
......@@ -44,7 +46,7 @@ class BandStructureNormalizer(Normalizer):
valence_band_maximum = scc.energy_reference_highest_occupied
# In order to resolve the special points and the reciprocal cell,
# we need informatoin about the system.
# we need information about the system.
system = scc.single_configuration_calculation_to_system_ref
for band in scc.section_k_band:
......@@ -52,6 +54,7 @@ class BandStructureNormalizer(Normalizer):
self.add_reciprocal_cell(band, system)
self.add_brillouin_zone(band)
self.add_band_gaps(band, valence_band_maximum)
self.add_path_labels(band, system)
def add_reciprocal_cell(self, band: section_k_band, system: section_system):
"""A reciprocal cell for this calculation. If the original unit cell is
......@@ -231,3 +234,189 @@ class BandStructureNormalizer(Normalizer):
# Add section to band
band.m_add_sub_section(section_k_band.section_band_gap, gap)
def add_path_labels(self, band: section_k_band, system: section_system) -> None:
"""Adds special high symmmetry point labels to the band path. Only k
points that land on the special points defined by Setyawan/Curtarolo
are automatically labeled.
"""
# If labels are already set by the parser dot nothing.
for segment in band.section_k_band_segment:
labels = segment.band_segm_labels
if labels is not None:
self.logger.info("Existing band segment labels detected, skipping label detection.")
return
# Try to get the required data. Fail if not found.
try:
cell = system.lattice_vectors.to("angstrom").magnitude
reciprocal_cell_trans = band.reciprocal_cell.magnitude.T
bravais_lattice = system.section_symmetry[0].bravais_lattice
except Exception:
self.logger.info("Could not resolve path labels as required information is missing.")
return
# Find special points for this lattice. If an error occurs, the labels
# are simply not written.
try:
special_points = self.get_special_points(bravais_lattice, cell)
except Exception as e:
self.logger.warning("Could not resolve high-symmetry points for the given simulation cell.", exception=e)
return
# Form a contiguous array of k points for faster operations
special_point_labels = list(special_points.keys())
special_k_points = np.empty((len(special_points), 3))
for i, kpt in enumerate(special_points.values()):
special_k_points[i, :] = kpt
special_k_points_cartesian = np.dot(special_k_points, reciprocal_cell_trans)
# Match tolerance in 1/m. Taken from the VASP parser.
eps = config.normalize.k_space_precision
# Try to find matches for the special points. We only attempt to match
# points at the start and end of a segment. Any labels set by the
# parser are overridden, because one cannot ascertain that those labels
# are consistent across codes.
for segment in band.section_k_band_segment:
seg_k_points = segment.band_k_points
if seg_k_points is None:
self.logger.info("Could not resolve band path as k points are missing.")
return
start_point_cartesian = np.dot(segment.band_k_points[0], reciprocal_cell_trans)
end_point_cartesian = np.dot(segment.band_k_points[-1], reciprocal_cell_trans)
# Calculate distance in cartesian space
start_index = atomutils.find_match(start_point_cartesian, special_k_points_cartesian, eps)
end_index = atomutils.find_match(end_point_cartesian, special_k_points_cartesian, eps)
if start_index is None:
start_label = ""
else:
start_label = special_point_labels[start_index]
if end_index is None:
end_label = ""
else:
end_label = special_point_labels[end_index]
segment.band_segm_labels = [start_label, end_label]
def get_special_points(self, bravais_lattice, 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
bravais_lattice: str
bravais lattice in Pearson notation.
cell: 3x3 ndarray
Unit cell.
eps: float
Tolerance for cell-check.
"""
# Special points that do not depend on lattice parameters. TODO: A lot
# of the bravais lattice are missing from this implementation that is
# copied from the VASP parser.
special_points = {
'cP': {
'Γ': [0, 0, 0],
'M': [1 / 2, 1 / 2, 0],
'R': [1 / 2, 1 / 2, 1 / 2],
'X': [0, 1 / 2, 0]
},
'cF': {
'Γ': [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]
},
'cI': {
'Γ': [0, 0, 0],
'H': [1 / 2, -1 / 2, 1 / 2],
'P': [1 / 4, 1 / 4, 1 / 4],
'N': [0, 0, 1 / 2]
},
'tP': {
'Γ': [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]
},
'oP': {
'Γ': [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]
},
'hP': {
'Γ': [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]
}
}
cellpar = ase.geometry.cell_to_cellpar(cell=cell)
abc = cellpar[:3]
angles = cellpar[3:] / 180 * pi
a, b, c = abc
alpha, _, gamma = angles
# Check that the unit cells are as in the Setyawana-Curtarolo paper:
if bravais_lattice == 'cP':
assert abc.ptp() < eps and abs(angles - pi / 2).max() < eps
elif bravais_lattice == 'cF':
assert abc.ptp() < eps and abs(angles - pi / 3).max() < eps
elif bravais_lattice == 'cI':
angle = np.arccos(-1 / 3)
assert abc.ptp() < eps and abs(angles - angle).max() < eps
elif bravais_lattice == 'tP':
assert abs(a - b) < eps and abs(angles - pi / 2).max() < eps
elif bravais_lattice == 'oP':
assert abs(angles - pi / 2).max() < eps
elif bravais_lattice == 'hP':
assert abs(a - b) < eps
assert abs(gamma - pi / 3 * 2) < eps
assert abs(angles[:2] - pi / 2).max() < eps
elif bravais_lattice == 'mP':
sin_alpha = np.sin(alpha)
cos_alpha = np.cos(alpha)
assert c >= a and c >= b
assert alpha < pi / 2
assert alpha < pi / 2
assert (np.abs(angles[1:] - pi / 2) < eps).all()
eta = (1 - b * cos_alpha / c) / (2 * sin_alpha**2)
nu = 1 / 2 - eta * c * 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]
}
return special_points[bravais_lattice]
......@@ -17,7 +17,6 @@ import json
from nomad.metainfo.encyclopedia import (
Calculation,
Properties,
Material,
)
from nomad.parsing.legacy import Backend
from nomad.metainfo import Section
......@@ -33,53 +32,60 @@ class PropertiesNormalizer():
self.logger = logger
def electronic_band_structure(self, properties: Properties, calc_type: str, material_type: str, context: Context, sec_system: Section) -> None:
"""Tries to resolve a reference to a representative band structure.
Will not store any data if it cannot be resolved unambiguously.
"""Tries to resolve a reference to a representative electronic band
structure. Will loop through the available band in reversed order until
a valid band is found.
"""
# Band structure data is extracted only from bulk calculations. This is
# because for now we need the reciprocal cell of the primitive cell
# that is only available from the symmetry analysis. Once the
# reciprocal cell is directly reported with the band structure this
# restriction can go away.
if calc_type != Calculation.calculation_type.type.single_point or material_type != Material.material_type.type.bulk:
return
try:
representative_scc = context.representative_scc
bands = representative_scc.section_k_band
if bands is None:
return
# Try to find a band structure.
representative_scc = context.representative_scc
bands = representative_scc.section_k_band
if bands is None or len(bands) == 0:
representative_band = None
for band in reversed(bands):
kind = band.band_structure_kind
if kind != "vibrational":
valid = True
for segment in band.section_k_band_segment:
energies = segment.band_energies
k_points = segment.band_k_points
labels = segment.band_segm_labels
if energies is None or k_points is None or labels is None:
valid = False
if valid:
representative_band = band
break
except Exception:
return
if representative_band is not None:
properties.electronic_band_structure = representative_band
# Loop over bands
representative_band = None
for band in bands:
# Check segments
segments = band.section_k_band_segment
if segments is None or len(segments) == 0:
continue
# Check data in segments
valid_band = True
for segment in segments:
try:
segment.band_k_points
segment.band_energies
except Exception:
valid_band = False
break
def electronic_dos(self, properties: Properties, context: Context) -> None:
"""Tries to resolve a reference to a representative electonic density
of states. Will not store any data if it cannot be resolved
unambiguously.
"""
try:
representative_scc = context.representative_scc
doses = representative_scc.section_dos
if doses is None:
return
# Save a valid band as representative, if multiple encountered log
# a warning and do not save anything.
if valid_band:
if valid_band is None:
representative_band = band
else:
self.logger.warn("Could not unambiguously select data to display for band structure.")
return
representative_dos = None
for dos in reversed(doses):
kind = dos.dos_kind
energies = dos.dos_energies
values = dos.dos_values
if kind != "vibrational" and energies is not None and values is not None:
representative_dos = dos
break
if representative_band is not None:
properties.electronic_band_structure = representative_band
except Exception:
return
if representative_dos is not None:
properties.electronic_dos = representative_dos
def elastic_constants_matrix(self) -> None:
pass
......@@ -135,19 +141,52 @@ class PropertiesNormalizer():
return
representative_phonon_band = None
for band in bands:
if band.band_structure_kind == "vibrational":
if representative_phonon_band is None:
for band in reversed(bands):
kind = band.band_structure_kind
if kind == "vibrational":
valid = True
for segment in band.section_k_band_segment:
energies = segment.band_energies
k_points = segment.band_k_points
labels = segment.band_segm_labels
if energies is None or k_points is None or labels is None or "?" in labels:
valid = False
if valid:
representative_phonon_band = band
else:
self.logger("Could not unambiguously select data to display for a phonon band structure.")
return
break
except Exception:
return
if representative_phonon_band is not None:
properties.phonon_band_structure = representative_phonon_band
def phonon_dos(self, properties: Properties, context: Context) -> None:
"""Tries to resolve a reference to a representative phonon density of
states. Will not store any data if it cannot be resolved unambiguously.
"""
try:
representative_scc = context.representative_scc
doses = representative_scc.section_dos
if doses is None:
return
representative_phonon_dos = None
for dos in doses:
kind = dos.dos_kind
energies = dos.dos_energies
values = dos.dos_values
if kind == "vibrational" and energies is not None and values is not None:
if representative_phonon_dos is None:
representative_phonon_dos = dos
else:
self.logger("Could not unambiguously select data to display for phonon density of states.")
return
except Exception:
return
if representative_phonon_dos is not None:
properties.phonon_dos = representative_phonon_dos
def energies(self, properties: Properties, gcd: int, representative_scc: Section) -> None:
energy_dict = {}
if representative_scc is not None:
......@@ -182,9 +221,11 @@ class PropertiesNormalizer():
# Save metainfo
self.electronic_band_structure(properties, calc_type, material_type, context, sec_system)
self.electronic_dos(properties, context)
self.energies(properties, gcd, representative_scc)
# Phonon calculations have a specific set of properties to extract
if context.calc_type == Calculation.calculation_type.type.phonon_calculation:
self.thermodynamical_properties(properties)
self.phonon_band_structure(properties, context)
self.phonon_dos(properties, context)
......@@ -115,13 +115,18 @@ class SystemBasedNormalizer(Normalizer, metaclass=ABCMeta):
except Exception:
frame_seqs = []
# If no frame sequences detected, try to find scc
# If no frame sequences detected, try to find valid scc by looping all
# available in reverse order until a valid one is found.
if len(frame_seqs) == 0:
try:
sccs = self.section_run.section_single_configuration_calculation
scc_idx = -1
scc = sccs[scc_idx]
system = scc.single_configuration_calculation_to_system_ref
for scc in reversed(sccs):
idx = scc.m_parent_index
isys = scc.single_configuration_calculation_to_system_ref
if isys is not None:
scc_idx = idx
system = isys
break
if system is None:
sccs = []
except Exception:
......
This diff is collapsed.
0.000000000 -0.9140833100
0.000000000 0.7717749073
0.4330295261 -0.9140833100
0.4330295261 0.7717749073
0.9633802174 -0.9140833100
0.9633802174 0.7717749073
1.575776446 -0.9140833100
1.575776446 0.7717749073
1.881974561 -0.9140833100
1.881974561 0.7717749073
2.098489324 -0.9140833100
2.098489324 0.7717749073
Distance between is = 1 (Si), ia = 1 and
is = 1 (Si), ia = 1 : 0.000000000
is = 1 (Si), ia = 2 : 4.442710321
Distance between is = 1 (Si), ia = 2 and
is = 1 (Si), ia = 1 : 4.442710321
is = 1 (Si), ia = 2 : 0.000000000
8 : nkpt
10 : nstsv
1 0.000000000 0.000000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.2497043491 2.000000000
2 0.1905948812 2.000000000
3 0.1905948817 2.000000000
4 0.1905948817 2.000000000
5 0.2828036167 0.000000000
6 0.2828036167 0.000000000
7 0.2828036167 0.000000000
8 0.3082264419 0.000000000
9 0.4718777237 0.000000000
10 0.4744426734 0.000000000
2 0.2500000000 0.000000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.2202526954 2.000000000
2 0.4518078475E-01 2.000000000
3 0.1625370114 2.000000000
4 0.1625370118 2.000000000
5 0.2608986433 0.000000000
6 0.3173979562 0.000000000
7 0.3173979563 0.000000000
8 0.4391479690 0.000000000
9 0.4535013189 0.000000000
10 0.4535013189 0.000000000
3 0.5000000000 0.000000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.1636113809 2.000000000
2 -0.6705647201E-01 2.000000000
3 0.1464645238 2.000000000
4 0.1464645242 2.000000000
5 0.2425344126 0.000000000
6 0.3110259881 0.000000000
7 0.3110259882 0.000000000
8 0.4656584713 0.000000000
9 0.5772057897 0.000000000
10 0.5772057897 0.000000000
4 0.2500000000 0.2500000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.2098782040 2.000000000
2 0.6064473086E-01 2.000000000
3 0.1206981399 2.000000000
4 0.1206981401 2.000000000
5 0.2281393790 0.000000000
6 0.3013954263 0.000000000
7 0.3999423456 0.000000000
8 0.3999423457 0.000000000
9 0.4822779571 0.000000000
10 0.5187767035 0.000000000
5 0.5000000000 0.2500000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.1495475026 2.000000000
2 -0.4478041538E-01 2.000000000
3 0.5924573236E-01 2.000000000
4 0.1091678230 2.000000000
5 0.2421724898 0.000000000
6 0.3503493565 0.000000000
7 0.3981897514 0.000000000
8 0.4058414463 0.000000000
9 0.5601507885 0.000000000
10 0.5757806497 0.000000000
6 0.7500000000 0.2500000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.1744052011 2.000000000
2 -0.1152956982E-01 2.000000000
3 0.5102065019E-01 2.000000000
4 0.1416083986 2.000000000
5 0.2768689325 0.000000000
6 0.3643197265 0.000000000
7 0.3756534945 0.000000000
8 0.4115824549 0.000000000
9 0.4256175595 0.000000000
10 0.5399567625 0.000000000
7 0.5000000000 0.5000000000 0.000000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.9721130198E-01 2.000000000
2 -0.9721130189E-01 2.000000000
3 0.8528953629E-01 2.000000000
4 0.8528953651E-01 2.000000000
5 0.2119290288 0.000000000
6 0.2119290289 0.000000000
7 0.5572021961 0.000000000
8 0.5572021962 0.000000000
9 0.5922517179 0.000000000
10 0.5922517179 0.000000000
8 0.7500000000 0.5000000000 0.2500000000 : k-point, vkl
(state, eigenvalue and occupancy below)
1 -0.9102030155E-01 2.000000000
2 -0.9102030144E-01 2.000000000
3 0.4734528698E-01 2.000000000
4 0.4734528716E-01 2.000000000
5 0.3437660776 0.000000000
6 0.3437660777 0.000000000
7 0.3711661333 0.000000000
8 0.3711661335 0.000000000
9 0.5758927616 0.000000000
10 0.5758927617 0.000000000